Back to index

glibc  2.9
mpsqrt.c
Go to the documentation of this file.
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 02111-1307, 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 /* Multi-Precision square root function subroutine for precision p >= 4.    */
00030 /* The relative error is bounded by 3.501*r**(1-p), where r=2**24.          */
00031 /*                                                                          */
00032 /****************************************************************************/
00033 #include "endian.h"
00034 #include "mpa.h"
00035 
00036 /****************************************************************************/
00037 /* Multi-Precision square root function subroutine for precision p >= 4.    */
00038 /* The relative error is bounded by 3.501*r**(1-p), 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 multi-precision 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 }