sysdeps/ieee754/dbl-64/mpsqrt.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001, 2011 Free Software Foundation
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU Lesser General Public License as published by
   8  * the Free Software Foundation; either version 2.1 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public License
  17  * along with this program; if not, see <http://www.gnu.org/licenses/>.
  18  */
  19 /****************************************************************************/
  20 /*  MODULE_NAME:mpsqrt.c                                                    */
  21 /*                                                                          */
  22 /*  FUNCTION:mpsqrt                                                         */
  23 /*           fastiroot                                                      */
  24 /*                                                                          */
  25 /* FILES NEEDED:endian.h mpa.h mpsqrt.h                                     */
  26 /*              mpa.c                                                       */
  27 /* Multi-Precision square root function subroutine for precision p >= 4.    */
  28 /* The relative error is bounded by 3.501*r**(1-p), where r=2**24.          */
  29 /*                                                                          */
  30 /****************************************************************************/
  31 #include "endian.h"
  32 #include "mpa.h"
  33
  34 #ifndef SECTION
  35 # define SECTION
  36 #endif
  37
  38 #include "mpsqrt.h"
  39
  40 /****************************************************************************/
  41 /* Multi-Precision square root function subroutine for precision p >= 4.    */
  42 /* The relative error is bounded by 3.501*r**(1-p), where r=2**24.          */
  43 /* Routine receives two pointers to  Multi Precision numbers:               */
  44 /* x (left argument) and y (next argument). Routine also receives precision */
  45 /* p as integer. Routine computes sqrt(*x) and stores result in *y          */
  46 /****************************************************************************/
  47
  48 static double fastiroot(double);
  49
  50 void
  51 SECTION
  52 __mpsqrt(mp_no *x, mp_no *y, int p) {
  53   int i,m,ey;
  54   double dx,dy;
  55   mp_no
  56     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,
  57                    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,
  58                    0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}},
  59     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,
  60                    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,
  61                    0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
  62   mp_no mpxn,mpz,mpu,mpt1,mpt2;
  63
  64   /* Prepare multi-precision 1/2 and 3/2 */
  65   mphalf.e  =0;  mphalf.d[0]  =ONE;  mphalf.d[1]  =HALFRAD;
  66   mp3halfs.e=1;  mp3halfs.d[0]=ONE;  mp3halfs.d[1]=ONE;  mp3halfs.d[2]=HALFRAD;
  67
  68   ey=EX/2;     __cpy(x,&mpxn,p);    mpxn.e -= (ey+ey);
  69   __mp_dbl(&mpxn,&dx,p);   dy=fastiroot(dx);    __dbl_mp(dy,&mpu,p);
  70   __mul(&mpxn,&mphalf,&mpz,p);
  71
  72   m=__mpsqrt_mp[p];
  73   for (i=0; i<m; i++) {
  74     __mul(&mpu,&mpu,&mpt1,p);
  75     __mul(&mpt1,&mpz,&mpt2,p);
  76     __sub(&mp3halfs,&mpt2,&mpt1,p);
  77     __mul(&mpu,&mpt1,&mpt2,p);
  78     __cpy(&mpt2,&mpu,p);
  79   }
  80   __mul(&mpxn,&mpu,y,p);  EY += ey;
  81
  82   return;
  83 }
  84
  85 /***********************************************************/
  86 /* Compute a double precision approximation for 1/sqrt(x)  */
  87 /* with the relative error bounded by 2**-51.              */
  88 /***********************************************************/
  89 static double
  90 SECTION
  91 fastiroot(double x) {
  92   union {int i[2]; double d;} p,q;
  93   double y,z, t;
  94   int n;
  95   static const double c0 = 0.99674, c1 = -0.53380, c2 = 0.45472, c3 = -0.21553;
  96
  97   p.d = x;
  98   p.i[HIGH_HALF] = (p.i[HIGH_HALF] & 0x3FFFFFFF ) | 0x3FE00000 ;
  99   q.d = x;
 100   y = p.d;
 101   z = y -1.0;
 102   n = (q.i[HIGH_HALF] - p.i[HIGH_HALF])>>1;
 103   z = ((c3*z + c2)*z + c1)*z + c0;            /* 2**-7         */
 104   z = z*(1.5 - 0.5*y*z*z);                    /* 2**-14        */
 105   p.d = z*(1.5 - 0.5*y*z*z);                  /* 2**-28        */
 106   p.i[HIGH_HALF] -= n;
 107   t = x*p.d;
 108   return p.d*(1.5 - 0.5*p.d*t);
 109 }