glibc
2.9

00001 00002 /* 00003 * IBM Accurate Mathematical Library 00004 * written by International Business Machines Corp. 00005 * Copyright (C) 2001 Free Software Foundation 00006 * 00007 * This program is free software; you can redistribute it and/or modify 00008 * it under the terms of the GNU Lesser General Public License as published by 00009 * the Free Software Foundation; either version 2.1 of the License, or 00010 * (at your option) any later version. 00011 * 00012 * This program is distributed in the hope that it will be useful, 00013 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 * GNU Lesser General Public License for more details. 00016 * 00017 * You should have received a copy of the GNU Lesser General Public License 00018 * along with this program; if not, write to the Free Software 00019 * Foundation, Inc., 59 Temple Place  Suite 330, Boston, MA 021111307, USA. 00020 */ 00021 /****************************************************************************/ 00022 /* MODULE_NAME:mpsqrt.c */ 00023 /* */ 00024 /* FUNCTION:mpsqrt */ 00025 /* fastiroot */ 00026 /* */ 00027 /* FILES NEEDED:endian.h mpa.h mpsqrt.h */ 00028 /* mpa.c */ 00029 /* MultiPrecision square root function subroutine for precision p >= 4. */ 00030 /* The relative error is bounded by 3.501*r**(1p), where r=2**24. */ 00031 /* */ 00032 /****************************************************************************/ 00033 #include "endian.h" 00034 #include "mpa.h" 00035 00036 /****************************************************************************/ 00037 /* MultiPrecision square root function subroutine for precision p >= 4. */ 00038 /* The relative error is bounded by 3.501*r**(1p), where r=2**24. */ 00039 /* Routine receives two pointers to Multi Precision numbers: */ 00040 /* x (left argument) and y (next argument). Routine also receives precision */ 00041 /* p as integer. Routine computes sqrt(*x) and stores result in *y */ 00042 /****************************************************************************/ 00043 00044 double fastiroot(double); 00045 00046 void __mpsqrt(mp_no *x, mp_no *y, int p) { 00047 #include "mpsqrt.h" 00048 00049 int i,m,ex,ey; 00050 double dx,dy; 00051 mp_no 00052 mphalf = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 00053 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 00054 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}, 00055 mp3halfs = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 00056 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 00057 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; 00058 mp_no mpxn,mpz,mpu,mpt1,mpt2; 00059 00060 /* Prepare multiprecision 1/2 and 3/2 */ 00061 mphalf.e =0; mphalf.d[0] =ONE; mphalf.d[1] =HALFRAD; 00062 mp3halfs.e=1; mp3halfs.d[0]=ONE; mp3halfs.d[1]=ONE; mp3halfs.d[2]=HALFRAD; 00063 00064 ex=EX; ey=EX/2; __cpy(x,&mpxn,p); mpxn.e = (ey+ey); 00065 __mp_dbl(&mpxn,&dx,p); dy=fastiroot(dx); __dbl_mp(dy,&mpu,p); 00066 __mul(&mpxn,&mphalf,&mpz,p); 00067 00068 m=mp[p]; 00069 for (i=0; i<m; i++) { 00070 __mul(&mpu,&mpu,&mpt1,p); 00071 __mul(&mpt1,&mpz,&mpt2,p); 00072 __sub(&mp3halfs,&mpt2,&mpt1,p); 00073 __mul(&mpu,&mpt1,&mpt2,p); 00074 __cpy(&mpt2,&mpu,p); 00075 } 00076 __mul(&mpxn,&mpu,y,p); EY += ey; 00077 00078 return; 00079 } 00080 00081 /***********************************************************/ 00082 /* Compute a double precision approximation for 1/sqrt(x) */ 00083 /* with the relative error bounded by 2**51. */ 00084 /***********************************************************/ 00085 double fastiroot(double x) { 00086 union {int i[2]; double d;} p,q; 00087 double y,z, t; 00088 int n; 00089 static const double c0 = 0.99674, c1 = 0.53380, c2 = 0.45472, c3 = 0.21553; 00090 00091 p.d = x; 00092 p.i[HIGH_HALF] = (p.i[HIGH_HALF] & 0x3FFFFFFF )  0x3FE00000 ; 00093 q.d = x; 00094 y = p.d; 00095 z = y 1.0; 00096 n = (q.i[HIGH_HALF]  p.i[HIGH_HALF])>>1; 00097 z = ((c3*z + c2)*z + c1)*z + c0; /* 2**7 */ 00098 z = z*(1.5  0.5*y*z*z); /* 2**14 */ 00099 p.d = z*(1.5  0.5*y*z*z); /* 2**28 */ 00100 p.i[HIGH_HALF] = n; 00101 t = x*p.d; 00102 return p.d*(1.5  0.5*p.d*t); 00103 }