src/libs/mapm/mapm_pow.c

   1
   2 /*
   3  *  M_APM  -  mapm_pow.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: mapm_pow.c,v 1.10 2007/12/03 01:46:07 mike Exp $
  23  *
  24  *      This file contains the POW function.
  25  *
  26  *      $Log: mapm_pow.c,v $
  27  *      Revision 1.10  2007/12/03 01:46:07  mike
  28  *      Update license
  29  *
  30  *      Revision 1.9  2002/11/05 23:39:42  mike
  31  *      use new set_to_zero call
  32  *
  33  *      Revision 1.8  2002/11/03 22:20:59  mike
  34  *      Updated function parameters to use the modern style
  35  *
  36  *      Revision 1.7  2001/07/16 19:24:26  mike
  37  *      add function M_free_all_pow
  38  *
  39  *      Revision 1.6  2000/09/05 22:15:03  mike
  40  *      minor tweak
  41  *
  42  *      Revision 1.5  2000/08/22 21:22:29  mike
  43  *      if parameter yy is an integer, call the more
  44  *      efficient _integer_pow function
  45  *
  46  *      Revision 1.4  2000/08/22 20:42:08  mike
  47  *      compute more digits in the log calculation
  48  *
  49  *      Revision 1.3  2000/05/24 20:08:21  mike
  50  *      update some comments
  51  *
  52  *      Revision 1.2  2000/05/23 23:20:11  mike
  53  *      return 1 when input is 0^0.
  54  *
  55  *      Revision 1.1  2000/05/18 22:10:43  mike
  56  *      Initial revision
  57  */
  58
  59 #include "m_apm_lc.h"
  60
  61 static  M_APM   M_last_xx_input;
  62 static  M_APM   M_last_xx_log;
  63 static  int     M_last_log_digits;
  64 static  int     M_size_flag = 0;
  65
  66 /****************************************************************************/
  67 void    M_free_all_pow()
  68 {
  69 if (M_size_flag != 0)
  70   {
  71    m_apm_free(M_last_xx_input);
  72    m_apm_free(M_last_xx_log);
  73    M_size_flag = 0;
  74   }
  75 }
  76 /****************************************************************************/
  77 /*
  78         Calculate the POW function by calling EXP :
  79
  80                   Y      A
  81                  X   =  e    where A = Y * log(X)
  82 */
  83 void    m_apm_pow(M_APM rr, int places, M_APM xx, M_APM yy)
  84 {
  85 int     iflag, pflag;
  86 char    sbuf[64];
  87 M_APM   tmp8, tmp9;
  88
  89 /* if yy == 0, return 1 */
  90
  91 if (yy->m_apm_sign == 0)
  92   {
  93    m_apm_copy(rr, MM_One);
  94    return;
  95   }
  96
  97 /* if xx == 0, return 0 */
  98
  99 if (xx->m_apm_sign == 0)
 100   {
 101    M_set_to_zero(rr);
 102    return;
 103   }
 104
 105 if (M_size_flag == 0)       /* init locals on first call */
 106   {
 107    M_size_flag       = M_get_sizeof_int();
 108    M_last_log_digits = 0;
 109    M_last_xx_input   = m_apm_init();
 110    M_last_xx_log     = m_apm_init();
 111   }
 112
 113 /*
 114  *  if 'yy' is a small enough integer, call the more
 115  *  efficient _integer_pow function.
 116  */
 117
 118 if (m_apm_is_integer(yy))
 119   {
 120    iflag = FALSE;
 121
 122    if (M_size_flag == 2)            /* 16 bit compilers */
 123      {
 124       if (yy->m_apm_exponent <= 4)
 125         iflag = TRUE;
 126      }
 127    else                             /* >= 32 bit compilers */
 128      {
 129       if (yy->m_apm_exponent <= 7)
 130         iflag = TRUE;
 131      }
 132
 133    if (iflag)
 134      {
 135       m_apm_to_integer_string(sbuf, yy);
 136       m_apm_integer_pow(rr, places, xx, atoi(sbuf));
 137       return;
 138      }
 139   }
 140
 141 tmp8 = M_get_stack_var();
 142 tmp9 = M_get_stack_var();
 143
 144 /*
 145  *    If parameter 'X' is the same this call as it
 146  *    was the previous call, re-use the saved log
 147  *    calculation from last time.
 148  */
 149
 150 pflag = FALSE;
 151
 152 if (M_last_log_digits >= places)
 153   {
 154    if (m_apm_compare(xx, M_last_xx_input) == 0)
 155      pflag = TRUE;
 156   }
 157
 158 if (pflag)
 159   {
 160    m_apm_round(tmp9, (places + 8), M_last_xx_log);
 161   }
 162 else
 163   {
 164    m_apm_log(tmp9, (places + 8), xx);
 165
 166    M_last_log_digits = places + 2;
 167
 168    /* save the 'X' input value and the log calculation */
 169
 170    m_apm_copy(M_last_xx_input, xx);
 171    m_apm_copy(M_last_xx_log, tmp9);
 172   }
 173
 174 m_apm_multiply(tmp8, tmp9, yy);
 175 m_apm_exp(rr, places, tmp8);
 176 M_restore_stack(2);                    /* restore the 2 locals we used here */
 177 }
 178 /****************************************************************************/