Archive for the ‘Computer Algebra’ Category

A readable account of TeX’s math layout algorithm

Saturday, December 1st, 2007

Via Ars Mathematica, I came across a paper describing TeX’s math layout algorithm and offering an implementation of it in Standard ML. Sweet! I only wish I had come across this paper (it is about 10 years old) when I was daydreaming about a WYSIWYG gui interface that worked for multiple CASes. If I ever get back around to it, this will be really useful.

Once more into the fray

Monday, June 18th, 2007

I was recently invited to help rework BasicCAS into a full-fledged Mathematica parser for the SymPy project, a challenge which I have accepted. The idea is that a Mathematica parser (along with a working backend, a much harder goal to meet, I’d say) would allow SymPy to benefit from the large body of free Mathematica code floating around out there. Part of the reason I’m taking this on is because I’d not like to see SymPy turn into a melange like Sage– the idea of tying together disparate CASes and mathematical libraries with a glue language is aesthetically unappealing to me, and pragmatically, what’s the point of dividing development effort between two Python-based metaCASes? A well-defined global syntax is the biggest preventative measure I can think of.

So now I have some fun goals for the summer, to round out my list of not-so-fun goals: figure out a grammar for Mathematica’s syntax (this may be, literally, impossible: I have heard it said that the only grammar Mathematica follows is the ex post facto grammar imposed by its parser’s behavior); write a parser for that grammar, test it; track down the inevitable incompatibilities, attempt to improve the grammar; recurse. It’s going to be a bitch, but fun, hopefully!

Aldor (free it!)

Saturday, August 26th, 2006

The programming language Aldor, aka A#, seems to be a well-kept secret. It started off as the compiler for Axiom, Axiom-xl, but has since developed into an interesting standalone language in its own right. It reminds me of Scheme, C/C++, Python, and Perl 6 to varying degrees. In particular (with respect to being similar to Perl), where Perl has Parrot and Java has the JVM, Aldor has FOAM, the First Order Abstract Machine:

A major part of the tex2html_wrap_inline302 compiler is concerned with producing optimized intermediate code, or Foam code. “Foam” is an acronym for “First Order Abstract Machine.” The abstract machine is first order in the sense that it does not treat its types as values.

Foam is designed to contain only those concepts which can have an efficient realization in both Lisp and C. For example it is not possible to take an address of a variable because that would be inefficient in Lisp (a closure would be created). Nor are dynamic type tests allowed, as that would be inefficient in C. We have been asked how the lack of address arithmetic limits the potential performance of compiled tex2html_wrap_inline302 vs hand-coded C which uses pointers to traverse arrays in inner loops. It is our experience that this is a minor concern on current architectures with optimizing compilers.

Foam is not restricted to the precise intersection of C and Lisp. Some aspects are handled by support libraries. Big integer arithmetic is assumed as part of Foam, and this is provided as a library for C. Also the memory model differs from both C and Lisp in some details: garbage collection is assumed (this is a run time support library in C) and it is possible to make an explicit request to free storage (in Lisp this is ignored).

A Foam program is comprised of a flat sequence of commands. Foam types have various sizes and uses. or example, “Char” is a text character whereas “Byte” is a character sized integer, “DFlo” is a double precision floating point, “Ptr” can point to an array, record, arbitrary sized integer, etc. Reference instructions contain the kind of reference and the position, e.g., “Loc 3” refers to the third local variable of the current function and “RElt 7 x 2” indicates the 2nd field of the record x, using the 7th layout format. Foam operations consist of instructions, such as “If b n,” which indicates that if b is true then proceed to label n, and builtin operations, e.g., “HIntLT a b” is a half-word-integer less-than comparison. The builtin operations are type specific and conversion operations are generally provided. A detailed description of Foam is given elsewhere [26].

The abstract machine does not support asharp types directly and relies on the code generator to produce appropriate calls to create and maintain types. This has the advantage that one can use these calls to add new representations of types to the system. These representations may be written in tex2html_wrap_inline302 itself, or some other language. This is used in order to interface with the Axiom system, and may be extended to other object/type systems, such as CLOS [21] and C++.

Also, since Aldor evolved out of Axiom, where one might want to coerce types like Polynomial Integers to Polynomial Reals on the fly, or reorder the types in a hierarchy, it supports functions and types as first class objects. This makes for the most interesting features of the language.

