contrib/gcc/ginclude/math-3300.h

   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