bigfunc.c

随机数算法

      b = (2^.5 + 1)/(2^.5 - 1)
*/

  sum.degree = 0;
  if( !get_space( &sum)) 
  {
    printf("no space left, calc_ln_coef\n");
    return 0;
  }

/*  compute bconst first  */

  int_to_float( 2, &r);
  square_root( &r, &root2);
  one(&temp);
  add( &root2, &temp, &r);
  subtract( &root2, &temp, &bconst);
  divide( &r, &bconst, &bconst);

/*  now compute r and r^2 */

  square_root( &root2, &r);
  subtract( &r, &temp, &root2);
  add( &r, &temp, &r);
  divide( &root2, &r, &r);
  multiply( &r, &r, &rsqrd);

  for( i=1; i<=maxdegree; i += 2)
  {
    int_to_float( i, &temp);
    divide( &r, &temp, &temp);
    multi_dup( cheb[i], &term);
    for( j=1; j<=i; j+=2)
    {
      iptr = Address( term) + j;
      multiply( &temp, iptr, iptr);
    }
    if( !multi_add( term, sum, &sum)) break;
    free_space( &term);
    multiply( &rsqrd, &r, &r);  /*  r^2k+1  */
  }
  if( i< maxdegree) free_space( &term);
  else i = maxdegree;

/*  multiply all coefficients by 4/ln(2)  */

  one( &temp);
  temp.expnt += 2;
  divide( &temp, &ln2, &temp);
  for( j=1; j<=i; j+=2)
  {
    iptr = Address(sum) + j;
    multiply( &temp, iptr, iptr);
  }
  multi_dup( sum, log2coef);
  free_space(&sum);
  return(i);
}

/*  evaluate simple polynomial.
	Input:  MULTIPOLY coefficients, pointer to x, pointer to y
	Output:  y = F(x)
	y can equal x
*/

void polyeval( MULTIPOLY coef, FLOAT *x, FLOAT *y)
{
	INDEX 	i;
	FLOAT	sum, *cof;
	
	null( &sum);
	for( i=coef.degree; i>0; i--)
	{
		cof = Address( coef) + i;
		add( cof, &sum, &sum);
		multiply( x, &sum, &sum);
	}
	cof = Address( coef);
	add( cof, &sum, y);
}

/*  compute 2^x for x in the range -1 ... 1.
	Most inputs will be in range +/- .5 ... 1 but this routine could handle
	unnormalized inputs.
	Enter with pointer to input and storage for output
	Returns y = 2^x ( works in place, both pointers can be the same)
*/

void twoexp( FLOAT *x, FLOAT *y)
{
	polyeval( twoxcoef, x, y);
}

/*  compute cos( x) for x in range +/- PI/2
	As above, works in place.
	This is a *core* routine, no range checking!
*/

void corecos(FLOAT *x, FLOAT *y)
{
	FLOAT	x2;
	
	divide( x, &P2, &x2);
	multiply( &x2, &x2, &x2);
	polyeval( coscoef, &x2, y);
}

/*  convert a float to a long.  Overflow is max
	possible result.
*/

int	float_to_int( FLOAT *f)
{
	FLOAT	dummy, intprt;
	int		value;
	
	if( f->expnt < 1) return 0;
	if( f->expnt > 31)
	{
		if( f->mntsa.e[MS_MNTSA] & SIGN_BIT)
			return SIGN_BIT;
		return ~SIGN_BIT;
	}
	split( f, &intprt, &dummy);
	if( intprt.mntsa.e[MS_MNTSA] & SIGN_BIT)
	{
		negate( &intprt);
		value = -(intprt.mntsa.e[MS_MNTSA] >> ( 31 - intprt.expnt));
	}
	else	value = intprt.mntsa.e[MS_MNTSA] >> ( 31 - intprt.expnt);
	return value;
}

/*  compute log base 2 of FLOAT.
    range of input for valid result is 1 <= t <= 2.
    See page 52 of "Chebyshev Methods in Numerical Approximation", 
    M. A. Snyder, Prentice-Hall, 1966 for details
    input pointer can equal output pointer and it will still work.
*/

void logbase2( FLOAT *t, FLOAT *y)
{
  INDEX  i;
  FLOAT  x, temp, root2, sum, xsqrd, *cofptr;

/*  construct x = (t - 2^.5)/(t + 2^.5) * bconst  */

  int_to_float( 2, &temp);
  square_root( &temp, &root2);
  subtract( t, &root2, &temp);
  add( t, &root2, &x);
  divide( &temp, &x, &x);
  multiply( &x, &bconst, &x);
  null( &sum);

/*  use x^2 to do polynomial expansion since only odd powers used  */

  multiply( &x, &x, &xsqrd);
  for( i=log2coef.degree; i>1; i -= 2)
  {
    cofptr = Address( log2coef) + i;
    add( cofptr, &sum, &sum);
    multiply( &xsqrd, &sum, &sum);
  }

/*  add on last term to mak it all odd powers */

  cofptr = Address( log2coef) + 1;
  add( cofptr, &sum, &sum);
  multiply( &x, &sum, &sum);

/*  add final 0.5  */

  one( &temp);
  temp.expnt--;
  add( &temp, &sum, y);
}

/*  compute e^x for any x.  |x| > 2^32/ln(2) will overflow
	and return max possible value and 0.
	Otherwise returns y = exp(x) and 1.
	works in place.
*/

