"Fossies" - the Fresh Open Source Software Archive

Member "mapm_4.9.5a/cpp_demo.cpp" (21 Feb 2010, 15320 Bytes) of package /linux/misc/old/mapm-4.9.5a.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "cpp_demo.cpp" see the Fossies "Dox" file reference documentation.

    1 
    2 /*
    3  *   $Id: cpp_demo.cpp,v 1.20 2007/12/03 01:18:31 mike Exp $
    4  *
    5  *   $Log: cpp_demo.cpp,v $
    6  *   Revision 1.20  2007/12/03 01:18:31  mike
    7  *   no changes
    8  *
    9  *   Revision 1.19  2004/05/31 22:11:27  mike
   10  *   add MOD operator demo
   11  *
   12  *   Revision 1.18  2003/05/14 21:21:43  mike
   13  *   *** empty log message ***
   14  *
   15  *   Revision 1.17  2003/05/06 22:27:32  mike
   16  *   increase array size
   17  *
   18  *   Revision 1.16  2003/05/06 21:56:39  mike
   19  *   add demo and check of library version functions
   20  *
   21  *   Revision 1.15  2002/11/03 23:34:36  mike
   22  *   added ipow_nr demo to show the new function
   23  *
   24  *   Revision 1.14  2001/08/28 18:41:16  mike
   25  *   add another fixed point format
   26  *
   27  *   Revision 1.13  2001/08/27 23:12:07  mike
   28  *   add demo of fixed point formatting functions
   29  *
   30  *   Revision 1.12  2001/08/07 19:49:06  mike
   31  *   add another PI example using atan, log, exp
   32  *
   33  *   Revision 1.11  2001/08/06 20:43:50  mike
   34  *   add some PI examples using asin and acos
   35  *
   36  *   Revision 1.10  2001/07/15 21:21:31  mike
   37  *   add is_odd, is_even examples
   38  *
   39  *   Revision 1.9  2001/07/15 21:14:45  mike
   40  *   add GCD and LCM example
   41  *
   42  *   Revision 1.8  2001/03/25 21:28:37  mike
   43  *   add new floor and ceil calls
   44  *
   45  *   Revision 1.7  2000/08/21 23:42:52  mike
   46  *   add is_integer demo
   47  *
   48  *   Revision 1.6  2000/06/20 21:40:15  mike
   49  *   added local versions of cbrt and inverse hyperbolics since
   50  *   not all C libraries have these functions yet.
   51  *
   52  *   Revision 1.5  2000/04/08 18:53:18  mike
   53  *   added more to Syntax Demo
   54  *
   55  *   Revision 1.4  2000/04/06 23:20:41  mike
   56  *   fix commented out code
   57  *
   58  *   Revision 1.3  2000/04/06 23:19:49  mike
   59  *   added SyntaxDemo function
   60  *
   61  *   Revision 1.2  2000/04/05 21:37:47  mike
   62  *   added a few more minor demos
   63  *
   64  *   Revision 1.1  2000/04/04 19:50:00  mike
   65  *   Initial revision
   66  *
   67  *   A simple file to demonstrate the use of the C++
   68  *   M_APM wrapper class.
   69  *
   70  *   Orion Sky Lawlor, olawlor@acm.org, 4/2/2000
   71  */
   72 
   73 #include <stdio.h>
   74 #include <stdlib.h>
   75 #include <string.h>
   76 #include <math.h>
   77 #include "m_apm.h"
   78 
   79 extern   double   cbrt_local(double);
   80 extern   double   asinh_local(double);
   81 extern   double   acosh_local(double);
   82 extern   double   atanh_local(double);
   83 
   84 
   85 void    library_version_check()
   86 {
   87 char    version_info[80];
   88 
   89 m_apm_lib_version(version_info);
   90 
   91 fprintf(stdout,"\n%s\n\n",version_info);
   92 
   93 /*
   94  *  check MAPM version info from m_apm.h to what was compiled into 
   95  *  the library. This check could be used in your application to
   96  *  confirm that you are using the correct include file (m_apm.h).
   97  */
   98 
   99 if (strcmp(version_info, MAPM_LIB_VERSION) != 0)
  100   {
  101    fprintf(stdout,
  102       "ERROR: MAPM Library compiled with different include file\n");
  103 
  104    fprintf(stdout, "Include file version : %s \n", MAPM_LIB_SHORT_VERSION);
  105    fprintf(stdout, "Library file version : %s \n", 
  106             m_apm_lib_short_version(version_info));
  107    exit(10);
  108   }
  109 }
  110 
  111 
  112 // Compute the factorial of the integer n
  113 
  114 MAPM Factorial(MAPM n)
  115 {
  116     MAPM ii;
  117     MAPM product=1;
  118     for (ii=2; ii <= n; ii++)
  119         product *= ii;
  120     return product;
  121 }
  122 
  123 
  124 void printMAPM(const char *caption,MAPM m)
  125 {
  126     char mBuf[1000];
  127     m.toString(mBuf,13);
  128     printf("%s%s\n",caption,mBuf);
  129 }
  130 
  131 void compare(const char *caption,double dVal,MAPM mVal)
  132 {
  133     char dBuf[100],mBuf[1000];
  134     mVal.toFixPtString(mBuf,13);
  135     printf(caption,dVal,mBuf);
  136     sprintf(dBuf,"%.13f",dVal);
  137     if (0 != strncmp(dBuf,mBuf,14))
  138       {
  139        fprintf(stderr,"ERROR!  %s != %s \n",dBuf,mBuf);
  140        exit(1);
  141       }
  142 }
  143 
  144 
  145 void testMAPMcppRoutines(double v);
  146 void SyntaxDemo(void);
  147 
  148 int main(void)
  149 {
  150     int n=200;
  151     library_version_check();
  152     MAPM c=Factorial(n);             /* use function from above */
  153     printMAPM("Factorial is ",c);
  154     m_apm_cpp_precision(25);         /* default is 30 */
  155     testMAPMcppRoutines(0.3);
  156     testMAPMcppRoutines(-0.6);
  157     testMAPMcppRoutines(1.0);
  158     testMAPMcppRoutines(5.4);
  159     SyntaxDemo();
  160     return 0;
  161 }
  162 
  163 void testMAPMcppRoutines(double x)
  164 {
  165     char fBuf[256];
  166     MAPM m=x;
  167     printf("-------------------------\nFor x=%f:\n",x);
  168 #define fmt "%.13f = %s\n"
  169     compare("x+0       "fmt,x+0,m+0);
  170     compare("x+1       "fmt,x+1,m+1);
  171     compare("3+x       "fmt,3+x,3+m);
  172     compare("2+(-x)    "fmt,2+(-x),2+(-m));
  173     compare("x-3.2     "fmt,x-3.2,m-3.2);
  174     compare("x-3.2 (as string) "fmt,x-3.2,m-"3.2");
  175     compare("3.2+(-x)  "fmt,3.2+(-x),3.2+(-m));
  176     compare("1.5*x     "fmt,1.5*x,1.5*m);
  177     compare("2/x       "fmt,2/x,2/m);
  178     int y=(int)(6+10*x+0.5);
  179     MAPM n=6+10*m;
  180     compare("y=6+10*x  "fmt,y,n);
  181     compare("y div 2   "fmt,y/2,n.div(2));
  182     compare("y mod 2   "fmt,y%2,n.rem(2));
  183     compare("y div -3  "fmt,y/-3,n.div(-3));
  184     compare("y mod -3  "fmt,y%-3,n.rem(-3));
  185     compare("-y div 3  "fmt,(-y)/-3,(-n).div(-3));
  186     compare("-y mod 3  "fmt,(-y)%-3,(-n).rem(-3));
  187     
  188     compare("fabs(x)   "fmt,fabs(x),fabs(m));
  189     compare("sin(x)    "fmt,sin(x),sin(m));
  190     compare("x.sin()   "fmt,sin(x),m.sin());
  191     compare("cos(x)    "fmt,cos(x),cos(m));
  192     compare("tan(x)    "fmt,tan(x),tan(m));
  193     compare("cbrt(x)   "fmt,cbrt_local(x),cbrt(m));
  194     compare("sinh(x)   "fmt,sinh(x),sinh(m));
  195     compare("cosh(x)   "fmt,cosh(x),cosh(m));
  196     compare("tanh(x)   "fmt,tanh(x),tanh(m));
  197     compare("asinh(x)  "fmt,asinh_local(x),asinh(m));
  198 
  199     if (fabs(m) < 1)
  200       {
  201            compare("asin(x)   "fmt,asin(x),asin(m));
  202        compare("acos(x)   "fmt,acos(x),acos(m));
  203        compare("atan(x)   "fmt,atan(x),atan(m));
  204        compare("atan2(x,-0.5)"fmt,atan2(x,-0.5),atan2(m,-0.5));
  205            compare("atanh(x)  "fmt,atanh_local(x),atanh(m));
  206       }
  207 
  208     if (m >= 1)
  209            compare("acosh(x)  "fmt,acosh_local(x),acosh(m));
  210 
  211     if (m > 0)
  212       {
  213        compare("sqrt(x)   "fmt,sqrt(x),sqrt(m));
  214        compare("log(x)    "fmt,log(x),log(m));
  215        compare("pow(x,-2.93)"fmt,pow(x,-2.93),pow(m,-2.93));
  216        compare("x.pow(-2.93)"fmt,pow(x,-2.93),m.pow(-2.93));
  217       }
  218 
  219     compare("x.pow(2)  "fmt,pow(x,2),m.ipow(2));
  220     compare("x.pow(-3) "fmt,pow(x,-3),m.ipow(-3));
  221     compare("exp(x)    "fmt,exp(x),exp(m));
  222 
  223     MAPM r1, r2;
  224     r1 = get_random();
  225     r2 = get_random();
  226     printMAPM("A random MAPM on [0,1) is ",r1);
  227     printMAPM("Another random MAPM on [0,1) is ",r2);
  228 
  229     printMAPM("m=",m);
  230     printMAPM("m++=",m++);
  231     printMAPM("now m is ",m);
  232     printMAPM("++m=",++m);
  233     printMAPM("now m is ",m);
  234 
  235     MAPM a,b,c;
  236     m.sincos(a,b);
  237     printMAPM("sincos [sin part]=",a);
  238     printMAPM("sincos [cos part]=",b);
  239 
  240     a = 4.409 * m;
  241     c = ceil(a);
  242     a = floor(a);                     /* make sure it's an integer */
  243     b = factorial(a);                 /* use library function      */
  244     b.toIntegerString(fBuf);    
  245 
  246     printMAPM("c = ",c);
  247     printMAPM("a = ",a);
  248     printf("factorial a = %s\n",fBuf);
  249 }
  250 
  251 
  252 
  253 void    SyntaxDemo()
  254 {
  255 #define ALL_DIGITS (-1)
  256 int       ii, jj, kk;
  257 static  char  obuf1[1024];
  258 static  char  obuf2[1024];
  259 static  char  obuf3[1024];
  260 static  char  obuf4[1024];
  261 
  262 MAPM a;               //a is uninitialized (value is undefined)
  263 MAPM b = "4.5e12";    //b is 4.5 trillion
  264 MAPM c = 3;           //c is now three
  265 MAPM d = 6.1;         //d is now six point one
  266 MAPM e = c;           //e is now three
  267 
  268 m_apm_cpp_precision(40);      // set minimum precision level
  269 
  270 fprintf(stdout, "----------------------------\n");
  271 fprintf(stdout, "Start of \'SyntaxDemo\' \n\n");
  272 
  273 b.toFixPtString(obuf1, ALL_DIGITS);
  274 c.toFixPtString(obuf2, ALL_DIGITS);
  275 d.toFixPtString(obuf3, ALL_DIGITS);
  276 e.toFixPtString(obuf4, ALL_DIGITS);
  277 fprintf(stdout, "b = [%s]    c = [%s]   d = [%s]   e = [%s]\n",
  278                     obuf1,      obuf2,      obuf3,    obuf4);
  279 
  280 e = -c * (23.9022 - d);
  281 
  282 e.toFixPtString(obuf4, ALL_DIGITS);
  283 fprintf(stdout, "e = [%s]\n", obuf4);
  284 
  285 
  286 a = sqrt(b);
  287 MAPM f;
  288 f = cbrt(e);
  289 
  290 a.toString(obuf1, ALL_DIGITS);
  291 f.toString(obuf2, 15);
  292 
  293 fprintf(stdout, "a = [%s]\n", obuf1);
  294 fprintf(stdout, "f = [%s]\n", obuf2);
  295 
  296 MAPM array2d[8][4];
  297 
  298 for (ii=0; ii < 8; ii++)
  299   {
  300    for (jj=0; jj < 4; jj++)
  301      {
  302       array2d[ii][jj] = 1000 * ii + 10 * jj;
  303      }
  304   }
  305 
  306 
  307 for (ii=0; ii < 8; ii++)
  308   {
  309    for (jj=0; jj < 4; jj++)
  310      {
  311       array2d[ii][jj].toIntegerString(obuf1);
  312       printf("[%d][%d] = %5s   ",ii,jj,obuf1);
  313      }
  314 
  315    printf("\n");
  316   }
  317 
  318 a = -23.93874;
  319 b = a.neg();
  320 
  321 if (b > a)
  322   printf("\'b\' > \'a\' \n");
  323 else
  324   printf("ERROR: \'a\' > \'b\' \n");
  325 
  326 a = 2;
  327 b = a.ipow(31);
  328 b.toIntegerString(obuf1);
  329 printf("2 ^ 31 = [%s] \n",obuf1);
  330 
  331 b = a.ipow(136);
  332 b.toIntegerString(obuf1);
  333 b.toString(obuf2,25);
  334 printf("2 ^ 136 = [%s] \n",obuf1);
  335 printf("2 ^ 136 = [%s] \n",obuf2);
  336 
  337 ii = b.exponent();
  338 jj = b.significant_digits();
  339 kk = b.is_integer();
  340 printf("exponent of \'b\' = %d   digits of \'b\' = %d   INTEGER flag = %d\n",
  341                           ii,jj,kk);
  342 a = -3.12;
  343 b = a.ipow(-5);
  344 b.toString(obuf1, 25);
  345 printf("-3.12 ^ -5 = [%s] \n",obuf1);
  346 
  347 
  348 //
  349 // ipow_nr is integer power function where the exponent
  350 // must be an integer >= 0.
  351 // this function will NOT do any rounding so it is ideal
  352 // for integer math
  353 //
  354 
  355 a = "13";
  356 b = a.ipow_nr(17);
  357 b.toIntegerString(obuf1);
  358 printf("13 ^ 17 = [%s] \n",obuf1);
  359 
  360 
  361 a = "81";
  362 b = a.ipow_nr(32);
  363 b.toIntegerString(obuf1);
  364 printf("81 ^ 32 = [%s] \n",obuf1);
  365 
  366 
  367 //   try out mod operator
  368 
  369 a = "3978381489565057925092524867";
  370 b = "53828048319274";
  371 
  372 f = a % b;
  373 
  374 a.toIntegerString(obuf1);
  375 b.toIntegerString(obuf2);
  376 f.toIntegerString(obuf3);
  377 printf("%s MOD %s = %s\n",obuf1,obuf2,obuf3);
  378 
  379 a %= b;
  380 
  381 a.toIntegerString(obuf1);
  382 printf("a = %s \n",obuf1);
  383 
  384 
  385 MAPM gcd1, lcm1;
  386 
  387 a = "8";
  388 a *= 3 * 7 * 19 * 53;
  389 
  390 b = "16";           // GCD = 8 * 3 * 19 = 456
  391 b *= 3 * 5 * 19;
  392 b *= 23 * 137;
  393 
  394 gcd1 = gcd(a, b);
  395 lcm1 = lcm(a, b);
  396 
  397 a.toIntegerString(obuf1);
  398 b.toIntegerString(obuf2);
  399 gcd1.toIntegerString(obuf3);
  400 lcm1.toIntegerString(obuf4);
  401 
  402 ii = gcd1.is_even();          // return boolean 
  403 jj = lcm1.is_odd();
  404 
  405 printf("gcd (%s, %s) = %s  \t\t ii = %d\n",obuf1,obuf2,obuf3,ii);
  406 printf("lcm (%s, %s) = %s  \t jj = %d\n",obuf1,obuf2,obuf4,jj);
  407 
  408 
  409 c = "1.0e+9131";
  410 d = log10(c);
  411 
  412 c.toString(obuf1, 10);
  413 d.toString(obuf2, 10);
  414 printf("log10( %s ) = %s \n",obuf1,obuf2);
  415 
  416 
  417 c = "6292734.982374929347e-817";
  418 d = log10(c);
  419 
  420 c.toString(obuf1, 20);
  421 d.toString(obuf2, 20);
  422 printf("log10( %s ) = %s \n",obuf1,obuf2);
  423 
  424 f  = 5.92734;
  425 ii = f.sign();
  426 f.toString(obuf1, ALL_DIGITS);
  427 printf("sign( %s ) = %d \n",obuf1,ii);
  428 
  429 f  = 0;
  430 ii = f.sign();
  431 f.toString(obuf1, ALL_DIGITS);
  432 printf("sign( %s ) = %d \n",obuf1,ii);
  433 
  434 f  = "-4.987234892347e-12073";
  435 ii = f.sign();
  436 f.toString(obuf1, ALL_DIGITS);
  437 printf("sign( %s ) = %d \n",obuf1,ii);
  438 
  439 
  440 e = sqrt(fabs(f)); 
  441 d = e.round(14);
  442 
  443 e.toString(obuf1, ALL_DIGITS);
  444 d.toString(obuf2, ALL_DIGITS);
  445 
  446 printf("e = [%s] \n",obuf1);
  447 printf("d = [%s] \n",obuf2);
  448 
  449 MAPM quot, remn, numer, denom;
  450 
  451 numer = "987598206958190375";
  452 denom = "289041345897";
  453 
  454 numer.toIntegerString(obuf1);
  455 denom.toIntegerString(obuf2);
  456 
  457 printf("numer = [%s] \n",obuf1);
  458 printf("denom = [%s] \n",obuf2);
  459 
  460 
  461 quot = numer.div(denom);
  462 quot.toIntegerString(obuf1);
  463 printf("quot  = [%s] \n",obuf1);
  464 
  465 remn = numer.rem(denom);
  466 
  467 remn.toIntegerString(obuf1);
  468 printf("remn  = [%s] \n",obuf1);
  469 
  470 
  471 numer.integer_div_rem(denom,quot,remn);
  472 quot.toIntegerString(obuf1);
  473 remn.toIntegerString(obuf2);
  474 printf("quot  = [%s] \n",obuf1);
  475 printf("remn  = [%s] \n",obuf2);
  476 
  477 
  478 /*
  479  *    show 2 typical and 3 not-so-typical
  480  *    methods to calculate PI.
  481  */
  482 
  483 MAPM  pi1, pi2, pi3, pi4, pi5;
  484 char  *obuf;
  485 
  486 obuf = obuf1;
  487 
  488 printf("\n");
  489 m_apm_cpp_precision(62);
  490 
  491 pi1 = acos("-1");
  492 pi2 = "2" * asin("1");
  493 pi3 = "6" * acos(sqrt("3") / -2) / 5;
  494 pi4 = "12" * asin((sqrt("3") - 1) * sqrt("2") * "0.25");
  495 pi5 = "12" * atan("2" + exp(log("3") / 2)) / 5;
  496 
  497 pi1.toFixPtString(obuf, 60);   printf("PI1 = [%s] \n",obuf); 
  498 pi2.toFixPtString(obuf, 60);   printf("PI2 = [%s] \n",obuf); 
  499 pi3.toFixPtString(obuf, 60);   printf("PI3 = [%s] \n",obuf); 
  500 pi4.toFixPtString(obuf, 60);   printf("PI4 = [%s] \n",obuf); 
  501 pi5.toFixPtString(obuf, 60);   printf("PI5 = [%s] \n",obuf); 
  502 
  503 
  504 
  505 MAPM  u, v, w, x, y, z;  
  506 
  507 m_apm_cpp_precision(50);
  508 
  509 printf("\n");
  510 
  511 x = 9.34231;
  512 y = -21;
  513 z = "-8.982349249829824921479824924792347921E-17";
  514 
  515 ii = x.is_integer();
  516 jj = y.is_integer();
  517 
  518 if (ii)
  519   printf("x is an integer  ");
  520 else
  521   printf("x is not an integer  ");
  522 
  523 if (jj)
  524   printf("y is an integer  ");
  525 else
  526   printf("y is not an integer  ");
  527 
  528 printf("\n");
  529 
  530 w = (82.30421 + sin(x / "21.11") / exp(y * 0.0726426)) * "4.32917E-2" / z;
  531 v = "3.742416" * log(-w);
  532 u = (2 + sqrt(v) + cbrt(v)) / 0.76;
  533 
  534 printf("\n");
  535 x.toString(obuf, 50);   printf("x = [%s] \n",obuf); 
  536 y.toString(obuf, 50);   printf("y = [%s] \n",obuf); 
  537 z.toString(obuf, 50);   printf("z = [%s] \n",obuf); 
  538 w.toString(obuf, 50);   printf("w = [%s] \n",obuf); 
  539 v.toString(obuf, 50);   printf("v = [%s] \n",obuf); 
  540 u.toString(obuf, 50);   printf("u = [%s] \n",obuf); 
  541 
  542 
  543 /*
  544  *   demo the fixed point formatting functions
  545  */
  546 
  547 printf("\n");
  548 
  549 a = u * "2.86521094E7";
  550 a.toString(obuf, 50);            printf("a = [%s] \n",obuf); 
  551 
  552 // convert 'a' to fixed point with 50 decimal places
  553 
  554 a.toFixPtString(obuf, 50);       printf("a = [%s] \n",obuf); 
  555 
  556 
  557 // convert 'a' to fixed point with 50 decimal places
  558 // change radix to comma, there is no separator char.
  559 
  560 a.toFixPtStringEx(obuf, 50, ',', 0, 0);   printf("a = [%s] \n",obuf); 
  561 
  562 
  563 // convert 'a' to fixed point with 50 decimal places
  564 // change radix to comma, use space as separator every 3 digits
  565 
  566 a.toFixPtStringEx(obuf, 50, ',', ' ', 3);   printf("a = [%s] \n",obuf); 
  567 
  568 
  569 // convert 'a' to fixed point with 40 decimal places
  570 // change radix to period, use comma as separator every 4 digits
  571 
  572 a.toFixPtStringEx(obuf, 40, '.', ',', 4);   printf("a = [%s] \n",obuf); 
  573 
  574 
  575 // convert 'a' to fixed point with 40 decimal places
  576 // change radix to period, use comma as separator every 5 digits
  577 // this function returns a char pointer to the string which
  578 // the caller is responsible for 'free'ing.
  579 
  580 char   *dest;
  581 char   radix = '.';
  582 char   separator = ',';
  583 int    separator_count = 5;
  584 int    decimal_places = 40;
  585 
  586 dest = a.toFixPtStringExp(decimal_places, radix, separator, separator_count);   
  587 printf("a = [%s] \n",dest); 
  588 free(dest);
  589 
  590 
  591 // convert 'a' to fixed point with 0 decimal places (round to integer)
  592 // set radix to dollar sign with no separator.
  593 // note that '$' will not be visible since we are using 
  594 // 0 decimal places (i.e., an integer format)
  595 
  596 radix = '$';
  597 decimal_places = 0;
  598 
  599 dest = a.toFixPtStringExp(decimal_places, radix, 0, 0);   
  600 printf("a = [%s] \n",dest); 
  601 free(dest);
  602 
  603 }
  604 
  605 /****************************************************************************/
  606 
  607 /*
  608  *      in case some compilers dont have a native cube root or
  609  *      inverse hyperbolic functions, we'll just make our own.
  610  *
  611  *      note that we are not doing any decent error checking on
  612  *      the inputs, these functions are here just to support 
  613  *      this program.
  614  */
  615 
  616 double  cbrt_local(double x)
  617 {
  618 if (x == 0.0)
  619   return(0.0);
  620 else
  621   {
  622    if (x < 0.0)
  623      return(-pow(-x, 0.333333333333333333));
  624    else
  625      return(pow(x, 0.333333333333333333));
  626   }
  627 }
  628 
  629 /****************************************************************************/
  630 
  631 double  asinh_local(double x)
  632 {
  633 return(log(x + sqrt(x * x + 1.0)));
  634 }
  635 
  636 /****************************************************************************/
  637 
  638 double  acosh_local(double x)
  639 {
  640 if (x >= 1.0)
  641   return(log(x + sqrt(x * x - 1.0)));
  642 else
  643   return(0.0);
  644 }
  645 
  646 /****************************************************************************/
  647 
  648 double  atanh_local(double x)
  649 {
  650 if (fabs(x) < 1.0)
  651   return(0.5 * log((1.0 + x) / (1.0 - x)));
  652 else
  653   return(0.0);
  654 }
  655 
  656 /****************************************************************************/
  657