bigfunc.c

随机数算法

/****************************************************************
*										        *
*		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 = -256
pi = 
E +3.14159265358979323846264338327950288419716939937\
	510582097494459230781640628585258
pi.expnt = 2
pi.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