src/libs/mapm/mapmrsin.c

   1
   2 /*
   3  *  M_APM  -  mapmrsin.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: mapmrsin.c,v 1.8 2007/12/03 01:57:00 mike Exp $
  23  *
  24  *      This file contains the basic series expansion functions for
  25  *      the SIN / COS functions.
  26  *
  27  *      $Log: mapmrsin.c,v $
  28  *      Revision 1.8  2007/12/03 01:57:00  mike
  29  *      Update license
  30  *
  31  *      Revision 1.7  2003/06/08 18:22:09  mike
  32  *      optimize the raw sin algorithm
  33  *
  34  *      Revision 1.6  2002/11/03 21:58:27  mike
  35  *      Updated function parameters to use the modern style
  36  *
  37  *      Revision 1.5  2001/07/10 22:14:43  mike
  38  *      optimize raw_sin & cos by using fewer digits
  39  *      as subsequent terms get smaller
  40  *
  41  *      Revision 1.4  2000/03/30 21:53:48  mike
  42  *      change compare to terminate series expansion using ints instead
  43  *      of MAPM numbers
  44  *
  45  *      Revision 1.3  1999/06/20 16:23:10  mike
  46  *      changed local static variables to MAPM stack variables
  47  *
  48  *      Revision 1.2  1999/05/12 21:06:36  mike
  49  *      changed global var names
  50  *
  51  *      Revision 1.1  1999/05/10 20:56:31  mike
  52  *      Initial revision
  53  */
  54
  55 #include "m_apm_lc.h"
  56
  57 /****************************************************************************/
  58 /*
  59                                   x^3     x^5     x^7     x^9
  60                  sin(x)  =  x  -  ---  +  ---  -  ---  +  ---  ...
  61                                    3!      5!      7!      9!
  62 */
  63 void    M_raw_sin(M_APM rr, int places, M_APM xx)
  64 {
  65 M_APM   sum, term, tmp2, tmp7, tmp8;
  66 int     tolerance, flag, local_precision, dplaces;
  67 long    m1, m2;
  68
  69 sum  = M_get_stack_var();
  70 term = M_get_stack_var();
  71 tmp2 = M_get_stack_var();
  72 tmp7 = M_get_stack_var();
  73 tmp8 = M_get_stack_var();
  74
  75 m_apm_copy(sum, xx);
  76 m_apm_copy(term, xx);
  77 m_apm_multiply(tmp8, xx, xx);
  78 m_apm_round(tmp2, (places + 6), tmp8);
  79
  80 dplaces   = (places + 8) - xx->m_apm_exponent;
  81 tolerance = xx->m_apm_exponent - (places + 4);
  82
  83 m1   = 2L;
  84 flag = 0;
  85
  86 while (TRUE)
  87   {
  88    m_apm_multiply(tmp8, term, tmp2);
  89
  90    if ((tmp8->m_apm_exponent < tolerance) || (tmp8->m_apm_sign == 0))
  91      break;
  92
  93    local_precision = dplaces + term->m_apm_exponent;
  94
  95    if (local_precision < 20)
  96      local_precision = 20;
  97
  98    m2 = m1 * (m1 + 1);
  99    m_apm_set_long(tmp7, m2);
 100
 101    m_apm_divide(term, local_precision, tmp8, tmp7);
 102
 103    if (flag == 0)
 104      {
 105       m_apm_subtract(tmp7, sum, term);
 106       m_apm_copy(sum, tmp7);
 107      }
 108    else
 109      {
 110       m_apm_add(tmp7, sum, term);
 111       m_apm_copy(sum, tmp7);
 112      }
 113
 114    m1 += 2;
 115    flag = 1 - flag;
 116   }
 117
 118 m_apm_round(rr, places, sum);
 119 M_restore_stack(5);
 120 }
 121 /****************************************************************************/
 122 /*
 123                                   x^2     x^4     x^6     x^8
 124                  cos(x)  =  1  -  ---  +  ---  -  ---  +  ---  ...
 125                                    2!      4!      6!      8!
 126 */
 127 void    M_raw_cos(M_APM rr, int places, M_APM xx)
 128 {
 129 M_APM   sum, term, tmp7, tmp8, tmp9;
 130 int     tolerance, flag, local_precision, prev_exp;
 131 long    m1, m2;
 132
 133 sum  = M_get_stack_var();
 134 term = M_get_stack_var();
 135 tmp7 = M_get_stack_var();
 136 tmp8 = M_get_stack_var();
 137 tmp9 = M_get_stack_var();
 138
 139 m_apm_copy(sum, MM_One);
 140 m_apm_copy(term, MM_One);
 141
 142 m_apm_multiply(tmp8, xx, xx);
 143 m_apm_round(tmp9, (places + 6), tmp8);
 144
 145 local_precision = places + 8;
 146 tolerance       = -(places + 4);
 147 prev_exp        = 0;
 148
 149 m1   = 1L;
 150 flag = 0;
 151
 152 while (TRUE)
 153   {
 154    m2 = m1 * (m1 + 1);
 155    m_apm_set_long(tmp7, m2);
 156
 157    m_apm_multiply(tmp8, term, tmp9);
 158    m_apm_divide(term, local_precision, tmp8, tmp7);
 159
 160    if (flag == 0)
 161      {
 162       m_apm_subtract(tmp7, sum, term);
 163       m_apm_copy(sum, tmp7);
 164      }
 165    else
 166      {
 167       m_apm_add(tmp7, sum, term);
 168       m_apm_copy(sum, tmp7);
 169      }
 170
 171    if ((term->m_apm_exponent < tolerance) || (term->m_apm_sign == 0))
 172      break;
 173
 174    if (m1 != 1L)
 175      {
 176       local_precision = local_precision + term->m_apm_exponent - prev_exp;
 177
 178       if (local_precision < 20)
 179         local_precision = 20;
 180      }
 181
 182    prev_exp = term->m_apm_exponent;
 183
 184    m1 += 2;
 185    flag = 1 - flag;
 186   }
 187
 188 m_apm_round(rr, places, sum);
 189 M_restore_stack(5);
 190 }
 191 /****************************************************************************/