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);
}

/*
	calculate log (1 + x) with the following series:

              x
	y = -----      ( |y| < 1 )
	    x + 2


            [ 1 + y ]                 y^3     y^5     y^7
	log [-------]  =  2 * [ y  +  ---  +  ---  +  ---  ... ] 
            [ 1 - y ]                  3       5       7 

*/
void	M_log_near_1(M_APM rr, int places, M_APM xx)
{
M_APM   tmp0, tmp1, tmp2, tmpS, term;
int	tolerance, dplaces, local_precision;
long    m1;

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

tolerance = xx->m_apm_exponent - (places + 6);
dplaces   = (places + 12) - xx->m_apm_exponent;

m_apm_add(tmp0, xx, MM_Two);
m_apm_divide(tmpS, (dplaces + 6), xx, tmp0);

m_apm_copy(term, tmpS);
m_apm_multiply(tmp0, tmpS, tmpS);
m_apm_round(tmp2, (dplaces + 6), tmp0);

m1 = 3L;

while (TRUE)
  {
   m_apm_multiply(tmp0, term, tmp2);

   if ((tmp0->m_apm_exponent < tolerance) || (tmp0->m_apm_sign == 0))
     break;

   local_precision = dplaces + tmp0->m_apm_exponent;

   if (local_precision < 20)
     local_precision = 20;

   m_apm_set_long(tmp1, m1);
   m_apm_round(term, local_precision, tmp0);
   m_apm_divide(tmp0, local_precision, term, tmp1);
   m_apm_add(tmp1, tmpS, tmp0);
   m_apm_copy(tmpS, tmp1);
   m1 += 2;
  }

m_apm_multiply(tmp0, MM_Two, tmpS);
m_apm_round(rr, places, tmp0);

M_restore_stack(5);                    /* restore the 5 locals we used here */
}

void	m_apm_divide(M_APM rr, int places, M_APM aa, M_APM bb)
{
M_APM   tmp0, tmp1;
int     sn, nexp, dplaces;

sn = aa->m_apm_sign * bb->m_apm_sign;

if (sn == 0)                  /* one number is zero, result is zero */
  {
   if (bb->m_apm_sign == 0)
     {
      M_apm_log_error_msg(M_APM_RETURN, 
                          "Warning! ... \'m_apm_divide\', Divide by 0");
     }

   M_set_to_zero(rr);
   return;
  }

/*
 *    Use the original 'Knuth' method for smaller divides. On the
 *    author's system, this was the *approx* break even point before
 *    the reciprocal method used below became faster.
 */

if (places < 250)
  {
   M_apm_sdivide(rr, places, aa, bb);
   return;
  }

/* mimic the decimal place behavior of the original divide */

nexp = aa->m_apm_exponent - bb->m_apm_exponent;

if (nexp > 0)
  dplaces = nexp + places;
else
  dplaces = places;

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

m_apm_reciprocal(tmp0, (dplaces + 8), bb);
m_apm_multiply(tmp1, tmp0, aa);
m_apm_round(rr, dplaces, tmp1);

M_restore_stack(2);
}

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

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

m_apm_multiply(tmp9, x, MM_5x_125R);        /* 1 / (5*5*5) */
M_raw_sin(tmp8, (places + 6), tmp9);
M_5x_do_it(tmp9, (places + 4), tmp8);
M_5x_do_it(tmp8, (places + 4), tmp9);
M_5x_do_it(r, places, tmp8);

M_restore_stack(2);
}

/*
        Calls the LOG function. The formula used is :

        log10(x)  =  A * log(x) where A = log  (e)  [0.43429448190325...]
                                             10
*/
void	m_apm_log10(M_APM rr, int places, M_APM aa)
{
int     dplaces;
M_APM   tmp8, tmp9;

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

dplaces = places + 4;
M_check_log_places(dplaces + 45);

m_apm_log(tmp9, dplaces, aa);
m_apm_multiply(tmp8, tmp9, MM_lc_log10R);
m_apm_round(rr, places, tmp8);
M_restore_stack(2);                    /* restore the 2 locals we used here */
}

void	m_apm_tan(M_APM r, int places, M_APM a)
{
	M_APM	tmps, tmpc, tmp0;

	tmps = M_get_stack_var();
	tmpc = M_get_stack_var();
	tmp0 = M_get_stack_var();

	m_apm_sin_cos(tmps, tmpc, (places + 4), a);

	/* tan(x) = sin(x) / cos(x) */

	m_apm_divide(tmp0, (places + 4), tmps, tmpc);
	m_apm_round(r, places, tmp0);
	M_restore_stack(3);
}