home *** CD-ROM | disk | FTP | other *** search
/ Tricks of the Mac Game Programming Gurus / TricksOfTheMacGameProgrammingGurus.iso / CodeWarrior Lite / Metrowerks C⁄C++ Lite / Libraries / MacOS 68K / MathLib68K / MathLib68K Sources / fp68k_glue.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-22  |  22.2 KB  |  1,168 lines  |  [TEXT/MMCC]
  1. /*
  2.  *    fp.h glue functions ...
  3.  */
  4.  
  5. #ifndef __FP__
  6. #include <fp.h>
  7. #endif
  8.  
  9. #ifndef __SANE__
  10. #include <SANE.h>
  11. #endif
  12.  
  13. #ifndef __STRINGS__
  14. #include <Strings.h>
  15. #endif
  16.  
  17. #ifndef _ERRNO
  18. #include <errno.h>
  19. #endif
  20.  
  21. #ifndef _STDLIB
  22. #include <stdlib.h>
  23. #endif
  24.  
  25. #define    MINVAL        (*(long double *)(&__C__))        /*    2^-33    =    1.1641532182693481450E-10 */
  26.  
  27. /*
  28.  *    Local functions and globals ...
  29.  */
  30.  
  31. #if __MC68881__
  32.  
  33. /* Trigonometric functions */
  34. static long double _cos(long double:__FP0):__FP0 = { 0xF200,0x001D };
  35. static long double _sin(long double:__FP0):__FP0 = { 0xF200,0x000E };
  36. static long double _tan(long double:__FP0):__FP0 = { 0xF200,0x000F };
  37. static long double _acos(long double:__FP0):__FP0 = { 0xF200,0x001C };
  38. static long double _asin(long double:__FP0):__FP0 = { 0xF200,0x000C };
  39. static long double _atan(long double:__FP0):__FP0 = { 0xF200,0x000A };
  40. /* Hyperbolic functions */
  41. static long double _cosh(long double:__FP0):__FP0 = { 0xF200,0x0019 };
  42. static long double _sinh(long double:__FP0):__FP0 = { 0xF200,0x0002 };
  43. static long double _tanh(long double:__FP0):__FP0 = { 0xF200,0x0009 };
  44. /* Exponential functions */
  45. static long double _exp(long double:__FP0):__FP0 = { 0xF200,0x0010 };
  46. static long double _ldexp(long double:__FP0,long:__D0):__FP0 = { 0xF200,0x4026 };
  47. static long double _log(long double:__FP0):__FP0 = { 0xF200,0x0014 };
  48. static long double _log10(long double:__FP0):__FP0 = { 0xF200,0x0015 };
  49. static long double fgetexp(long double:__FP0):__FP0 = { 0xF200,0x001E };
  50. static long double fgetman(long double:__FP0):__FP0 = { 0xF200,0x001F };
  51. /* Power and absolute value functions */
  52. static long double _fabs(long double:__FP0):__FP0 = { 0xF200,0x0018 };
  53. static long double _sqrt(long double:__FP0):__FP0 = { 0xF200,0x0004 };
  54. /* Nearest integer functions */
  55. static long double fint(long double:__FP0):__FP0 = { 0xF200,0x0001 };
  56. static long double fintrz(long double:__FP0):__FP0 = { 0xF200,0x0003 };
  57.  
  58. static long    FSetRound(long:__D0) = { 0xF200,0x9000 };
  59. static long FGetRound(void):__D0 = { 0xF200,0xB000 };
  60. static void FRoundUp(void)
  61. {
  62.     long l;
  63.  
  64.     l=FGetRound(); l=(l&0xFFFFFFCF)|0x30; FSetRound(l);
  65. }
  66.  
  67. static void FRoundDown(void)
  68. {
  69.     long l;
  70.  
  71.     l=FGetRound(); l=(l&0xFFFFFFCF)|0x20; FSetRound(l);
  72. }
  73. /* Remainder functions */
  74. static long double _fmod(long double:__FP0,long double:__FP1):__FP0 = { 0xF200,0x0421 };
  75.  
  76. static short __C__[6]             =    { 0x3fde,0x0000,0x8000,0x0000,0x0000,0x0000 };
  77. static short __inf__[6]            =    { 0x7FFF,0x0000,0x0000,0x0000,0x0000,0x0000 };
  78.  
  79. #else
  80.  
  81. static short __C__[5]            =    { 0x3fde,0x8000,0x0000,0x0000,0x0000 };
  82. static short __inf__[5]            =    { 0x7FFF,0x0000,0x0000,0x0000,0x0000 };
  83.  
  84. static pascal void MyRemainder(extended80 *,extended80 *) = { 0x3F3C,0x000C,0xA9EB };
  85.  
  86. #endif
  87.  
  88. static void MyCopySign(long double *xp,register long double *yp)
  89. {
  90.     *(short *)yp&=0x7FFF; if(*(short *)xp<0) *(short *)yp|=0x8000;
  91. }
  92.  
  93. /*******************************************************************************
  94. *                                Constants                                     *
  95. *******************************************************************************/
  96.  
  97. const double_t pi = 3.14159265358979323846;
  98. static const double_t log2_10 = 3.321928094887362348;
  99.  
  100. /*******************************************************************************
  101. *                            Trigonometric functions                           *
  102. *******************************************************************************/
  103.  
  104. double_t cos ( double_t x )
  105. {
  106. #if __MC68881__
  107.     return(_cos(x));
  108. #else    /* SANE */
  109.     extended80 x80=x;
  110.  
  111.     Cos(&x80);
  112.  
  113.     return(x80);
  114. #endif
  115. }
  116.  
  117. double_t sin ( double_t x )
  118. {
  119. #if __MC68881__
  120.     return(_sin(x));
  121. #else    /* SANE */
  122.     extended80 x80=x;
  123.  
  124.     Sin(&x80);
  125.  
  126.     return(x80);
  127. #endif
  128. }
  129.  
  130. double_t tan ( double_t x )
  131. {
  132. #if __MC68881__
  133.     return(_tan(x));
  134. #else    /* SANE */
  135.     extended80 x80=x;
  136.  
  137.     Tan(&x80);
  138.  
  139.     return(x80);
  140. #endif
  141. }
  142.  
  143. double_t acos ( double_t x )
  144. {
  145. #if __MC68881__
  146.     if(x<-1 || x>1) { errno=EDOM; return(0); }
  147.     return(_acos(x));
  148. #else    /* SANE */
  149.     if(x<-1 || x>1) { errno=EDOM; return(0); }
  150.     return(2*atan(sqrt((1-x)/(1+x))));
  151. #endif
  152. }
  153.  
  154. double_t asin ( double_t x )
  155. {
  156. #if __MC68881__
  157.     if(x<-1 || x>1) { errno=EDOM; return(0); }
  158.     return(_asin(x));
  159. #else    /* SANE */
  160.     extended80 y;
  161.  
  162.     if (x<-1 || x>1) { errno=EDOM; return(0); }
  163.     y=x; Abs(&y); if (y<=MINVAL) return(x);
  164.     if (y>0.5) { y=1-y; y=2*y-y*y; } else y=1-y*y; 
  165.  
  166.     return(atan(x/sqrt(y)));
  167. #endif
  168. }
  169.  
  170. double_t atan ( double_t x )
  171. {
  172. #if __MC68881__
  173.     return(_atan(x));
  174. #else    /* SANE */
  175.     extended80 x80=x;
  176.  
  177.     Atan(&x80);
  178.  
  179.     return(x80);
  180. #endif
  181. }
  182.  
  183. double_t atan2 ( double_t y, double_t x )
  184. {
  185.     double_t z;
  186.  
  187.     if (x==-0.0) x=0.0;
  188.     z=atan(y/x);
  189.     if(x<0) { if(y<0) z-=pi; else z+=pi; }
  190.  
  191.     return(z);
  192. }
  193.  
  194. /*******************************************************************************
  195. *                              Hyperbolic functions                            *
  196. *******************************************************************************/
  197.  
  198. double_t cosh ( double_t x )
  199. {
  200. #if __MC68881__
  201.     return(_cosh(x));
  202. #else    /* SANE */
  203.     extended80 y=x;
  204.  
  205.     y=0.5*exp(fabs(y));
  206.  
  207.     return(y+0.25/y);
  208. #endif
  209. }
  210.  
  211. double_t sinh ( double_t x )
  212. {
  213. #if __MC68881__
  214.     return(_sinh(x));
  215. #else    /* SANE */
  216.     extended80 y, x80=x;
  217.  
  218.     y=x; Abs(&y); if(y<=MINVAL) return(x);
  219.     Exp1(&y); y+=y/(1+y); y*=0.5;
  220.     MyCopySign(&x80,&y);
  221.  
  222.     return(y);
  223. #endif
  224. }
  225.  
  226. double_t tanh ( double_t x )
  227. {
  228. #if __MC68881__
  229.     return(_tanh(x));
  230. #else    /* SANE */
  231.     extended80 y,x80=x;
  232.  
  233.     y=x; Abs(&y);
  234.     if(y>MINVAL) { y*=-2; Exp1(&y); y=-y/(2+y); }
  235.     MyCopySign(&x80,&y);
  236.  
  237.     return(y);
  238. #endif
  239. }
  240.  
  241. //  asinh, acosh, atanh are taken from Apple Numerics Manual, 2nd Ed., p74.
  242.  
  243. double_t asinh(double_t x)
  244. {
  245.     double_t y;
  246.     
  247.     y = fabs(x);
  248.     if (y > MINVAL) y = log1p(y + y/(1.0/y + sqrt(1 + 1/(y*y))));
  249.     MyCopySign(&x,&y);
  250.     return(y);
  251. }
  252.  
  253. double_t acosh(double_t x)
  254. {
  255.     double_t y;
  256.     
  257.     y = sqrt(x - 1.0);
  258.     return log1p(y*(y + sqrt(x + 1.0)));
  259. }
  260.  
  261. double_t atanh(double_t x)
  262. {
  263.     double_t y;
  264.     
  265.     y = fabs(x);
  266.     if (y > 1.164e-10) y = 0.5*log1p(2.0*(y/(1.0 - y)));
  267.     MyCopySign(&x,&y);
  268.     return(y);
  269. }
  270.  
  271. /*******************************************************************************
  272. *                              Exponential functions                           *
  273. *******************************************************************************/
  274.  
  275. double_t exp ( double_t x )
  276. {
  277. #if __MC68881__
  278.     return(_exp(x));
  279. #else    /* SANE */
  280.     extended80 x80=x;
  281.  
  282.     Exp(&x80);
  283.  
  284.     return(x80);
  285. #endif
  286. }
  287.  
  288. double_t expm1  ( double_t x )
  289. {
  290.     extended80 x80;
  291. #if __MC68881__
  292.     extended96 x96=x;
  293.     x96tox80(&x96,&x80);
  294. #else    /* SANE */
  295.     x80=x;
  296. #endif
  297.  
  298.     Exp1(&x80);
  299.  
  300. #if __MC68881__
  301.     x80tox96(&x80,&x96); return(x96);
  302. #else    /* SANE */
  303.     return(x80);
  304. #endif
  305. }
  306.  
  307. double_t exp2  ( double_t x )
  308. {
  309.     extended80 x80;
  310. #if __MC68881__
  311.     extended96 x96=x;
  312.     x96tox80(&x96,&x80);
  313. #else    /* SANE */
  314.     x80=x;
  315. #endif
  316.  
  317.     Exp2(&x80);
  318.  
  319. #if __MC68881__
  320.     x80tox96(&x80,&x96); return(x96);
  321. #else    /* SANE */
  322.     return(x80);
  323. #endif
  324. }
  325.  
  326. double_t frexp ( double_t x, int *exponent )
  327. {
  328. #if __MC68881__
  329.     if(x==0) { *exponent=0; return(0); }
  330.     *exponent=(int)fgetexp(x)+1; return(fgetman(x)/2);
  331. #else    /* SANE */
  332.     long n;
  333.     extended80 x80=x, y;
  334.  
  335.     if(x == 0) { *exponent=0; return(0); }
  336.     y=fabs(x); Log2(&y); n=y; y-=n; y=pow(2,y);
  337.     if(y>=1) { y=y/2; n++; } else if(y<0.5 && y!=0.0) { y+=y; n--; }
  338.     *exponent=n;
  339.     MyCopySign(&x80,&y);
  340.  
  341.     return(y);
  342. #endif
  343. }
  344.  
  345. double_t ldexp ( double_t x, int n )
  346. {
  347. #if __MC68881__
  348.     return(_ldexp(x,n));
  349. #else    /* SANE */
  350.     extended x80=x;
  351.     short sn=n;
  352.  
  353.     Scalb(&sn,&x80);
  354.  
  355.     return(x80);
  356. #endif
  357. }
  358.  
  359. double_t log ( double_t x )
  360. {
  361. #if __MC68881__
  362.     if(x<0) errno=EDOM;
  363.     return(_log(x));
  364. #else    /* SANE */
  365.     extended80 x80=x;
  366.  
  367.     if(x<0) errno=EDOM;
  368.     Ln(&x80);
  369.  
  370.     return(x80);
  371. #endif
  372. }
  373.  
  374. double_t log2 ( double_t x )
  375. {
  376.     extended80 x80;
  377. #if __MC68881__
  378.     extended96 x96=x;
  379.     x96tox80(&x96,&x80);
  380. #else    /* SANE */
  381.     x80=x;
  382. #endif
  383.  
  384.     Log2(&x80);
  385.  
  386. #if __MC68881__
  387.     x80tox96(&x80,&x96); return(x96);
  388. #else    /* SANE */
  389.     return(x80);
  390. #endif
  391. }
  392.  
  393. double_t log1p ( double_t x )
  394. {
  395.     extended80 x80;
  396. #if __MC68881__
  397.     extended96 x96=x;
  398.     x96tox80(&x96,&x80);
  399. #else    /* SANE */
  400.     x80=x;
  401. #endif
  402.  
  403.     Ln1(&x80);
  404.  
  405. #if __MC68881__
  406.     x80tox96(&x80,&x96); return(x96);
  407. #else    /* SANE */
  408.     return(x80);
  409. #endif
  410. }
  411.  
  412. double_t log10 ( double_t x )
  413. {
  414. #if __MC68881__
  415.     if(x<0) errno=EDOM;
  416.     return(_log10(x));
  417. #else    /* SANE */
  418.     extended80 x80=x;
  419.     if(x<0) errno=EDOM;
  420.     Log2(&x80); x80/=log2_10;
  421.     return(x80);
  422. #endif
  423. }
  424.  
  425. double_t logb ( double_t x )
  426. {
  427.     extended80 x80;
  428. #if __MC68881__
  429.     extended96 x96=x;
  430.     x96tox80(&x96,&x80);
  431. #else    /* SANE */
  432.     x80=x;
  433. #endif
  434.  
  435.     Logb(&x80);
  436.  
  437. #if __MC68881__
  438.     x80tox96(&x80,&x96); return(x96);
  439. #else    /* SANE */
  440.     return(x80);
  441. #endif
  442. }
  443.  
  444. long double modfl ( long double x, long double *iptrl )
  445. {
  446.     return(modf(x, iptrl));
  447. }
  448.  
  449. double_t    modf  ( double_t x, double_t *iptr )
  450. {
  451. #if __MC68881__
  452.     *iptr=fintrz(x); return(x-*iptr);
  453. #else    /* SANE */
  454.     extended80    x80=x;
  455.  
  456.     Tint(&x80); *iptr=x80;
  457.  
  458.     return(x-x80);
  459. #endif
  460. }
  461.  
  462. float       modff ( float x, float *iptrf )
  463. {
  464.     double_t x_flt, result;
  465.  
  466.     result  = modf(x, &x_flt);
  467.     *iptrf = x_flt;
  468.     return(result);
  469. }
  470.  
  471. double_t scalb ( double_t x, long int n )
  472. {
  473.     short scale;
  474.     extended80 x80;
  475. #if __MC68881__
  476.     extended96 x96=x;
  477.     x96tox80(&x96,&x80);
  478. #else    /* SANE */
  479.     x80=x;
  480. #endif
  481.  
  482.     if (n > 32767) scale =  32767;
  483.     else if (n < -32768) scale = -32768;
  484.     else scale = n;
  485.  
  486.     Scalb(&scale,&x80);
  487.  
  488. #if __MC68881__
  489.     x80tox96(&x80,&x96); return(x96);
  490. #else
  491.     return(x80);
  492. #endif
  493. }
  494.  
  495. /*******************************************************************************
  496. *                     Power and absolute value functions                       *
  497. *******************************************************************************/
  498.  
  499. double_t fabs ( double_t x )
  500. {
  501. #if __MC68881__
  502.     return(_fabs(x));
  503. #else    /* SANE */
  504.     extended80 x80=x;
  505.  
  506.     Abs(&x80); return(x80);
  507. #endif
  508. }
  509.  
  510. double_t hypot ( double_t x, double_t y )
  511. {
  512.     long int fpclassx = __fpclassify(x);
  513.     long int fpclassy = __fpclassify(y);
  514.  
  515.     if (fpclassx == FP_SNAN || fpclassx == FP_QNAN ||
  516.         fpclassy == FP_SNAN || fpclassy == FP_QNAN)
  517.         return(x + y);
  518.  
  519.     return(sqrt(x*x + y*y));
  520. }
  521.  
  522. double_t pow   ( double_t x, double_t y )
  523. {
  524.     extended80    x80,y80;
  525.     double_t iptr;
  526. #if __MC68881__
  527.     extended96    x96=x, y96=y;
  528.     x96tox80(&x96,&x80);
  529.     x96tox80(&y96,&y80);
  530. #else    /* SANE */
  531.     x80=x;y80=y;
  532. #endif
  533.  
  534.     if(x==0) { if(y<=0) errno=EDOM; return(0); }
  535.     if(y==0) return(1);
  536.     if((x<0) && modf(y,&iptr)) { errno=EDOM; return(0);}
  537.  
  538.     Power(&y80,&x80);
  539.  
  540. #if __MC68881__
  541.     x80tox96(&x80,&x96); return(x96);
  542. #else    /* SANE */
  543.     return(x80);
  544. #endif
  545. }
  546.  
  547. double_t sqrt  ( double_t x )
  548. {
  549. #if __MC68881__
  550.     if (x < 0) { errno=EDOM; return(0); }
  551.     return(_sqrt(x));
  552. #else    /* SANE */
  553.     extended80 x80=x;
  554.  
  555.     if (x < 0) { errno=EDOM; return(0); }
  556.     Sqrt(&x80); return(x80);
  557. #endif
  558. }
  559.  
  560. /*******************************************************************************
  561. *                        Gamma and Error functions                             *
  562. *******************************************************************************/
  563.  
  564. #if __MC68881__
  565. extern double_t __FPUerf(double_t);
  566. extern double_t __FPUerfc(double_t);
  567. extern double_t __FPUgamma(double_t);
  568. extern double_t __FPUlgamma(double_t);
  569. #else
  570. extern double_t __nonFPUerf(double_t);
  571. extern double_t __nonFPUerfc(double_t);
  572. extern double_t __nonFPUgamma(double_t);
  573. extern double_t __nonFPUlgamma(double_t);
  574. #endif
  575.  
  576. double_t erf  ( double_t x )
  577. {
  578. #if __MC68881__
  579.     return(__FPUerf(x));
  580. #else
  581.     return(__nonFPUerf(x));
  582. #endif
  583. }
  584.  
  585. double_t erfc ( double_t x )
  586. {
  587. #if __MC68881__
  588.     return(__FPUerfc(x));
  589. #else
  590.     return(__nonFPUerfc(x));
  591. #endif
  592. }
  593.  
  594. double_t gamma ( double_t x )
  595. {
  596. #if __MC68881__
  597.     return(__FPUgamma(x));
  598. #else
  599.     return(__nonFPUgamma(x));
  600. #endif
  601. }
  602.  
  603. double_t lgamma ( double_t x )
  604. {
  605. #if __MC68881__
  606.     return(__FPUlgamma(x));
  607. #else
  608.     return(__nonFPUlgamma(x));
  609. #endif
  610. }
  611.  
  612. /*******************************************************************************
  613. *                        Nearest integer functions                             *
  614. *******************************************************************************/
  615.  
  616. double_t ceil  ( double_t x )
  617. {
  618. #if __MC68881__
  619.     long round;
  620.     double_t temp;
  621.  
  622.     round=FGetRound(); FRoundUp();
  623.     temp=fint(x);
  624.     FSetRound(round); return(temp);
  625. #else    /* SANE */
  626.     Environment    env_old,env_new;
  627.     extended80 x80=x;
  628.  
  629.     GetEnvironment(&env_old);
  630.     env_new=0x2000;        /*    Round Upward */
  631.     SetEnvironment(&env_new);
  632.  
  633.     Rint(&x80);
  634.  
  635.     SetEnvironment(&env_old);
  636.     return(x80);
  637. #endif
  638. }
  639.  
  640. double_t floor ( double_t x )
  641. {
  642. #if __MC68881__
  643.     long round;
  644.     double_t temp;
  645.  
  646.     round=FGetRound(); FRoundDown();
  647.     temp=fint(x);
  648.     FSetRound(round); return(temp);
  649. #else    /* SANE */
  650.     Environment    env_old,env_new;
  651.     extended80 x80=x;
  652.  
  653.     GetEnvironment(&env_old);
  654.     env_new=0x4000;                /*    Round Downward */
  655.     SetEnvironment(&env_new);
  656.  
  657.     Rint(&x80);
  658.  
  659.     SetEnvironment(&env_old);
  660.     return(x80);
  661. #endif
  662. }
  663.  
  664. double_t rint ( double_t x )
  665. {
  666.     extended80 x80;
  667. #if __MC68881__
  668.     extended96 x96=x;
  669.     x96tox80(&x96,&x80);
  670. #else    /* SANE */
  671.     x80=x;
  672. #endif
  673.  
  674.     Rint(&x80);
  675.  
  676. #if __MC68881__
  677.     x80tox96(&x80,&x96); return(x96);
  678. #else    /* SANE */
  679.     return(x80);
  680. #endif
  681. }
  682.  
  683.  
  684. double_t nearbyint ( double_t x )
  685. {
  686.     return(rint(x));
  687. }
  688.  
  689. long int rinttol ( double_t x )
  690. {
  691.     return((long)rint(x));
  692. }
  693.  
  694. double_t round ( double_t x )
  695. {
  696.     short env=0;
  697.     double_t temp,y;
  698.  
  699.     SetEnvironment((Environment *)&env);
  700.     temp=rint(x);
  701.     GetEnvironment((Environment *)&env);
  702.     if (env & (((short)16) << 8)) {
  703.         y=0.5+fabs(x);
  704.         MyCopySign(&x,&y);
  705.         temp=floor(y);
  706.     }
  707.     return(temp);
  708. }
  709.  
  710. long int roundtol ( double_t x )
  711. {
  712.     return((long)round(x));
  713. }
  714.  
  715. int      trunc ( double_t x )
  716. {
  717.     return((int)x);
  718. }
  719.  
  720. /*******************************************************************************
  721. *                            Remainder functions                               *
  722. *******************************************************************************/
  723.  
  724. double_t fmod ( double_t x, double_t y )
  725. {
  726. #if __MC68881__
  727.     return((_fmod)(x,y));
  728. #else    /* SANE */
  729.     extended80 z,x80=x,y80=y;
  730.  
  731.     Abs(&y80); z=x80; MyRemainder(&y80,&z);
  732.     if(x80>0 && z<0) z+=y80;
  733.     else if(x80<0 && z>0) z-=y80;
  734.     return(z);
  735. #endif
  736. }
  737.  
  738. double_t remainder ( double_t x, double_t y )
  739. {
  740.     short quotient;
  741.     extended80 x80, y80;
  742. #if __MC68881__
  743.     extended96 x96=x, y96=y;
  744.     x96tox80(&x96,&x80);
  745.     x96tox80(&y96,&y80);
  746. #else    /* SANE */
  747.     x80=x; y80=y;
  748. #endif
  749.  
  750.     Remainder(&y80,&x80,"ient);
  751.  
  752. #if __MC68881__
  753.     x80tox96(&x80,&x96); return(x96);
  754. #else    /* SANE */
  755.     return(x80);
  756. #endif
  757. }
  758.  
  759. double_t remquo    ( double_t x, double_t y, int *quo )
  760. {
  761.     short quotient;
  762.     extended80 x80, y80;
  763. #if __MC68881__
  764.     extended96 x96=x, y96=y;
  765.     x96tox80(&x96,&x80);
  766.     x96tox80(&y96,&y80);
  767. #else    /* SANE */
  768.     x80=x; y80=y;
  769. #endif
  770.  
  771.     Remainder(&y80,&x80,"ient);
  772.  
  773.     /* Fill in the quo parameter */
  774.     if (quotient>=0)        /* return only lower 7 quotient bits (????) */
  775.         *quo = ((long)quotient) & 0x7fL;
  776.     else
  777.         *quo = -((-((long)quotient)) & 0x7fL);
  778.  
  779. #if __MC68881__
  780.     x80tox96(&x80,&x96); return(x96);
  781. #else    /* SANE */
  782.     return(x80);
  783. #endif
  784. }
  785.  
  786. /*******************************************************************************
  787. *                             Auxiliary functions                              *
  788. *******************************************************************************/
  789.  
  790. double_t copysign ( double_t x, double_t y )
  791. {
  792.     MyCopySign(&y,&x); return(x);
  793. }
  794.  
  795. long double nanl ( const char *tagp )
  796. {
  797.     return(NaN(strtol(tagp, NULL, 0)));
  798. }
  799.  
  800. double      nan  ( const char *tagp )
  801. {
  802.     return(nanl(tagp));
  803. }
  804.  
  805. float       nanf ( const char *tagp )
  806. {
  807.     return(nanl(tagp));
  808. }
  809.  
  810. long double nextafterl ( long double x, long double y )
  811. {
  812.     extended80 x80, y80;
  813. #if __MC68881__
  814.     extended x96=x, y96=y;
  815.     x96tox80(&x96,&x80);
  816.     x96tox80(&y96,&y80);
  817. #else    /* SANE */
  818.     x80=x; y80=y;
  819. #endif
  820.  
  821.     NextExtended(&x80,&y80);
  822.  
  823. #if __MC68881__
  824.     x80tox96(&y80,&y96); return(y96);
  825. #else    /* SANE */
  826.     return(y80);
  827. #endif
  828. }
  829.  
  830. double      nextafterd ( double x, double y )
  831. {
  832.     double temp=x;
  833.  
  834.     NextDouble(&temp,&y);
  835.  
  836.     return(temp);
  837. }
  838.  
  839. float       nextafterf ( float x, float y )
  840. {
  841.     float temp=x;
  842.  
  843.     NextFloat(&temp,&y);
  844.  
  845.     return(temp);
  846. }
  847.  
  848. /*******************************************************************************
  849. *                      Max, Min and Positive Difference                        *
  850. *******************************************************************************/
  851.  
  852. double_t fdim ( double_t x, double_t y )
  853. {
  854.     if (__isnan(x)) return(x);
  855.     if (__isnan(y)) return(y);
  856.     if (x > y)  return(x-y);  else return(+0.0);
  857. }
  858.  
  859. double_t fmax ( double_t x, double_t y )
  860. {
  861.     if (__isnan(x) && __isnan(y)) return(x);
  862.     if (__isnan(x)) return(y);
  863.     else if (__isnan(y)) return(x);
  864.  
  865.     return((x>=y) ? x : y);
  866. }
  867.  
  868. double_t fmin ( double_t x, double_t y )
  869. {
  870.     if (__isnan(x) && __isnan(y)) return(x);
  871.     if (__isnan(x)) return(y);
  872.     else if (__isnan(y)) return(x);
  873.  
  874.     return((x<=y) ? x : y);
  875. }
  876.  
  877. /*******************************************************************************
  878. *                              Internal prototypes                             *
  879. *******************************************************************************/
  880.  
  881. long int __fpclassify  ( long double x )
  882. {
  883.     NumClass fpclass;
  884.     extended80 x80;
  885. #if __MC68881__
  886.     extended96 x96=x;
  887.     x96tox80(&x96,&x80);
  888. #else    /* SANE */
  889.     x80=x;
  890. #endif
  891.  
  892.     ClassExtended(&x80,&fpclass);
  893.     if (fpclass >= 0)
  894.         return(fpclass-1);        /* SANE is one off from fp */
  895.     else
  896.         return(-(fpclass+1));    /* SANE is one off from fp */
  897. }
  898.  
  899. long int __fpclassifyd ( double x )
  900. {
  901.     NumClass fpclass;
  902.     double temp=x;
  903.  
  904.     ClassDouble(&temp,&fpclass);
  905.     if (fpclass >= 0)
  906.         return(fpclass-1);        /* SANE is one off from fp */
  907.     else
  908.         return(-(fpclass+1));    /* SANE is one off from fp */
  909. }
  910.  
  911. long int __fpclassifyf ( float x )
  912. {
  913.     NumClass fpclass;
  914.     float temp=x;
  915.  
  916.     ClassFloat(&temp,&fpclass);
  917.     if (fpclass >= 0)
  918.         return(fpclass-1);        /* SANE is one off from fp */
  919.     else
  920.         return(-(fpclass+1));    /* SANE is one off from fp */
  921. }
  922.  
  923. long int __isnormal  ( long double x )
  924. {
  925.     return(__fpclassify(x) == FP_NORMAL);
  926. }
  927.  
  928. long int __isnormald ( double x )
  929. {
  930.     return(__fpclassifyd(x) == FP_NORMAL);
  931. }
  932.  
  933. long int __isnormalf ( float x )
  934. {
  935.     return(__fpclassifyf(x) == FP_NORMAL);
  936. }
  937.  
  938. long int __isfinite  ( long double x )
  939. {
  940.     long int fpclass = __fpclassify(x);
  941.  
  942.     return ((fpclass==FP_ZERO) ||
  943.             (fpclass==FP_NORMAL) ||
  944.             (fpclass==FP_SUBNORMAL));
  945. }
  946.  
  947. long int __isfinited ( double x )
  948. {
  949.     long int fpclass = __fpclassifyd(x);
  950.  
  951.     return ((fpclass==FP_ZERO) ||
  952.             (fpclass==FP_NORMAL) ||
  953.             (fpclass==FP_SUBNORMAL));
  954. }
  955.  
  956. long int __isfinitef ( float x )
  957. {
  958.     long int fpclass = __fpclassifyf(x);
  959.  
  960.     return ((fpclass==FP_ZERO) ||
  961.             (fpclass==FP_NORMAL) ||
  962.             (fpclass==FP_SUBNORMAL));
  963. }
  964.  
  965. long int __isnan  ( long double x )
  966. {
  967.     long int fpclass = __fpclassify(x);
  968.  
  969.     return((fpclass==FP_SNAN) || (fpclass==FP_QNAN));
  970. }
  971.  
  972. long int __isnand ( double x )
  973. {
  974.     long int fpclass = __fpclassifyd(x);
  975.  
  976.     return((fpclass==FP_SNAN) || (fpclass==FP_QNAN));
  977. }
  978.  
  979. long int __isnanf ( float x )
  980. {
  981.     long int fpclass = __fpclassifyf(x);
  982.  
  983.     return((fpclass==FP_SNAN) || (fpclass==FP_QNAN));
  984. }
  985.  
  986. long int __signbit  ( long double x )
  987. {
  988.     extended80 x80;
  989. #if __MC68881__
  990.     extended x96=x;
  991.     x96tox80(&x96,&x80);
  992. #else
  993.     x80=x;
  994. #endif
  995.  
  996.     return(SignNum(x80));
  997. }
  998.  
  999. long int __signbitd ( double x )
  1000. {
  1001.     return(__signbit(x));
  1002. }
  1003.  
  1004. long int __signbitf ( float x )
  1005. {
  1006.     return(__signbit(x));
  1007. }
  1008.  
  1009. double_t __inf ( void )
  1010. {
  1011.     return(*(double_t *)__inf__);
  1012. }
  1013.  
  1014. /*******************************************************************************
  1015. *                              Non NCEG extensions                             *
  1016. *******************************************************************************/
  1017.  
  1018. /*******************************************************************************
  1019. *                              Financial functions                             *
  1020. *******************************************************************************/
  1021.  
  1022. double_t compound ( double_t rate, double_t periods )
  1023. {
  1024.     extended80 rate80, periods80, result80;
  1025. #if __MC68881__
  1026.     extended96 rate96=rate, periods96=periods, result96;
  1027.     x96tox80(&rate96,&rate80);
  1028.     x96tox80(&periods96,&periods80);
  1029. #else
  1030.     rate80=rate, periods80=periods;
  1031. #endif
  1032.  
  1033.     Compound(&rate80,&periods80,&result80);
  1034.  
  1035. #if __MC68881__
  1036.     x80tox96(&result80,&result96); return(result96);
  1037. #else
  1038.     return(result80);
  1039. #endif
  1040. }
  1041.  
  1042. double_t annuity ( double_t rate, double_t periods )
  1043. {
  1044.     extended80 rate80, periods80, result80;
  1045. #if __MC68881__
  1046.     extended96 rate96=rate, periods96=periods, result96;
  1047.     x96tox80(&rate96,&rate80);
  1048.     x96tox80(&periods96,&periods80);
  1049. #else
  1050.     rate80=rate, periods80=periods;
  1051. #endif
  1052.  
  1053.     Annuity(&rate80, &periods80, &result80);
  1054.  
  1055. #if __MC68881__
  1056.     x80tox96(&result80,&result96); return(result96);
  1057. #else
  1058.     return(result80);
  1059. #endif
  1060. }
  1061.  
  1062. /*******************************************************************************
  1063. *                              Random function                                 *
  1064. *******************************************************************************/
  1065.  
  1066. double_t randomx ( double_t *x )
  1067. {
  1068.     extended80 x80;
  1069. #if __MC68881__
  1070.     extended96 x96=*x;
  1071.     x96tox80(&x96,&x80);
  1072. #else
  1073.     x80=*x;
  1074. #endif
  1075.  
  1076.     RandomX(&x80);
  1077.  
  1078. #if __MC68881__
  1079.     x80tox96(&x80,&x96); return(x96);
  1080. #else
  1081.     return(x80);
  1082. #endif
  1083. }
  1084.  
  1085. /*******************************************************************************
  1086. *                              Relational operator                             *
  1087. *******************************************************************************/
  1088.  
  1089. relop relation ( double_t x, double_t y )
  1090. {
  1091.     extended80 x80, y80;
  1092. #if __MC68881__
  1093.     extended96 x96=x, y96=y;
  1094.     x96tox80(&x96,&x80);
  1095.     x96tox80(&y96,&y80);
  1096. #else
  1097.     x80=x, y80=y;
  1098. #endif
  1099.  
  1100.     return(Relation(&y80, &x80));
  1101. }
  1102.  
  1103. /*******************************************************************************
  1104. *                         Binary to decimal conversions                        *
  1105. *******************************************************************************/
  1106.  
  1107. void num2dec ( const decform *f, double_t x, decimal *d )
  1108. {
  1109. #if __MC68881__
  1110.     extended80 x80;
  1111.     extended96 x96=x;
  1112.     x96tox80(&x96,&x80);
  1113. #else
  1114.     extended80 x80=x;
  1115. #endif    
  1116.  
  1117.     Num2Dec((DecForm *)f, &x80, (Decimal *)d);
  1118. }
  1119.  
  1120. double_t dec2num ( const decimal *d )
  1121. {
  1122.     extended80 x80;
  1123. #if __MC68881__
  1124.     extended96 x96;
  1125. #endif
  1126.  
  1127.     Dec2Num((Decimal *)d,&x80);
  1128.  
  1129. #if __MC68881__
  1130.     x80tox96(&x80,&x96); return(x96);
  1131. #else
  1132.     return(x80);
  1133. #endif
  1134. }
  1135.  
  1136. void dec2str ( const decform *f, const decimal *d, char *s )
  1137. {
  1138.     Dec2Str((DecForm *)f,(Decimal *)d,(StringPtr)s);
  1139.     p2cstr((unsigned char *)s);
  1140. }
  1141.  
  1142. void str2dec ( const char *s, short *ix, decimal *d, short *vp )
  1143. {
  1144.     *vp = *vp << 8;    /* since pascal puts the boolean byte in the upper byte */
  1145.     CStr2Dec((char *)s,ix,(Decimal *)d,(Boolean *)vp);
  1146.     *vp = *vp >> 8;    /* since pascal puts the boolean byte in the upper byte */
  1147. }
  1148.  
  1149. double dec2d ( const decimal *d )
  1150. {
  1151.     return((double)dec2num(d));
  1152. }
  1153.  
  1154. float dec2f ( const decimal *d )
  1155. {
  1156.     return((float)dec2num(d));
  1157. }
  1158.  
  1159. short int dec2s ( const decimal *d )
  1160. {
  1161.     return((short)dec2num(d));
  1162. }
  1163.  
  1164. long int dec2l  ( const decimal *d )
  1165. {
  1166.     return((long)dec2num(d));
  1167. }
  1168.