Then there are several libraries which make Aldor look like a particularly good system for implementing computer algebra algorithms on: Algebra and Sumit, by Manuel Bronstein, support basic CAS facilities and the solution of first order linear differential/difference equations, respectively. There’s also an older BasicMath library from NAG which seems to support a rather respectible subset of CAS facilities, e.g. subresultants. Of course, considering its origins, it’s also probably trivial to use Axiom from Aldor.

Maybe the biggest reason that Aldor isn’t as well known as it could be is that it is not free– binaries are available for non-commercial use, but the source is not, even though Aldor is not being commercially distributed. Anyhow, in this day and age, what profit is there in commercial general purpose languages? Although it would be more convenient to use Aldor, I suspect most researchers/implementors of CASes (the two classes of people I see most likely to be attracted by Aldor’s particular strengths) will be more willing to reinvent several wheels than depend on a potentially capricious licensor. Too bad, but some Axiom developers are gathering a petition together– go sign it!

More CAS stuff

Friday, December 16th, 2005

I just finished grading the last two assignments for this semester: it took me more than 5 hours, but now I’m done until next year. On the other hand, I’m broke, because I was so busy studying for finals and doing research that I didn’t turn in my timesheets, so I’ve missed *3* paychecks!

Now I’m starting to look back into BasicCAS. I printed out the code a few days ago so I can refamiliarize myself with the product of my genius. I’ve also been reading Computing with Real Numbers with the intention of implementing an exact arithmetic background for BasicCAS. Not because I need to– there are a couple of C libraries out there I could call from Python– but because I’d like to learn the concepts behind it. I’ve read a significant portion of the paper, and I still don’t see how the structure they’re building (a system for manipulating digit streams based on linear fractional transformations) will lead to an exact arithmetic system. It’s a cliff-hanger!

I’m going to implement the arithmetic system in pure Python first, and then write it up in C so when BasicCAS is used in C-based Pythons there is the option for increased speed. I also intend to rework the parser to make it more Mathematica compatible.

Other than that, over the break I’m going to be working a lot on the DCA, since I finally have time. In fact, that’s what I’m off to do right now. It looks like the high pass channel may contain nothing but noise, so we may drop it from consideration in the algorithm.

PyPI package basicCAS

Sunday, July 24th, 2005

I spent this morning figuring out how to use Distutils to 1)compile python extensions/non-pure modules, 2) install modules/packages. Then I packaged up basicCAS, version 1.0, and made my first submission to the PyPI library.

Although it is certainly usable— e.g. Didier is using it in a polynomial class he’s writing—, the first version of basicCAS is not all that wonderful. The eventual goal is full compatibility with Mathematica’s text mode parser, but I have a ways to go on that— even after reading the Mathematica Book, and appendices, it isn’t clear to me exactly what’s going on sometimes. So it looks like I’m going to have to play around with Mathematica for a while to get a feel before I start rewriting the parser. The things that will be– need to be– different in the next version of basicCAS are are: 1) precedence and associativity are respected, 2)all forms of character input taken by Mathematica are accepted, 3) the Fullform constructions returned are identical to Mathematica’s, and 4) the error reporting system as much as possible, reflects Mathematica’s. All of these are tall orders.

Fateman online CA lecture notes

Sunday, July 17th, 2005

I ran into this site maybe a year ago, but lost the link. It is a collection of materials from a class on Computer Algebra that Fateman taught at UCB. I looked through the additional material, but not yet the powerpoint presentations that make up the course– judging from the quality of those, the syllabus, and Fateman’s reputation as one of the prominent figures in CA, this material should provide a good foundation in CA.

Math parser update

Wednesday, June 22nd, 2005

I found a massive error in the parsing of factors: as soon as implied multiplication occurred, any explicit multiplication and division (using * or /) caused an error. That was silly, and easily fixed. I also fixed the problem with unary minuses and pluses. I think I’m going to reimplement the whole thing anyhow, to fit much more closely to the Mathematica grammar (which I can’t seem to find anywhere; even MockMMA doesn’t come with an explicit grammar), so that it could be used for the same purpose MockMMA has in LISP: to allow the implementation of a free system which can use code written for Mathematica. I tried looking at the MockMMA code, but Common LISP is alien to me. I haven’t uploaded the update yet; I’m at school right now.

Mathematica-style parser

Monday, June 20th, 2005

