sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2014 Free Software Foundation, Inc.
   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: uroot.c                                              */
  21 /*                                                                   */
  22 /* FUNCTION:    usqrt                                                */
  23 /*                                                                   */
  24 /* FILES NEEDED: dla.h endian.h mydefs.h uroot.h                     */
  25 /*               uroot.tbl                                           */
  26 /*                                                                   */
  27 /* An ultimate sqrt routine. Given an IEEE double machine number x   */
  28 /* it computes the correctly rounded (to nearest) value of square    */
  29 /* root of x.                                                        */
  30 /* Assumption: Machine arithmetic operations are performed in        */
  31 /* round to nearest mode of IEEE 754 standard.                       */
  32 /*                                                                   */
  33 /*********************************************************************/
  34
  35 #include <math_private.h>
  36
  37 typedef union {int64_t i[2]; long double x; double d[2]; } mynumber;
  38
  39 static const double
  40   t512 = 0x1p512,
  41   tm256 = 0x1p-256,
  42   two54 = 0x1p54,       /* 0x4350000000000000 */
  43   twom54 = 0x1p-54;     /* 0x3C90000000000000 */
  44
  45 /*********************************************************************/
  46 /* An ultimate sqrt routine. Given an IEEE double machine number x   */
  47 /* it computes the correctly rounded (to nearest) value of square    */
  48 /* root of x.                                                        */
  49 /*********************************************************************/
  50 long double __ieee754_sqrtl(long double x)
  51 {
  52   static const long double big = 134217728.0, big1 = 134217729.0;
  53   long double t,s,i;
  54   mynumber a,c;
  55   uint64_t k, l;
  56   int64_t m, n;
  57   double d;
  58
  59   a.x=x;
  60   k=a.i[0] & INT64_C(0x7fffffffffffffff);
  61   /*----------------- 2^-1022  <= | x |< 2^1024  -----------------*/
  62   if (k>INT64_C(0x000fffff00000000) && k<INT64_C(0x7ff0000000000000)) {
  63     if (x < 0) return (big1-big1)/(big-big);
  64     l = (k&INT64_C(0x001fffffffffffff))|INT64_C(0x3fe0000000000000);
  65     if ((a.i[1] & INT64_C(0x7fffffffffffffff)) != 0) {
  66       n = (int64_t) ((l - k) * 2) >> 53;
  67       m = (a.i[1] >> 52) & 0x7ff;
  68       if (m == 0) {
  69         a.d[1] *= two54;
  70         m = ((a.i[1] >> 52) & 0x7ff) - 54;
  71       }
  72       m += n;
  73       if (m > 0)
  74         a.i[1] = (a.i[1] & INT64_C(0x800fffffffffffff)) | (m << 52);
  75       else if (m <= -54) {
  76         a.i[1] &= INT64_C(0x8000000000000000);
  77       } else {
  78         m += 54;
  79         a.i[1] = (a.i[1] & INT64_C(0x800fffffffffffff)) | (m << 52);
  80         a.d[1] *= twom54;
  81       }
  82     }
  83     a.i[0] = l;
  84     s = a.x;
  85     d = __ieee754_sqrt (a.d[0]);
  86     c.i[0] = INT64_C(0x2000000000000000)+((k&INT64_C(0x7fe0000000000000))>>1);
  87     c.i[1] = 0;
  88     i = d;
  89     t = 0.5L * (i + s / i);
  90     i = 0.5L * (t + s / t);
  91     return c.x * i;
  92   }
  93   else {
  94     if (k>=INT64_C(0x7ff0000000000000)) {
  95       if (a.i[0] == INT64_C(0xfff0000000000000))
  96         return (big1-big1)/(big-big); /* sqrt (-Inf) = NaN.  */
  97       return x; /* sqrt (NaN) = NaN, sqrt (+Inf) = +Inf.  */
  98     }
  99     if (x == 0) return x;
 100     if (x < 0) return (big1-big1)/(big-big);
 101     return tm256*__ieee754_sqrtl(x*t512);
 102   }
 103 }
 104 strong_alias (__ieee754_sqrtl, __sqrtl_finite)