home *** CD-ROM | disk | FTP | other *** search
/ Il CD di internet / CD.iso / SOURCE / D / LIBC / LIBC-4.6 / LIBC-4 / libc-linux / sysdeps / linux / i386 / math / pow.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-01-09  |  1.9 KB  |  78 lines
  1. /* Copyright (C) 1993  Hongjiu Lu
  2. This file is part of the Linux C Library.
  3.  
  4. The Linux C Library is free software; you can redistribute it and/or
  5. modify it under the terms of the GNU Library General Public License as
  6. published by the Free Software Foundation; either version 2 of the
  7. License, or (at your option) any later version.
  8.  
  9. The Linux C Library is distributed in the hope that it will be useful,
  10. but WITHOUT ANY WARRANTY; without even the implied warranty of
  11. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12. Library General Public License for more details. */
  13.  
  14. #include <errno.h>
  15. #include <math.h>
  16. #include <fpu_control.h>
  17.  
  18.     
  19. double pow (x,y)      
  20. double x,y;
  21. {
  22.   int negative;
  23.   __volatile__ unsigned short cw, saved_cw;    /* 80387 control word */
  24.  
  25.   if (y == 0.0) return 1.0;
  26.  
  27.   if (x == 0.0) return (y < 0.0) ? __infnan(EDOM) : 0.0;
  28.  
  29.   if (y == 1.0) return x;
  30.  
  31.   if (x < 0.0) {
  32.     long long tmp;
  33.  
  34.     tmp = (long long) y; 
  35.  
  36.     /* Even or odd */
  37.     negative = tmp & 1;
  38.  
  39.     /* That is only valid if |y| < 2^64. */
  40.     if (y != (double) tmp) {
  41.        return __infnan( EDOM);
  42.     }
  43.  
  44.     x = -x;
  45.     } else {
  46.     negative = 0;
  47.     }
  48.  
  49.  /*
  50.   * Inline assembler implements the following algorithm:
  51.   *
  52.   * 1 - x' = y *log2(x)
  53.   * 2 - find x1, x2 where x' = x1 + x2 and |x2| <= 0.5
  54.   * 3 - x = 2^x2 scaled by x1
  55.   */
  56.  
  57.   __asm__ __volatile__ ("fnstcw %0" : "=m" (cw) : );
  58.   saved_cw = cw;
  59.  
  60.   cw &= 0xf3ff;    /* force round to nearest */
  61.   cw |= 0x003f; /* turn off exceptions (see ANSI/ISO 9989 section 7.5.1) */
  62.  
  63.   __asm__ __volatile__ ("fldcw %0" : : "m" (cw));
  64.  
  65.   __asm__ __volatile__ ("fyl2x;fstl %2;frndint;fstl %%st(2);fsubrp;f2xm1;"
  66.             "fld1;faddp;fscale"
  67.                   : "=t" (x) : "0" (x), "u" (y));
  68.  
  69.   /* Restore the control word */
  70.   __asm__ __volatile__ ("fldcw %0" : : "m" (saved_cw));
  71.  
  72.   /* Check for results that can not be represented (kludgy): */
  73.   if (__isinf(x) || __isnan(x))
  74.     errno = ERANGE;
  75.  
  76.   return (negative) ? -x : x;
  77. }
  78.