/*****************************************************************										        **		Add on some more complex functions.  These are mostly used to find	**	different constants to be used in various expansions of trig and exp		**	functions.  Need to combine with polynomials of floats to create all terms	**	of a full expansion using Chebyshev polynomials and reducing to standard	**	polynomials.									**					Author = Mike Rosing				**					 date   = May 1, 2000				**											*****************************************************************/#include "bigfloat.h"#include "multipoly.h"void split( FLOAT *x, FLOAT *intprt, FLOAT *frac);extern RAMDATA ram_block[];MULTIPOLY	twoxcoef;		/*  2^x expansion series coefficients  */MULTIPOLY	coscoef;		/*  cos(x)  expansion series coefficients  */MULTIPOLY       log2coef;               /*  log_2(x) expansion series coefficients */FLOAT		P2;			/*  PI/2  */FLOAT		ln2;			/*  ln(2)	*/FLOAT           bconst;                 /*  needed in log_2(x)  *//*  set a float to a constant 1  */void one( FLOAT *x){	null( x);	x->expnt = 1;	x->mntsa.e[MS_MNTSA] = 0x40000000;}/*  compute pi to 250 bits  or so.  	Uses formula arcsin(1/2) = pi/6 = pi/2 - 1 - sum(		1*3*5*...*(2k-1))/(2^3k (2k+1) k!)*/void calcpi( FLOAT *pi){	FLOAT tn, constant;	int i;		null( pi);	int_to_float( 1, &tn);	tn.expnt = -2;	int_to_float( 3, &constant);	divide( &tn, &constant, &tn);	add( pi, &tn, pi);	for( i=2; i<123; i++)	{		int_to_float( 2*i-1, &constant);		multiply(&constant, &constant, &constant);		multiply( &constant, &tn, &tn);		int_to_float( 2*i+1, &constant);		divide( &tn, &constant, &tn);		int_to_float( i, &constant);		divide( &tn, &constant, &tn);		tn.expnt -= 3;		add( &tn, pi, pi);	}	int_to_float( 1, &constant);	add( &constant, pi, pi);	int_to_float( 3, &constant);	multiply( &constant, pi, pi);}/*  Output of above routine is:tn.expnt = -256pi = E +3.14159265358979323846264338327950288419716939937\	510582097494459230781640628585258pi.expnt = 2pi.mntsa.e[] = {	 0x1d89cd8c,  0x105df53,  0x4533e63a,  0x94812704,  	 0xc06e0e68,  0x62633145,  0x10b4611a,  0x6487ed51                      }*//*  compute ln(2) to 256 bits.	ln(2) = 2 sum( 1/( (2k+1) * 3^(2k+1) ))	which converges in about 70+ terms.*/void calcln2( FLOAT *ln2){	int	k, epsilon, startxp;	FLOAT	tk, constant, nine;		null( ln2);	int_to_float( 3, &constant);	reciprical( &constant, &tk);	// gives me t0	startxp = tk.expnt;	multiply( &tk, &tk, &nine);	// additional 3^2k term	epsilon = 1;	k = 0;		while( epsilon > -256)	{		add( &tk, ln2, ln2);		int_to_float( 2*k+1, &constant);		multiply( &constant, &tk, &tk);		int_to_float( 2*k+3, &constant);		divide( &tk, &constant, &tk);		multiply( &nine, &tk, &tk);		epsilon = tk.expnt - startxp;		k++;	}	ln2->expnt++;		// final multiply by 2}/*  compute y = x^k	where k is a signed integer in range +/- 2^31 and x, y are FLOAT.	Copied from complex version.	Returns 1 if ok, 0 if x = 0 and k<0*/int intpwr( FLOAT *x, int k, FLOAT *y){	int signflag, n;	FLOAT z, t;	/*  initialize Knuth's algorithm A pg 442 V2  */	copy( x, &z);	if( k<0)	{		signflag = -1;		n = -k;	}	else	{		signflag = 0;		n = k;	}	int_to_float( 1, &t);	while( n)	{		if( n & 1) multiply( &t, &z, &t);		multiply( &z, &z, &z);		n >>= 1;	}	if( signflag) return reciprical( &t, y);	copy( &t, y);	return 1;}/*  This routine computes bessel functions for real arguments	for nth order to accuracy of 256 bits.  Accuracy is easy to 	change, assuming storage chages.  Purpose is coefficients	for exp and trig functions approximated by chebyshev	polynomials.	Enter with type = +1 for In(x) and type = -1 for Jn(x),		FLOAT pointing to x, integer n	Returns with Float y = Jn(x) or In(x).	uses |n|*/void bessel( int type, int n, FLOAT *x, FLOAT *y){	FLOAT	z2, z4, constant, t1, sum;	int		startxp, epsilon, j, k;		if( n<0) n = -n;	copy( x, &z2);	z2.expnt--;		// divide input by 2	multiply( &z2, &z2, &z4);		// z^2/4	if( type < 0) negate( &z4);	int_to_float( 1, &t1);	j = n;			// compute 1/n!	while( j>1)	{		int_to_float( j, &constant);		multiply( &constant, &t1, &t1);		j--;	}	reciprical( &t1, &sum);	startxp = sum.expnt;	copy( &sum, &t1);	/*  when next term exponent is 256 less than starting	exponent, we have more bits of accuracy.	j and k increment by 1 each time to create factorial	terms 1/k! and 1/(n+j)!*/	j = n;	k = 0;	epsilon = 1;	while( epsilon > -250)	{		j++;		k++;		int_to_float( j, &constant);		divide( &t1, &constant, &t1);		multiply( &z4, &t1, &t1);		int_to_float( k, &constant);		divide( &t1, &constant, &t1);		add( &t1, &sum, &sum);		epsilon = t1.expnt - startxp;	}	intpwr( &z2, n, &t1);	multiply( &t1, &sum, y);}/* create table of chebyshev polynomials.	Enter with maximum degree desired and array large	enough to hold the MULTIPOLY data pointers.	Returns array of pointers filled and maximum	degree successfully created.		T(n,x) = 2xT(n-1, x) - T(n-2, x)*/int gen_chebyshev( MULTIPOLY *chebary, int degree){	int 		numgen;	FLOAT	*fptr;	MULTIPOLY twox, temp;		if( degree < 0) return 0;	twox.degree = 1;	if( !get_space( &twox)) return 0;	fptr = Address(twox);	null( fptr);	fptr++;	int_to_float( 2, fptr);	numgen = 0;	chebary[0].degree = 0;	if( !get_space( chebary))  goto chebdie;	fptr = Address( chebary[0]);	int_to_float( 1, fptr);	numgen++;	if (degree < 1) return 1;	chebary[1].degree = 1;	if( !get_space( &chebary[1])) goto chebdie;	fptr = Address( chebary[1]);	null( fptr);	fptr++;	int_to_float( 1, fptr);	numgen++;	while( numgen <= degree)	{		chebary[numgen].degree = numgen;		if( !get_space( &chebary[numgen])) break;		multi_mul( twox, chebary[numgen - 1], &temp);		multi_sub( temp, chebary[numgen - 2], &chebary[numgen]);		free_space( &temp);		numgen++;	}chebdie:	free_space( &twox);	return numgen;}/*  compute coefficients of 2^x using Chebyshev polynomials.		2^x  =  I(0, ln(2)) + sum{ 2*I(n, ln(2))*T(n, x)} n = 1... oo		Enter with pointer to table of cheybshev polynomials, 	maximum degree of approximation and pointer to where	you want result.	Returns with power series in x of above, and max degree	actually computed.*/int calc_2x_coef( MULTIPOLY *cheb,  int maxdegree, 				MULTIPOLY *twoxcoef){	INDEX		i, j;	FLOAT		ibesl, *iptr;	MULTIPOLY	tnterm, sum;		char		test[32];		calcln2( &ln2);		sum.degree = 0;	if( !get_space( &sum))	{		printf(" no space left, calc_2x_coef \n");		return 0;	}	iptr = Address(sum);	bessel (+1, 0, &ln2, iptr);	for( i=1; i<=maxdegree; i++)	{		bessel( +1, i, &ln2, &ibesl);		ibesl.expnt++;		multi_dup( cheb[i], &tnterm);		for( j = (i&1); j<=i; j += 2)		{			iptr = Address( tnterm) + j;			multiply( &ibesl, iptr, iptr);		}		if( !multi_add( tnterm, sum, &sum)) break;		free_space( &tnterm);	}	if( i< maxdegree)  free_space( &tnterm);	else i = maxdegree;	multi_dup( sum, twoxcoef);	free_space( &sum);	return (i);}/*  compute coefficients of cos(x*PI/2) using Chebyshev polynomials.		cos(x*P2)  =   J(0, P2) + 2*sum{(-1)^n*J(2n, P2)*T(2n+1, x)} n = 1... oo		Enter with pointer to table of cheybshev polynomials, 	maximum degree of approximation and pointer to where	you want result.	Returns with power series in x^2 of above, and max degree	actually computed.  cos( x*PI/2) = sum{ a_j * (x^2)^j}.		NOTE:  All formulas in all the books I found are WRONG.  None have	the (-1)^n factor, but they act like it's there when constructing	the polynomial.*/int calc_cos_coef( MULTIPOLY *cheb,  int maxdegree, 				MULTIPOLY *coscoef){	INDEX		i, j, k;	FLOAT		jbesl, *jptr, *kptr;	MULTIPOLY	tnterm, sum;		char		test[32];	calcpi( &P2);	P2.expnt--;	sum.degree = 0;	if( !get_space( &sum))	{		printf(" no space left, calc_cos_coef \n");		return 0;	}	jptr = Address(sum);	bessel(-1, 0, &P2, jptr);	for( i=1; i<=maxdegree/2; i++)	{		k = 2*i;		bessel( -1, k, &P2, &jbesl);		jbesl.expnt++;		if( i&1) negate( &jbesl);		multi_dup( cheb[k], &tnterm);		jptr = Address( tnterm);		for( j = 0; j<=k; j += 2)		{			multiply( &jbesl, jptr, jptr);			jptr += 2;		}		if( !multi_add( tnterm, sum, &sum)) break;		free_space( &tnterm);	}	if( i< maxdegree/2)  	{		free_space( &tnterm);		i = 2*i;	}	else i = maxdegree;/*  convert to polynomial in x^2  */	j = sum.degree/2;	coscoef->degree = j;	if( !get_space(coscoef))	{		printf(" no space for cosine.\n");		return 0;	}	for( k=0; k<=j; k++)	{		kptr = AddressOf( coscoef) + k;		jptr = Address( sum) + 2*k;		copy( jptr, kptr);	}	free_space( &sum);	return (i);}/*  compute coefficients for logarithm base 2 expansion.    Assumes ln(2) already computed!    Enter with pointer to table of chebyshev polynomials,    maximum degree (should be odd) and    pointer to where you want result.    creates bconst for use in main routine.    Value returned is maximum degree actually computed.*/int calc_ln_coef( MULTIPOLY *cheb, int maxdegree, MULTIPOLY *log2coef){  INDEX i, j;  FLOAT *iptr, r, temp, root2, rsqrd;  MULTIPOLY sum, term;/*  log_2(t) = 0.5 + 4/ln(2) * sum( r^(2k+1)/(2k+1) * T_(2k+1)(x))      r = (2^.25 - 1)/(2^.25 + 1)      x = (t - 2^.5)/( t + 2^.5) * b