Symbolic computing module
Contents
Contents
Abstract
- This project consist in a module that will let work with symbolic math on scilab. There is a module that connect scilab with Maxima through an intermediary program but this way is very slow, a faster way is to create a dynamic library of Maxima with ECL, but a native tool is a better choice.
A lighter way is using GiNac library, the integration with scilab will be easier, and there is more programmers on C++ than Lisp, so more people would enhance the tool. This approach is the one I prefer.
Student time line
See here for the Google Summer of Code planning.
Problems to solve
Ideas
Implemetation
Scope
At this moment the scope has a scope with named variables, and soon will have temporary variables, with a way to clean unused variables.
Types of variables
A symbolic object is a tlist with a sub-type, a reference to the internal symbolic scope. ie: a = tlist(sym","type","ref,"symbol", "a"]; b = tlist(sym","type","ref,"symbol", "b",]; a+b // is like tlist(sym", "type", "ref, "ex", "#temporal#0")
Sub-types are at this moment:
- symbol: Symbolic variable
- ex: Expression based on symbols.
- number: Real, complex and rational numbers.
- constant: Constant predefined values, like pi, e, ...
Integrals
We can define an integral in GiNaC, and work with it without solve the integral, and call to an expression the eval_integ method and that it will try to solve all the integrals inside the expression.
I will define two functions, one to define an integral, and other one to solve integrals.
sym_integral it will define an integral on an interval.
sym_solve_integral it will try to symbolic solve all the integrals on an expression.
Note that evalf it will numerical solve an integral.
'''Symbolic Integration Tutorial''' '''Symbolic integration: the stormy decade'''
Lists
GiNaC lead us to define list with the comma operator, but in scilab there is no way to overload the comma operand, we can use the list datatype, but that will increase data access time, and it will be created a sym_list" function, that create a symbolic list, and sym_list2sci and sym_sci2list.
Matrices
We can represent a matrix in two ways on GiNaC. One is a matrix class that has several functions to create matrix froms vectos, diagonal matrix, unit matrix, symbolic matrx and more, the other ways is using indexed objects and construct a two dimension object. We will use here the first one method, in that way we will define this functions:
- sym_matrix(rows, cols)
- sym_matrix(rows, cols, list)
- sym_matrix_diag(list)
- sym_matrix_sub(matrix, row, col, nrow, ncol)
- sym_matrix_reduced(matrix, row, col)
Limits
'''On computing limits in a Symbolic Manipulation System'''
Solvers
- Lineal solver:
Given a system Ax=B, where A is a Matrix nxn and B a matrix nx1, the solution of this system is given by sym_solver_lineal
Combinatorics
Polynomials
Legendre Polynomial
Well, there is several ways to generate the legendre Polynomial Legendre Polynomial:
- With an contour integral.
- With a sums (R. Schmied).
- Wit a special form of a Jacobi Polynomial.
Chebyshev Polynomial
There is several way to calculate this polynomial Efficient Computation of Chebyshev Polynomials in Computer Algebra
- As the determinant of a special constructed matrix.
- With a generating function expanded on a taylr series.
- With the Rodrigues representation.
- Recurrence equation.
- Differential equation.
- Series Representation.
- Divide and conquer.
Jacobi Polynomial
Laguerre Polynomial
Associated polynomials
Several polynomials has associated polynomials, that are defined in most of cases in terms of the main polynomial, some recurrency law, or other
TODO
- Define the integration between similar datatypes, like lists, polynomials, integers.
- free some pointers!!! (Done)
- numeric to scilab float
- constants
- built-in functions of GiNaC
- relations.
- integrals
- matrices
- Determinant (Done)
- Trace (Done)
- Characteristic polynomial (Done)
- rank (Done)
- inverse (Done)
- solve (Done)
- Indexed elements.
- indexes
- tensors
- limits
TODO Rosetta
This list (checklist) is constructed based on The Rosetta Stone for Computer Algebra Systems (I think that some of them don't apply to the module, but i will give them a try)
- (DONE) substitution (sym_subs)
- (DONE) list (sym_list_new)
- (DONE) matrix (sym_matrix_new)
- equation
- (DONE) list element (sym_list_item_get, sym_list_item_set)
- (DONE) matrix element (sym_matrix_item_get, sym_matrix_item_set)
- (DONE) list length (sym_list_size)
- (DONE) prepend element to list (sym_list_prepend)
- (DONE) append element to list (sym_list_append)
- (DONE) append two list(sym_list_concatenation)
- (DONE) matrix column didmension(sym_matrx_size, sym_matrix_size_rows, sym_matrix_size_cols)
- (DONE) list to column vector(sym_matrix_new(sym_list_size(l),1,l)
- (DONE) column vector to list
- true, false, and, or, not, equal, not equal, if-then-else
- generate list
- complex looping
- function definition
- assignment
- clear vars and fun
- define functions from symbol
- function with undefine number of arguments
- sum lists
- pattern matching
- pi, e, i, +/-infinite, sqrt, power
- (DONE) euler's constant, log, atan, factorial
- (DONE) legendre polynomial, (sym_poly_legendre)
- (DONE) chebyshev polynomial (sym_poly_chebyshev_T, sym_poly_chebyshev_U)
- fibonnaci, elliptic integral
- (DONE) gamma function, psi function
- cos integral, bessel function
- hypergeometric function
- dirac delta
- (DONE) unit step
- define abs(x) as piecewise (checklist??)
- assumptions, remove assumption
- simplification of expressions
- (DONE)unknown functions
- (DONE) numerical evaluation (sym_evalf)
- module operation
- Solve e =_ 0 mod m for x
- Put over common denominator
- expand
- roots of polynomials
- noncommutative multiplication
- solve pair of equations
- Decrease/increase angles in trigonometric functions
- Gröbner basis
- Factorization of e over i
- real part, compex part
- pass from complex to rectangular form
- (DONE) matrix add, mul and transpose
- (DONE) solve Ax=B
- Sum
- Prodc
- Limit
- (DONE) Taylor/Laurent Series
- (DONE) Differentiate
- Integrate
- Laplace transform, Inverse Laplacec transform
- solve ode with initial condition
- Define the differential operator L = Dx + I and apply it to sin x
- 2D plot of two separate curves overlayed
- Simple 3D plotting
GiNaC
- Some features are missing on GiNaC:
Integrals
- There is already a basic::eval_integ, but there is no a definition of integrals of function on the definitions, so, a new function_options will be added to handle integrals functions, in order to define the integrals of built-in functions.
- ginac/function.pl: this perl code generate function.h and funcction.cpp we need to add here the options to definition of integral.
- basic::eval_integ is a virtual method that defined in basic just walk on operands and call its own eval_integ if its redefined on the class of the operand. In integral class is redefined and its return for now just the integral of polynomials, in that part will reside the integrations method.
Well, GiNaC has a really nice architecture, i was reading the source code and thinking in the best way to deal with the integration capability. Basically this is how it work now:
There is an integral class, that act as an algebraic object, is defined with 4 parameteres, the function, the bounds of the integration range, and the integration variable.
This class implements the eval_integ method (it's a virtual method defined on basic class):
* Check if the function is expanded in sums terms if it doesn't then call the expand method, which expand the function in sums terms and call the eval_integ method for each term, and also check for constant factor of the terms.
- * If the function is already expanded checks if its equal to the integration variable and return the integral of it, also check if it's a power on the integration variable and return its integral (just compute polynomial integrals).
There is an adaptivesimpson method when the numeric evaluation is needed.
I need to add then:
- Support to indefinite integrals.
- A wrapper method that calls this indefinite integral, and define this method in basic, in that way all the entities could define its own integration method.
- Define when is computed the integral, i think that it's not really cleaver to try to compute it when the wrapper method is called, because the integral could be used in some computation that reduce its complexity or even make it go away (Fundamental theorem of calculus).
- Define what kind of integration algorithm it will be used (even when this could be a complete GSoC, there are really good implementation (and not very long) of Risch algorithm (or part of it). I will try my best in this part :).adaptivesimpson
What is already done:
There is a eval_integ that is called when we want to compute the integral, it can return an integral object in case that it can't find the integral or an expression if it can, I can use this method in the same way and implement here the integration algorithm.
The integral class could be extended to store indefinite integrals, or create a uintegral new class could be better to split the indefinite and definite integrals.
Define a integrate virtual method on basic class, with two signatures, for defined and indefinite integrals.
Define a integrate_function option to the function class in which is defined the integral of the function defined at the library.
I think that a c++ code could be like this:
symbols x("x"), a("a"), b("b"); ex f = 2*x; cout << f.integrate(x) << endl; // This print a symbolic indefinite integral but it doesn't computer cout << f.integrate(x,1,2) << endl; // This print a symbolic definite integral but it doesn't computer cout << f.integrate(x,a,b) << endl; // This print a symbolic definite integral but it doesn't computer cout << f.integrate(x).eval_integ() << endl; // This compute the symbolic indefinite integral. cout << f.integrate(x,1,2).eval_integ() << endl; // This compute the symbolic definite integral. cout << f.integrate(x,a,b).eval_integ() << endl; // This compute the symbolic definite integral. cout << f.integrate(x).evalf() << endl; // This should return an error. cout << f.integrate(x,1,2).evalf() << endl; // Compute the numerical integration. cout << f.integrate(x,a,b).evalf() << endl; // Compute the symbolic integration.
- 2009/07/09:
I decide not to change the basic or ex classes for sake of simplicity, add a new iintegral class basically copied from integral, and i will take that at a start point to the algorithm.
Right now i can define function like this on GiNaC:
static ex log_iinteg(const ex & x, unsigned iinteg_param) { GINAC_ASSERT(iinteg_param==0); // int {log(x)} -> x * (log (x) - 1) return x * (log (x) - _ex1); } ... REGISTER_FUNCTION(log, eval_func(log_eval). evalf_func(log_evalf). derivative_func(log_deriv). iintegral_func(log_iinteg). series_func(log_series). real_part_func(log_real_part). imag_part_func(log_imag_part). latex_name("\\ln"));
With that now i can compute this:
#include "ginac/ginac.h" using namespace std; using namespace GiNaC; int main() { symbol x("x"); ex f,g,h; f = 2*x; cout << f << endl; g = iintegral(x,f); cout << g << endl; h = g.eval_integ(); cout << h << endl; f = 2*sin(x) + exp(x) + log(x); cout << f << endl; g = iintegral(x,f); cout << g << endl; h = g.eval_integ(); cout << h << endl; }
Which return:
jcardona@terminus:~/stuff/personal/symbolic/spikes/Integrals$ ./test 2*x iintegral(x,2*x) x^2 exp(x)+log(x)+2*sin(x) iintegral(x,exp(x)+log(x)+2*sin(x)) exp(x)-2*cos(x)+(-1+log(x))*x
Limits
This has to be implemented from scratch.
Roots Of Polynomials
Content
Calculus
define_integral
integral_eval
Polynomials
Combinatorics
Usage
(I will try to write a preliminary howto of the module)
Symbols
The basic symbolic datatype is the symbol, it's created with the function sym_symbol("a"), its arguments is a string with the name of the symbolic data to be created:
-->a=sym_symbol("a") a = a -->b=sym_symbol("b") b = b
If you want to define several symbol at once use symbols that receive a variable number of strings arguments and create all the symbols at once, and return all the symbols:
-->[x,y,z] = symbols("x", "y", "z") z = z y = y x = x
We can construct other expression in base of the symbols previously defined:
-->a+b ans = a+b -->sin(a) ans = sin(a)
Numbers
There is three way to define a number:
-->sym_number(1) ans = 1.0 -->sym_number(1.2) ans = 1.1999999999999999556
We can see that the second number has an machine error by the float representation of the machine, because scilab pass to the symbolic number a number that aproch to 1.2, this could be improved by the second way of number definition:
-->sym_number("1") ans = 1 -->sym_number("1.2") ans = 1.2
Passing a string to sym_number we can create more accurate numbers, a third way is defining a fraction:
-->sym_number(1,2) ans = 1/2 -->sym_number(12,10) ans = 6/5
Evaluation of a expression is made by sym_evalf:
-->sym_evalf(sym_number(1,3)) ans = 0.33333333333333333334
And we can define an arbitrary accurate of the floar number with sym_set_digits:
-->sym_set_digits(4) ans = 4. -->sym_evalf(sym_number(1,3)) ans = 0.33333334 -->sym_set_digits(40) ans = 40. -->sym_evalf(sym_number(1,3)) ans = 0.33333333333333333333333333333333333333333333333333333333336 -->sym_set_digits(100) ans = 100. -->sym_evalf(sym_number(1,3)) ans = 0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333334 -->sym_set_digits(200) ans = 200. -->sym_evalf(sym_number(1,3)) ans = 0.333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333 333333333333333333333333333333333333333333333333333333333333335 -->sym_evalf(cos(sym_number(1,3))) ans = 0.944956946314737664388284007675880607845852699565140737677645733750099562196500364824428158805698556592705504904288646216725285458873269654325997662961 78553927751040668266258167782046572205997367707266399322702533
Lists
A symbolic list is created with a sym_list_new and you can add items with sym_list_append sym_list_prepend, get and set items with sym_list_item_get and sym_list_item_set, remove the last and first element with sym_list_remove_first and sym_list_remove_last, and clear with sym_list_clear. ie:
-->[a,b,c,d,e] = symbols("a","b","c","d","e"); -->l = sym_list_new() l = {} -->sym_list_append(l,a) ans = {a} -->sym_list_append(l,b) ans = {b} -->l=sym_list_append(l,a) l = {a} -->l=sym_list_append(l,b) l = {a,b} -->l=sym_list_prepend(l,c) l = {c,a,b} -->l=sym_list_prepend(l,d) l = {d,c,a,b} -->l1 = sym_list_new() l1 = {} -->l=sym_list_prepend(l,l1) l = {{},d,c,a,b} -->l=sym_list_item_set(l,0,e) l = {e,d,c,a,b} -->sym_list_item_get(l,2) ans = c -->l=sym_list_remove_first(l) l = {d,c,a,b} -->l=sym_list_remove_last(l) l = {d,c,a} -->l=sym_list_clear(l) l = {}
=== Matrices ===
-->[x,y,z] = symbols("x","y","z"); -->m = sym_matrix_symbolic("a", 2, 3) m = [[a00,a01,a02],[a10,a11,a12]] -->sym_matrix_transpose(m) ans = [[a00,a10],[a01,a11],[a02,a12]] -->sym_matrix_size_rows(m) ans = 2 -->sym_matrix_size_cols(m) ans = 3 -->l = sym_list_new();l=sym_list_append(l,x);l=sym_list_append(l,y);l=sym_list_append(l,z);l=sym_list_append(l,x); -->l l = {x,y,z,x} -->m1 = sym_matrix_diag(l) m1 = [[x,0,0,0],[0,y,0,0],[0,0,z,0],[0,0,0,x]] -->m1 = sym_matrix_new(2,2,l) m1 = [[x,y],[z,x]] -->m1 = sym_matrix_new(4,1,l) m1 = [[x],[y],[z],[x]] -->m1 = sym_matrix_new(1,4,l) m1 = [[x,y,z,x]] -->m1 = sym_matrix_new(2,2,l) m1 = [[x,y],[z,x]] -->sym_matrix_rank(m1) ans = 2 -->sym_matrix_determinant(m1,x) !--error 999 Error: Incorrect number of parameters. -->sym_matrix_determinant(m1) ans = x^2-y*z -->sym_matrix_charpoly(m1,x) ans = -y*z -->a = sym_matrix_symbolic("a", 2, 2) a = [[a00,a01],[a10,a11]] -->sym_matrix_inverse(a) ans = [[a11*(a11*a00-a01*a10)^(-1),-a01*(a11*a00-a01*a10)^(-1)],[-(a11*a00-a01*a10)^(-1)*a10,a00*(a11*a00-a01*a10)^(-1)]] -->b = sym_matrix_symbolic("b", 2, 1) b = [[b0],[b1]] -->x = sym_matrix_symbolic("x", 2, 1) x = [[x0],[x1]] -->sym_matrix_solve(a,b,x) ans = [[(a11*b0-a01*b1)*(a11*a00-a01*a10)^(-1)],[-(b0*a10-a00*b1)*(a11*a00-a01*a10)^(-1)]]
Polynomials
-->x=symbols("x"); -->T = sym_poly_chebyshev_T(sym_number("10"),x) T = Tn(10,x) -->U = sym_poly_chebyshev_U(sym_number("10"),x) U = Un(10,x) -->sym_expand(T) ans = -1-1280*x^8+50*x^2+1120*x^6-400*x^4+512*x^10 -->sym_expand(U) ans = -1-2304*x^8+60*x^2+1792*x^6-560*x^4+1024*x^10 -->sym_diff(T,x) ans = 10*Un(9,x) -->sym_expand(sym_diff(T,x)) ans = 100*x-1600*x^3+5120*x^9-10240*x^7+6720*x^5
Constants
Expressions
Functions
Derivatives
Derivates are calculated with sym_diff:
-->a = sym_symbol("a") ; b = sym_symbol("b"); -->sym_diff((((a*b)^2-b^a)/(a^b))*(1-a*b^(1-a*b)), a) ans = -(b^a-(a*b)^(2.0))*(a^b)^(-1)*(-1.0+b^(1.0-a*b)*a)*a^(-1)*b+(-(2.0)*a*b^(2.0)+b^a*log(b))*(a^b)^(-1)*(-1.0+b^(1.0-a*b)*a)+(b^(1.0-a*b)-b^(1.0-a*b)*a*log (b)*b)*(b^a-(a*b)^(2.0))*(a^b)^(-1) -->c = (((a*b)^2-b^a)/(a^b))*(1-a*b^(1-a*b)); -->sym_diff(c,a) ans = -(b^a-(a*b)^(2.0))*(a^b)^(-1)*(-1.0+b^(1.0-a*b)*a)*a^(-1)*b+(-(2.0)*a*b^(2.0)+b^a*log(b))*(a^b)^(-1)*(-1.0+b^(1.0-a*b)*a)+(b^(1.0-a*b)-b^(1.0-a*b)*a*log (b)*b)*(b^a-(a*b)^(2.0))*(a^b)^(-1)
Integrals
There are two basic functions (more comming) to handle integrals. The first is to define a integral:
-->x = symbols("x"); -->y = sin(x); -->sym_integral(y,x,0,1) !--error 999 Error: Not a tlist at 3 argument -->sym_integral(y,x,sym_number(0),sym_number(1)) ans = integral(x,0.0,1.0,sin(x))
And the way to solve an integral is usign sym_eval_integral:
-->g = sym_integral(y,x,sym_number(0),sym_number(1)) g = integral(x,0.0,1.0,sin(x)) -->sym_eval_integral(g) ans = integral(x,0.0,1.0,sin(x)) -->y = x + x^2 +x ^3; -->g = sym_integral(y,x,sym_number(0),sym_number(1)) g = integral(x,0.0,1.0,x^(2.0)+x^(3.0)+x) -->sym_eval_integral(g) ans = 1.0833333333333333334
sadly, ginac only solve polynomials integrals yet!, i will add some methods to the library using some paper that i'm collecting from ACM Digital Library.
Command List
- sym_symbol
- sym_number
- sym_
- sym_add
- sym_sub
- sym_mul
- sym_res
- sym_pow
- sym_sin
- sym_cos
- sym_tan
- sym_sinh
- sym_cosh
- sym_tanh
- sym_diff
- sym_evalf
- sym_set_digits
- sym_subs
- sym_scope_clean
- sym_scope_list
- sym_scope_get
- sym_Pi
- sym_Catalan
- sym_Euler
- sym_integral
- sym_integral_eval
- sym_list_new
- sym_list_append
- sym_list_prepend
- sym_list_item_get
- sym_list_item_set
- sym_list_remove_first
- sym_list_remove_last
- sym_list_clear
- sym_list_size
- sym_matrix_new
- sym_matrix_diag
- sym_matrix_unit
- sym_matrix_symbolic
- sym_matrix_sub
- sym_matrix_reduced
- sym_matrix_item_get
- sym_matrix_item_set
- sym_matrix_eval
- sym_matrix_determinant
- sym_matrix_trace
- sym_matrix_charpoly
- sym_matrix_rank
Benchmarks
I will try to put here some benchmark that i made in order to increase the performance of some operation, all the benchmark were tested on my Laptop (Dell 1420 with Intel Core2Duo T7300 2.00Ghz 4MBCache 2GBRam):
Define a symbol
This test measure the time that takes to create a new symbol given that there is already n symbols defined, is a constant time (the resolution of tic-toc is not enough):
Define n symbols
This test measure the time that takes to create n new symbols given that the scope is clean:
Clean scope
This test measure the time that takes to clean all the scope when it has n symbols defined:
Matrix symbolic creation
Basically the creation of a symbolic square matrix should has a complexity of O(n^2).
Concatenation of lists
This test measure the time that takes to concatenate two lists of size n:
Factorial
This test measure the time that takes to calculate a factorial(n):
Chebyshev polynomial
This test measure the time that takes to calculate the expansion of sym_poly_chebyshev_T(n,x):
Legendre polynomial
This test measure the time that takes to calculate the expansion of sym_poly_legendre(n,x):
Jacobi polynomial
This test measure the time that takes to calculate the expansion of sym_poly_jacobi(n,a,b,x) with a and b defined as symbols:
This test measure the time that takes to calculate the expansion of sym_poly_jacobi(n,a,b,x) with a=1 and b defined as symbol:
This test measure the time that takes to calculate the expansion of sym_poly_jacobi(n,a,b,x) with b=1 and a defined as symbol:
This test measure the time that takes to calculate the expansion of sym_poly_jacobi(n,a,b,x) with b=1 and a=1:
Comparisons
Chebyshev Polynomial First Kind
The built-in function on Mathematica give us this:
After increase the performance on scilab creating a new function that return the representation of a expression without store it at the tlist we get this:
Unit test
Original Proposal
Title: Symbolic computing module Student: Jorge Eduardo Cardona Abstract: This project consists in a module that will let work with symbolic math on Scilab. There is a module that connect Scilab with Maxima through an intermediary program but this way is very slow, a faster way is to create a dynamic library of Maxima with ECL, but a native tool is a better choice. A lighter way is using GiNac library, the integration with Scilab will be easier, and there is more programmers on C++ than Lisp, so more people would enhance the tool. This approach is the one I prefer. Content:
Detailed description
There is a module created by Jean-Francois Magni that consist on a intermediary code, that create a connection between scilab and Maxima, through TCP/IP sockets. A better way to connect SciLab with Maxima could be creating a dynamic library (.lib or .dll) with ECL, and then linking this library on the module code.
Axiom is other CAS writed also on Lisp, could be used in the same way that Maxima.
GiNac is a C++ Library that work with symbolic math, maybe this could be the best approach to have a native tool, and have a better performance on it, this tool it has been developed through 10 years, and would be the more easily to maintain in the future.
With GiNac there is a better chance that the code could be maintained by several coders, and the integration just require to define a datatype on scilab, could be a mlist(['sym','name','value'],'symbol',s); and then create wrappers for the operations and the function defined on GiNaC, like expand, derivate, lsolve, evalf, determinant, and more.
GiNaC is Not a Cas, is a library, so, wee have to build a CAS system based on the library, a initial approach is creating basics datatypes that match those of GiNaC, or even more like a sym_real, sym_integer (GiNaC, has only a general Symbol, so we could create on top of it a SymbolReal, SymbolInteger, SymbolRational, SymbolIrrational), but we use the GiNaC differentiation, integration, and more.
As it say on the GiNaC tutorial, it has some disadvantages with other tools like Maxima, this is time of development, in 30 years of development of others tool the debug and enhance process is deeper than the 10 years of GiNaC, but the time that has GSoC, it's just proper to a ligther solution, and in other hand we can see this as something good, at the end the CAS system it will be in absolute control of Scilab, so, it could enhance or add feautures in pro of scilabs's goals.
We have then several approach that will be studied:
- A library with ECL of Maxima.
- A library with ECL of Axiom.
- A native tool from scratch.
- A native tool using GiNaC (I preffer this one).
- A bunch of .sce files with wrappers of a 'sym' with mlists.
Action Plan
- Decide which way will be used, based on what i have searched the best way is with GiNaC, and i will assume this for the next plan.
- Define the features that a will be enable on scilab, like which datatypes, function and operators, and look which match with GiNaC directly. If we want something like a Symbolic_Integer or Symbolic_Irrational, that doesn't exists on GiNaC, but exists a General Symbol, we create then a new class that inheritance from Symbol with this capabilities and overload the necessary operators on c++.
- Define unit tests.
Create like a dictionary (could be a simple std::map<std::string, ginac::symbol>) to storage the symbolic variables, this could be avoided if we add a new data type to scilab, at the code, and used the way in which scilab use this variables. At this moment i'm using a tlist to create 'sym' objects, and then overload operators, but if i pass this variables to some module functions, the code has to recreate the GiNaC variable and that add extra time, we should be able to store the variables on some data structure, and when a function is called ask for that variable.
- With the mechanism to save and recover the symbolic variables , either with a buffer, o recreating at every function call, we build the basic wrappers from the c code at sci_gateway to the cpp code with GiNaC basic functions.
- Add features that doesn't exists directly on GiNaC.
Preliminary work
I'm already work on this, and following the topic at https://scilab.gitlab.io/legacy_wiki/Scilab_Module_Architecture I build a initial 'sym' module in which i'm learning to handle with the passing of arguments and jumping to a c++ code. Looking the way that Scilab stored variables, i think the best way to stored the symbolic data, that belong to GiNaC could be a static map living on that c++ code, with some basic functions like sym_register, sym_ask. In this way when you want to find the derivative you write:
x=sym_register("x"); //the x returned is a tlist, with 'sym' as type y=sym_register("y"); sym_check("x"); y=sym_diff(x+y^2-2*x,x); // As it return a tlist, this will print a not very pretty representation. sym_print(y); // this should print a 1.
With the sym_register("x") we call a function that return a tlist, and insert a pair<string, symbol> at the static map, so it will be ready for a future use.
The x+y^2-2*x create a new tlist, that act as a tree with some overloading like this (it could be implemented at the c++ code, and ):
function x = sym(s)
- x = tlist(['sym','name','value'],'symbol',s);
endfunction
function x = sym_print(s)
- if ((s.name == "+")|(s.name == "-")|(s.name == "*")|(s.name == "/")|(s.name == "^"))
- x = "(" + sym_walker(s.value(1)) + s.name + sym_walker(s.value(2)) + ")";
- x = string(s.value);
- x = "wtf(" + s.name + ")";
endfunction
// Basic arithmetic operations deff("x=%sym_a_sym(a,b)",'x = tlist([sym,name,value],+,list(a,b));'); deff("x=%sym_m_sym(a,b)",'x = tlist([sym,name,value],*,list(a,b));'); deff("x=%sym_s_sym(a,b)",'x = tlist([sym,name,value],-,list(a,b));'); deff("x=%sym_r_sym(a,b)",'x = tlist([sym,name,value],/,list(a,b));'); deff("x=%sym_p_sym(a,b)",'x = tlist([sym,name,value],^,list(a,b));');
// Operations with constant deff("x=%s_a_sym(a,b)",'x = tlist([sym,name,value],constant,a) + b;'); deff("x=%s_m_sym(a,b)",'x = tlist([sym,name,value],constant,a) * b;'); deff("x=%s_s_sym(a,b)",'x = tlist([sym,name,value],constant,a) - b;'); deff("x=%s_r_sym(a,b)",'x = tlist([sym,name,value],constant,a) / b;'); deff("x=%s_p_sym(a,b)",'x = tlist([sym,name,value],constant,a) ^ b;');
deff("x=%sym_a_s(a,b)",'x = a + tlist([sym,name,value],constant,b);'); deff("x=%sym_m_s(a,b)",'x = a * tlist([sym,name,value],constant,b);'); deff("x=%sym_s_s(a,b)",'x = a - tlist([sym,name,value],constant,b);'); deff("x=%sym_r_s(a,b)",'x = a / tlist([sym,name,value],constant,b);'); deff("x=%sym_p_s(a,b)",'x = a ^ tlist([sym,name,value],constant,b);');
So, we pass this "tree" in a tlist to the sym_diff.
At this time the module already has the sym_registry and sym_check functions, the first one receive a string and then call a sym_register function that create a new symbol and put it an a map<string, symbol>, and return a tlist with type sym. The second receive a string and ask to the map if it has a symbol with that name, if it does return 1 if not return 0; I can build then the symbolic scope easily with this map.
With sym_diff we call a function that take two tlist as arguments, the second one should be a simple sym (not a tree), and capture the first tlist, and start to replace while gets down at the tree and changing symbolic string with symbols at the scope map, and at the end it has an GiNaC expression.
Look at the map handler:
#include "sym_register.hxx" static map<string, symbol> scope; int sym_register(const char* name) { pair<map<string,symbol>::iterator,bool> ret; ret = scope.insert(pair<string, symbol>(string(name),symbol(name))); if (ret.second) return 1; else return 0; } int sym_check(const char* name) { if (scope.find(name) == scope.end()) return 0; return 1; } ****** int sci_sym_register(char *fname,unsigned long fname_len) { int m, n, sp; int numbers[] = {0,1,2,3}; char **name = NULL; char *sym_name = NULL; char *sym_symbol_name = NULL; int number_of_elements; int actual_postion; CheckRhs(1,1); CheckLhs(0,1); GetRhsVar(1, MATRIX_OF_STRING_DATATYPE,&m,&n,&name); /* Register symbol at the scope */ sym_register(name[0]); sym_name = strdup("sym"); sym_symbol_name = strdup("symbol"); actual_postion = Rhs+1; CreateVar(actual_postion,TYPED_LIST_DATATYPE, &numbers[3], &numbers[1], &sp ); CreateListVarFromPtr(actual_postion, 1,MATRIX_OF_STRING_DATATYPE, &numbers[1], &numbers[1], &sym_name); CreateListVarFromPtr(actual_postion, 2,MATRIX_OF_STRING_DATATYPE, &numbers[1], &numbers[1], &sym_symbol_name); CreateListVarFromPtr(actual_postion, 3,MATRIX_OF_STRING_DATATYPE, &numbers[1], &numbers[1], name); LhsVar(1) = Rhs+1; C2F(putlhsvar)(); // Free strings!! return TRUE; } ****** int sci_sym_check(char *fname,unsigned long fname_len) { int n, m; char **name = NULL; int actual_position; int one=1; int *ret = NULL; ret = (int*) MALLOC(sizeof(int)*1); CheckRhs(1,1); CheckLhs(0,1); GetRhsVar(1, MATRIX_OF_STRING_DATATYPE,&m,&n,&name); ret[0] = sym_check(name[0]); actual_position = Rhs+1; CreateVarFromPtr(actual_position,MATRIX_OF_INTEGER_DATATYPE, &one, &one, &ret ); LhsVar(1) = actual_position; C2F(putlhsvar)(); return TRUE; }
Better ideas will come with more time, but this spike show that it's no hard to implement the complete system, and the GSoC period will be right for it, even when i'm not interact directly with GiNaC, is just to storage expressions, even that could be the only scope, and share the same map, with, expressions, symbols, matrix, numerics.
Well, having all the interaction part defined and correcty implemented is just a matter of use the features of GiNaC, and also enhnace some of them.
Future
- I'm a full-time user of free software, and an engineering student, I will deal with software like Scilab for a long time, I can use part of my time to maintain the tool on the future. And in my future plans I will be able to have the time for it.
About you
* I'm a student of electronics engineering at Universidad Pontificia Bolivariana, with experience on tools like scilab, scicos, matlab, simulink, labview, openmodelica, and several python libraries for simulation, numerical math.
* I have experience in programming with C and C++, like 5 years now.
* I study but i'm also work on a company half-time (until the next week, i need more time to do my undergraduate thesis) and I'm been involved in a couple of software development project, using XP methodology.
* In my engineer career i received great math and physics background, and personally i really enjoy all the modeling and simulation theory (I actually expend a lot of money on amazon buying books of theory of modeling and simulation of discrete systems )
* I knew about scilab like 4 years ago, but i'm start using it like 2 and half, and for the last year i used a lot because my undergraduate work project.
* I'm really interested on this project, i found it like a opportunity also to maybe publish a paper about the construction of a CAS and how it works for a Mathematical meeting at August at Cali.
Original Comments
I end a prototype, that store all expressions on the scope, and overload basic arithmetic operations, and create a sym_diff to do this:
-->b=sym_symbol("b"); -->a=sym_symbol("a"); -->sym_diff(a,b); -->c=sym_diff(a,b); -->c.repr ans = 0 -->c=sym_diff(a,a); -->c.repr ans = 1 -->c=sym_diff(a^2,a); -->c.repr ans = (2.0)*a -->c=sym_diff(a^2+b^a-2*(a-b)/(b^(2-a)),a); -->c.repr ans = log(b)*b^a-log(b)*(b^(2.0-a))^(-1)*(-(2.0)*b+(2.0)*a)-(2.0)*(b^(2.0-a))^(-1)+(2.0)*a -->
I'm storing all expression at the scope, even intermediary expressions, because i can find a way to identify when someone type just a+b or c=a+b, and at the first case, the expression exists only for a while, but the other case is a more long-lived value, and even if i delete the c variable of the scilab scope, i should do the same with the expression at the symbolic scope, i was thinking on a clean_symbolic_scope function to do that. But it's a good start.