mapm_exp.c

任意精度的数学库

/* 
 *  M_APM  -  mapm_exp.c
 *
 *  Copyright (C) 1999   Michael C. Ring
 *
 *  Permission to use, copy, and distribute this software and its
 *  documentation for any purpose with or without fee is hereby granted, 
 *  provided that the above copyright notice appear in all copies and 
 *  that both that copyright notice and this permission notice appear 
 *  in supporting documentation.
 *
 *  Permission to modify the software is granted, but not the right to
 *  distribute the modified code.  Modifications are to be distributed 
 *  as patches to released version.
 *  
 *  This software is provided "as is" without express or implied warranty.
 */

/*
 *      $Id: mapm_exp.c,v 1.11 1999/09/18 01:27:40 mike Exp $
 *
 *      This file contains the POW and EXP functions.
 *
 *      $Log: mapm_exp.c,v $
 *      Revision 1.11  1999/09/18 01:27:40  mike
 *      if X is 0 on the pow function, return 0 right away
 *
 *      Revision 1.10  1999/06/19 20:54:07  mike
 *      changed local static MAPM to stack variables
 *
 *      Revision 1.9  1999/06/01 22:37:44  mike
 *      adjust decimal places passed to raw function
 *
 *      Revision 1.8  1999/06/01 01:44:03  mike
 *      change threshold from 1000 to 100 for 65536 divisor
 *
 *      Revision 1.7  1999/06/01 01:03:31  mike
 *      vary 'q' instead of checking input against +/- 10 and +/- 40
 *
 *      Revision 1.6  1999/05/15 01:54:27  mike
 *      add check for number of decimal places
 *
 *      Revision 1.5  1999/05/15 01:09:49  mike
 *      minor tweak to POW decimal places
 *
 *      Revision 1.4  1999/05/13 00:14:00  mike
 *      added more comments
 *
 *      Revision 1.3  1999/05/12 23:39:05  mike
 *      change #places passed to sub functions
 *
 *      Revision 1.2  1999/05/10 21:35:13  mike
 *      added some comments
 *
 *      Revision 1.1  1999/05/10 20:56:31  mike
 *      Initial revision
 */

#include "m_apm_lc.h"

