gcc/ginclude/math-68881.h

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