PySymbolic design proposal

by Pearu Peterson <pearu@cens.ioc.ee>

$Revision: 1.4 $

$Date: 2000/10/24 15:20:06 $

• PySymbolic home page is http://cens.ioc.ee/projects/pysymbolic/
• Mailing list is pysymbolic-users@cens.ioc.ee

This document is OBSOLUTE. Read PySymbolic Users Guide instead.

Contents

• Introduction
• Getting started
• Developers guide
  • Overview of PySymbolic
  • Parsing Python strings
  • Manipulating ast-lists
  • Operations with Symbolic objects
  • Extensions
    • Differentiation
    • Simplification
    • Wish list for extensions
• Links

Introduction

PySymbolic - Doing Symbolics in Python - is a Python package for symbolic mathematical calculation. One of the design features of PySymbolic is that mathematical notations can be represented using Python native syntax. Although Python was not meant to manipulate its objects symbolically, its clear and rich language design makes it possible to extend Python to the world of symbolic calculations of mathematics. However, it is clear that Pythons fixed set of operators will limit PySymbolics ability to provide the level of completeness of mathematical notations as it is in programs like Mathematica, Maple, etc. But this rather pessimistic opinion of mine may turn to be wrong (I certainly hope so!). It is also clear that in short term PySymbolic cannot compete with available symbolic manipulation tools except in its low cost (0 cents) and its availability in Python programs (this was my main point for starting this project).

A very important factor for this project to be successful, is the support from Python community. I hope that there will be enough interested users and contributors to reach the critical mass and to get Python to crunch symbols (TM).

Getting started

To use PySymbolic, you must import Symbolic module:

>>> from Symbolic import *

Module Symbolic provides class Symbolic and various functions such as simplify, diff, etc. To manipulate objects symbolically, you must define them as Symbolic instances. For example,

>>> x = Symbolic("x")

This Symbolic instance x can be used to generate new and more complex Symbolic objects by just performing operations with it. For example,

>>> y = 5*x*x + 4*x/3
>>> y
Symbolic("(5 * x) * x + (4 * x) / 3")

In general, the constructor Symbolic(<obj>) converts arbitrary Python object to a Symbolic instance. Typically, a Symbolic object is created from a string. For example,

>>> z = Symbolic("(3*y-2)*(2*x-x/2)")
>>> print z
(3 * x - 2) * (2 * x - x / 2)

Note that a Symbolic object may contain arbitrary (but valid) Python statements. For example,

>>> c = Symbolic('if x==2:\n\ty=3\nelse:\n\tz=1\n')
>>> print c
if x == 2: 
    y = 3
else: 
    z = 1

An application for the full support of Python language syntax may be symbolic differentiation of functions defined in Python, for example. However, in the following I will concentrate more on simple Python statements that make more sense for symbolic mathematical calculations. No need to mention that Python environment itself provides a rich variety of tools for programming tasks.

As already seen above, you can perform various operations with Symbolic objects. In expressions typical operations are addition, subtraction, multiplication, division, raising to a power, etc. For example,

>>> 1*x/y+2*z+0
Symbolic("(x / ((5 * x) * x + (4 * x) / 3) + 2 * ((3 * y - 2) * (2 * x - x / 2)))")

While constructing the new Symbolic object for 1*x/y+2*z+0, only most obvious simplifications are carried out. To get a deeper simplification level, use simplify(<obj>) function. For example,

>>> f = simplify(1*x/y+2*z+0)
Symbolic("-6 * x + x / (4 * x / 3 + 5 * x**2) + 9 * x * y")

