C++ M_get_stack_var函数代码示例

void	m_apm_sin_cos(M_APM sinv, M_APM cosv, int places, M_APM aa)
{
	M_APM	tmp5, tmp6, tmp7;

	tmp5 = M_get_stack_var();
	tmp6 = M_get_stack_var();
	tmp7 = M_get_stack_var();

	M_limit_angle_to_pi(tmp5, (places + 6), aa);
	M_4x_cos(tmp7, (places + 6), tmp5);

	/*
	 *   compute sin(x) = sqrt(1 - cos(x) ^ 2).
	 *
	 *   note that the sign of 'sin' will always be positive after the
	 *   sqrt call. we need to adjust the sign based on what quadrant
	 *   the original angle is in.
	 */

	M_cos_to_sin(tmp6, (places + 6), tmp7);
	if (tmp6->m_apm_sign != 0)
		tmp6->m_apm_sign = tmp5->m_apm_sign;

	m_apm_round(sinv, places, tmp6);
	m_apm_round(cosv, places, tmp7);
	M_restore_stack(3);
}

/*
 *      arccosh(x) == log [ x + sqrt(x^2 - 1) ]
 *
 *      x >= 1.0
 */
void	m_apm_arccosh(M_APM rr, int places, M_APM aa)
{
M_APM	tmp1, tmp2;
int     ii;

ii = m_apm_compare(aa, MM_One);

if (ii == -1)       /* x < 1 */
  {
   M_apm_log_error_msg(M_APM_RETURN, "\'m_apm_arccosh\', Argument < 1");
   M_set_to_zero(rr);
   return;
  }

tmp1 = M_get_stack_var();
tmp2 = M_get_stack_var();

m_apm_multiply(tmp1, aa, aa);
m_apm_subtract(tmp2, tmp1, MM_One);
m_apm_sqrt(tmp1, (places + 6), tmp2);
m_apm_add(tmp2, aa, tmp1);
m_apm_log(rr, places, tmp2);

M_restore_stack(2);
}

/*
        Calculate arctan using the identity :

                                      x
        arctan (x) == arcsin [ --------------- ]
                                sqrt(1 + x^2)

*/
void	m_apm_arctan(M_APM rr, int places, M_APM xx)
{
M_APM   tmp8, tmp9;

if (xx->m_apm_sign == 0)			/* input == 0 ?? */
  {
   M_set_to_zero(rr);
   return;
  }

if (xx->m_apm_exponent <= -4)			/* input close to 0 ?? */
  {
   M_arctan_near_0(rr, places, xx);
   return;
  }

if (xx->m_apm_exponent >= 4)			/* large input */
  {
   M_arctan_large_input(rr, places, xx);
   return;
  }

tmp8 = M_get_stack_var();
tmp9 = M_get_stack_var();

m_apm_multiply(tmp9, xx, xx);
m_apm_add(tmp8, tmp9, MM_One);
m_apm_sqrt(tmp9, (places + 6), tmp8);
m_apm_divide(tmp8, (places + 6), xx, tmp9);
m_apm_arcsin(rr, places, tmp8);
M_restore_stack(2);
}

