src/libs/mapm/mapmcbrt.c

   1
   2 /*
   3  *  M_APM  -  mapmcbrt.c
   4  *
   5  *  Copyright (C) 2000 - 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: mapmcbrt.c,v 1.8 2007/12/03 01:50:32 mike Exp $
  23  *
  24  *      This file contains the CBRT (cube root) function.
  25  *
  26  *      $Log: mapmcbrt.c,v $
  27  *      Revision 1.8  2007/12/03 01:50:32  mike
  28  *      Update license
  29  *
  30  *      Revision 1.7  2003/05/05 18:17:38  mike
  31  *      simplify some logic
  32  *
  33  *      Revision 1.6  2003/04/16 16:55:58  mike
  34  *      use new faster algorithm which finds 1 / cbrt(N)
  35  *
  36  *      Revision 1.5  2002/11/03 21:34:34  mike
  37  *      Updated function parameters to use the modern style
  38  *
  39  *      Revision 1.4  2000/10/30 16:42:22  mike
  40  *      minor speed optimization
  41  *
  42  *      Revision 1.3  2000/07/11 18:03:39  mike
  43  *      make better estimate for initial precision
  44  *
  45  *      Revision 1.2  2000/04/08 18:34:35  mike
  46  *      added some more comments
  47  *
  48  *      Revision 1.1  2000/04/03 17:58:04  mike
  49  *      Initial revision
  50  */
  51
  52 #include "m_apm_lc.h"
  53
  54 /****************************************************************************/
  55 void    m_apm_cbrt(M_APM rr, int places, M_APM aa)
  56 {
  57 M_APM   last_x, guess, tmpN, tmp7, tmp8, tmp9;
  58 int     ii, nexp, bflag, tolerance, maxp, local_precision;
  59
  60 /* result is 0 if input is 0 */
  61
  62 if (aa->m_apm_sign == 0)
  63   {
  64    M_set_to_zero(rr);
  65    return;
  66   }
  67
  68 last_x = M_get_stack_var();
  69 guess  = M_get_stack_var();
  70 tmpN   = M_get_stack_var();
  71 tmp7   = M_get_stack_var();
  72 tmp8   = M_get_stack_var();
  73 tmp9   = M_get_stack_var();
  74
  75 /* compute the cube root of the positive number, we'll fix the sign later */
  76
  77 m_apm_absolute_value(tmpN, aa);
  78
  79 /*
  80     normalize the input number (make the exponent near 0) so
  81     the 'guess' function will not over/under flow on large
  82     magnitude exponents.
  83 */
  84
  85 nexp = aa->m_apm_exponent / 3;
  86 tmpN->m_apm_exponent -= 3 * nexp;
  87
  88 M_get_cbrt_guess(guess, tmpN);
  89
  90 tolerance       = places + 4;
  91 maxp            = places + 16;
  92 bflag           = FALSE;
  93 local_precision = 14;
  94
  95 m_apm_negate(last_x, MM_Ten);
  96
  97 /*   Use the following iteration to calculate 1 / cbrt(N) :
  98
  99                                  4
 100          X     =  [ 4 * X - N * X ] / 3
 101           n+1
 102 */
 103
 104 ii = 0;
 105
 106 while (TRUE)
 107   {
 108    m_apm_multiply(tmp8, guess, guess);
 109    m_apm_multiply(tmp7, tmp8, tmp8);
 110    m_apm_round(tmp8, local_precision, tmp7);
 111    m_apm_multiply(tmp9, tmpN, tmp8);
 112
 113    m_apm_multiply(tmp8, MM_Four, guess);
 114    m_apm_subtract(tmp7, tmp8, tmp9);
 115    m_apm_divide(guess, local_precision, tmp7, MM_Three);
 116
 117    if (bflag)
 118      break;
 119
 120    /* force at least 2 iterations so 'last_x' has valid data */
 121
 122    if (ii != 0)
 123      {
 124       m_apm_subtract(tmp8, guess, last_x);
 125
 126       if (tmp8->m_apm_sign == 0)
 127         break;
 128
 129       if ((-4 * tmp8->m_apm_exponent) > tolerance)
 130         bflag = TRUE;
 131      }
 132
 133    local_precision *= 2;
 134
 135    if (local_precision > maxp)
 136      local_precision = maxp;
 137
 138    m_apm_copy(last_x, guess);
 139    ii = 1;
 140   }
 141
 142 /* final cbrt = N * guess ^ 2 */
 143
 144 m_apm_multiply(tmp9, guess, guess);
 145 m_apm_multiply(tmp8, tmp9, tmpN);
 146 m_apm_round(rr, places, tmp8);
 147
 148 rr->m_apm_exponent += nexp;
 149 rr->m_apm_sign = aa->m_apm_sign;
 150 M_restore_stack(6);
 151 }
 152 /****************************************************************************/
 153