int big_exp( FLOAT *x, FLOAT *y)
{
	FLOAT	z, xp;
	long		xpnt;
	INDEX	i;
	
/*  convert to base 2  */

	divide( x, &ln2, &z);
	
/*  check range is possible to do  */

	if (z.expnt > 32)
	{
		if( x->mntsa.e[MS_MNTSA] & SIGN_BIT)
			null(y);
		else
		{
			OPLOOP(i) y->mntsa.e[MS_MNTSA] = ~0;
			y->mntsa.e[MS_MNTSA] >>= 1;
			y->expnt = ~0 >> 1;
		}
		return 0;
	}
	
/*  we can perform operation, send z mods 2 to core */

	split(&z, &xp, &z);
	xpnt = float_to_int( &xp);
	twoexp( &z, y);
	
/*  next add xpnt to exponent of y  */

	y->expnt +=  xpnt;
	return 1;
}

/*  split a FLOAT into its integer and fractional parts  */

void split( FLOAT *x, FLOAT *intprt, FLOAT *frac)
{
	FLOAT	fracpart;
	INDEX	i, signflag;
	ELEMENT	mask, xpchk;

/*  if number < 1, return just a fraction  */
	
	if( x->expnt <= 0)
	{
		null( intprt);
		copy( x, frac);
		return;
	}

/*  zero out 31 bits of ms word, then one block at a time */

	copy( x, &fracpart);
	signflag = 0;
	if( fracpart.mntsa.e[MS_MNTSA] & SIGN_BIT)
	{
		negate( &fracpart);
		signflag = 1;
	}
	i = MS_MNTSA;
	xpchk = fracpart.expnt;
	if( xpchk > 31)
	{
		fracpart.mntsa.e[i] = 0;
		i--;
		xpchk -= 31;
	}
	while( (xpchk > 32) && ( i>0 ))
	{
		fracpart.mntsa.e[i] = 0;
		i--;
		xpchk -= 32;
	}

/*  next zero out the remaining integer bits in a left over word  */

	if( i != MS_MNTSA) mask = ~0UL >> xpchk;
	else 	mask = ( ~0UL >> (xpchk+1));
	fracpart.mntsa.e[i] &= mask;
	if( signflag) negate( &fracpart);
	normal( &fracpart);
	subtract( x, &fracpart, intprt);
	copy( &fracpart, frac);
}

/*  compute cosine(x) for any x.
	x values larger than PI*2^200 will be in gross error, so watch out!
	works in place, returns y = cos(x)
*/

void cosine( FLOAT *x, FLOAT *y)
{
	FLOAT	z, PI, dummy, PI3;
	int		cmpr;
	
/*  create 2*PI  */

	copy( &P2, &PI);
	PI.expnt += 2;

/*  check range of input and force modulo 2PI operation  */

	cmpr = compare( x, &PI);
	if( cmpr > 0 )
	{
		divide( x, &PI, &z);
		split( &z, &dummy, &z);
		multiply( &PI, &z, &z);
	}
	else copy( x, &z);
	if( z.mntsa.e[MS_MNTSA] & SIGN_BIT) negate( &z);

/*  z is now in range 0...2PI.  Now convert to range of core cos */

	cmpr = compare( &z, &P2);
	if( cmpr <= 0)
	{
		corecos( &z, y);
		return;
	}
	PI.expnt--;
	add( &PI, &P2, &PI3);	// 3 PI/2
	cmpr = compare( &z, &PI3);
	if( cmpr > 0)
	{
		PI.expnt++;
		subtract( &PI, &z, &z);	// 2PI - x
		corecos( &z, y);
		return;
	}
	subtract( &PI, &z, &z);	// PI - x
	corecos( &z, y);
	negate( y);
}

/*  compute sine(x) for any x.
	same as cosine, jus move arguments around.
	works in place, 
	returns y = sin(x)
*/

void sine( FLOAT *x, FLOAT *y)
{
	FLOAT	z, PI, dummy;
	int		cmpr, signflag;
	
/*  create 2*PI and reduce x modulo 2PI signed  */

	copy( &P2, &PI);
	PI.expnt += 2;
	cmpr = compare( x, &PI);
	if( cmpr > 0)
	{
		divide( x, &PI, &z);
		split( &z, &dummy, &z);
		multiply( &PI, &z, &z);
	}
	else copy( x, &z);
	if( z.mntsa.e[MS_MNTSA] & SIGN_BIT)
	{
		negate( &z);
		signflag = 1;
	}
	else signflag = 0;

/*  z is no in range 0... 2*PI.
	if bigger than PI, flip sign of result and fold back to 0..PI,
	then subtract PI/2 to put into corecos range.
*/
	PI.expnt--;
	cmpr = compare( &z, &PI);
	if( cmpr > 0)
	{
		signflag ^= 1;
		subtract( &z, &PI, &z);
	}
	subtract( &z, &P2, &z);
	corecos( &z, y);
	if( signflag) negate(y);
}

/*  compute natural logarithm of big FLOAT value.
    values <= 0 return 0 and null result.
    Returns value to desired location.  Pointers can
    be the same and this will work.  Returns 1 if
    calculation valid.
*/

int big_ln( FLOAT *x, FLOAT *y)
{
  FLOAT  t, kn;
  ELEMENT n;

/*  check if input negative or zero  */

  if( (x->mntsa.e[MS_MNTSA] & SIGN_BIT) || iszero(x)) 
  {
    null( y);
    return 0;
  }

/*  split out fractional part in range 1..2 and feed to
    approximation routine  */

  n = x->expnt - 1;
  copy( x, &t);
  t.expnt = 1; 
  logbase2( &t, &t);
  int_to_float( n, &kn);
  add( &kn, &t, &t);
  multiply( &ln2, &t, y);
}