home *** CD-ROM | disk | FTP | other *** search
/ Carousel Volume 2 #1 / carousel.iso / mactosh / lang / thinkc30.sit / Math.c < prev    next >
Text File  |  1989-01-16  |  8KB  |  613 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.     Only one version of each function is defined by this library.
  12.     Whether error checking is done is determined by how this file
  13.     is compiled.  Error checking is ENABLED when the following
  14.     #define is COMMENTED OUT.
  15.     
  16. */
  17.  
  18. /*#define _NOERRORCHECK_*/
  19. #include "math.h"
  20. #include "sane.h"
  21.  
  22. /*  useful constants  */
  23. static double Zero = 0.0;
  24. static double One = 1.0;
  25. static double Two = 2.0;
  26. static double MinusOne = -1.0;
  27. static double MinusTwo = -2.0;
  28. static double PointFive = 0.5;
  29. static double PointTwoFive = 0.25;
  30. static double Pi = PI;
  31. static double Pi2 = PI2;
  32. static double Log2Ten = 3.321928094887362348;
  33.  
  34. static short _Max[] = { 0x7FFE, 0x7FFF, 0xFFFF, 0xFFFF, 0xFFFF };
  35. static short _MinusMax[] = { 0xFFFE, 0x7FFF, 0xFFFF, 0xFFFF, 0xFFFF };
  36. #define Max                (* (double *) _Max)
  37. #define MinusMax        (* (double *) _MinusMax)
  38.  
  39. /*  seed for pseudo-random number generator  */
  40. static unsigned long seed = 1;
  41.  
  42.  
  43. /*  environment word masks  */
  44. #define ROUND            0x6000
  45. #define ROUND_UP        0x2000
  46. #define ROUND_DOWN        0x4000
  47. #define ERROR            0x0F00
  48.  
  49.  
  50. /* ---------- error checking ---------- */
  51.  
  52. #ifdef _ERRORCHECK_
  53.  
  54. #define DomainCheck(test, result)        if (test) {                    \
  55.                                             errno = EDOM;            \
  56.                                             return(result);            \
  57.                                         }
  58. #define RangeCheck(target)                if (GetState() & ERROR) {    \
  59.                                             errno = ERANGE;            \
  60.                                             target = Max;            \
  61.                                         }
  62.  
  63. #else _ERRORCHECK_
  64.  
  65. #define ClearExceptions()
  66. #define DomainCheck(test, result)
  67. #define RangeCheck(target)
  68.  
  69. #endif _ERRORCHECK_
  70.  
  71. /* ---------- environment ---------- */
  72.  
  73.  
  74. /*
  75.  *  GetState - query environment word
  76.  *
  77.  */
  78.  
  79. static
  80. GetState()
  81. {
  82.     int state;
  83.     
  84.     GetEnvironment(&state);
  85.     return(state);
  86. }
  87.  
  88.  
  89. /*
  90.  *  SetState - define environment word
  91.  *
  92.  */
  93.  
  94. static
  95. SetState(state)
  96. {
  97.     SetEnvironment(&state);
  98. }
  99.  
  100.  
  101. /*
  102.  *  ClearExceptions - clear error conditions
  103.  *
  104.  *  This must be called if RangeCheck is to be used.
  105.  *
  106.  */
  107.  
  108. #ifdef _ERRORCHECK_
  109. static
  110. ClearExceptions()
  111. {
  112.     int state;
  113.     
  114.     GetEnvironment(&state);
  115.     state &= ~ERROR;
  116.     SetEnvironment(&state);
  117. }
  118. #endif _ERRORCHECK_
  119.  
  120.  
  121. /*
  122.  *  SetRound - define rounding direction
  123.  *
  124.  *  The previous environment word is returned.  The rounding direction can
  125.  *  be restored by passing this value to SetState.
  126.  *
  127.  */
  128.  
  129. static
  130. SetRound(direction)
  131. {
  132.     int state;
  133.  
  134.     GetEnvironment(&state);
  135.     SetState((state & ~ROUND) | direction);
  136.     return(state);
  137. }
  138.  
  139.  
  140. /*
  141.  *  xfersign - transfer sign from one floating number to another
  142.  *
  143.  */
  144.  
  145. static
  146. xfersign(x, yp)
  147. double x, *yp;
  148. {
  149.     asm {
  150.         movea.l    yp,a0
  151.         bclr    #7,(a0)
  152.         tst.w    x
  153.         bpl.s    @1
  154.         bset    #7,(a0)
  155. @1    }
  156. }
  157.  
  158.  
  159. /* ---------- math functions (alphabetically) ---------- */
  160.  
  161.  
  162. /*
  163.  *  abs - absolute value of an integer
  164.  *
  165.  */
  166.  
  167. int abs(x)
  168. int x;
  169. {
  170.     return(x < 0 ? -x : x);
  171. }
  172.  
  173.  
  174. /*
  175.  *  acos - inverse circular cosine
  176.  *
  177.  */
  178.  
  179. double acos(x)
  180. double x;
  181. {
  182.     DomainCheck(x > One || x < MinusOne, Zero);
  183.     if (x == MinusOne)
  184.         return(Pi);
  185.     return(Two * atan(sqrt((One - x) / (One + x))));
  186. }
  187.  
  188.  
  189. /*
  190.  *  asin - inverse circular sine
  191.  *
  192.  */
  193.  
  194. double asin(x)
  195. double x;
  196. {
  197.     double y = fabs(x);
  198.  
  199.     DomainCheck(y > One, Zero);
  200.     if (y == One) {
  201.         y = Pi2;
  202.         xfersign(x, &y);
  203.         return(y);
  204.     }
  205.     if (y > PointFive) {
  206.         y = One - y;
  207.         y = Two * y - y * y;
  208.     }
  209.     else
  210.         y = One - y * y;
  211.     return(atan(x / sqrt(y)));
  212. }
  213.  
  214.  
  215. /*
  216.  *  atan - inverse circular tangent
  217.  *
  218.  */
  219.  
  220. double atan(x)
  221. double x;
  222. {
  223.     elems68k(&x, FOATANX);
  224.     return(x);
  225. }
  226.  
  227.  
  228. /*
  229.  *  atan2 - inverse circular tangent (y/x)
  230.  *
  231.  */
  232.  
  233. double atan2(y, x)
  234. double y, x;
  235. {
  236.     double z;
  237.  
  238.     if (x == 0) {
  239.         DomainCheck(y == 0, Zero);
  240.         z = Pi2;
  241.     }
  242.     else {
  243.         z = atan(fabs(y / x));
  244.         if (x < 0)
  245.             z = Pi - z;
  246.     }
  247.     xfersign(y, &z);
  248.     return(z);
  249. }
  250.  
  251.  
  252. /*
  253.  *  ceil - round up to an integer
  254.  *
  255.  */
  256.  
  257. double ceil(x)
  258. double x;
  259. {
  260.     int state = SetRound(ROUND_UP);
  261.     fp68k(&x, FORTI);
  262.     SetState(state);
  263.     return(x);
  264. }
  265.  
  266.  
  267. /*
  268.  *  cos - circular cosine
  269.  *
  270.  */
  271.  
  272. double cos(x)
  273. double x;
  274. {    
  275.     elems68k(&x, FOCOSX);
  276.     return(x);
  277. }
  278.  
  279.  
  280. /*
  281.  *  cosh - hyperbolic cosine
  282.  *
  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.  
  296. /*
  297.  *  exp - exponential function
  298.  *
  299.  */
  300.  
  301. double exp(x)
  302. double x;
  303. {
  304.     ClearExceptions();
  305.     elems68k(&x, FOEXPX);
  306.     RangeCheck(x);
  307.     return(x);
  308. }
  309.  
  310.  
  311. /*
  312.  *  fabs - absolute value of a floating number
  313.  *
  314.  */
  315.  
  316. double fabs(x)
  317. double x;
  318. {
  319.     fp68k(&x, FOABS);
  320.     return(x);
  321. }
  322.  
  323.  
  324. /*
  325.  *  floor - round down to an integer
  326.  *
  327.  */
  328.  
  329. double floor(x)
  330. double x;
  331. {
  332.     int state = SetRound(ROUND_DOWN);
  333.     
  334.     fp68k(&x, FORTI);
  335.     SetState(state);
  336.     return(x);
  337. }
  338.  
  339.  
  340. /*
  341.  *  fmod - remainder function
  342.  *
  343.  *  This computes a value z, with the same sign as x, such that for some
  344.  *  integer k, k*y + z == x.
  345.  *
  346.  */
  347.  
  348. double fmod(x, y)
  349. double x, y;
  350. {
  351.     double z = x;
  352.  
  353.     fp68k(&y, FOABS);
  354.     fp68k(&y, &z, FOREM);
  355.     
  356.     if (x > 0 && z < 0)
  357.         z += y;
  358.     else if (x < 0 && z > 0)
  359.         z -= y;
  360.     return(z);
  361. }
  362.  
  363.  
  364. /*
  365.  *  frexp - split floating number into fraction/exponent
  366.  *
  367.  *  This computes a value z, where 0.5 <= fabs(z) < 1.0, and an integer n such
  368.  *  that z*(2^n) == x.
  369.  *
  370.  */
  371.  
  372. double frexp(x, nptr)
  373. double x;
  374. register int *nptr;
  375. {
  376.     double y = fabs(x), z = Two;
  377.     
  378.     if (y == 0) {
  379.         *nptr = 0;
  380.         return(Zero);
  381.     }
  382.     
  383.     elems68k(&y, FOLOG2X);
  384.     
  385.     y -= *nptr = y;
  386.     
  387.     elems68k(&y, &z, FOPWRY);
  388.     
  389.     if (z >= One) {
  390.         z *= PointFive;
  391.         ++*nptr;
  392.     }
  393.     else if (z < PointFive) {
  394.         z += z;
  395.         --*nptr;
  396.     }
  397.     xfersign(x, &z);
  398.     return(z);
  399. }
  400.  
  401.  
  402. /*
  403.  *  labs - absolute value of a long integer
  404.  *
  405.  */
  406.  
  407. long labs(x)
  408. long x;
  409. {
  410.     return(x < 0 ? -x : x);
  411. }
  412.  
  413.  
  414. /*
  415.  *  ldexp - combine fraction/exponent into a floating number
  416.  *
  417.  */
  418.  
  419. double ldexp(x, n)
  420. double x;
  421. int n;
  422. {
  423.     fp68k(&n, &x, FOSCALB);
  424.     return(x);
  425. }
  426.  
  427.  
  428. /*
  429.  *  log - natural logarithm
  430.  *
  431.  */
  432.  
  433. double log(x)
  434. double x;
  435. {
  436.     DomainCheck(x <= 0, MinusMax);
  437.     elems68k(&x, FOLNX);
  438.     return(x);
  439. }
  440.  
  441.  
  442. /*
  443.  *  log10 - logarithm base 10
  444.  *
  445.  */
  446.  
  447. double log10(x)
  448. double x;
  449. {
  450.     DomainCheck(x <= 0, MinusMax);
  451.     elems68k(&x, FOLOG2X);    /* LOG2 is much faster than LN */
  452.     return(x / Log2Ten);
  453. }
  454.  
  455.  
  456. /*
  457.  *  modf - split a floating number into fraction/integer
  458.  *
  459.  */
  460.  
  461. double modf(x, nptr)
  462. double x;
  463. register double *nptr;
  464. {
  465.     *nptr = x;
  466.     fp68k(nptr, FOTTI);
  467.     return(x - *nptr);
  468. }
  469.  
  470.  
  471. /*
  472.  *  pow - power function (exponentiation)
  473.  *
  474.  */
  475.  
  476. double pow(x, y)
  477. double x, y;
  478. {
  479.     int odd = 0;
  480.     double z;
  481.     
  482.     ClearExceptions();
  483.     if (x == 0) {
  484.         if (y <= 0) {
  485. #ifdef _ERRORCHECK_
  486.             errno = EDOM;
  487. #endif
  488.             return (MinusMax);
  489.         }
  490.         return(Zero);
  491.     }
  492.     if (y == 0)
  493.         return(One);
  494.     if (x < 0) {
  495.         if (modf(y, &y)) {
  496. #ifdef _ERRORCHECK_
  497.             errno = EDOM;
  498. #endif
  499.             return (MinusMax);
  500.         }
  501.         x = -x;
  502.         odd = fmod(y, Two);
  503.     }
  504.     elems68k(&y, &x, FOPWRY);
  505.     
  506.     RangeCheck(x);
  507.     return(odd ? -x : x);
  508. }
  509.  
  510. /*
  511.  *  rand - pseudo-random number generator (ANSI C standard)
  512.  *
  513.  */
  514.  
  515. int rand()
  516. {
  517.     seed = seed * 1103515245 + 12345;
  518.     asm {
  519.         move.w    seed,d0            ;  high word of long
  520.         andi.w    #0x7FFF,d0        ;  remove high bit
  521.     }
  522. }
  523.  
  524.  
  525. /*
  526.  *  sin - circular sine
  527.  *
  528.  */
  529.  
  530. double sin(x)
  531. double x;
  532. {
  533.     elems68k(&x, FOSINX);
  534.     return(x);
  535. }
  536.  
  537.  
  538. /*
  539.  *  sinh - hyperbolic sine
  540.  *
  541.  */
  542.  
  543. double sinh(x)
  544. double x;
  545. {
  546.     double y = fabs(x);
  547.     
  548.     ClearExceptions();
  549.     elems68k(&y, FOEXP1X);
  550.     y += y / (y + One);
  551.     y *= PointFive;
  552.     RangeCheck(y);
  553.     xfersign(x, &y);
  554.     return(y);
  555. }
  556.  
  557.  
  558. /*
  559.  *  sqrt - square root
  560.  *
  561.  */
  562.  
  563. double sqrt(x)
  564. double x;
  565. {
  566.     DomainCheck(x < 0, Zero);
  567.     fp68k(&x, FOSQRT);
  568.     return(x);
  569. }
  570.  
  571.  
  572. /*
  573.  *  srand - seed pseudo-random number generator
  574.  *
  575.  */
  576.  
  577. void srand(n)
  578. unsigned n;
  579. {
  580.     seed = n;
  581. }
  582.  
  583.  
  584. /*
  585.  *  tan - circular tangent
  586.  *
  587.  */
  588.  
  589. double tan(x)
  590. double x;
  591. {
  592.     ClearExceptions();
  593.     elems68k(&x, FOTANX);
  594.     RangeCheck(x);
  595.     return(x);
  596. }
  597.  
  598.  
  599. /*
  600.  *  tanh - hyperbolic tangent
  601.  *
  602.  */
  603.  
  604. double tanh(x)
  605. double x;
  606. {
  607.     double y = MinusTwo * fabs(x);
  608.     elems68k(&y, FOEXP1X);
  609.     y = -y / (y + Two);
  610.     xfersign(x, &y);
  611.     return(y);
  612. }
  613.