home *** CD-ROM | disk | FTP | other *** search
/ OpenStep (Enterprise) / OpenStepENTCD.toast / OEDEV / GNUSRC.Z / math-3300.h < prev    next >
C/C++ Source or Header  |  1993-09-28  |  8KB  |  462 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. #include <errno.h>
  23.  
  24. #undef HUGE_VAL
  25. #define HUGE_VAL                            \
  26. ({                                    \
  27.   double huge_val;                            \
  28.                                     \
  29.   __asm ("fmove%.d %#0x7ff0000000000000,%0"    /* Infinity */        \
  30.      : "=f" (huge_val)                        \
  31.      : /* no inputs */);                        \
  32.   huge_val;                                \
  33. })
  34.  
  35. __inline static const double sin (double x)
  36. {
  37.   double value;
  38.  
  39.   __asm ("fsin%.x %1,%0"
  40.      : "=f" (value)
  41.      : "f" (x));
  42.   return value;
  43. }
  44.  
  45. __inline static const double cos (double x)
  46. {
  47.   double value;
  48.  
  49.   __asm ("fcos%.x %1,%0"
  50.      : "=f" (value)
  51.      : "f" (x));
  52.   return value;
  53. }
  54.  
  55. __inline static const double tan (double x)
  56. {
  57.   double value;
  58.  
  59.   __asm ("ftan%.x %1,%0"
  60.      : "=f" (value)
  61.      : "f" (x));
  62.   return value;
  63. }
  64.  
  65. __inline static const double asin (double x)
  66. {
  67.   double value;
  68.  
  69.   __asm ("fasin%.x %1,%0"
  70.      : "=f" (value)
  71.      : "f" (x));
  72.   return value;
  73. }
  74.  
  75. __inline static const double acos (double x)
  76. {
  77.   double value;
  78.  
  79.   __asm ("facos%.x %1,%0"
  80.      : "=f" (value)
  81.      : "f" (x));
  82.   return value;
  83. }
  84.  
  85. __inline static const double atan (double x)
  86. {
  87.   double value;
  88.  
  89.   __asm ("fatan%.x %1,%0"
  90.      : "=f" (value)
  91.      : "f" (x));
  92.   return value;
  93. }
  94.  
  95. __inline static const double atan2 (double y, double x)
  96. {
  97.   double pi, pi_over_2;
  98.  
  99.   __asm ("fmovecr%.x %#0,%0"        /* extended precision pi */
  100.      : "=f" (pi)
  101.      : /* no inputs */ );
  102.   __asm ("fscale%.b %#-1,%0"        /* no loss of accuracy */
  103.      : "=f" (pi_over_2)
  104.      : "0" (pi));
  105.   if (x > 0)
  106.     {
  107.       if (y > 0)
  108.     {
  109.       if (x > y)
  110.         return atan (y / x);
  111.       else
  112.         return pi_over_2 - atan (x / y);
  113.     }
  114.       else
  115.     {
  116.       if (x > -y)
  117.         return atan (y / x);
  118.       else
  119.         return - pi_over_2 - atan (x / y);
  120.     }
  121.     }
  122.   else
  123.     {
  124.       if (y > 0)
  125.     {
  126.       if (-x > y)
  127.         return pi + atan (y / x);
  128.       else
  129.         return pi_over_2 - atan (x / y);
  130.     }
  131.       else
  132.     {
  133.       if (-x > -y)
  134.         return - pi + atan (y / x);
  135.       else if (y < 0)
  136.         return - pi_over_2 - atan (x / y);
  137.       else
  138.         {
  139.           double value;
  140.  
  141.           errno = EDOM;
  142.           __asm ("fmove%.d %#0x7fffffffffffffff,%0"     /* quiet NaN */
  143.              : "=f" (value)
  144.              : /* no inputs */);
  145.           return value;
  146.         }
  147.     }
  148.     }
  149. }
  150.  
  151. __inline static const double sinh (double x)
  152. {
  153.   double value;
  154.  
  155.   __asm ("fsinh%.x %1,%0"
  156.      : "=f" (value)
  157.      : "f" (x));
  158.   return value;
  159. }
  160.  
  161. __inline static const double cosh (double x)
  162. {
  163.   double value;
  164.  
  165.   __asm ("fcosh%.x %1,%0"
  166.      : "=f" (value)
  167.      : "f" (x));
  168.   return value;
  169. }
  170.  
  171. __inline static const double tanh (double x)
  172. {
  173.   double value;
  174.  
  175.   __asm ("ftanh%.x %1,%0"
  176.      : "=f" (value)
  177.      : "f" (x));
  178.   return value;
  179. }
  180.  
  181. __inline static const double atanh (double x)
  182. {
  183.   double value;
  184.  
  185.   __asm ("fatanh%.x %1,%0"
  186.      : "=f" (value)
  187.      : "f" (x));
  188.   return value;
  189. }
  190.  
  191. __inline static const double exp (double x)
  192. {
  193.   double value;
  194.  
  195.   __asm ("fetox%.x %1,%0"
  196.      : "=f" (value)
  197.      : "f" (x));
  198.   return value;
  199. }
  200.  
  201. __inline static const double expm1 (double x)
  202. {
  203.   double value;
  204.  
  205.   __asm ("fetoxm1%.x %1,%0"
  206.      : "=f" (value)
  207.      : "f" (x));
  208.   return value;
  209. }
  210.  
  211. __inline static const double log (double x)
  212. {
  213.   double value;
  214.  
  215.   __asm ("flogn%.x %1,%0"
  216.      : "=f" (value)
  217.      : "f" (x));
  218.   return value;
  219. }
  220.  
  221. __inline static const double log1p (double x)
  222. {
  223.   double value;
  224.  
  225.   __asm ("flognp1%.x %1,%0"
  226.      : "=f" (value)
  227.      : "f" (x));
  228.   return value;
  229. }
  230.  
  231. __inline static const double log10 (double x)
  232. {
  233.   double value;
  234.  
  235.   __asm ("flog10%.x %1,%0"
  236.      : "=f" (value)
  237.      : "f" (x));
  238.   return value;
  239. }
  240.  
  241. __inline static const double sqrt (double x)
  242. {
  243.   double value;
  244.  
  245.   __asm ("fsqrt%.x %1,%0"
  246.      : "=f" (value)
  247.      : "f" (x));
  248.   return value;
  249. }
  250.  
  251. __inline static const double pow (const double x, const double y)
  252. {
  253.   if (x > 0)
  254.     return exp (y * log (x));
  255.   else if (x == 0)
  256.     {
  257.       if (y > 0)
  258.     return 0.0;
  259.       else
  260.     {
  261.       double value;
  262.  
  263.       errno = EDOM;
  264.       __asm ("fmove%.d %#0x7fffffffffffffff,%0"        /* quiet NaN */
  265.          : "=f" (value)
  266.          : /* no inputs */);
  267.       return value;
  268.     }
  269.     }
  270.   else
  271.     {
  272.       double temp;
  273.  
  274.       __asm ("fintrz%.x %1,%0"
  275.          : "=f" (temp)            /* integer-valued float */
  276.          : "f" (y));
  277.       if (y == temp)
  278.         {
  279.       int i = (int) y;
  280.       
  281.       if (i & 1 == 0)            /* even */
  282.         return exp (y * log (x));
  283.       else
  284.         return - exp (y * log (x));
  285.         }
  286.       else
  287.         {
  288.       double value;
  289.  
  290.       errno = EDOM;
  291.       __asm ("fmove%.d %#0x7fffffffffffffff,%0"        /* quiet NaN */
  292.          : "=f" (value)
  293.          : /* no inputs */);
  294.       return value;
  295.         }
  296.     }
  297. }
  298.  
  299. __inline static const double fabs (double x)
  300. {
  301.   double value;
  302.  
  303.   __asm ("fabs%.x %1,%0"
  304.      : "=f" (value)
  305.      : "f" (x));
  306.   return value;
  307. }
  308.  
  309. __inline static const double ceil (double x)
  310. {
  311.   int rounding_mode, round_up;
  312.   double value;
  313.  
  314.   __asm volatile ("fmove%.l %%fpcr,%0"
  315.           : "=dm" (rounding_mode)
  316.           : /* no inputs */ );
  317.   round_up = rounding_mode | 0x30;
  318.   __asm volatile ("fmove%.l %0,%%fpcr"
  319.           : /* no outputs */
  320.           : "dmi" (round_up));
  321.   __asm volatile ("fint%.x %1,%0"
  322.           : "=f" (value)
  323.           : "f" (x));
  324.   __asm volatile ("fmove%.l %0,%%fpcr"
  325.           : /* no outputs */
  326.           : "dmi" (rounding_mode));
  327.   return value;
  328. }
  329.  
  330. __inline static const double floor (double x)
  331. {
  332.   int rounding_mode, round_down;
  333.   double value;
  334.  
  335.   __asm volatile ("fmove%.l %%fpcr,%0"
  336.           : "=dm" (rounding_mode)
  337.           : /* no inputs */ );
  338.   round_down = (rounding_mode & ~0x10)
  339.         | 0x20;
  340.   __asm volatile ("fmove%.l %0,%%fpcr"
  341.           : /* no outputs */
  342.           : "dmi" (round_down));
  343.   __asm volatile ("fint%.x %1,%0"
  344.           : "=f" (value)
  345.           : "f" (x));
  346.   __asm volatile ("fmove%.l %0,%%fpcr"
  347.           : /* no outputs */
  348.           : "dmi" (rounding_mode));
  349.   return value;
  350. }
  351.  
  352. __inline static const double rint (double x)
  353. {
  354.   int rounding_mode, round_nearest;
  355.   double value;
  356.  
  357.   __asm volatile ("fmove%.l %%fpcr,%0"
  358.           : "=dm" (rounding_mode)
  359.           : /* no inputs */ );
  360.   round_nearest = rounding_mode & ~0x30;
  361.   __asm volatile ("fmove%.l %0,%%fpcr"
  362.           : /* no outputs */
  363.           : "dmi" (round_nearest));
  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 fmod (double x, double y)
  374. {
  375.   double value;
  376.  
  377.   __asm ("fmod%.x %2,%0"
  378.      : "=f" (value)
  379.      : "0" (x),
  380.        "f" (y));
  381.   return value;
  382. }
  383.  
  384. __inline static const double drem (double x, double y)
  385. {
  386.   double value;
  387.  
  388.   __asm ("frem%.x %2,%0"
  389.      : "=f" (value)
  390.      : "0" (x),
  391.        "f" (y));
  392.   return value;
  393. }
  394.  
  395. __inline static const double scalb (double x, int n)
  396. {
  397.   double value;
  398.  
  399.   __asm ("fscale%.l %2,%0"
  400.      : "=f" (value)
  401.      : "0" (x),
  402.        "dmi" (n));
  403.   return value;
  404. }
  405.  
  406. __inline static double logb (double x)
  407. {
  408.   double exponent;
  409.  
  410.   __asm ("fgetexp%.x %1,%0"
  411.      : "=f" (exponent)
  412.      : "f" (x));
  413.   return exponent;
  414. }
  415.  
  416. __inline static const double ldexp (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 frexp (double x, int *exp)
  428. {
  429.   double float_exponent;
  430.   int int_exponent;
  431.   double mantissa;
  432.  
  433.   __asm ("fgetexp%.x %1,%0"
  434.      : "=f" (float_exponent)     /* integer-valued float */
  435.      : "f" (x));
  436.   int_exponent = (int) float_exponent;
  437.   __asm ("fgetman%.x %1,%0"
  438.      : "=f" (mantissa)        /* 1.0 <= mantissa < 2.0 */
  439.      : "f" (x));
  440.   if (mantissa != 0)
  441.     {
  442.       __asm ("fscale%.b %#-1,%0"
  443.          : "=f" (mantissa)        /* mantissa /= 2.0 */
  444.          : "0" (mantissa));
  445.       int_exponent += 1;
  446.     }
  447.   *exp = int_exponent;
  448.   return mantissa;
  449. }
  450.  
  451. __inline static double modf (double x, double *ip)
  452. {
  453.   double temp;
  454.  
  455.   __asm ("fintrz%.x %1,%0"
  456.      : "=f" (temp)            /* integer-valued float */
  457.      : "f" (x));
  458.   *ip = temp;
  459.   return x - temp;
  460. }
  461.  
  462.