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