src/libs/mapm/mapmsqrt.c

   1
   2 /*
   3  *  M_APM  -  mapmsqrt.c
   4  *
   5  *  Copyright (C) 1999 - 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: mapmsqrt.c,v 1.19 2007/12/03 01:57:31 mike Exp $
  23  *
  24  *      This file contains the SQRT function.
  25  *
  26  *      $Log: mapmsqrt.c,v $
  27  *      Revision 1.19  2007/12/03 01:57:31  mike
  28  *      Update license
  29  *
  30  *      Revision 1.18  2003/07/21 20:39:00  mike
  31  *      Modify error messages to be in a consistent format.
  32  *
  33  *      Revision 1.17  2003/05/07 16:36:04  mike
  34  *      simplify 'nexp' logic
  35  *
  36  *      Revision 1.16  2003/03/31 21:50:14  mike
  37  *      call generic error handling function
  38  *
  39  *      Revision 1.15  2003/03/11 21:29:00  mike
  40  *      round an intermediate result for faster runtime.
  41  *
  42  *      Revision 1.14  2002/11/03 22:00:46  mike
  43  *      Updated function parameters to use the modern style
  44  *
  45  *      Revision 1.13  2001/07/10 22:50:31  mike
  46  *      minor optimization
  47  *
  48  *      Revision 1.12  2000/09/26 18:32:04  mike
  49  *      use new algorithm which only uses multiply and subtract
  50  *      (avoids the slower version which used division)
  51  *
  52  *      Revision 1.11  2000/07/11 17:56:22  mike
  53  *      make better estimate for initial precision
  54  *
  55  *      Revision 1.10  1999/07/21 02:48:45  mike
  56  *      added some comments
  57  *
  58  *      Revision 1.9  1999/07/19 00:25:44  mike
  59  *      adjust local precision again
  60  *
  61  *      Revision 1.8  1999/07/19 00:09:41  mike
  62  *      adjust local precision during loop
  63  *
  64  *      Revision 1.7  1999/07/18 22:57:08  mike
  65  *      change to dynamically changing local precision and
  66  *      change tolerance checks using integers
  67  *
  68  *      Revision 1.6  1999/06/19 21:18:00  mike
  69  *      changed local static variables to MAPM stack variables
  70  *
  71  *      Revision 1.5  1999/05/31 01:40:39  mike
  72  *      minor update to normalizing the exponent
  73  *
  74  *      Revision 1.4  1999/05/31 01:21:41  mike
  75  *      optimize for large exponents
  76  *
  77  *      Revision 1.3  1999/05/12 20:59:35  mike
  78  *      use a better 'guess' function
  79  *
  80  *      Revision 1.2  1999/05/10 21:15:26  mike
  81  *      added some comments
  82  *
  83  *      Revision 1.1  1999/05/10 20:56:31  mike
  84  *      Initial revision
  85  */
  86
  87 #include "m_apm_lc.h"
  88
  89 /****************************************************************************/
  90 void    m_apm_sqrt(M_APM rr, int places, M_APM aa)
  91 {
  92 M_APM   last_x, guess, tmpN, tmp7, tmp8, tmp9;
  93 int     ii, bflag, nexp, tolerance, dplaces;
  94
  95 if (aa->m_apm_sign <= 0)
  96   {
  97    if (aa->m_apm_sign == -1)
  98      {
  99       M_apm_log_error_msg(M_APM_RETURN, "\'m_apm_sqrt\', Negative argument");
 100      }
 101
 102    M_set_to_zero(rr);
 103    return;
 104   }
 105
 106 last_x = M_get_stack_var();
 107 guess  = M_get_stack_var();
 108 tmpN   = M_get_stack_var();
 109 tmp7   = M_get_stack_var();
 110 tmp8   = M_get_stack_var();
 111 tmp9   = M_get_stack_var();
 112
 113 m_apm_copy(tmpN, aa);
 114
 115 /*
 116     normalize the input number (make the exponent near 0) so
 117     the 'guess' function will not over/under flow on large
 118     magnitude exponents.
 119 */
 120
 121 nexp = aa->m_apm_exponent / 2;
 122 tmpN->m_apm_exponent -= 2 * nexp;
 123
 124 M_get_sqrt_guess(guess, tmpN);    /* actually gets 1/sqrt guess */
 125
 126 tolerance = places + 4;
 127 dplaces   = places + 16;
 128 bflag     = FALSE;
 129
 130 m_apm_negate(last_x, MM_Ten);
 131
 132 /*   Use the following iteration to calculate 1 / sqrt(N) :
 133
 134          X    =  0.5 * X * [ 3 - N * X^2 ]
 135           n+1
 136 */
 137
 138 ii = 0;
 139
 140 while (TRUE)
 141   {
 142    m_apm_multiply(tmp9, tmpN, guess);
 143    m_apm_multiply(tmp8, tmp9, guess);
 144    m_apm_round(tmp7, dplaces, tmp8);
 145    m_apm_subtract(tmp9, MM_Three, tmp7);
 146    m_apm_multiply(tmp8, tmp9, guess);
 147    m_apm_multiply(tmp9, tmp8, MM_0_5);
 148
 149    if (bflag)
 150      break;
 151
 152    m_apm_round(guess, dplaces, tmp9);
 153
 154    /* force at least 2 iterations so 'last_x' has valid data */
 155
 156    if (ii != 0)
 157      {
 158       m_apm_subtract(tmp7, guess, last_x);
 159
 160       if (tmp7->m_apm_sign == 0)
 161         break;
 162
 163       /*
 164        *   if we are within a factor of 4 on the error term,
 165        *   we will be accurate enough after the *next* iteration
 166        *   is complete.  (note that the sign of the exponent on
 167        *   the error term will be a negative number).
 168        */
 169
 170       if ((-4 * tmp7->m_apm_exponent) > tolerance)
 171         bflag = TRUE;
 172      }
 173
 174    m_apm_copy(last_x, guess);
 175    ii++;
 176   }
 177
 178 /*
 179  *    multiply by the starting number to get the final
 180  *    sqrt and then adjust the exponent since we found
 181  *    the sqrt of the normalized number.
 182  */
 183
 184 m_apm_multiply(tmp8, tmp9, tmpN);
 185 m_apm_round(rr, places, tmp8);
 186 rr->m_apm_exponent += nexp;
 187
 188 M_restore_stack(6);
 189 }
 190 /****************************************************************************/