Back to index

glibc  2.9
halfulp.c
Go to the documentation of this file.
00001 /*
00002  * IBM Accurate Mathematical Library
00003  * written by International Business Machines Corp.
00004  * Copyright (C) 2001, 2005 Free Software Foundation
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU Lesser General Public License as published by
00008  * the Free Software Foundation; either version 2.1 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU Lesser General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00019  */
00020 /************************************************************************/
00021 /*                                                                      */
00022 /* MODULE_NAME:halfulp.c                                                */
00023 /*                                                                      */
00024 /*  FUNCTIONS:halfulp                                                   */
00025 /*  FILES NEEDED: mydefs.h dla.h endian.h                               */
00026 /*                uroot.c                                               */
00027 /*                                                                      */
00028 /*Routine halfulp(double x, double y) computes x^y where result does    */
00029 /*not need rounding. If the result is closer to 0 than can be           */
00030 /*represented it returns 0.                                             */
00031 /*     In the following cases the function does not compute anything    */
00032 /*and returns a negative number:                                        */
00033 /*1. if the result needs rounding,                                      */
00034 /*2. if y is outside the interval [0,  2^20-1],                         */
00035 /*3. if x can be represented by  x=2**n for some integer n.             */
00036 /************************************************************************/
00037 
00038 #include "endian.h"
00039 #include "mydefs.h"
00040 #include "dla.h"
00041 #include "math_private.h"
00042 
00043 double __ieee754_sqrt(double x);
00044 
00045 static const int4 tab54[32] = {
00046    262143, 11585, 1782, 511, 210, 107, 63, 42,
00047        30,    22,   17,  14,  12,  10,  9,  7,
00048         7,     6,    5,   5,   5,   4,  4,  4,
00049         3,     3,    3,   3,   3,   3,  3,  3 };
00050 
00051 
00052 double __halfulp(double x, double y)
00053 {
00054   mynumber v;
00055   double z,u,uu,j1,j2,j3,j4,j5;
00056   int4 k,l,m,n;
00057   if (y <= 0) {               /*if power is negative or zero */
00058     v.x = y;
00059     if (v.i[LOW_HALF] != 0) return -10.0;
00060     v.x = x;
00061     if (v.i[LOW_HALF] != 0) return -10.0;
00062     if ((v.i[HIGH_HALF]&0x000fffff) != 0) return -10;   /* if x =2 ^ n */
00063     k = ((v.i[HIGH_HALF]&0x7fffffff)>>20)-1023;         /* find this n */
00064     z = (double) k;
00065     return (z*y == -1075.0)?0: -10.0;
00066   }
00067                               /* if y > 0  */
00068   v.x = y;
00069     if (v.i[LOW_HALF] != 0) return -10.0;
00070 
00071   v.x=x;
00072                               /*  case where x = 2**n for some integer n */
00073   if (((v.i[HIGH_HALF]&0x000fffff)|v.i[LOW_HALF]) == 0) {
00074     k=(v.i[HIGH_HALF]>>20)-1023;
00075     return (((double) k)*y == -1075.0)?0:-10.0;
00076   }
00077 
00078   v.x = y;
00079   k = v.i[HIGH_HALF];
00080   m = k<<12;
00081   l = 0;
00082   while (m)
00083     {m = m<<1; l++; }
00084   n = (k&0x000fffff)|0x00100000;
00085   n = n>>(20-l);                       /*   n is the odd integer of y    */
00086   k = ((k>>20) -1023)-l;               /*   y = n*2**k                   */
00087   if (k>5) return -10.0;
00088   if (k>0) for (;k>0;k--) n *= 2;
00089   if (n > 34) return -10.0;
00090   k = -k;
00091   if (k>5) return -10.0;
00092 
00093                             /*   now treat x        */
00094   while (k>0) {
00095     z = __ieee754_sqrt(x);
00096     EMULV(z,z,u,uu,j1,j2,j3,j4,j5);
00097     if (((u-x)+uu) != 0) break;
00098     x = z;
00099     k--;
00100  }
00101   if (k) return -10.0;
00102 
00103   /* it is impossible that n == 2,  so the mantissa of x must be short  */
00104 
00105   v.x = x;
00106   if (v.i[LOW_HALF]) return -10.0;
00107   k = v.i[HIGH_HALF];
00108   m = k<<12;
00109   l = 0;
00110   while (m) {m = m<<1; l++; }
00111   m = (k&0x000fffff)|0x00100000;
00112   m = m>>(20-l);                       /*   m is the odd integer of x    */
00113 
00114             /*   now check whether the length of m**n is at most 54 bits */
00115 
00116   if  (m > tab54[n-3]) return -10.0;
00117 
00118              /* yes, it is - now compute x**n by simple multiplications  */
00119 
00120   u = x;
00121   for (k=1;k<n;k++) u = u*x;
00122   return u;
00123 }