[GiNaC-devel] proposed integration algorithm

[Jose Antonio Garcia Peiro]

Hi

I am very interested in GiNaC, but I also need symbolic integration,
and GiNaC only compute basic integrals. I wrote a simple recursive
routine that applies a basic integration rules, matching the input
function.
It is a good way to implement integration algoritm?

Thanks


int(c,x) = c*x
int(x,x) = x^2/2
int(f+g,x) = int(f,x) + int(g,x)
int(c*f,x) = c*int(f,x)
int( f^c*1,x) = f^(c+1)/(c+1) or log(f)
int( f^c*f',x) = f^(c+1)/(c+1) or log(f)

some examples:

cout << "int(1+x,x) = " << myinteg( 1 + x, x ) << endl;
                 // 1/2*x^2+x
cout << "int(3*x^2+x+6,x) = " << myinteg( 3*pow(x,2)+x+6, x ) << endl;
       // 1/2*x^2+x^3+6*x
cout << "int((3*x+5)^2*3,x) = " << myinteg( pow(3*x+5,2)*3, x ) <<
endl;       // 1/3*(5+3*x)^3
cout << "int((3*x+5)^(-1)*3,x) = " << myinteg( pow(3*x+5,-1)*3, x ) <<
endl;   // log(5+3*x)
cout << "int((3*x+5)^b*3,x) = " << myinteg( pow(3*x+5,b)*3, x ) <<
endl;       // (1+b)^(-1)*(5+3*x)^(1+b)


the code (with debug messages):
ex myinteg( const ex &f, const symbol &x ){
	
	exmap repls;

	cout << endl << "f(x) = " << f << endl;
	//int(c,x)		c*x
	if( !has( f, x ) ){
		cout << "match int(c,x)" << endl;
		return f*x;
	}

	//int(x,x)		x^2/2
	if( f.is_equal(x) ){
		cout << "match int(x,x)" << endl;	
		return pow(x,2)/2;
	}	

	//int(f+g,x)	int(f,x)+int(g,x)
	if( f.match( wild(0)+wild(1), repls ) ){
		cout << "match int(f+g,x)" << endl;
		cout << "wild(0)+wild(1) = " << repls << endl;
		ex f, g;
		for( exmap::const_iterator i=repls.begin(); i!=repls.end(); ++i){
			if(i->first.is_equal(wild(0)) ){
				f = i->second;
			}else if(i->first.is_equal(wild(1))){
				g = i->second;
			}
		}
		cout << "f = " << f << endl;
		cout << "g = " << g << endl;
		return myinteg( f, x ) + myinteg( g, x );
	}

	//int(c*f,x)		c*int(f,x)
	if( f.match( wild(0)*wild(1), repls ) ){
		cout << "match int(c*f,x)" << endl;
		cout << "wild(0)*wild(1) = " << repls << endl;
		ex c, f;
		for( exmap::const_iterator i=repls.begin(); i!=repls.end(); ++i){
			if(i->first.is_equal(wild(0)) ){
				c = i->second;
			}else if(i->first.is_equal(wild(1))){
				f = i->second;
			}
		}
		if( !has( c, x ) && has( f, x )){
			return c*myinteg( f, x );
		}if( has( c, x ) && !has( f, x )){
			return f*myinteg( c, x );
		}
	}


	//int( pow( f, c ), x )				pow( f, c+1 ) * pow( c+1, -1 ) if c != -1,
 or ln( f ) if c == -1
	if( f.match( pow(wild(0),wild(1)), repls ) ){
		cout << "match int( pow( f, c ), x )" << endl;
		cout << "pow(wild(0),wild(1)) = " << repls << endl;

		ex f, c;
		for( exmap::const_iterator i=repls.begin(); i!=repls.end(); ++i){
			if(i->first.is_equal(wild(0)) ){
				f = i->second;
			}else if(i->first.is_equal(wild(1))){
				c = i->second;
			}
		}
		cout << "f = " << f << endl;
		cout << "c = " << c << endl;
	
		if( has( f, x ) && !has( c, x ) ){
			ex dft = f.diff(x,1);
			ex df = 1;
			cout << "df/dft = " << df/dft << endl;
			if( !has(df/dft,x) ){
				cout << "!has(df/dft,x)" << endl;
		 		if( c.is_equal(-1) ){
					return log( f ) * (df/dft);
				}else{
					cout << "df = " << df << endl;
					cout << "dft = " << dft << endl;
					cout << "(df/dft) = " << (df/dft) << endl;
					cout << "(1/dft) = " << (1/dft) << endl;
					cout << "pow(dft,-1) = " << pow(dft,-1) << endl;
					cout << "f = " << f << endl;
					cout << "c = " << c << endl;
					cout << "(pow( f, c+1 )/( c+1 )) = " << (pow( f, c+1 )/( c+1 )) << endl;
					cout << "pow( f, c+1 )*pow( c+1, -1) = " << pow( f, c+1 )*pow(
c+1, -1) << endl;
					cout << "(pow( f, c+1 )/( c+1 ))*(df/dft) = " << (pow( f, c+1 )/(
c+1 ))*(df/dft) << endl;
					cout << "pow( f, c+1 )*pow( c+1, -1) = " << pow( f, c+1 )*pow(
c+1, -1)*df*pow( dft, -1) << endl;
					return (pow( f, c+1 )/( c+1 )) * (df/dft);
				}
			}
		}

	}
	

	
	//int( pow( f, c ) * f', x )				pow( f, c+1 ) * pow( c+1, -1 ) if c
!= -1,  or ln( f ) if c == -1
	if( f.match( wild(0)*pow(wild(1),wild(2)), repls ) ){
		cout << "match int( pow( f, c ) * f', x )" << endl;
		cout << "wild(0)*pow(wild(1),wild(2)) = " << repls << endl;
	 	ex df, f, c;
		for( exmap::const_iterator i=repls.begin(); i!=repls.end(); ++i){
			if(i->first.is_equal(wild(0)) ){
				df = i->second;
			}else if(i->first.is_equal(wild(1))){
				f = i->second;
			}else if(i->first.is_equal(wild(2))){
				c = i->second;
			}
		}

		cout << "f' = " << df << endl;
		cout << "f = " << f << endl;
		cout << "c = " << c << endl;
	
		if( has( f, x ) && !has( c, x ) ){
			ex dft = f.diff(x,1);
			cout << "df/dft = " << df/dft << endl;
			if( !has(df/dft,x) ){
				cout << "!has(df/dft,x)" << endl;
		 		if( c.is_equal(-1) ){
					return log( f ) * (df/dft);
				}else{
					return pow( f, c+1 )*pow( c+1, -1) * (df/dft);
				}
			}
		}
	}
	

	return 0;
}