home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / c / math.c < prev    next >
C/C++ Source or Header  |  1994-03-04  |  12KB  |  520 lines
  1. /*****************************************************************************
  2.  *                                         *
  3.  * MATH.C  Mathematical subroutines for use with Lattice or Microsoft-         *
  4.  *       compiler.                                 *
  5.  *       By Max R. D^Arsteler                             *
  6.  *                                         *
  7.  *****************************************************************************
  8.  */
  9.  
  10. /*****************************************************************************
  11.  *  This file contains the combined source files of the most
  12.  *  important mathe- matical formulas written for the Lattice C Compiler.
  13.  *  It makes no use of the 8087.  Some of the functions make use of the
  14.  *  internal binary representation of double precision numbers.  As the
  15.  *  internal data representation may vary between different compilers,
  16.  *  watch out for errors e.g.  in the power function.  The function values
  17.  *  are calculated through polynominal or rational approximations.  The
  18.  *  results are for allmost all functions precise for the first 9 to 10
  19.  *  decimal digits.  For still higher precision or checking of values use
  20.  *  functions from file mathf.c (not available here), which contains
  21.  *  Taylor's expansions series.  As this requires many iterations with
  22.  *  multiplications and divisions, the calculation time is too long for use
  23.  *  in regular programs.  In order to use these functions, a program has to
  24.  *  include the header file "math.h".  Furthermore, either the compiled
  25.  *  version (math.obj) or a library with the compiled functions has to be
  26.  *  linked to the program.  As exemples see included testprograms
  27.  *  "tstsin.c" and "tstpow.c" To compile:  "mcb tstsin <CR>", to link:"link
  28.  *  c+%1+math.obj, %1, ,mc <CR>".  I recommend to compile the source of
  29.  *  each function seperately with your personal compiler and use the MS-DOS
  30.  *  librarian to build a library of just the functions you need.  Some of
  31.  *  the more esotheric functions used on "UNIX" are such as Bessel
  32.  *  functions or the gamma function are not implemented.  Evidently, this
  33.  *  is not a professional implementation.
  34.  *
  35.  *   Rockville, 11/27/83      Max R. D^Arsteler
  36.  *                   12405 Village Sq. Terrace
  37.  *                   Rockville Md. 20852    H 301 984-1168
  38.  ******************************************************************************/
  39.  
  40. /* SQRT.C
  41.  * Calculates square root
  42.  * Precision: 11.05 E-10
  43.  * Range: 0 to e308
  44.  * Author: Max R. D^Arsteler
  45.  */
  46. double sqrt(x)
  47. double x;
  48. {
  49.     typedef union {           /* Makes use of internal data  */
  50.         double d;          /* representation */
  51.         unsigned u[4];
  52.     } DBL;
  53.     double y, re, p, q;
  54.     unsigned ex;        /* exponens*/
  55.     DBL *xp;
  56.  
  57.     xp = (DBL *)&x;
  58.     ex = xp->u[3] & ~0100017; /* save exponens */
  59.     re = ex & 020 ? 1.4142135623730950488 : 1.0;
  60.     ex = ex - 037740 >> 1;
  61.     ex &= ~0100017;
  62.     xp->u[3] &= ~0177760;      /* erase exponens and sign*/
  63.     xp->u[3] |= 037740;      /* arrange for mantissa in range 0.5 to 1.0 */
  64.     if (xp->d < .7071067812) {  /* multiply by sqrt(2) if mantissa < 0.7 */
  65.         xp->d *= 1.4142135623730950488;
  66.         if (re > 1.0) re = 1.189207115002721;
  67.         else re = 1.0 / 1.18920711500271;
  68.     }
  69.     p =           .54525387389085 +   /* Polynomial approximation */
  70.         xp->d * (13.65944682358639 +   /* from: COMPUTER APPROXIMATIONS*/
  71.         xp->d * (27.090122712655   +   /* Hart, J.F. et al. 1968 */
  72.         xp->d *   6.43107724139784));
  73.     q =          4.17003771413707 +
  74.         xp->d * (24.84175210765124 +
  75.         xp->d * (17.71411083016702 +
  76.         xp->d ));
  77.     xp->d = p / q;
  78.     xp->u[3] = xp->u[3] + ex & ~0100000;
  79.     xp->d *= re;
  80.     return(xp->d);
  81. }
  82.  
  83.  
  84. /* EXP2.C
  85.  * Calculates exponens to base 2 using polynomial approximations
  86.  * Precision: 10E-10
  87.  * Range:
  88.  * Author: Max R. D^Arsteler
  89.  */
  90. double exp2(x)
  91. double x;
  92. {
  93.     typedef union {
  94.         double d;
  95.         unsigned u[4];
  96.     } DBL;
  97.     double y, x2, p, q, re;
  98.     int ix;
  99.     DBL *xp, *yp;
  100.  
  101.     xp = (DBL *)&x;
  102.     y = 0.0;
  103.     yp = (DBL *)&y;
  104.     if(xp->d > 1023.0 || xp->d < -1023.0) return (1E307);
  105.     ix = (int) xp->d;
  106.     yp->u[3] += ix + 1023 << 4;
  107.     yp->u[3] &= ~0100017;
  108.     if ((xp->d -= (double) ix) == 0.0) return(yp->d);
  109.     if (xp->d < 0.0) {
  110.           yp->u[3] -= 1 << 4;
  111.           yp->u[3] &= ~0100017;
  112.           xp->d++;
  113.     }
  114.     if (xp->d >= 0.5) {  /* adjust to range 0-0.5 */
  115.         xp->d -= 0.5;
  116.         re = 1.41421356237309504880;
  117.     }
  118.     else re = 1.0;
  119.     x2 = xp->d * xp->d;
  120.     p = xp->d * (7.2152891511493 +
  121.            x2 *  0.0576900723731);
  122.     q =        20.8189237930062 +
  123.            x2 ;
  124.     xp->d = (q + p) / (q - p);
  125.     yp->d *= re * xp->d;
  126.     return(yp->d);
  127. }
  128.  
  129.  
  130. /* POW.C
  131.  * Calculates y^x
  132.  * Precision:
  133.  * Range: o to big for y, +/-1023 for x
  134.  * Author: Max R. D^Arsteler, 10/2/83
  135.  */
  136.  
  137. double pow(y,x)
  138. double y, x;
  139. {
  140.     typedef union {
  141.         double d;
  142.         unsigned u[4];
  143.     } DBL;
  144.     double z, w, p, p2, q, re;
  145.     unsigned ex;        /* exponens*/
  146.     int iz;
  147.     DBL *yp, *zp, *wp;
  148.  
  149.     yp = (DBL *)&y;
  150.     if (yp->d <= 0.0) y = -y;
  151.     z = 0.0;
  152.     zp = (DBL *)&z;
  153.     zp->u[3] = yp->u[3] & ~0100017; /* save exponens */
  154.     iz = (zp->u[3] >> 4)-1023;
  155.     if ((yp->d - zp->d) == 0.0)
  156.         z = (double)iz;
  157.     else {
  158.         yp->u[3] -= ++iz << 4; /* arrange for range 0.5 to 0.99999999999 */
  159.         yp->d *= 1.4142135623730950488;  /* shift for 1/sqrt(2) to sqrt(2) */
  160.         p = (yp->d - 1.0) / (yp->d + 1.0);
  161.         p2 = p * p;
  162.         z = p  * (2.000000000046727 +    /* Polynomial approximation */
  163.             p2 * (0.666666635059382 +    /* from: COMPUTER APPROXIMATIONS*/
  164.             p2 * (0.4000059794795   +    /* Hart, J.F. et al. 1968 */
  165.             p2 * (0.28525381498     +
  166.             p2 *  0.2376245609 ))));
  167.         z = z * 1.442695040888634 + (double)iz    - 0.5;
  168.     }
  169.     z *= x;
  170.     w = 0.0;
  171.     wp = (DBL *)&w;
  172.     if(zp->d > 1023.0 || zp->d < -1023.0) return (1E307);
  173.     iz = (int) zp->d;
  174.     wp->u[3] += iz + 1023 << 4;
  175.     wp->u[3] &= ~0100017;
  176.     if ((zp->d -= (double) iz) == 0.0) return(wp->d);
  177.     while (zp->d < 0.0) {
  178.           wp->u[3] -= 1 << 4;
  179.           wp->u[3] &= ~0100017;
  180.           zp->d++;
  181.     }
  182.     if (zp->d >= 0.5) {  /* adjust to range 0-0.5 */
  183.         zp->d -= 0.5;
  184.         re = 1.41421356237309504880;
  185.     }
  186.     else re = 1.0;
  187.     p2 = zp->d * zp->d;
  188.     p = zp->d * (7.2152891511493 +
  189.            p2 *  0.0576900723731);
  190.     q =        20.8189237930062 +
  191.            p2 ;
  192.     zp->d = re * wp->d * (q + p) / (q - p);
  193.     return(zp->d);
  194. }
  195.  
  196.  
  197. /* FMOD.C
  198.  * Returns the number f such that x = i*y + f. i is an integer, and
  199.  * 0 <= f < y.
  200.  */
  201. double fmod(x, y)
  202. double x, y;
  203. {
  204.     double zint, z;
  205.     int i;
  206.  
  207.     z = x / y;
  208.     zint = 0.0;
  209.     while (z > 32768.0) {
  210.         zint += 32768.0;
  211.         z -= 32768;
  212.     }
  213.     while (z < -32768.0) {
  214.         zint -= 32768.0;
  215.         z += 32768.0;
  216.     }
  217.     i = (int) z;
  218.     zint += (double) i;
  219.     return( x - zint * y);
  220. }
  221.  
  222. /*ASIN.C
  223.  *Calculates arcsin(x)
  224.  *Range: 0 <= x <= 1
  225.  *Precision: +/- .000,000,02
  226.  *Header: math.h
  227.  *Author: Max R. D^Arsteler
  228.  */
  229. extern double sqrt();
  230.  
  231. double asin(x)
  232. double x;
  233.  
  234. {
  235.     double y;
  236.     int sign;
  237.  
  238.     if (x > 1.0 || x < -1.0) exit(1);
  239.     sign = 0;
  240.     if (x < 0) {
  241.          sign = 1;
  242.          x = -x;
  243.     }
  244.     y = ((((((-.0012624911    * x
  245.          + .0066700901) * x
  246.          - .0170881256) * x
  247.          + .0308918810) * x
  248.          - .0501743046) * x
  249.          + .0889789874) * x
  250.          - .2145988016) * x
  251.          +1.5707963050;
  252.     y = 1.57079632679 - sqrt(1.0 - x) * y;
  253.     if (sign) y = -y;
  254.     return(y);
  255. }
  256.  
  257. /* SIN.C
  258.  * Calculates sin(x), angle x must be in rad.
  259.  * Range: -pi/2 <= x <= pi/2
  260.  * Precision: +/- .000,000,005
  261.  * Header: math.h
  262.  * Author: Max R. D^Arsteler
  263.  */
  264.  
  265. double sin(x)
  266. double x;
  267. {
  268.     double xi, y, q, q2;
  269.     int sign;
  270.  
  271.     xi = x; sign = 1;
  272.     while (xi < -1.57079632679489661923) xi += 6.28318530717958647692;
  273.     while (xi > 4.71238898038468985769) xi -= 6.28318530717958647692;
  274.     if (xi > 1.57079632679489661923) {
  275.         xi -= 3.141592653589793238462643;
  276.         sign = -1;
  277.     }
  278.     q = xi / 1.57079632679; q2 = q * q;
  279.     y = ((((.00015148419  * q2
  280.           - .00467376557) * q2
  281.           + .07968967928) * q2
  282.           - .64596371106) * q2
  283.           +1.57079631847) * q;
  284.     return(sign < 0? -y : y);
  285. }
  286.  
  287.  
  288. /* LOG.C
  289.  * Calculates natural logarithmus
  290.  * Precision: 11.56 E-10
  291.  * Range: 0 to e308
  292.  * Author: Max R. D^Arsteler
  293.  */
  294. double log(x)
  295. double x;
  296. {
  297.     typedef union {
  298.         double d;
  299.         unsigned u[4];
  300.     } DBL;
  301.     double y, z, z2, p;
  302.     unsigned ex;        /* exponens*/
  303.     int ix;
  304.     DBL *xp, *yp;
  305.  
  306.     xp = (DBL *)&x;
  307.     if (xp->d <= 0.0) return(y);
  308.     y = 0.0;
  309.     yp = (DBL *)&y;
  310.     yp->u[3] = xp->u[3] & ~0100017; /* save exponens */
  311.     ix = (yp->u[3] >> 4)-1023;
  312.     if ((xp->d - yp->d) == 0.0) return( .693147180559945 * (double)ix);
  313.     xp->u[3] -= ++ix << 4; /* arrange for range 0.5 to 0.99999999999 */
  314.     xp->d *= 1.4142135623730950488;  /* shift for 1/sqrt(2) to sqrt(2) */
  315.     z = (xp->d - 1.0) / (xp->d + 1.0);
  316.     z2 = z * z;
  317.     y = z  * (2.000000000046727 +    /* Polynomial approximation */
  318.         z2 * (0.666666635059382 +    /* from: COMPUTER APPROXIMATIONS*/
  319.         z2 * (0.4000059794795   +    /* Hart, J.F. et al. 1968 */
  320.         z2 * (0.28525381498     +
  321.         z2 *  0.2376245609 ))));
  322.     y = y + .693147180559945 * ((double)ix    - 0.5);
  323.     return(yp->d);
  324. }
  325.  
  326. /* ATAN.C
  327.  * Calculates arctan(x)
  328.  * Range: -infinite <= x <= infinite (Output -pi/2 to +pi/2)
  329.  * Precision: +/- .000,000,04
  330.  * Header: math.h
  331.  * Author: Max R. D^Arsteler 9/15/83
  332.  */
  333.  
  334. double atan(x)
  335. double x;
  336. {
  337.     double xi, q, q2, y;
  338.     int sign;
  339.  
  340.     xi = (x < 0. ? -x : x);
  341.     q = (xi - 1.0) / (xi + 1.0); q2 = q * q;
  342.     y = ((((((( - .0040540580 * q2
  343.              + .0218612286) * q2
  344.              - .0559098861) * q2
  345.              + .0964200441) * q2
  346.              - .1390853351) * q2
  347.              + .1994653599) * q2
  348.              - .3332985605) * q2
  349.              + .9999993329) * q + 0.785398163397;
  350.     return(x < 0. ? -y: y);
  351. }
  352.  
  353.  
  354. /* FLOOR.C
  355.  * Returns largest integer not greater than x
  356.  * Author: Max R. D^Arsteler 9/26/83
  357.  */
  358. double floor(x)
  359. double x;
  360. {
  361.     double y;
  362.     int ix;
  363.  
  364.     y = 0.0;
  365.     while (x >= 32768.0) {
  366.         y += 32768.0;
  367.         x -= 32768.0;
  368.     }
  369.     while (x <= -32768.0) {
  370.         y -= 32768.0;
  371.         x += 32768.0;
  372.     }
  373.     if (x > 0.0) ix = (int) x;
  374.     else ix = (int)(x - 0.9999999999);
  375.     return( y + (double) ix);
  376. }
  377.  
  378.  
  379. /* CEIL.C
  380.  * Returns smallest integer not less than x
  381.  * Author: Max R. D^Arsteler 9/26/83
  382.  */
  383. double ceil(x)
  384. double x;
  385. {
  386.     double y;
  387.     int ix;
  388.  
  389.     y = 0.0;
  390.     while (x >= 32768.0) {
  391.         y += 32768.0;
  392.         x -= 32768.0;
  393.     }
  394.     while (x <= -32768.0) {
  395.         y -= 32768.0;
  396.         x += 32768.0;
  397.     }
  398.     if (x > 0.0) ix = (int) (x + 0.999999999999999);
  399.     else ix = (int) x;
  400.     return( y + (double) ix);
  401. }
  402.  
  403. /* EXP.C
  404.  * Calculates exponens of x to base e
  405.  * Range: +/- exp(88)
  406.  * Precision: +/- .000,000,000,1
  407.  * Author: Max R. D^Arsteler, 9/20/83
  408.  */
  409.  double exp(xi)
  410.  double xi;
  411.  {
  412.     double y, ex, px, nn, ds, in;
  413.  
  414.     if (xi > 88.0) return(1.7014117331926443e38);
  415.     if (xi < -88.0) return(0.0);
  416.     ex = 1.0;
  417.     while( xi > 1.0) {
  418.         ex *= 2.718281828459; /* const. e */
  419.         xi--;
  420.     }
  421.     while( xi < -1.0) {
  422.         ex *= .367879441171;  /* 1/e */
  423.         xi++;
  424.     }
  425. /* Slow, but more precise Taylor expansion series */
  426.     y = ds = 1.0; nn = 0.0;
  427.     while ((ds < 0.0 ? -ds : ds) > 0.000000000001) {
  428.         px = xi/++nn;          /* above precision required */
  429.         ds *= px;
  430.         y += ds;
  431.     }
  432.        y *= ex;
  433.  
  434. /*  Chebyshev polynomials: fast, but less precise then expected!
  435.  *    xi = -xi;
  436.  *    y = (((((.0000006906  * xi
  437.  *        +.0000054302) * xi
  438.  *        +.0001715620) * xi
  439.  *        +.0025913712) * xi
  440.  *        +.0312575832) * xi
  441.  *        +.2499986842) * xi
  442.  *        +1.0;
  443.  *     y = ex / (y * y * y * y);
  444.  */
  445.     return(y);
  446. }
  447.  
  448.  
  449. /* LOG10.C
  450.  * Approximation for logarithm of basis 10
  451.  * Range 0 < x < 1e+38
  452.  * Precision: +/- 0.000,000,1
  453.  * Header: math.h
  454.  * Author: Max R. D^Arsteler 9/15/83
  455.  */
  456.  
  457. /* Method of Chebyshev polynomials */
  458. /* C. Hastings, jr. 1955 */
  459.  
  460. double log10(x)
  461. double x;
  462. {
  463.     double xi, y, q, q2;
  464.     int ex;
  465.  
  466.     if (x <= 0.0) return(0.); /* Error!! */
  467.     ex = 0.0; xi = x;
  468.     while (xi < 1.0 ) {
  469.         xi *= 10.0;
  470.         ex--;
  471.     }
  472.     while (xi > 10.0) {
  473.         xi *= 0.1;
  474.         ex++;
  475.     }
  476.     q = (xi - 3.16227766) / (xi + 3.16227766); q2 = q * q;
  477.     y = ((((.191337714 * q2
  478.           + .094376476) * q2
  479.           + .177522071) * q2
  480.           + .289335524) * q2
  481.           + .868591718) * q + .5;
  482.     y += (double) ex;
  483.     return(y);
  484. }
  485.  
  486. /* PW10.C
  487.  * Calculates 10 power x
  488.  * Range: 0 <= x <= 1
  489.  * Precision: +/- 0.000,000,005
  490.  * Header: math.lib
  491.  * Author: Max R. D^Arsteler 9/15/83
  492.  */
  493.  
  494. double pw10(x)
  495. double x;
  496. {
  497.     double xi,ex,y;
  498.  
  499.     if (x > 38.0) return(1.7014117331926443e38);
  500.     if (x < -38.0) return(0.0);
  501.     xi = x; ex = 1.0;
  502.     while (xi > 1.0) {
  503.         ex *= 10.0;
  504.         xi--;
  505.     }
  506.     while (xi < 0.0) {
  507.         ex *= 0.1;
  508.         xi++;
  509.     }
  510.     y = ((((((.00093264267    * xi
  511.         + .00255491796) * xi
  512.         + .01742111988) * xi
  513.         + .07295173666) * xi
  514.         + .25439357484) * xi
  515.         + .66273088429) * xi
  516.         +1.15129277603) * xi + 1.0;
  517.     y *= y;
  518.     return(y*ex);
  519. }
  520.