home *** CD-ROM | disk | FTP | other *** search
/ Carousel Volume 2 #1 / carousel.iso / mactosh / lang / thinkc30.sit / MathHybrid.c < prev    next >
Text File  |  1989-01-16  |  10KB  |  711 lines
  1.  
  2. /*
  3.  
  4.     Math library for LightspeedC
  5.     
  6.     (C) Copyright 1986 THINK Technologies.  All rights reserved.
  7.     
  8.     For details, refer to Harbison & Steele's "C: A Reference Manual",
  9.     Chapter 11.
  10.     
  11.     Two versions of each function are defined by this library:  an
  12.     error-checking version (e.g. "sin") and a non-error-checking
  13.     version (e.g. "_sin").
  14.     
  15.     The non-underscore names can be made to refer to the non-error-
  16.     checking functions by #defining _NOERRORCHECK_ before #including
  17.     "math.h".  Doing so in THIS file suppresses the definitions of
  18.     the error-checking versions of the functions altogether.
  19.  
  20. */
  21.  
  22. /*#define _NOERRORCHECK_*/
  23. #include "math.h"
  24. #include "sane.h"
  25.  
  26. /*  useful constants  */
  27. static double Zero = 0.0;
  28. static double One = 1.0;
  29. static double Two = 2.0;
  30. static double MinusOne = -1.0;
  31. static double MinusTwo = -2.0;
  32. static double PointFive = 0.5;
  33. static double PointTwoFive = 0.25;
  34. static double Pi = PI;
  35. static double Pi2 = PI2;
  36. static double Log2Ten = 3.321928094887362348;
  37.  
  38. /* the 80-bit extended "_x80" is used as the destination to convert the 
  39. arguments to SANE doubles. */
  40. static Extended80 _x80;
  41.  
  42. /* the extra word of zeros is for 68881 compatibility. */
  43. static short _Max[] = { 0x7FFE, 0x0000, 0x7FFF, 0xFFFF, 0xFFFF, 0xFFFF };
  44. static short _MinusMax[] = { 0xFFFE, 0x0000, 0x7FFF, 0xFFFF, 0xFFFF, 0xFFFF };
  45. #define Max                (* (double *) _Max)
  46. #define MinusMax        (* (double *) _MinusMax)
  47.  
  48. /* routines for conversion between 80- and 96-bit formats */
  49.  
  50. void x80tox96(x80, x96)
  51. register Extended80 *x80;
  52. register Extended96 *x96;
  53. {
  54.     (*x96).exponent = (*x80).exponent;
  55.     (*x96).reserved = 0;
  56.     (*x96).mantissa = (*x80).mantissa;
  57. }
  58.  
  59. void x96tox80(x96, x80)
  60. register Extended96 *x96;
  61. register Extended80 *x80;
  62.  
  63. {
  64.     (*x80).exponent = (*x96).exponent;
  65.     (*x80).mantissa = (*x96).mantissa;
  66. }
  67.  
  68. #define _x80tox96(x80, x96) x80tox96(&x80, &x96)
  69. #define _x96tox80(x96, x80) x96tox80(&x96, &x80)
  70.  
  71. #define _to80(x)  _x96tox80(x, _x80);
  72. #define _to96(x)  _x80tox96(_x80, x);
  73.  
  74. /*  seed for pseudo-random number generator  */
  75. static unsigned long seed = 1;
  76.  
  77.  
  78. /*  environment word masks  */
  79. #define ROUND            0x6000
  80. #define ROUND_UP        0x2000
  81. #define ROUND_DOWN        0x4000
  82. #define ERROR            0x0F00
  83.  
  84.  
  85. /* ---------- error checking ---------- */
  86.  
  87. #ifdef _ERRORCHECK_
  88.  
  89. #define DomainCheck(test, result)        if (test) {                    \
  90.                                             errno = EDOM;            \
  91.                                             return(result);            \
  92.                                         }
  93. #define RangeCheck(target)                if (GetState() & ERROR) {    \
  94.                                             errno = ERANGE;            \
  95.                                             target = Max;            \
  96.                                         }
  97.  
  98. #else _ERRORCHECK_
  99.  
  100. #define ClearExceptions()
  101. #define DomainCheck(test, result)
  102. #define RangeCheck(target)
  103.  
  104. #endif _ERRORCHECK_
  105.  
  106. /*
  107.  *  GetState - query environment word
  108.  *
  109.  */
  110.  
  111. static
  112. GetState()
  113. {
  114.     int state;
  115.     
  116.     GetEnvironment(&state);
  117.     return(state);
  118. }
  119.  
  120.  
  121. /*
  122.  *  SetState - define environment word
  123.  *
  124.  */
  125.  
  126. static
  127. SetState(state)
  128. {
  129.     SetEnvironment(&state);
  130. }
  131.  
  132.  
  133. /*
  134.  *  ClearExceptions - clear error conditions
  135.  *
  136.  *  This must be called if RangeCheck is to be used.
  137.  *
  138.  */
  139.  
  140. #ifdef _ERRORCHECK_
  141. static
  142. ClearExceptions()
  143. {
  144.     int state;
  145.     
  146.     GetEnvironment(&state);
  147.     state &= ~ERROR;
  148.     SetEnvironment(&state);
  149. }
  150. #endif _ERRORCHECK_
  151.  
  152.  
  153. /*
  154.  *  SetRound - define rounding direction
  155.  *
  156.  *  The previous environment word is returned.  The rounding direction can
  157.  *  be restored by passing this value to SetState.
  158.  *
  159.  */
  160.  
  161. static
  162. SetRound(direction)
  163. {
  164.     int state;
  165.  
  166.     GetEnvironment(&state);
  167.     SetState((state & ~ROUND) | direction);
  168.     return(state);
  169. }
  170.  
  171.  
  172. /*
  173.  *  xfersign - transfer sign from one floating number to another
  174.  *
  175.  */
  176.  
  177. static
  178. xfersign(x, yp)
  179. double x, *yp;
  180. {
  181.     asm {
  182.         movea.l    yp,a0
  183.         bclr    #7,(a0)
  184.         tst.w    x
  185.         bpl.s    @1
  186.         bset    #7,(a0)
  187. @1    }
  188. }
  189.  
  190.  
  191. /* ---------- math functions (alphabetically) ---------- */
  192.  
  193.  
  194. /*
  195.  *  abs - absolute value of an integer
  196.  *
  197.  */
  198.  
  199. int abs(x)
  200. int x;
  201. {
  202.     return(x < 0 ? -x : x);
  203. }
  204.  
  205. /* the real functions are here, in fact, when _ERRORCHECK_ is undefined,
  206. they will be called directly.*/
  207.  
  208. /* some pre-declarations..*/
  209. double _exp();
  210. double _fabs();
  211.  
  212. double _acos(x)
  213. double x;
  214. {
  215.     return(Two * atan(sqrt((One - x) / (One + x))));
  216. }
  217.  
  218. double _asin(x)
  219. double x;
  220. {
  221.     double y = _fabs(x);
  222.  
  223.     DomainCheck(y > One, Zero);
  224.     if (y == One) {
  225.         y = Pi2;
  226.         xfersign(x, &y);
  227.         return(y);
  228.     }
  229.     if (y > PointFive) {
  230.         y = One - y;
  231.         y = Two * y - y * y;
  232.     }
  233.     else
  234.         y = One - y * y;
  235.     return(atan(x / sqrt(y)));
  236. }
  237.  
  238. double _atan(x)
  239. double x;
  240. {
  241.     _to80(x);
  242.     elems68k(&_x80, FOATANX);
  243.     _to96(x);
  244.     return(x);
  245. }
  246.  
  247. double _atan2(y, x)
  248. double y, x;
  249. {
  250.     double z;
  251.  
  252.     if (x == 0) {
  253.         DomainCheck(y == 0, Zero);
  254.         z = Pi2;
  255.     }
  256.     else {
  257.         z = atan(fabs(y / x));
  258.         if (x < 0)
  259.             z = Pi - z;
  260.     }
  261.     xfersign(y, &z);
  262.     return(z);
  263. }
  264.  
  265. double _ceil(x)
  266. double x;
  267. {
  268.     register int state = SetRound(ROUND_UP);
  269.     _to80(x);
  270.     fp68k(&_x80, FORTI);
  271.     SetState(state);
  272.     _to96(x);
  273.     return(x);
  274. }
  275.  
  276. double _cos(x)
  277. double x;
  278. {    
  279.     _to80(x)
  280.     elems68k(&_x80, FOCOSX);
  281.     _to96(x)
  282.     return(x);
  283. }
  284.  
  285. double _cosh(x)
  286. double x;
  287. {
  288.     ClearExceptions();
  289.     x = PointFive * _exp(_fabs(x));
  290.     x += PointTwoFive / x;
  291.     RangeCheck(x);
  292.     return(x);
  293. }
  294.  
  295. double _exp(x)
  296. double x;
  297. {
  298.     _to80(x);
  299.     elems68k(&_x80, FOEXPX);
  300.     _to96(x);
  301.     return(x);
  302. }
  303.  
  304. double _fabs(x)
  305. double x;
  306. {
  307.     _to80(x);
  308.     fp68k(&_x80, FOABS);
  309.     _to96(x);
  310.     return(x);
  311. }
  312.  
  313. double _floor(x)
  314. double x;
  315. {
  316.     register int state = SetRound(ROUND_DOWN);
  317.     
  318.     _to80(x);
  319.     fp68k(&_x80, FORTI);
  320.     _to96(x);
  321.     SetState(state);
  322.     return(x);
  323. }
  324.  
  325. double _fmod(x, y)
  326. double x, y;
  327. {
  328.     double z;
  329.     Extended80 tz;
  330.     _x96tox80(x, tz);
  331.     _to80(y);
  332.     fp68k(&_x80, FOABS);
  333.     fp68k(&_x80, &tz, FOREM);
  334.     _to96(y);
  335.     _x80tox96(tz, z);
  336.     
  337.     if (x > 0 && z < 0)
  338.         z += y;
  339.     else if (x < 0 && z > 0)
  340.         z -= y;
  341.     return(z);
  342. }
  343.  
  344. double _frexp(x, nptr)
  345. double x;
  346. register int *nptr;
  347. {
  348.     double y = fabs(x), z = Two;
  349.     Extended80 ty, tz;
  350.     
  351.     if (y == 0) {
  352.         *nptr = 0;
  353.         return(Zero);
  354.     }
  355.     
  356.     _x96tox80(y, ty);
  357.     elems68k(&ty, FOLOG2X);
  358.     _x80tox96(ty, y);
  359.     
  360.     y -= *nptr = y;
  361.     
  362.     _x96tox80(y, ty);
  363.     _x96tox80(z, tz);
  364.     elems68k(&ty, &tz, FOPWRY);
  365.     _x80tox96(ty, y);
  366.     _x80tox96(tz, z);
  367.     
  368.     if (z >= One) {
  369.         z *= PointFive;
  370.         ++*nptr;
  371.     }
  372.     else if (z < PointFive) {
  373.         z += z;
  374.         --*nptr;
  375.     }
  376.     xfersign(x, &z);
  377.     return(z);
  378. }
  379.  
  380. double _ldexp(x, n)
  381. double x;
  382. int n;
  383. {
  384.     _to80(x);
  385.     fp68k(&n, &_x80, FOSCALB);
  386.     _to96(x);
  387.     return(x);
  388. }
  389.  
  390. double _log(x)
  391. double x;
  392. {
  393.     _to80(x);
  394.     elems68k(&_x80, FOLNX);
  395.     _to96(x);
  396.     return(x);
  397. }
  398.  
  399. double _log10(x)
  400. double x;
  401. {
  402.     _to80(x);
  403.     elems68k(&_x80, FOLOG2X);    /* LOG2 is much faster than LN */
  404.     _to96(x);
  405.     return(x/ Log2Ten);
  406. }
  407.  
  408. double _modf(x, nptr)
  409. double x;
  410. register double *nptr;
  411. {
  412.     _to80(x);
  413.     fp68k(&_x80, FOTTI);
  414.     x80tox96(&_x80, nptr);
  415.     return(x - *nptr);
  416. }
  417.  
  418. double _pow(x, y)
  419. double x, y;
  420. {
  421.     double tx, ty;
  422.     _x96tox80(y, ty);
  423.     _x96tox80(x, tx);
  424.     elems68k(&ty, &tx, FOPWRY);
  425. /*    _x80tox96(ty, y);    */
  426.     _x80tox96(tx, x);
  427.     return(x);
  428. }
  429.  
  430. double _sin(x)
  431. double x;
  432. {
  433.     _to80(x);
  434.     elems68k(&_x80, FOSINX);
  435.     _to96(x);
  436.     return(x);
  437. }
  438.  
  439. double _sinh(x)
  440. double x;
  441. {
  442.     double y = _fabs(x);
  443.     
  444.     _to80(y);
  445.     elems68k(&_x80, FOEXP1X);
  446.     _to96(y);
  447.     y += y / (y + One);
  448.     y *= PointFive;
  449.  
  450.     xfersign(x, &y);
  451.     return(y);
  452. }
  453.  
  454. double _sqrt(x)
  455. double x;
  456. {
  457.     _to80(x);
  458.     fp68k(&_x80, FOSQRT);
  459.     _to96(x);
  460.     return(x);
  461. }
  462.  
  463. double _tan(x)
  464. double x;
  465. {
  466.     _to80(x);
  467.     elems68k(&_x80, FOTANX);
  468.     _to96(x);
  469.     return(x);
  470. }
  471.  
  472. double _tanh(x)
  473. double x;
  474. {
  475.     double y = MinusTwo * _fabs(x);
  476.     _to80(y);
  477.     elems68k(&_x80, FOEXP1X);
  478.     _to96(y);
  479.     y = -y / (y + Two);
  480.     xfersign(x, &y);
  481.     return(y);
  482. }
  483.  
  484. long labs(x)
  485. long x;
  486. {
  487.     return(x < 0 ? x: -x);
  488. }
  489.  
  490. int rand()
  491. {
  492.     seed = seed * 1103515245 + 12345;
  493.     asm {
  494.         move.w    seed,d0            ;  high word of long
  495.         andi.w    #0x7FFF,d0        ;  remove high bit
  496.     }
  497. }
  498.  
  499. void srand(n)
  500. unsigned n;
  501. {
  502.     seed = n;
  503. }
  504.  
  505.  
  506. #ifdef _ERRORCHECK_
  507. double acos(x)
  508. double x;
  509. {
  510.     DomainCheck(x > One || x < MinusOne, Zero);
  511.     if (x == MinusOne)
  512.         return(Pi);
  513.     return(_acos(x));
  514. }
  515.  
  516. double asin(x)
  517. double x;
  518. {
  519.     return(_asin(x));
  520. }
  521.  
  522.  
  523. double atan(x)
  524. double x;
  525. {
  526.     return(_atan(x));
  527. }
  528.  
  529.  
  530. double atan2(y, x)
  531. double y, x;
  532. {
  533.     return(_atan2(y, x));
  534. }
  535.  
  536. double ceil(x)
  537. double x;
  538. {
  539.     return(_ceil(x));
  540. }
  541.  
  542. double cos(x)
  543. double x;
  544. {    
  545.     return(_cos(x));
  546. }
  547.  
  548. double cosh(x)
  549. double x;
  550. {
  551.     return(_cosh(x));
  552. }
  553.  
  554. double exp(x)
  555. double x;
  556. {
  557.     ClearExceptions();
  558.     x = _exp(x);
  559.     RangeCheck(x);
  560.     return(x);
  561. }
  562.  
  563. double fabs(x)
  564. double x;
  565. {
  566.     return(_fabs(x));
  567. }
  568.  
  569. double floor(x)
  570. double x;
  571. {
  572.     return(_floor(x));
  573. }
  574.  
  575. double fmod(x, y)
  576. double x, y;
  577. {
  578.     return(_fmod(x, y));
  579. }
  580.  
  581. double frexp(x, nptr)
  582. double x;
  583. register int *nptr;
  584. {
  585.     return(_frexp(x, nptr));
  586. }
  587.  
  588.  
  589. /*
  590.  *  ldexp - combine fraction/exponent into a floating number
  591.  *
  592.  */
  593.  
  594. double ldexp(x, n)
  595. double x;
  596. int n;
  597. {
  598.     return(_ldexp(x, n));
  599. }
  600.  
  601.  
  602. double log(x)
  603. double x;
  604. {
  605.     DomainCheck(x <= 0, MinusMax);
  606.     return(_log(x));
  607. }
  608.  
  609. double log10(x)
  610. double x;
  611. {
  612.     DomainCheck(x <= 0, MinusMax);
  613.     return(_log10(x));
  614. }
  615.  
  616.  
  617. double modf(x, nptr)
  618. double x, *nptr;
  619. {
  620.     return(_modf(x, nptr));
  621. }
  622.  
  623. double pow(x, y)
  624. double x, y;
  625. {
  626.     int odd = 0;
  627.     Extended80 tx, ty;
  628.     
  629.     ClearExceptions();
  630.     if (x == 0) {
  631.         if (y <= 0) {
  632. #ifdef _ERRORCHECK_
  633.             errno = EDOM;
  634. #endif
  635.             return (MinusMax);
  636.         }
  637.         return(Zero);
  638.     }
  639.  
  640.     if (y == 0)
  641.         return(One);
  642.     if (x < 0) {
  643.         if (_modf(y, &y)) {
  644. #ifdef _ERRORCHECK_
  645.             errno = EDOM;
  646. #endif
  647.             return (MinusMax);
  648.         }
  649.         x = -x;
  650.         odd = _fmod(y, Two);
  651.     }
  652.  
  653.     x = _pow(x, y);
  654.     RangeCheck(x);
  655.     return(odd ? -x : x);
  656. }
  657.  
  658. double sin(x)
  659. double x;
  660. {
  661.     return(_sin(x));
  662. }
  663.  
  664.  
  665. double sinh(x)
  666. double x;
  667. {
  668.     double y;
  669.     
  670.     ClearExceptions();
  671.     y = _sinh(x);
  672.     RangeCheck(y);
  673.     return(y);
  674. }
  675.  
  676. double sqrt(x)
  677. double x;
  678. {
  679.     DomainCheck(x < 0, Zero);
  680.     return(_sqrt(x));
  681. }
  682.  
  683.  
  684. /*
  685.  *  tan - circular tangent
  686.  *
  687.  */
  688.  
  689. double tan(x)
  690. double x;
  691. {
  692.     ClearExceptions();
  693.     x = _tan(x);
  694.     RangeCheck(x);
  695.     return(x);
  696. }
  697.  
  698.  
  699. /*
  700.  *  tanh - hyperbolic tangent
  701.  *
  702.  */
  703.  
  704. double tanh(x)
  705. double x;
  706. {
  707.     return(_tanh(x));
  708. }
  709.  
  710. #endif
  711.