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 / powl.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-26  |  2.0 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 "mathl.h"
  16. #include <fpu_control.h>
  17. #include <fp.h>
  18.     
  19. long double powl (long double x, long double y)      
  20. {
  21.   int negative;
  22.   __volatile__ unsigned short cw, saved_cw;    /* 80387 control word */
  23.  
  24.   if (y == 0.0L) return 1.0L;
  25.  
  26.   if (x == 0.0L) return (y < 0.0L) ? __infnanl(EDOM) : 0.0L;
  27.  
  28.   if (y == 1.0L) return x;
  29.  
  30.   if (x < 0.0L) {
  31.     long long tmp;
  32.  
  33.     tmp = (long long) y; 
  34.  
  35.     /* Even or odd */
  36.     negative = tmp & 1;
  37.  
  38.     /* That is only valid if |y| < 2^64. */
  39.     if (y != (double) tmp) {
  40.        return __infnanl( EDOM);
  41.     }
  42.  
  43.     x = -x;
  44.     } else {
  45.     negative = 0;
  46.     }
  47.  
  48.  /*
  49.   * Inline assembler implements the following algorithm:
  50.   *
  51.   * 1 - x' = y *log2(x)
  52.   * 2 - find x1, x2 where x' = x1 + x2 and |x2| <= 0.5
  53.   * 3 - x = 2^x2 scaled by x1
  54.   */
  55.  
  56.   __asm__ __volatile__ ("fnstcw %0" : "=m" (cw) : );
  57.   saved_cw = cw;
  58.  
  59.   cw &= 0xf3ff;    /* force round to nearest */
  60.   cw |= 0x003f; /* turn off exceptions (see ANSI/ISO 9989 section 7.5.1) */
  61.  
  62.   __asm__ __volatile__ ("fldcw %0" : : "m" (cw));
  63.  
  64.   __asm__ __volatile__ ("fyl2x;fstl %2;frndint;fstl %%st(2);"
  65.             "fsubrp;f2xm1;fld1;faddp;fscale"
  66.                   : "=t" (x) : "0" (x), "u" (y));
  67.  
  68.   /* Restore the control word */
  69.   __asm__ __volatile__ ("fldcw %0" : : "m" (saved_cw));
  70.  
  71.   /* Check for results that can not be represented (kludgy): */
  72.  
  73.   if (isinfl((long double)x) || isnanl((long double)x))
  74.     errno = ERANGE;
  75.  
  76.   return (negative) ? -x : x;
  77. }
  78.