void	M_log_solve_cubic(M_APM rr, int places, M_APM nn)
{
M_APM   tmp0, tmp1, tmp2, tmp3, guess;
int	ii, maxp, tolerance, local_precision;

guess = M_get_stack_var();
tmp0  = M_get_stack_var();
tmp1  = M_get_stack_var();
tmp2  = M_get_stack_var();
tmp3  = M_get_stack_var();

M_get_log_guess(guess, nn);

tolerance       = -(places + 4);
maxp            = places + 16;
local_precision = 18;

/*    Use the following iteration to solve for log :

                        exp(X) - N 
      X     =  X - 2 * ------------
       n+1              exp(X) + N 

   
      this is a cubically convergent algorithm 
      (each iteration yields 3X more digits)
*/

ii = 0;

while (TRUE)
  {
   m_apm_exp(tmp1, local_precision, guess);

   m_apm_subtract(tmp3, tmp1, nn);
   m_apm_add(tmp2, tmp1, nn);

   m_apm_divide(tmp1, local_precision, tmp3, tmp2);
   m_apm_multiply(tmp0, MM_Two, tmp1);
   m_apm_subtract(tmp3, guess, tmp0);

   if (ii != 0)
     {
      if (((3 * tmp0->m_apm_exponent) < tolerance) || (tmp0->m_apm_sign == 0))
        break;
     }

   m_apm_round(guess, local_precision, tmp3);

   local_precision *= 3;

   if (local_precision > maxp)
     local_precision = maxp;

   ii = 1;
  }

m_apm_round(rr, places, tmp3);
M_restore_stack(5);
}

void	M_4x_cos(M_APM r, int places, M_APM x)
{
M_APM   tmp8, tmp9;

tmp8 = M_get_stack_var();
tmp9 = M_get_stack_var();

/* 
 *  if  |x| >= 1.0   use multiple angle identity 4 times
 *  if  |x|  < 1.0   use multiple angle identity 3 times
 */

if (x->m_apm_exponent > 0)
  {
   m_apm_multiply(tmp9, x, MM_5x_256R);        /* 1 / (4*4*4*4) */
   M_raw_cos(tmp8, (places + 8), tmp9);
   M_4x_do_it(tmp9, (places + 8), tmp8);
   M_4x_do_it(tmp8, (places + 6), tmp9);
   M_4x_do_it(tmp9, (places + 4), tmp8);
   M_4x_do_it(r, places, tmp9);
  }
else
  {
   m_apm_multiply(tmp9, x, MM_5x_64R);         /* 1 / (4*4*4) */
   M_raw_cos(tmp8, (places + 6), tmp9);
   M_4x_do_it(tmp9, (places + 4), tmp8);
   M_4x_do_it(tmp8, (places + 4), tmp9);
   M_4x_do_it(r, places, tmp8);
  }

M_restore_stack(2);
}

/*
 *      arctanh(x) == 0.5 * log [ (1 + x) / (1 - x) ]
 *
 *      |x| < 1.0
 */
void	m_apm_arctanh(M_APM rr, int places, M_APM aa)
{
M_APM	tmp1, tmp2, tmp3;
int     ii, local_precision;

tmp1 = M_get_stack_var();

m_apm_absolute_value(tmp1, aa);

ii = m_apm_compare(tmp1, MM_One);

if (ii >= 0)       /* |x| >= 1.0 */
  {
   M_apm_log_error_msg(M_APM_RETURN, "\'m_apm_arctanh\', |Argument| >= 1");
   M_set_to_zero(rr);
   M_restore_stack(1);
   return;
  }

tmp2 = M_get_stack_var();
tmp3 = M_get_stack_var();

local_precision = places + 8;

m_apm_add(tmp1, MM_One, aa);
m_apm_subtract(tmp2, MM_One, aa);
m_apm_divide(tmp3, local_precision, tmp1, tmp2);
m_apm_log(tmp2, local_precision, tmp3);
m_apm_multiply(tmp1, tmp2, MM_0_5);
m_apm_round(rr, places, tmp1);

M_restore_stack(3);
}

/*
 *      arcsinh(x) == log [ x + sqrt(x^2 + 1) ]
 *
 *      also, use arcsinh(-x) == -arcsinh(x)
 */
void	m_apm_arcsinh(M_APM rr, int places, M_APM aa)
{
M_APM	tmp0, tmp1, tmp2;

/* result is 0 if input is 0 */

if (aa->m_apm_sign == 0)
  {
   M_set_to_zero(rr);
   return;
  }

tmp0 = M_get_stack_var();
tmp1 = M_get_stack_var();
tmp2 = M_get_stack_var();

m_apm_absolute_value(tmp0, aa);
m_apm_multiply(tmp1, tmp0, tmp0);
m_apm_add(tmp2, tmp1, MM_One);
m_apm_sqrt(tmp1, (places + 6), tmp2);
m_apm_add(tmp2, tmp0, tmp1);
m_apm_log(rr, places, tmp2);

rr->m_apm_sign = aa->m_apm_sign; 			  /* fix final sign */

M_restore_stack(3);
}