(Clearly, there is room for improvement, but that's the current behavior of vers. 0.14.)

Warning: Note that 3/4*x is 0 because integer division. You must use 3*x/4.

In PySymbolic you can differentiate expressions using diff(<expr>,<var>) function. For example,

>>> print simplify(diff(x**3/3-x**2+5*x,x))
5 - 2 * x + x**2

Developers guide

In this section the internals of PySymbolic package are described. I hope it will help possible contributors to get better understanding of the code. Current implementation of PySymbolic is certainly in alpha stage and it may change considerably in future.

Overview of PySymbolic

In PySymbolic one can distinguish three closely related tasks:

1. Creating Symbolic objects;
2. Performing operations with Symbolic objects;
3. Applying extension functions.

New Symbolic objects are created either directly with the constructor Symbolic(<string>) or as a result of operations with already defined Symbolic objects. In the former case, Python code in <string> is parsed and the result (ast-list, see below) is saved into its attribute ast. All manipulations with Symbolic objects are actually manipulations with their ast-lists. The so-called extension functions are functions that actually do relevant manipulations with Symbolic objects, e.g. simplify, differentiate, integrate, etc.

Parsing Python strings

Parsing of a Python string containing arbitrary Python code means constructing a lisp-like list - I will call it ast-list. An ast-list is a recursive list of lists where 1st element (I will call it ast-list key) is always a string and it defines how to interpret the following elements which can be either ast-lists or a single string. For example, Python string a+b is parsed to an ast-list ['arith_expr',['NAME','a'],['PLUS','+'],['NAME','b']].

Relevant functions are

• symparser.py2ast(<string>) - returns ast-list corresponding to <string>.
• symparser.ast2py(<astlist>) - returns Python string corresponding to <astlist>.

Complete specification of ast-lists is determined by Python grammar. Python grammar is defined in the file Grammar/Grammar located in the Python source tree.

PySymbolic parsing part (symparser.py) can be considered being complete. For better performance, it should eventually implemented in C.

Manipulating ast-lists

To perform symbolic operations in the level of ast-lists, use the following functions:

• auxfuncs.astmul(a,b) - return the product of ast-lists a and b. This corresponds to a*b.
• auxfuncs.astdiv(a,b) - return ast-list corresponding to a/b.
• auxfuncs.astadd(a,b) - return ast-list corresponding to a+b.
• auxfuncs.astsub(a,b) - return ast-list corresponding to a-b.
• auxfuncs.astpow(a,b) - return ast-list corresponding to a**b.
• auxfuncs.astneg(a,b) - return ast-list corresponding to -a.
• auxfuncs.makefunc(f,l) - return ast-list corresponding to f(l1,..,lk). Here f is an ast-list representing function name and l is a k-list of ast-lists representing functions arguments.
• (This list of functions can be expanded in future).
• The functions above perform elementary simplifications that include
  • 1*a and a*1are returned as a, 0*a and a*0 as 0;
  • a/1 as a, 0/b as 0, a/a as 1;
  • 0+a and a+0 as a, n+m as k=n+m if n,m are numbers;
  • 0-a as -a, a-0 as a, n-m as k=n-m if n,m are numbers;
  • a**0 and 1**a as 1, a**1 as a;
  • -(a+b) as -a-b, -(-a) as a, -(a) as (-a), -(n*a) as -n*a, -(-1*a) as a, -(a*b) as -1*a*b, -(a**b) as -a**b;
  • (This list of simplifications can be extended or reduced - any suggestions? E.g. how to interpret 0/0, a/0?).
  • Need simplification for some special function arguments: ln(e), ln(1), exp(0), sin(n*2*pi), etc. Or should this be in simplification extension?

To create elementary ast-lists, use functions

• auxfuncs.astNAME(s)- return ast-list corresponding to symbol s. s is string.
• auxfuncs.astNUM(n) - return ast-list corresponding to number n. n is number or string (`n`).
• auxfuncs.astSTR(s) - return ast-list corresponding to string s. s is string.

and the dictionary

• auxfuncs.astm[k] - is ast-list corresponding to symbol k. Here k is a string from a list ['~', '+=', '|', '}', '{', '-=', '/=', '^', ']', '[', '**', '<=', '<>', '&=', '|=', '==', '**=', '*=', '>', '<', '=', ':', ';', '>>=', '.', '/', ',', '-', '*', '+', '(', ')', '&', '<<=', '^=', '%', '%=', '<<', '>>', '>='].

The following test functions on ast-lists are available:

• auxfuncs.iszero(a) - returns 1 if ast-list a corresponds to 0;
• auxfuncs.isone(a) - if a corresponds to 1;
• auxfuncs.isratterm(a) - if a corresponds to rational number n/m;
• auxfuncs.isatom(a) - if a corresponds to either symbol, number, string, or atom. Atom is either list, tuple, dictionary, or `` statement;
• auxfuncs.isnegative(a) - if a corresponds to expression starting with `-' sign;
• auxfuncs.ispower(a) - if a corresponds to power c**p;
• auxfuncs.isinverse(a) - if a corresponds to power c**-1;
• auxfuncs.isintpower(a) - if a corresponds to integer power c**n;
• auxfuncs.isratpower(a) - if a corresponds to rational power c**(n/m);
• auxfuncs.isfunc(a) - if a corresponds to function f(...);
• auxfuncs.istimes(a) - if a corresponds to term n*c;
• auxfuncs.israttimes(a) - if a corresponds to term n/m*c;
• (Should these functions be renamed to astis...?)

Some auxiliary functions that can be useful:

• auxfuncs.gcd(n,m) - returns greatest common divisor of integers n and m;
• auxfuncs.factor_int(n) - returns factors and the corresponding powers of an integer n. For example. factor_int(200) returns [2,5],[3,2] because 200==2**3*5**2;
• auxfuncs.addcoefs(n,m,k,l) - returns the terms of collected rational sum of n/m and k/l. Also performs normalization of returned rational number, for example, addcoefs(2,3,5,6) returns 3,2 because 2/3+5/6==3/2;
• (for performance, these functions should be implemented in C or be used from some specialized number theory C library);
• auxfuncs.atom_ize(a) - if a is not an atom then returns (a), otherwise a;
• auxfuncs.join_ize(l,j) - returns a list of ast-lists [l1,j,l2,j,...,j,lk] where l is k-list of ast-lists and j==astm[s] for some s. This function is useful for constructing composed ast-lists, e.g. ['artih_expr']+join_ize([['NAME','a'],['NAME','b']],'+') corresponds to a+b;
• auxfuncs.factor_ize(a,s) - returns ast-list corresponding to '%sa'%s, where s is in ['+','-','~'];
• auxfuncs.getpower(a) - returns ast-lists n,c such that a==c**n;
• auxfuncs.getcoefs(a) - returns ast-lists n,m,c such that a==n/m*c;

Operations with Symbolic objects

Symbolic object has the following methods:

• __init__(obj=None,ast=None) - the constructor. Argument ast is used for constructing Symbolic object from a raw ast-list if obj==None;
• __str__(), __repr__() - return proper Python strings;
• __call__(*args,**kws) - return new Symbolic object representing function self(<args>,<keywords>);
• __pos__() - return +self;
• __neg__() - return -self;
• __add__(obj) - return self+obj, where obj is first transformed to Symbolic object;
• __sub__(obj) - return self-obj;
• __radd__(obj) - return obj+self;
• __rsub__(obj) - return obj-self;
• __mul__(obj) - return self*obj;
• __div__(obj) - return self/obj;
• __rmul__(obj) - return obj*self;
• __rdiv__(obj) - return obj/self;
• __pow__(obj) - return self**obj;
• __rpow__(obj) - return obj**self;
• issimpleexpr() - returns 1 if self is simple expression (needs revision);

In addition, Symbolic has the following extension methods:

• simplify() - return self after simplifications;
• diff(obj) - return derivative of self in respect to obj which must be ast-list with a key 'NAME';

Extensions

Currently there are only two extensions that are partially implemented: simplification and differentiation. I am not completely happy with the current simplification procedures as its many tasks should actually be implemented elsewhere, e.g. in functions like expand, collect, etc.

Differentiation module, however, needs no rewrite, just extension and bug fixes.

For better performance all extension functions should be eventually implemented in ast-list level, or better yet, in C.

Current experience shows that extension functions are, in general, recursive.

Differentiation

Differentiation routines locate in file differentiator.py and it provides a function:

astdifferentiate(a,v)

It returns an ast-list containing derivative of a in respect to v. (a,v are ast-lists, of course.) It behaves as follows (let's denote diff=astdifferentiate):

• diff(s,x) - returns 0, where s is number or a symbol not equal to x;
• diff(x,x) - returns 1;
• diff((a),x) - returns (diff(a,x));
• diff(-a,x) - returns -diff(a,x);
• diff(a+b,x) - returns diff(a,x)+diff(b,x);
• diff(a*b,x) - returns diff(a,x)*b+a*diff(b,x);
• diff(a/s,x) - returns diff(a,x)/s;
• diff(a/b,x) - returns (diff(a)*b-a*diff(b))/b**2;
• diff(a**s,x) - returns s*a**(s-1)*diff(a,x);
• diff(s**a,x) - returns s**a*ln(s)*diff(a,x);
• diff(a**b,x) - returns a**b*(diff(b,x)*ln(a)-diff(a,x)*b/a)/(b*ln(a))**2;
• diff(f(),x) - returns 0;
• diff(f(g(x)),x) - returns f_0(g(x))*g_0(x)(or should function derivatives denoted in some other way, e.g. D(f,0)(g(x))*D(g,0)(x)?);
• diff(f(g1(x),...,gk(x)),x) - returns f_0(g1(x),...)*g1_0(x)+...+f_km1(g1(x),...)*g_k(x);
• diff((a1,...,ak),x) - returns (diff(a1,x),...,diff(ak,x));
• Differentiation of elementary functions are described in the dictionary math_function_derivatives:
  • diff(ln(a),x) and diff(log(a),x)- return 1/a*diff(a,x);
  • diff(exp(a),x) - returns exp(a)*diff(a,x);
  • diff(log10(a),x) - returns ln(10)/a*diff(a,x);
  • diff(sin(a),x) - returns cos(a)*diff(a,x);
  • diff(cos(a),x) - returns -sin(a)*diff(a,x);
  • diff(tan(a),x) - returns 1/cos(a)**2*diff(a,x);
  • diff(cot(a),x) - returns -1/sin(a)**2*diff(a,x);
  • diff(sinh(a),x) - returns cosh(a)*diff(a,x);
  • diff(cosh(a),x) - returns sinh(a)*diff(a,x);
  • diff(tanh(a),x) - returns 1/cosh(a)**2*diff(a,x);
  • diff(coth(a),x) - returns -1/sinh(a)**2*diff(a,x);
  • diff(sqrt(a),x) - returns 1/(2*sqrt(a))*diff(a,x);
  • (other?)
• in all other cases an exception is raised.

Simplification

Simple is always a very relative concept and so there should be several simplification functions that simplify expressions in different ways, often in very controversial ways. For example, in some cases it is useful to expand powers so that lots of terms will cancel, in other cases, however, it gives little effect, expressions may become even more complicated.

On the other hand, simplification functions should try to make expressions as simple as possible without users intervention. It probably requires different strategies in different parts of an expression. It is not clear to me, how to make this in the most optimal way. Should it try all possible combinations at all recurrence levels and then tries to find the most compact combination? What is the measure for simplicity? How to optimize this? Caching?

Ok, to tackle these problems, some basic simplification functions are required. I would propose the following functions:

• expand(expr,*opts) - performs expansion that includes opening parenthesis, optionally expanding integer powers if power is less than, say, 10;
• collect(expr,*opts) - collects terms (what else?);
• factor(expr,*opts) - should this be in polynomial extension?
• assume(...) - a mechanism to say that x is positive, for example;
• (what else?)
• simplify(expr,*opts) - tries hard to simplify an expression, possibly using combinations of the functions above;

Wish list for extensions

For completeness, PySymbolic should contain the following extensions (in random order):

• Factorization of polynomials, both over finite and infinite fields;
• Limits;
• Integration;
• Special functions;
• Sums;
• Series;
• Sets;
• Combinatorics, Discrete mathematics;
• Group theory;
• Number theory;
• Root finding, solving equations;
• Linear Algebra;
• Input/output tools (TeX, MathML,C,Fortran,etc);
• (what else?)

Many of these disciplines are not simple and working with them requires special algorithms with different complexity levels. Some of them can be implemented directly in Python, and for some of them interfaces to freely available C libraries should be created.

Anyway, if you feel that you wish to contribute in any of these disciplines, please, let us know. This is to avoid parallel contributions.

Thanks, Pearu

Links

The following libraries should be considered interfaced to Python in order to use them in PySymbolic:

• NTL - a library for doing Number Theoryit contains very fast algorithm for factorization of polynomial;
• GiNaC- "it does not try to provide extensive algebraic capabilities and a simple programming language but instead accepts a given language (C++) and extends it by a set of algebraic capabilities";
• CASs - quite excensive list of CASs;
• Pari-GP- a software package for computer-aided number theory;
• LiDIA- A C++ Library For Computational Number Theory;
• YACAS- is an open source general purpose computer algebra system with a pleasant syntax;
• CoCoA- is a special-purpose system for doing Computations in Commutative Algebra;
• SIMATH- is a computer algebra system, especially for number theoretic purpose;
• SACLIB- is a library of C programs for computer algebra;
• Maxima- is a large computer algebra system for symbolic and numerical computations;
• Macaulay 2- is a new software system devoted to supporting research in algebraic geometry and commutative algebra;
• GAP- Groups, Algorithms and Programming.
• General Purpose Computer Algebra Systems - an overview;
• KANT/KASH- is a software package for mathematicians interested in algebraic number theory;
• GROEBNER- is a system for computing Groebner bases, is free, is C library
• (what else? comments?)

Other attempts to do symbolics in Python:

• Pythonica - a system for doing mathematics with Python;
• Pysimplex- provides some basic symbolic programming tools for constructing, solving and optimizing systems of linear equations and inequalities.