src/libs/mapm/mapm_lg3.c

   1
   2 /*
   3  *  M_APM  -  mapm_lg3.c
   4  *
   5  *  Copyright (C) 2003 - 2007   Michael C. Ring
   6  *
   7  *  Permission to use, copy, and distribute this software and its
   8  *  documentation for any purpose with or without fee is hereby granted,
   9  *  provided that the above copyright notice appear in all copies and
  10  *  that both that copyright notice and this permission notice appear
  11  *  in supporting documentation.
  12  *
  13  *  Permission to modify the software is granted. Permission to distribute
  14  *  the modified code is granted. Modifications are to be distributed by
  15  *  using the file 'license.txt' as a template to modify the file header.
  16  *  'license.txt' is available in the official MAPM distribution.
  17  *
  18  *  This software is provided "as is" without express or implied warranty.
  19  */
  20
  21 /*
  22  *      $Id: mapm_lg3.c,v 1.7 2007/12/03 01:42:59 mike Exp $
  23  *
  24  *      This file contains the function to compute log(2), log(10),
  25  *      and 1/log(10) to the desired precision using an AGM algorithm.
  26  *
  27  *      $Log: mapm_lg3.c,v $
  28  *      Revision 1.7  2007/12/03 01:42:59  mike
  29  *      Update license
  30  *
  31  *      Revision 1.6  2003/12/09 01:25:06  mike
  32  *      actually compute the first term of the AGM iteration instead
  33  *      of assuming the inputs a=1 and b=10^-N.
  34  *
  35  *      Revision 1.5  2003/12/04 03:19:16  mike
  36  *      rearrange logic in AGM to be more straight-forward
  37  *
  38  *      Revision 1.4  2003/05/01 22:04:37  mike
  39  *      rearrange some code
  40  *
  41  *      Revision 1.3  2003/05/01 21:58:31  mike
  42  *      remove math.h
  43  *
  44  *      Revision 1.2  2003/03/30 22:14:58  mike
  45  *      add comments
  46  *
  47  *      Revision 1.1  2003/03/30 21:18:04  mike
  48  *      Initial revision
  49  */
  50
  51 #include "m_apm_lc.h"
  52
  53 /*
  54  *  using the 'R' function (defined below) for 'N' decimal places :
  55  *
  56  *
  57  *                          -N             -N
  58  *  log(2)  =  R(1, 0.5 * 10  )  -  R(1, 10  )
  59  *
  60  *
  61  *                          -N             -N
  62  *  log(10) =  R(1, 0.1 * 10  )  -  R(1, 10  )
  63  *
  64  *
  65  *  In general:
  66  *
  67  *                    -N                -N
  68  *  log(x)  =  R(1, 10  / x)  -  R(1, 10  )
  69  *
  70  *
  71  *  I found this on a web site which went into considerable detail
  72  *  on the history of log(2). This formula is algebraically identical
  73  *  to the formula specified in J. Borwein and P. Borwein's book
  74  *  "PI and the AGM". (reference algorithm 7.2)
  75  */
  76
  77 /****************************************************************************/
  78 /*
  79  *      check if our local copy of log(2) & log(10) is precise
  80  *      enough for our purpose. if not, calculate them so it's
  81  *      as precise as desired, accurate to at least 'places'.
  82  */
  83 void    M_check_log_places(int places)
  84 {
  85 M_APM   tmp6, tmp7, tmp8, tmp9;
  86 int     dplaces;
  87
  88 dplaces = places + 4;
  89
  90 if (dplaces > MM_lc_log_digits)
  91   {
  92    MM_lc_log_digits = dplaces + 4;
  93
  94    tmp6 = M_get_stack_var();
  95    tmp7 = M_get_stack_var();
  96    tmp8 = M_get_stack_var();
  97    tmp9 = M_get_stack_var();
  98
  99    dplaces += 6 + (int)log10((double)places);
 100
 101    m_apm_copy(tmp7, MM_One);
 102    tmp7->m_apm_exponent = -places;
 103
 104    M_log_AGM_R_func(tmp8, dplaces, MM_One, tmp7);
 105
 106    m_apm_multiply(tmp6, tmp7, MM_0_5);
 107
 108    M_log_AGM_R_func(tmp9, dplaces, MM_One, tmp6);
 109
 110    m_apm_subtract(MM_lc_log2, tmp9, tmp8);               /* log(2) */
 111
 112    tmp7->m_apm_exponent -= 1;                            /* divide by 10 */
 113
 114    M_log_AGM_R_func(tmp9, dplaces, MM_One, tmp7);
 115
 116    m_apm_subtract(MM_lc_log10, tmp9, tmp8);              /* log(10) */
 117    m_apm_reciprocal(MM_lc_log10R, dplaces, MM_lc_log10); /* 1 / log(10) */
 118
 119    M_restore_stack(4);
 120   }
 121 }
 122 /****************************************************************************/
 123
 124 /*
 125  *      define a notation for a function 'R' :
 126  *
 127  *
 128  *
 129  *                                    1
 130  *      R (a0, b0)  =  ------------------------------
 131  *
 132  *                          ----
 133  *                           \
 134  *                            \     n-1      2    2
 135  *                      1  -   |   2    *  (a  - b )
 136  *                            /              n    n
 137  *                           /
 138  *                          ----
 139  *                         n >= 0
 140  *
 141  *
 142  *      where a, b are the classic AGM iteration :
 143  *
 144  *
 145  *      a    =  0.5 * (a  + b )
 146  *       n+1            n    n
 147  *
 148  *
 149  *      b    =  sqrt(a  * b )
 150  *       n+1          n    n
 151  *
 152  *
 153  *
 154  *      define a variable 'c' for more efficient computation :
 155  *
 156  *                                      2     2     2
 157  *      c    =  0.5 * (a  - b )    ,   c  =  a  -  b
 158  *       n+1            n    n          n     n     n
 159  *
 160  */
 161
 162 /****************************************************************************/
 163 void    M_log_AGM_R_func(M_APM rr, int places, M_APM aa, M_APM bb)
 164 {
 165 M_APM   tmp1, tmp2, tmp3, tmp4, tmpC2, sum, pow_2, tmpA0, tmpB0;
 166 int     tolerance, dplaces;
 167
 168 tmpA0 = M_get_stack_var();
 169 tmpB0 = M_get_stack_var();
 170 tmpC2 = M_get_stack_var();
 171 tmp1  = M_get_stack_var();
 172 tmp2  = M_get_stack_var();
 173 tmp3  = M_get_stack_var();
 174 tmp4  = M_get_stack_var();
 175 sum   = M_get_stack_var();
 176 pow_2 = M_get_stack_var();
 177
 178 tolerance = places + 8;
 179 dplaces   = places + 16;
 180
 181 m_apm_copy(tmpA0, aa);
 182 m_apm_copy(tmpB0, bb);
 183 m_apm_copy(pow_2, MM_0_5);
 184
 185 m_apm_multiply(tmp1, aa, aa);               /* 0.5 * [ a ^ 2 - b ^ 2 ] */
 186 m_apm_multiply(tmp2, bb, bb);
 187 m_apm_subtract(tmp3, tmp1, tmp2);
 188 m_apm_multiply(sum, MM_0_5, tmp3);
 189
 190 while (TRUE)
 191   {
 192    m_apm_subtract(tmp1, tmpA0, tmpB0);      /* C n+1 = 0.5 * [ An - Bn ] */
 193    m_apm_multiply(tmp4, MM_0_5, tmp1);      /* C n+1 */
 194    m_apm_multiply(tmpC2, tmp4, tmp4);       /* C n+1 ^ 2 */
 195
 196    /* do the AGM */
 197
 198    m_apm_add(tmp1, tmpA0, tmpB0);
 199    m_apm_multiply(tmp3, MM_0_5, tmp1);
 200
 201    m_apm_multiply(tmp2, tmpA0, tmpB0);
 202    m_apm_sqrt(tmpB0, dplaces, tmp2);
 203
 204    m_apm_round(tmpA0, dplaces, tmp3);
 205
 206    /* end AGM */
 207
 208    m_apm_multiply(tmp2, MM_Two, pow_2);
 209    m_apm_copy(pow_2, tmp2);
 210
 211    m_apm_multiply(tmp1, tmpC2, pow_2);
 212    m_apm_add(tmp3, sum, tmp1);
 213
 214    if ((tmp1->m_apm_sign == 0) ||
 215       ((-2 * tmp1->m_apm_exponent) > tolerance))
 216      break;
 217
 218    m_apm_round(sum, dplaces, tmp3);
 219   }
 220
 221 m_apm_subtract(tmp4, MM_One, tmp3);
 222 m_apm_reciprocal(rr, places, tmp4);
 223
 224 M_restore_stack(9);
 225 }
 226 /****************************************************************************/