src/libs/mapm/mapmasn0.c

   1
   2 /*
   3  *  M_APM  -  mapmasn0.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: mapmasn0.c,v 1.8 2007/12/03 01:49:49 mike Exp $
  23  *
  24  *      This file contains the 'ARC' family of functions; ARC-SIN,
  25  *      ARC-COS, ARC-TAN when the input arg is very close to 0 (zero).
  26  *
  27  *      $Log: mapmasn0.c,v $
  28  *      Revision 1.8  2007/12/03 01:49:49  mike
  29  *      Update license
  30  *
  31  *      Revision 1.7  2003/06/02 16:51:13  mike
  32  *      *** empty log message ***
  33  *
  34  *      Revision 1.6  2003/06/02 16:49:48  mike
  35  *      tweak the decimal places
  36  *
  37  *      Revision 1.5  2003/06/02 16:47:39  mike
  38  *      tweak arctan algorithm some more
  39  *
  40  *      Revision 1.4  2003/05/31 22:38:07  mike
  41  *      optimize arctan by using fewer digits as subsequent
  42  *      terms get smaller
  43  *
  44  *      Revision 1.3  2002/11/03 21:36:43  mike
  45  *      Updated function parameters to use the modern style
  46  *
  47  *      Revision 1.2  2000/12/02 20:11:37  mike
  48  *      add comments
  49  *
  50  *      Revision 1.1  2000/12/02 20:08:27  mike
  51  *      Initial revision
  52  */
  53
  54 #include "m_apm_lc.h"
  55
  56 /****************************************************************************/
  57 /*
  58         Calculate arcsin using the identity :
  59
  60                                       x
  61         arcsin (x) == arctan [ --------------- ]
  62                                 sqrt(1 - x^2)
  63
  64 */
  65 void    M_arcsin_near_0(M_APM rr, int places, M_APM aa)
  66 {
  67 M_APM   tmp5, tmp6;
  68
  69 tmp5 = M_get_stack_var();
  70 tmp6 = M_get_stack_var();
  71
  72 M_cos_to_sin(tmp5, (places + 8), aa);
  73 m_apm_divide(tmp6, (places + 8), aa, tmp5);
  74 M_arctan_near_0(rr, places, tmp6);
  75
  76 M_restore_stack(2);
  77 }
  78 /****************************************************************************/
  79 /*
  80         Calculate arccos using the identity :
  81
  82         arccos (x) == PI / 2 - arcsin (x)
  83
  84 */
  85 void    M_arccos_near_0(M_APM rr, int places, M_APM aa)
  86 {
  87 M_APM   tmp1, tmp2;
  88
  89 tmp1 = M_get_stack_var();
  90 tmp2 = M_get_stack_var();
  91
  92 M_check_PI_places(places);
  93 M_arcsin_near_0(tmp1, (places + 4), aa);
  94 m_apm_subtract(tmp2, MM_lc_HALF_PI, tmp1);
  95 m_apm_round(rr, places, tmp2);
  96
  97 M_restore_stack(2);
  98 }
  99 /****************************************************************************/
 100 /*
 101         calculate arctan (x) with the following series:
 102
 103                               x^3     x^5     x^7     x^9
 104         arctan (x)  =   x  -  ---  +  ---  -  ---  +  ---  ...
 105                                3       5       7       9
 106
 107 */
 108 void    M_arctan_near_0(M_APM rr, int places, M_APM aa)
 109 {
 110 M_APM   tmp0, tmp2, tmpR, tmpS, digit, term;
 111 int     tolerance, dplaces, local_precision;
 112 long    m1;
 113
 114 tmp0  = M_get_stack_var();
 115 tmp2  = M_get_stack_var();
 116 tmpR  = M_get_stack_var();
 117 tmpS  = M_get_stack_var();
 118 term  = M_get_stack_var();
 119 digit = M_get_stack_var();
 120
 121 tolerance = aa->m_apm_exponent - (places + 4);
 122 dplaces   = (places + 8) - aa->m_apm_exponent;
 123
 124 m_apm_copy(term, aa);
 125 m_apm_copy(tmpS, aa);
 126 m_apm_multiply(tmp0, aa, aa);
 127 m_apm_round(tmp2, (dplaces + 8), tmp0);
 128
 129 m1 = 1L;
 130
 131 while (TRUE)
 132   {
 133    /*
 134     *   do the subtraction term
 135     */
 136
 137    m_apm_multiply(tmp0, term, tmp2);
 138
 139    if ((tmp0->m_apm_exponent < tolerance) || (tmp0->m_apm_sign == 0))
 140      {
 141       m_apm_round(rr, places, tmpS);
 142       break;
 143      }
 144
 145    local_precision = dplaces + tmp0->m_apm_exponent;
 146
 147    if (local_precision < 20)
 148      local_precision = 20;
 149
 150    m1 += 2;
 151    m_apm_set_long(digit, m1);
 152    m_apm_round(term, local_precision, tmp0);
 153    m_apm_divide(tmp0, local_precision, term, digit);
 154    m_apm_subtract(tmpR, tmpS, tmp0);
 155
 156    /*
 157     *   do the addition term
 158     */
 159
 160    m_apm_multiply(tmp0, term, tmp2);
 161
 162    if ((tmp0->m_apm_exponent < tolerance) || (tmp0->m_apm_sign == 0))
 163      {
 164       m_apm_round(rr, places, tmpR);
 165       break;
 166      }
 167
 168    local_precision = dplaces + tmp0->m_apm_exponent;
 169
 170    if (local_precision < 20)
 171      local_precision = 20;
 172
 173    m1 += 2;
 174    m_apm_set_long(digit, m1);
 175    m_apm_round(term, local_precision, tmp0);
 176    m_apm_divide(tmp0, local_precision, term, digit);
 177    m_apm_add(tmpS, tmpR, tmp0);
 178   }
 179
 180 M_restore_stack(6);                    /* restore the 6 locals we used here */
 181 }
 182 /****************************************************************************/