I spent this weekend writing a parser for a Mathematica-style grammar. I posted the files (scanner, parser, and grammar). They should be pretty interesting for anyone who is interested in implementing a non-trivial parser from scratch. I based them on the example calculator parser in Programming Python, but the complexity of the grammar I used, compounded with the fact that it doesn’t evaluate anything, makes it different enough to be interesting as is. The only big problem I have is that I haven’t been able to figure out how to implement unary minuses; everything else is licked, but that’s getting me. I’m going to implement the equivalent of Mathematica’s FullForm command for printing out the trees returned by the parser; currently, I’m using a crude tree dump that’s pretty hard to read (and extremely long! for even short expressions), but contains every last detail. Also, I’m going to gut the current error handling system, based on the rather inadequate system in Programming Python, and put in a more informative one; currently, half the time when an expression is ill-formed, the parser crashes instead of signaling the location of the error as it should, and it never makes any suggestions about what the error could be. Another caveat is that, for some reason, when the parser crashes on one expression, it sometime crashes on everything else after it, so you have to create a new parser for each expression to be safe. Oh, and the parser (this is minor, I’m just too tired to fix right now) expects each expression to be on separate lines.

Example of parsing f[{x_,y_}] = {0-1*y, x} / Norm[{x,y}] ( the weird 0-1*y is because of the incapability of dealing with unary minuses).

BasicCAS Mathematica FullForm

Type: factor
Value: None
Children:
   Type: term
   Value: None
   Children:
      Type: ImmediateAssign
      Value: Type: factor
      Value: None
      Children:
         Type: term
         Value: None
         Children:
            Type: list
            Value: None
            Children:
               Type: sum
               Value: None
               Children:
                  Type: factor
                  Value: None
                  Children:
                     Type: term
                     Value: None
                     Children:
                        Type: integer
                        Value: 0
                        Children:

                     *

                  +
                  Type: factor
                  Value: None
                  Children:
                     Type: term
                     Value: None
                     Children:
                        Type: integer
                        Value: 1
                        Children:

                     *
                     Type: term
                     Value: None
                     Children:
                        Type: identifier
                        Value: y
                        Children:

                     *

                  -

               Type: factor
               Value: None
               Children:
                  Type: term
                  Value: None
                  Children:
                     Type: identifier
                     Value: x
                     Children:

                  *

         *
         Type: term
         Value: None
         Children:
            Type: funccall
            Value: None
            Children:
               Norm
               Type: factor
               Value: None
               Children:
                  Type: term
                  Value: None
                  Children:
                     Type: list
                     Value: None
                     Children:
                        Type: factor
                        Value: None
                        Children:
                           Type: term
                           Value: None
                           Children:
                              Type: identifier
                              Value: x
                              Children:

                           *

                        Type: factor
                        Value: None
                        Children:
                           Type: term
                           Value: None
                           Children:
                              Type: identifier
                              Value: y
                              Children:

                           *

                  *

         /

      Children:
         Type: funccall
         Value: None
         Children:
            f
            Type: factor
            Value: None
            Children:
               Type: term
               Value: None
               Children:
                  Type: list
                  Value: None
                  Children:
                     Type: factor
                     Value: None
                     Children:
                        Type: term
                        Value: None
                        Children:
                           Type: identifier
                           Value: x_
                           Children:

                        *

                     Type: factor
                     Value: None
                     Children:
                        Type: term
                        Value: None
                        Children:
                           Type: identifier
                           Value: y_
                           Children:

                        *

               *

   *
List[Times[-1, y, Power[Norm[List[x,y]], -1]], Times[x, Power[Norm[List[x, y]], -1]]]

Note that Mathematica performs some pre-processing before displaying the output (the zero disappears, for one, and the list takes the division inside). I plan on implementing this also, because among other things, I’m going to use GMP or a similar library for the numbers, and I want to keep rationals, reals, and integers distinct.

The next step is to implement a pattern based function definition system, like the one Mathematica uses. I love being able to say f[{x_, y_}] := {-y, x} /Norm[{x, y}] in Mathematica, so I want to be able to define functions in the same way. Or specify the value for a particular class of inputs, e.g. f[x_Integer] := Factorial[x] . Before starting the parser, this seemed like a difficult task, but now it seems like cakewalk. To implement this system, all I have to do is store the expression tree used to define the parameters to the function, and then match that tree against the parameters used to call the function. In the first example, the tree would be something like List[x, y] , so I’d match up the passed in parameters to this tree, binding x and y to the value of the expressions in the corresponding space, and running any tests on them (like _Integer); if the tree fit and the parameters passed the test, then I’d know what definition of the function to use.

After that, I need to make a namespace system. Then it should be a simple task to code in the simple interpretive programming utilities like For, Do, While, etc. Then I can think about the real hard stuff– the symbolic stuff– at leisure.

Comments and suggestions are more than welcome.