/*
 *     calculate the multiple angle identity for sin (5x)
 *
 *     sin (5x) == 16 * sin^5 (x)  -  20 * sin^3 (x)  +  5 * sin(x)  
 */
void	M_5x_do_it(M_APM rr, int places, M_APM xx)
{
M_APM   tmp0, tmp1, t2, t3, t5;

tmp0 = M_get_stack_var();
tmp1 = M_get_stack_var();
t2   = M_get_stack_var();
t3   = M_get_stack_var();
t5   = M_get_stack_var();

m_apm_multiply(tmp1, xx, xx);
m_apm_round(t2, (places + 4), tmp1);     /* x ^ 2 */

m_apm_multiply(tmp1, t2, xx);
m_apm_round(t3, (places + 4), tmp1);     /* x ^ 3 */

m_apm_multiply(t5, t2, t3);              /* x ^ 5 */

m_apm_multiply(tmp0, xx, MM_Five);
m_apm_multiply(tmp1, t5, MM_5x_Sixteen);
m_apm_add(t2, tmp0, tmp1);
m_apm_multiply(tmp1, t3, MM_5x_Twenty);
m_apm_subtract(tmp0, t2, tmp1);

m_apm_round(rr, places, tmp0);
M_restore_stack(5);
}

/*
 *      From Knuth, The Art of Computer Programming:
 *
 *      To compute GCD(u,v)
 *          
 *      A1:
 *	    if (v == 0)  return (u)
 *      A2:
 *          t = u mod v
 *	    u = v
 *	    v = t
 *	    goto A1
 */
void	m_apm_gcd_traditional(M_APM r, M_APM u, M_APM v)
{
M_APM   tmpD, tmpN, tmpU, tmpV;

tmpD = M_get_stack_var();
tmpN = M_get_stack_var();
tmpU = M_get_stack_var();
tmpV = M_get_stack_var();

m_apm_absolute_value(tmpU, u);
m_apm_absolute_value(tmpV, v);

while (TRUE)
  {
   if (tmpV->m_apm_sign == 0)
     break;

   m_apm_integer_div_rem(tmpD, tmpN, tmpU, tmpV);
   m_apm_copy(tmpU, tmpV);
   m_apm_copy(tmpV, tmpN);
  }

m_apm_copy(r, tmpU);
M_restore_stack(4);
}

/*
 *      find log(N)
 *
 *      if places < 360
 *         solve with cubically convergent algorithm above
 *
 *      else
 *
 *      let 'X' be 'close' to the solution   (we use ~110 decimal places)
 *
 *      let Y = N * exp(-X) - 1
 *
 *	then
 *
 *	log(N) = X + log(1 + Y)
 *
 *      since 'Y' will be small, we can use the efficient log_near_1 algorithm.
 *
 */
void	M_log_basic_iteration(M_APM rr, int places, M_APM nn)
{
M_APM   tmp0, tmp1, tmp2, tmpX;

if (places < 360)
  {
   M_log_solve_cubic(rr, places, nn);
  }
else
  {
   tmp0 = M_get_stack_var();
   tmp1 = M_get_stack_var();
   tmp2 = M_get_stack_var();
   tmpX = M_get_stack_var();
   
   M_log_solve_cubic(tmpX, 110, nn);
   
   m_apm_negate(tmp0, tmpX);
   m_apm_exp(tmp1, (places + 8), tmp0);
   m_apm_multiply(tmp2, tmp1, nn);
   m_apm_subtract(tmp1, tmp2, MM_One);
   
   M_log_near_1(tmp0, (places - 104), tmp1);
   
   m_apm_add(tmp1, tmpX, tmp0);
   m_apm_round(rr, places, tmp1);
   
   M_restore_stack(4);
  }
}