/****************************************************************************/
/*
	Calculate the POW function by calling EXP :

                  Y      A                 
                 X   =  e    where A = Y * log(X)
*/
void	m_apm_pow(rr,places,xx,yy)
M_APM	rr, xx, yy;
int	places;
{
M_APM   tmp8;
M_APM   tmp9;

if (xx->m_apm_sign == 0)
  {
   rr->m_apm_datalength = 1;
   rr->m_apm_sign       = 0;
   rr->m_apm_exponent   = 0;
   rr->m_apm_data[0]    = 0;
   return;
  }

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

/* if 'X' is 2 or 10, use the predefined constants */

if (m_apm_compare(xx, MM_Ten) == 0)
  {
   M_check_dec_places(M_POW, places);
   m_apm_round(tmp9, (places + 8), MM_LOG_10_BASE_E);
  }
else
  {
   if (m_apm_compare(xx, MM_Two) == 0)
     {
      M_check_dec_places(M_POW, places);
      m_apm_round(tmp9, (places + 8), MM_LOG_2_BASE_E);
     }
   else
     {
      m_apm_log(tmp9, (places + 4), xx);
     }
  }

m_apm_multiply(tmp8, tmp9, yy);
m_apm_exp(tmp9, (places + 4), tmp8);
m_apm_round(rr, places, tmp9);
M_restore_stack(2);                    /* restore the 2 locals we used here */
}
/****************************************************************************/
void	m_apm_exp(r,places,x)
M_APM	r, x;
int	places;
{
M_APM   tmp8;
M_APM   tmp9;
static  M_APM  MM_exp_256R;
static  M_APM  MM_exp_4096R;
static  M_APM  MM_exp_65536R;
static	int    firsttime1 = TRUE;
int	dplaces, nn, ii;

if (firsttime1)
  {
   firsttime1    = FALSE;
   MM_exp_256R   = m_apm_init();
   MM_exp_4096R  = m_apm_init();
   MM_exp_65536R = m_apm_init();

   m_apm_set_string(MM_exp_256R, "3.90625E-3");         /* 1 / 256   */
   m_apm_set_string(MM_exp_4096R, "2.44140625E-4");     /* 1 / 4096  */
   m_apm_set_string(MM_exp_65536R, "1.52587890625E-5"); /* 1 / 65536 */
  }

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

/*
    From David H. Bailey's MPFUN Fortran package :

    exp (t) =  (1 + r + r^2 / 2! + r^3 / 3! + r^4 / 4! ...) ^ q * 2 ^ n

    where q = 256, r = t' / q, t' = t - n Log(2) and where n is chosen so
    that -0.5 Log(2) < t' <= 0.5 Log(2).  Reducing t mod Log(2) and
    dividing by 256 insures that -0.001 < r <= 0.001, which accelerates
    convergence in the above series.
    
    In the MAPM implementation, set n = 0 (we need an efficient way to 
    compute 2 ^ n) but use a different 'q' :

    if |input| < 1        :  use q = 256
    if |input|  1 - 100   :  use q = 4096
    if |input| >= 100     :  use q = 65536
*/

if (x->m_apm_exponent >= 3)        /* >= 100 */
  {
   nn = 16;
   dplaces = places + 7;
   m_apm_multiply(tmp9, x, MM_exp_65536R);
  }
else
  {
   if (x->m_apm_exponent <= 0)     /* < 1 */
     {
      nn = 8;
      dplaces = places + 4;
      m_apm_multiply(tmp9, x, MM_exp_256R);
     }
   else                            /* 1 - 100 */
     {
      nn = 12;
      dplaces = places + 6;
      m_apm_multiply(tmp9, x, MM_exp_4096R);
     }
  }

/* perform the series expansion ... */

M_raw_exp(tmp8, dplaces, tmp9);

/* 
 *   raise result to the 256 / 4096 / 65536 power (q = 256 / 4096 / 65536) 
 *
 *   note : x ^ 256    == (((x ^ 2) ^ 2) ^ 2) ... 8  times
 *        : x ^ 4096   == (((x ^ 2) ^ 2) ^ 2) ... 12 times
 *        : x ^ 65536  == (((x ^ 2) ^ 2) ^ 2) ... 16 times
 */

for (ii=1; ii < nn; ii++)              /* do first N-1 multiplies */
  {
   m_apm_multiply(tmp9, tmp8, tmp8);
   m_apm_round(tmp8, dplaces, tmp9);
  }
				       /* then do the last one */
m_apm_multiply(tmp9, tmp8, tmp8);
m_apm_round(r, places, tmp9);
M_restore_stack(2);                    /* restore the 2 locals we used here */
}
/****************************************************************************/
/*
	calculate the exponential function using the following series :

                              x^2     x^3     x^4    x^5
	exp(x) == 1  +  x  +  ---  +  ---  +  ---  + ---  ...
                               2!      3!      4!     5!

*/
void	M_raw_exp(rr,places,xx)
M_APM	rr, xx;
int	places;
{
M_APM   tolerance;
M_APM   tmp0;
M_APM   digit;
M_APM   term;
int	flag,  local_precision;
long    m1;

tolerance = M_get_stack_var();
tmp0      = M_get_stack_var();
term      = M_get_stack_var();
digit     = M_get_stack_var();

local_precision = places + 4;

m_apm_copy(tolerance, MM_One);
tolerance->m_apm_exponent = -local_precision;

m_apm_add(rr, MM_One, xx);
m_apm_copy(term, xx);

flag = 0;
m1   = 2;

while (TRUE)
  {
   m_apm_set_long(digit, m1);
   m_apm_multiply(tmp0, term, xx);
   m_apm_divide(term, local_precision, tmp0, digit);
   m_apm_add(tmp0, rr, term);
   m_apm_copy(rr, tmp0);

   /*
    *  only check tolerance on even powers of 'x', guaranteed 
    *  to be positive so it saves an absolute value
    */

   if (flag == 0)
     {
      if (m_apm_compare(term, tolerance) <= 0)
        break;
     }

   m1++;
   flag = 1 - flag;
  }

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