home *** CD-ROM | disk | FTP | other *** search
/ Education Sampler 1992 [NeXTSTEP] / Education_1992_Sampler.iso / NeXT / GnuSource / cc-61.0.1 / cc / math-68881.h < prev    next >
C/C++ Source or Header  |  1991-05-02  |  9KB  |  483 lines
  1. /******************************************************************\
  2. *                                   *
  3. *  <math-68881.h>        last modified: 18 May 1989.       *
  4. *                                   *
  5. *  Copyright (C) 1989 by Matthew Self.                   *
  6. *  You may freely distribute verbatim copies of this software       *
  7. *  provided that this copyright notice is retained in all copies.  *
  8. *  You may distribute modifications to this software under the     *
  9. *  conditions above if you also clearly note such modifications    *
  10. *  with their author and date.                               *
  11. *                                   *
  12. *  Note:  errno is not set to EDOM when domain errors occur for    *
  13. *  most of these functions.  Rather, it is assumed that the       *
  14. *  68881's OPERR exception will be enabled and handled           *
  15. *  appropriately by the    operating system.  Similarly, overflow       *
  16. *  and underflow do not set errno to ERANGE.               *
  17. *                                   *
  18. *  Send bugs to Matthew Self (self@bayes.arc.nasa.gov).           *
  19. *                                   *
  20. \******************************************************************/
  21.  
  22. /* Changed by Richard Stallman: % inserted before a #.
  23.    New function `hypot' added.
  24.    Nans written in hex to avoid 0rnan.
  25.    December 1989, add parens around `&' in pow.
  26.    November 1990, added alternate definition of HUGE_VAL for Sun.  */
  27.  
  28. #include <errno.h>
  29.  
  30. #ifndef HUGE_VAL
  31. #ifdef __sun__
  32. /* The Sun assembler fails to handle the hex constant in the usual defn.  */
  33. #define HUGE_VAL                            \
  34. ({                                    \
  35.   static union { int i[2]; double d; } u = { {0x7ff00000, 0} };        \
  36.   u.d;                                    \
  37. })
  38. #else
  39. #define HUGE_VAL                            \
  40. ({                                    \
  41.   double huge_val;                            \
  42.                                     \
  43.   __asm ("fmove%.d %#0x7ff0000000000000,%0"    /* Infinity */        \
  44.      : "=f" (huge_val)                        \
  45.      : /* no inputs */);                        \
  46.   huge_val;                                \
  47. })
  48. #endif
  49. #endif
  50.  
  51. __inline static const double sin (double x)
  52. {
  53.   double value;
  54.  
  55.   __asm ("fsin%.x %1,%0"
  56.      : "=f" (value)
  57.      : "f" (x));
  58.   return value;
  59. }
  60.  
  61. __inline static const double cos (double x)
  62. {
  63.   double value;
  64.  
  65.   __asm ("fcos%.x %1,%0"
  66.      : "=f" (value)
  67.      : "f" (x));
  68.   return value;
  69. }
  70.  
  71. __inline static const double tan (double x)
  72. {
  73.   double value;
  74.  
  75.   __asm ("ftan%.x %1,%0"
  76.      : "=f" (value)
  77.      : "f" (x));
  78.   return value;
  79. }
  80.  
  81. __inline static const double asin (double x)
  82. {
  83.   double value;
  84.  
  85.   __asm ("fasin%.x %1,%0"
  86.      : "=f" (value)
  87.      : "f" (x));
  88.   return value;
  89. }
  90.  
  91. __inline static const double acos (double x)
  92. {
  93.   double value;
  94.  
  95.   __asm ("facos%.x %1,%0"
  96.      : "=f" (value)
  97.      : "f" (x));
  98.   return value;
  99. }
  100.  
  101. __inline static const double atan (double x)
  102. {
  103.   double value;
  104.  
  105.   __asm ("fatan%.x %1,%0"
  106.      : "=f" (value)
  107.      : "f" (x));
  108.   return value;
  109. }
  110.  
  111. __inline static const double atan2 (double y, double x)
  112. {
  113.   double pi, pi_over_2;
  114.  
  115.   __asm ("fmovecr%.x %#0,%0"        /* extended precision pi */
  116.      : "=f" (pi)
  117.      : /* no inputs */ );
  118.   __asm ("fscale%.b %#-1,%0"        /* no loss of accuracy */
  119.      : "=f" (pi_over_2)
  120.      : "0" (pi));
  121.   if (x > 0)
  122.     {
  123.       if (y > 0)
  124.     {
  125.       if (x > y)
  126.         return atan (y / x);
  127.       else
  128.         return pi_over_2 - atan (x / y);
  129.     }
  130.       else
  131.     {
  132.       if (x > -y)
  133.         return atan (y / x);
  134.       else
  135.         return - pi_over_2 - atan (x / y);
  136.     }
  137.     }
  138.   else
  139.     {
  140.       if (y > 0)
  141.     {
  142.       if (-x > y)
  143.         return pi + atan (y / x);
  144.       else
  145.         return pi_over_2 - atan (x / y);
  146.     }
  147.       else
  148.     {
  149.       if (-x > -y)
  150.         return - pi + atan (y / x);
  151.       else if (y < 0)
  152.         return - pi_over_2 - atan (x / y);
  153.       else
  154.         {
  155.           double value;
  156.  
  157.           errno = EDOM;
  158.           __asm ("fmove%.d %#0x7fffffffffffffff,%0"     /* quiet NaN */
  159.              : "=f" (value)
  160.              : /* no inputs */);
  161.           return value;
  162.         }
  163.     }
  164.     }
  165. }
  166.  
  167. __inline static const double sinh (double x)
  168. {
  169.   double value;
  170.  
  171.   __asm ("fsinh%.x %1,%0"
  172.      : "=f" (value)
  173.      : "f" (x));
  174.   return value;
  175. }
  176.  
  177. __inline static const double cosh (double x)
  178. {
  179.   double value;
  180.  
  181.   __asm ("fcosh%.x %1,%0"
  182.      : "=f" (value)
  183.      : "f" (x));
  184.   return value;
  185. }
  186.  
  187. __inline static const double tanh (double x)
  188. {
  189.   double value;
  190.  
  191.   __asm ("ftanh%.x %1,%0"
  192.      : "=f" (value)
  193.      : "f" (x));
  194.   return value;
  195. }
  196.  
  197. __inline static const double atanh (double x)
  198. {
  199.   double value;
  200.  
  201.   __asm ("fatanh%.x %1,%0"
  202.      : "=f" (value)
  203.      : "f" (x));
  204.   return value;
  205. }
  206.  
  207. __inline static const double exp (double x)
  208. {
  209.   double value;
  210.  
  211.   __asm ("fetox%.x %1,%0"
  212.      : "=f" (value)
  213.      : "f" (x));
  214.   return value;
  215. }
  216.  
  217. __inline static const double expm1 (double x)
  218. {
  219.   double value;
  220.  
  221.   __asm ("fetoxm1%.x %1,%0"
  222.      : "=f" (value)
  223.      : "f" (x));
  224.   return value;
  225. }
  226.  
  227. __inline static const double log (double x)
  228. {
  229.   double value;
  230.  
  231.   __asm ("flogn%.x %1,%0"
  232.      : "=f" (value)
  233.      : "f" (x));
  234.   return value;
  235. }
  236.  
  237. __inline static const double log1p (double x)
  238. {
  239.   double value;
  240.  
  241.   __asm ("flognp1%.x %1,%0"
  242.      : "=f" (value)
  243.      : "f" (x));
  244.   return value;
  245. }
  246.  
  247. __inline static const double log10 (double x)
  248. {
  249.   double value;
  250.  
  251.   __asm ("flog10%.x %1,%0"
  252.      : "=f" (value)
  253.      : "f" (x));
  254.   return value;
  255. }
  256.  
  257. __inline static const double sqrt (double x)
  258. {
  259.   double value;
  260.  
  261.   __asm ("fsqrt%.x %1,%0"
  262.      : "=f" (value)
  263.      : "f" (x));
  264.   return value;
  265. }
  266.  
  267. __inline static const double hypot (const double x, const double y)
  268. {
  269.   return sqrt (x*x + y*y);
  270. }
  271.  
  272. __inline static const double pow (const double x, const double y)
  273. {
  274.   if (x > 0)
  275.     return exp (y * log (x));
  276.   else if (x == 0)
  277.     {
  278.       if (y > 0)
  279.     return 0.0;
  280.       else
  281.     {
  282.       double value;
  283.  
  284.       errno = EDOM;
  285.       __asm ("fmove%.d %#07fffffffffffffff,%0"        /* quiet NaN */
  286.          : "=f" (value)
  287.          : /* no inputs */);
  288.       return value;
  289.     }
  290.     }
  291.   else
  292.     {
  293.       double temp;
  294.  
  295.       __asm ("fintrz%.x %1,%0"
  296.          : "=f" (temp)            /* integer-valued float */
  297.          : "f" (y));
  298.       if (y == temp)
  299.         {
  300.       int i = (int) y;
  301.       
  302.       if (i & 1 == 0)            /* even */
  303.         return exp (y * log (-x));
  304.       else
  305.         return - exp (y * log (-x));
  306.         }
  307.       else
  308.         {
  309.       double value;
  310.  
  311.       errno = EDOM;
  312.       __asm ("fmove%.d %#0x7fffffffffffffff,%0"        /* quiet NaN */
  313.          : "=f" (value)
  314.          : /* no inputs */);
  315.       return value;
  316.         }
  317.     }
  318. }
  319.  
  320. __inline static const double fabs (double x)
  321. {
  322.   double value;
  323.  
  324.   __asm ("fabs%.x %1,%0"
  325.      : "=f" (value)
  326.      : "f" (x));
  327.   return value;
  328. }
  329.  
  330. __inline static const double ceil (double x)
  331. {
  332.   int rounding_mode, round_up;
  333.   double value;
  334.  
  335.   __asm volatile ("fmove%.l fpcr,%0"
  336.           : "=dm" (rounding_mode)
  337.           : /* no inputs */ );
  338.   round_up = rounding_mode | 0x30;
  339.   __asm volatile ("fmove%.l %0,fpcr"
  340.           : /* no outputs */
  341.           : "dmi" (round_up));
  342.   __asm volatile ("fint%.x %1,%0"
  343.           : "=f" (value)
  344.           : "f" (x));
  345.   __asm volatile ("fmove%.l %0,fpcr"
  346.           : /* no outputs */
  347.           : "dmi" (rounding_mode));
  348.   return value;
  349. }
  350.  
  351. __inline static const double floor (double x)
  352. {
  353.   int rounding_mode, round_down;
  354.   double value;
  355.  
  356.   __asm volatile ("fmove%.l fpcr,%0"
  357.           : "=dm" (rounding_mode)
  358.           : /* no inputs */ );
  359.   round_down = (rounding_mode & ~0x10)
  360.         | 0x20;
  361.   __asm volatile ("fmove%.l %0,fpcr"
  362.           : /* no outputs */
  363.           : "dmi" (round_down));
  364.   __asm volatile ("fint%.x %1,%0"
  365.           : "=f" (value)
  366.           : "f" (x));
  367.   __asm volatile ("fmove%.l %0,fpcr"
  368.           : /* no outputs */
  369.           : "dmi" (rounding_mode));
  370.   return value;
  371. }
  372.  
  373. __inline static const double rint (double x)
  374. {
  375.   int rounding_mode, round_nearest;
  376.   double value;
  377.  
  378.   __asm volatile ("fmove%.l fpcr,%0"
  379.           : "=dm" (rounding_mode)
  380.           : /* no inputs */ );
  381.   round_nearest = rounding_mode & ~0x30;
  382.   __asm volatile ("fmove%.l %0,fpcr"
  383.           : /* no outputs */
  384.           : "dmi" (round_nearest));
  385.   __asm volatile ("fint%.x %1,%0"
  386.           : "=f" (value)
  387.           : "f" (x));
  388.   __asm volatile ("fmove%.l %0,fpcr"
  389.           : /* no outputs */
  390.           : "dmi" (rounding_mode));
  391.   return value;
  392. }
  393.  
  394. __inline static const double fmod (double x, double y)
  395. {
  396.   double value;
  397.  
  398.   __asm ("fmod%.x %2,%0"
  399.      : "=f" (value)
  400.      : "0" (x),
  401.        "f" (y));
  402.   return value;
  403. }
  404.  
  405. __inline static const double drem (double x, double y)
  406. {
  407.   double value;
  408.  
  409.   __asm ("frem%.x %2,%0"
  410.      : "=f" (value)
  411.      : "0" (x),
  412.        "f" (y));
  413.   return value;
  414. }
  415.  
  416. __inline static const double scalb (double x, int n)
  417. {
  418.   double value;
  419.  
  420.   __asm ("fscale%.l %2,%0"
  421.      : "=f" (value)
  422.      : "0" (x),
  423.        "dmi" (n));
  424.   return value;
  425. }
  426.  
  427. __inline static double logb (double x)
  428. {
  429.   double exponent;
  430.  
  431.   __asm ("fgetexp%.x %1,%0"
  432.      : "=f" (exponent)
  433.      : "f" (x));
  434.   return exponent;
  435. }
  436.  
  437. __inline static const double ldexp (double x, int n)
  438. {
  439.   double value;
  440.  
  441.   __asm ("fscale%.l %2,%0"
  442.      : "=f" (value)
  443.      : "0" (x),
  444.        "dmi" (n));
  445.   return value;
  446. }
  447.  
  448. __inline static double frexp (double x, int *exp)
  449. {
  450.   double float_exponent;
  451.   int int_exponent;
  452.   double mantissa;
  453.  
  454.   __asm ("fgetexp%.x %1,%0"
  455.      : "=f" (float_exponent)     /* integer-valued float */
  456.      : "f" (x));
  457.   int_exponent = (int) float_exponent;
  458.   __asm ("fgetman%.x %1,%0"
  459.      : "=f" (mantissa)        /* 1.0 <= mantissa < 2.0 */
  460.      : "f" (x));
  461.   if (mantissa != 0)
  462.     {
  463.       __asm ("fscale%.b %#-1,%0"
  464.          : "=f" (mantissa)        /* mantissa /= 2.0 */
  465.          : "0" (mantissa));
  466.       int_exponent += 1;
  467.     }
  468.   *exp = int_exponent;
  469.   return mantissa;
  470. }
  471.  
  472. __inline static double modf (double x, double *ip)
  473. {
  474.   double temp;
  475.  
  476.   __asm ("fintrz%.x %1,%0"
  477.      : "=f" (temp)            /* integer-valued float */
  478.      : "f" (x));
  479.   *ip = temp;
  480.   return x - temp;
  481. }
  482.  
  483.