/*
	compute  int *n  = round_to_nearest_int(a / log(2))
	         M_APM b = MAPM version of *n

        returns      0: OK
		 -1, 1: failure
*/
int	M_exp_compute_nn(int *n, M_APM b, M_APM a)
{
M_APM	tmp0, tmp1;
void	*vp;
char    *cp, sbuf[48];
int	kk;

*n   = 0;
vp   = NULL;
cp   = sbuf;
tmp0 = M_get_stack_var();
tmp1 = M_get_stack_var();

/* find 'n' and convert it to a normal C int            */
/* we just need an approx 1/log(2) for this calculation */

m_apm_multiply(tmp1, a, MM_exp_log2R);

/* round to the nearest int */

if (tmp1->m_apm_sign >= 0)
  {
   m_apm_add(tmp0, tmp1, MM_0_5);
   m_apm_floor(tmp1, tmp0);
  }
else
  {
   m_apm_subtract(tmp0, tmp1, MM_0_5);
   m_apm_ceil(tmp1, tmp0);
  }

kk = tmp1->m_apm_exponent;
if (kk >= 42)
  {
   if ((vp = (void *)MAPM_MALLOC((kk + 16) * sizeof(char))) == NULL)
     {
      /* fatal, this does not return */

      M_apm_log_error_msg(M_APM_FATAL, "\'M_exp_compute_nn\', Out of memory");
     }

   cp = (char *)vp;
  }

m_apm_to_integer_string(cp, tmp1);
*n = atoi(cp);

m_apm_set_long(b, (long)(*n));

kk = m_apm_compare(b, tmp1);

if (vp != NULL)
  MAPM_FREE(vp);

M_restore_stack(2);
return(kk);
}

/*
 *   compute  r = sqrt(1 - a ^ 2).
 */
void	M_cos_to_sin(M_APM r, int places, M_APM a)
{
M_APM	tmp1, tmp2;

tmp1 = M_get_stack_var();
tmp2 = M_get_stack_var();

m_apm_multiply(tmp1, a, a);
m_apm_subtract(tmp2, MM_One, tmp1);
m_apm_sqrt(r, places, tmp2);
M_restore_stack(2);
}

void	m_apm_lcm(M_APM r, M_APM u, M_APM v)
{
M_APM   tmpN, tmpG;

tmpN = M_get_stack_var();
tmpG = M_get_stack_var();

m_apm_multiply(tmpN, u, v);
m_apm_gcd(tmpG, u, v);
m_apm_integer_divide(r, tmpN, tmpG);

M_restore_stack(2);
}

/*
        Calculate arcsin using the identity :

                                      x
        arcsin (x) == arctan [ --------------- ]
                                sqrt(1 - x^2)

*/
void	M_arcsin_near_0(M_APM rr, int places, M_APM aa)
{
M_APM   tmp5, tmp6;

tmp5 = M_get_stack_var();
tmp6 = M_get_stack_var();

M_cos_to_sin(tmp5, (places + 8), aa);
m_apm_divide(tmp6, (places + 8), aa, tmp5);
M_arctan_near_0(rr, places, tmp6);

M_restore_stack(2);
}

/*
        Calculate arccos using the identity :

        arccos (x) == PI / 2 - arcsin (x)

*/
void	M_arccos_near_0(M_APM rr, int places, M_APM aa)
{
M_APM   tmp1, tmp2;

tmp1 = M_get_stack_var();
tmp2 = M_get_stack_var();

M_check_PI_places(places);
M_arcsin_near_0(tmp1, (places + 4), aa);
m_apm_subtract(tmp2, MM_lc_HALF_PI, tmp1);
m_apm_round(rr, places, tmp2);

M_restore_stack(2);
}