/*
 * IBM Accurate Mathematical Library
 * Copyright (c) International Business Machines Corp., 2001
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or 
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  
 */
/****************************************************************************/
/*  MODULE_NAME:mpsqrt.c                                                    */
/*                                                                          */
/*  FUNCTION:mpsqrt                                                         */
/*           fastiroot                                                      */
/*                                                                          */
/* FILES NEEDED:endian.h mpa.h mpsqrt.h                                     */
/*              mpa.c                                                       */
/* Multi-Precision square root function subroutine for precision p >= 4.    */
/* The relative error is bounded by 3.501*r**(1-p), where r=2**24.          */
/*                                                                          */
/****************************************************************************/
#include "endian.h"
#include "mpa.h"

/****************************************************************************/
/* Multi-Precision square root function subroutine for precision p >= 4.    */
/* The relative error is bounded by 3.501*r**(1-p), where r=2**24.          */
/* Routine receives two pointers to  Multi Precision numbers:               */
/* x (left argument) and y (next argument). Routine also receives precision */
/* p as integer. Routine computes sqrt(*x) and stores result in *y          */
/****************************************************************************/

double fastiroot(double);

void mpsqrt(mp_no *x, mp_no *y, int p) {
#include "mpsqrt.h"

  int i,m,ex,ey;
  double dx,dy;
  mp_no
    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,
                   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.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,},
    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,
                   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.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,};
  mp_no mpxn,mpz,mpu,mpt1,mpt2;

  /* Prepare multi-precision 1/2 and 3/2 */
  mphalf.e  =0;  mphalf.d[0]  =ONE;  mphalf.d[1]  =HALFRAD;
  mp3halfs.e=1;  mp3halfs.d[0]=ONE;  mp3halfs.d[1]=ONE;  mp3halfs.d[2]=HALFRAD;

  ex=EX;      ey=EX/2;     cpy(x,&mpxn,p);    mpxn.e -= (ey+ey);
  mp_dbl(&mpxn,&dx,p);   dy=fastiroot(dx);    dbl_mp(dy,&mpu,p);
  mul(&mpxn,&mphalf,&mpz,p);

  m=mp[p];
  for (i=0; i<m; i++) {
    mul(&mpu,&mpu,&mpt1,p);
    mul(&mpt1,&mpz,&mpt2,p);
    sub(&mp3halfs,&mpt2,&mpt1,p);
    mul(&mpu,&mpt1,&mpt2,p);
    cpy(&mpt2,&mpu,p);
  }
  mul(&mpxn,&mpu,y,p);  EY += ey;

  return;
}

/***********************************************************/
/* Compute a double precision approximation for 1/sqrt(x)  */
/* with the relative error bounded by 2**-51.              */
/***********************************************************/
double fastiroot(double x) {
  union {long i[2]; double d;} p,q;
  double y,z, t;
  long n;
  static const double c0 = 0.99674, c1 = -0.53380, c2 = 0.45472, c3 = -0.21553;
  
  p.d = x;
  p.i[HIGH_HALF] = (p.i[HIGH_HALF] & 0x3FFFFFFF ) | 0x3FE00000 ;
  q.d = x;
  y = p.d;
  z = y -1.0;
  n = (q.i[HIGH_HALF] - p.i[HIGH_HALF])>>1;
  z = ((c3*z + c2)*z + c1)*z + c0;            /* 2**-7         */
  z = z*(1.5 - 0.5*y*z*z);                    /* 2**-14        */
  p.d = z*(1.5 - 0.5*y*z*z);                  /* 2**-28        */
  p.i[HIGH_HALF] -= n;
  t = x*p.d;
  return p.d*(1.5 - 0.5*p.d*t);
}