src/libs/mapm/mapmgues.c

   1
   2 /*
   3  *  M_APM  -  mapmgues.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: mapmgues.c,v 1.20 2007/12/03 01:52:55 mike Exp $
  23  *
  24  *      This file contains the functions that generate the initial
  25  *      'guesses' for the sqrt, cbrt, log, arcsin, and arccos functions.
  26  *
  27  *      $Log: mapmgues.c,v $
  28  *      Revision 1.20  2007/12/03 01:52:55  mike
  29  *      Update license
  30  *
  31  *      Revision 1.19  2003/07/21 20:03:49  mike
  32  *      check for invalid inputs to _set_double
  33  *
  34  *      Revision 1.18  2003/05/01 21:58:45  mike
  35  *      remove math.h
  36  *
  37  *      Revision 1.17  2003/04/16 16:52:47  mike
  38  *      change cbrt guess to use reciprocal value for new cbrt algorithm
  39  *
  40  *      Revision 1.16  2003/04/11 14:18:13  mike
  41  *      add comments
  42  *
  43  *      Revision 1.15  2003/04/10 22:28:35  mike
  44  *      .
  45  *
  46  *      Revision 1.14  2003/04/09 21:33:19  mike
  47  *      induce same error for asin and acos
  48  *
  49  *      Revision 1.13  2003/04/09 20:11:38  mike
  50  *      induce error of 10 ^ -5 in log guess for known starting
  51  *      point in the iterative algorithm
  52  *
  53  *      Revision 1.12  2003/03/27 19:32:59  mike
  54  *      simplify log guess since caller guarantee's limited input magnitude
  55  *
  56  *      Revision 1.11  2002/11/03 21:45:53  mike
  57  *      Updated function parameters to use the modern style
  58  *
  59  *      Revision 1.10  2001/03/20 22:08:27  mike
  60  *      delete unneeded logic in asin guess
  61  *
  62  *      Revision 1.9  2000/09/26 17:05:11  mike
  63  *      guess 1/sqrt instead of sqrt due to new sqrt algorithm
  64  *
  65  *      Revision 1.8  2000/04/10 21:13:13  mike
  66  *      minor tweaks to log_guess
  67  *
  68  *      Revision 1.7  2000/04/03 17:25:45  mike
  69  *      added function to estimate the cube root (cbrt)
  70  *
  71  *      Revision 1.6  1999/07/18 01:57:35  mike
  72  *      adjust arc-sin guess for small exponents
  73  *
  74  *      Revision 1.5  1999/07/09 22:32:50  mike
  75  *      optimize some functions
  76  *
  77  *      Revision 1.4  1999/05/12 21:22:27  mike
  78  *      add more comments
  79  *
  80  *      Revision 1.3  1999/05/12 21:00:51  mike
  81  *      added new sqrt guess function
  82  *
  83  *      Revision 1.2  1999/05/11 02:10:12  mike
  84  *      added some comments
  85  *
  86  *      Revision 1.1  1999/05/10 20:56:31  mike
  87  *      Initial revision
  88  */
  89
  90 #include "m_apm_lc.h"
  91
  92 /****************************************************************************/
  93 void    M_get_sqrt_guess(M_APM r, M_APM a)
  94 {
  95 char    buf[48];
  96 double  dd;
  97
  98 m_apm_to_string(buf, 15, a);
  99 dd = atof(buf);                     /* sqrt algorithm actually finds 1/sqrt */
 100 m_apm_set_double(r, (1.0 / sqrt(dd)));
 101 }
 102 /****************************************************************************/
 103 /*
 104  *      for cbrt, log, asin, and acos we induce an error of 10 ^ -5.
 105  *      this enables the iterative routine to be more efficient
 106  *      by knowing exactly how accurate the initial guess is.
 107  *
 108  *      this also prevents some corner conditions where the iterative
 109  *      functions may terminate too soon.
 110  */
 111 /****************************************************************************/
 112 void    M_get_cbrt_guess(M_APM r, M_APM a)
 113 {
 114 char    buf[48];
 115 double  dd;
 116
 117 m_apm_to_string(buf, 15, a);
 118 dd = atof(buf);
 119 dd = log(dd) / 3.0;                 /* cbrt algorithm actually finds 1/cbrt */
 120 m_apm_set_double(r, (1.00001 / exp(dd)));
 121 }
 122 /****************************************************************************/
 123 void    M_get_log_guess(M_APM r, M_APM a)
 124 {
 125 char    buf[48];
 126 double  dd;
 127
 128 m_apm_to_string(buf, 15, a);
 129 dd = atof(buf);
 130 m_apm_set_double(r, (1.00001 * log(dd)));        /* induce error of 10 ^ -5 */
 131 }
 132 /****************************************************************************/
 133 /*
 134  *      the implementation of the asin & acos functions
 135  *      guarantee that 'a' is always < 0.85, so it is
 136  *      safe to multiply by a number > 1
 137  */
 138 void    M_get_asin_guess(M_APM r, M_APM a)
 139 {
 140 char    buf[48];
 141 double  dd;
 142
 143 m_apm_to_string(buf, 15, a);
 144 dd = atof(buf);
 145 m_apm_set_double(r, (1.00001 * asin(dd)));       /* induce error of 10 ^ -5 */
 146 }
 147 /****************************************************************************/
 148 void    M_get_acos_guess(M_APM r, M_APM a)
 149 {
 150 char    buf[48];
 151 double  dd;
 152
 153 m_apm_to_string(buf, 15, a);
 154 dd = atof(buf);
 155 m_apm_set_double(r, (1.00001 * acos(dd)));       /* induce error of 10 ^ -5 */
 156 }
 157 /****************************************************************************/
 158 /*
 159         convert a C 'double' into an M_APM value.
 160 */
 161 void    m_apm_set_double(M_APM atmp, double dd)
 162 {
 163 char    *cp, *p, *ps, buf[64];
 164
 165 if (dd == 0.0)                     /* special case for 0 exactly */
 166    M_set_to_zero(atmp);
 167 else
 168   {
 169    sprintf(buf, "%.14E", dd);
 170
 171    if ((cp = strstr(buf, "E")) == NULL)
 172      {
 173       M_apm_log_error_msg(M_APM_RETURN,
 174       "\'m_apm_set_double\', Invalid double input (likely a NAN or +/- INF)");
 175
 176       M_set_to_zero(atmp);
 177       return;
 178      }
 179
 180    if (atoi(cp + sizeof(char)) == 0)
 181      *cp = '\0';
 182
 183    p = cp;
 184
 185    while (TRUE)
 186      {
 187       p--;
 188       if (*p == '0' || *p == '.')
 189         *p = ' ';
 190       else
 191         break;
 192      }
 193
 194    ps = buf;
 195    p  = buf;
 196
 197    while (TRUE)
 198      {
 199       if ((*p = *ps) == '\0')
 200         break;
 201
 202       if (*ps++ != ' ')
 203         p++;
 204      }
 205
 206    m_apm_set_string(atmp, buf);
 207   }
 208 }
 209 /****************************************************************************/