newlib/libm/common/sqrtl.c

   1 /*
   2 (C) Copyright IBM Corp. 2009
   3
   4 All rights reserved.
   5
   6 Redistribution and use in source and binary forms, with or without
   7 modification, are permitted provided that the following conditions are met:
   8
   9 * Redistributions of source code must retain the above copyright notice,
  10 this list of conditions and the following disclaimer.
  11 * Redistributions in binary form must reproduce the above copyright
  12 notice, this list of conditions and the following disclaimer in the
  13 documentation and/or other materials provided with the distribution.
  14 * Neither the name of IBM nor the names of its contributors may be
  15 used to endorse or promote products derived from this software without
  16 specific prior written permission.
  17
  18 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  19 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  20 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  21 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  22 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  23 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  24 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  25 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  26 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  27 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  28 POSSIBILITY OF SUCH DAMAGE.
  29 */
  30
  31 #include <math.h>
  32 #include "local.h"
  33
  34 #ifdef _LDBL_EQ_DBL
  35 /* On platforms where long double is as wide as double.  */
  36 long double
  37 sqrtl (long double x)
  38 {
  39   return sqrt(x);
  40 }
  41
  42 #else
  43
  44   /* This code is based upon the version in the BSD math's library.
  45      That code is...
  46    *
  47    * Copyright (c) 2007 Steven G. Kargl
  48    * All rights reserved.
  49    *
  50    * Redistribution and use in source and binary forms, with or without
  51    * modification, are permitted provided that the following conditions
  52    * are met:
  53    * 1. Redistributions of source code must retain the above copyright
  54    *    notice unmodified, this list of conditions, and the following
  55    *    disclaimer.
  56    * 2. Redistributions in binary form must reproduce the above copyright
  57    *    notice, this list of conditions and the following disclaimer in the
  58    *    documentation and/or other materials provided with the distribution.
  59    *
  60    * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
  61    * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
  62    * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
  63    * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
  64    * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
  65    * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  66    * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  67    * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  68    * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  69    * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  70    */
  71 #include <float.h>
  72 #include "ieeefp.h"
  73
  74 #ifndef LDBL_NBIT
  75 #define LDBL_NBIT       0
  76 #endif
  77
  78 #ifndef LDBL_MAX_EXP
  79 #define LDBL_MAX_EXP    DBL_MAX_EXP
  80 #endif
  81
  82 /* Return (x + ulp) for normal positive x.  Assumes no overflow.  */
  83
  84 static inline long double
  85 inc (long double x)
  86 {
  87   union ieee_ext_u ux = { .extu_ld = x, };
  88
  89   if (++ux.extu_ext.ext_fracl == 0)
  90     {
  91       if (++ux.extu_ext.ext_frach == 0)
  92         {
  93           ux.extu_ext.ext_exp++;
  94           ux.extu_ext.ext_frach |= LDBL_NBIT;
  95         }
  96     }
  97
  98   return ux.extu_ld;
  99 }
 100
 101 /* Return (x - ulp) for normal positive x.  Assumes no underflow.  */
 102
 103 static inline long double
 104 dec (long double x)
 105 {
 106   union ieee_ext_u ux = { .extu_ld = x, };
 107
 108   if (ux.extu_ext.ext_fracl-- == 0)
 109     {
 110       if (ux.extu_ext.ext_frach-- == LDBL_NBIT)
 111         {
 112           ux.extu_ext.ext_exp--;
 113           ux.extu_ext.ext_frach |= LDBL_NBIT;
 114         }
 115     }
 116
 117   return ux.extu_ld;
 118 }
 119
 120 /* This is slow, but simple and portable.  */
 121
 122 long double
 123 sqrtl (long double x)
 124 {
 125   union ieee_ext_u ux = { .extu_ld = x, };
 126   int k;
 127   long double lo, xn;
 128
 129   /* If x = NaN, then sqrt(x) = NaN.  */
 130   /* If x = Inf, then sqrt(x) = Inf.  */
 131   /* If x = -Inf, then sqrt(x) = NaN.  */
 132   if (ux.extu_ext.ext_exp == LDBL_MAX_EXP * 2 - 1)
 133     return (x * x + x);
 134
 135   /* If x = +-0, then sqrt(x) = +-0.  */
 136   if (x == 0.0L || x == -0.0L)
 137     return x;
 138
 139   /* If x < 0, then raise invalid and return NaN.  */
 140   if (ux.extu_ext.ext_sign)
 141     return ((x - x) / (x - x));
 142
 143   if (ux.extu_ext.ext_exp == 0)
 144     {
 145       /* Adjust subnormal numbers.  */
 146       ux.extu_ld *= 0x1.0p514;
 147       k = -514;
 148     }
 149   else
 150     k = 0;
 151
 152   /* ux.extu_ld is a normal number, so break it into ux.extu_ld = e*2^n where
 153      ux.extu_ld = (2*e)*2^2k for odd n and ux.extu_ld = (4*e)*2^2k for even n.  */
 154
 155   if ((ux.extu_ext.ext_exp - EXT_EXP_BIAS) & 1)
 156     {
 157       /* n is even.  */
 158       k += ux.extu_ext.ext_exp - EXT_EXP_BIAS - 1; /* 2k = n - 2.  */
 159       ux.extu_ext.ext_exp = EXT_EXP_BIAS + 1;      /* ux.extu_ld in [2,4).  */
 160     }
 161   else
 162     {
 163       k += ux.extu_ext.ext_exp - EXT_EXP_BIAS;  /* 2k = n - 1.  */
 164       ux.extu_ext.ext_exp = EXT_EXP_BIAS;       /* ux.extu_ld in [1,2).  */
 165     }
 166
 167   /* Newton's iteration.
 168      Split ux.extu_ld into a high and low part to achieve additional precision.  */
 169
 170   xn = sqrt ((double) ux.extu_ld);      /* 53-bit estimate of sqrtl(x).  */
 171
 172 #if LDBL_MANT_DIG > 100
 173   xn = (xn + (ux.extu_ld / xn)) * 0.5;  /* 106-bit estimate.  */
 174 #endif
 175
 176   lo = ux.extu_ld;
 177   ux.extu_ext.ext_fracl = 0;            /* Zero out lower bits.  */
 178   lo = (lo - ux.extu_ld) / xn;          /* Low bits divided by xn.  */
 179   xn = xn + (ux.extu_ld / xn);          /* High portion of estimate.  */
 180   ux.extu_ld = xn + lo;                 /* Combine everything.  */
 181   ux.extu_ext.ext_exp += (k >> 1) - 1;
 182
 183   xn = x / ux.extu_ld;                  /* Chopped quotient (inexact?).  */
 184
 185   /* For simplicity we round to nearest.  */
 186   xn = inc (xn);                        /* xn = xn + ulp.  */
 187
 188   ux.extu_ld = ux.extu_ld + xn;         /* Chopped sum.  */
 189   ux.extu_ext.ext_exp--;
 190
 191   return ux.extu_ld;
 192 }
 193 #endif /* ! _LDBL_EQ_DBL */
 194
 195