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