lib/cbrtl.c

   1 /* Compute cubic root of long double value.
   2    Copyright (C) 2012-2020 Free Software Foundation, Inc.
   3    Cephes Math Library Release 2.2: January, 1991
   4    Copyright 1984, 1991 by Stephen L. Moshier
   5    Adapted for glibc October, 2001.
   6
   7    This program is free software: you can redistribute it and/or modify
   8    it under the terms of the GNU General Public License as published by
   9    the Free Software Foundation; either version 3 of the License, or
  10    (at your option) any later version.
  11
  12    This program is distributed in the hope that it will be useful,
  13    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15    GNU General Public License for more details.
  16
  17    You should have received a copy of the GNU General Public License
  18    along with this program.  If not, see <https://www.gnu.org/licenses/>.  */
  19
  20 #include <config.h>
  21
  22 /* Specification.  */
  23 #include <math.h>
  24
  25 #if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
  26
  27 long double
  28 cbrtl (long double x)
  29 {
  30   return cbrt (x);
  31 }
  32
  33 #else
  34
  35 /* Code based on glibc/sysdeps/ieee754/ldbl-128/s_cbrtl.c.  */
  36
  37 /*                                                      cbrtl.c
  38  *
  39  *      Cube root, long double precision
  40  *
  41  *
  42  *
  43  * SYNOPSIS:
  44  *
  45  * long double x, y, cbrtl();
  46  *
  47  * y = cbrtl( x );
  48  *
  49  *
  50  *
  51  * DESCRIPTION:
  52  *
  53  * Returns the cube root of the argument, which may be negative.
  54  *
  55  * Range reduction involves determining the power of 2 of
  56  * the argument.  A polynomial of degree 2 applied to the
  57  * mantissa, and multiplication by the cube root of 1, 2, or 4
  58  * approximates the root to within about 0.1%.  Then Newton's
  59  * iteration is used three times to converge to an accurate
  60  * result.
  61  *
  62  *
  63  *
  64  * ACCURACY:
  65  *
  66  *                      Relative error:
  67  * arithmetic   domain     # trials      peak         rms
  68  *    IEEE       -8,8       100000      1.3e-34     3.9e-35
  69  *    IEEE    exp(+-707)    100000      1.3e-34     4.3e-35
  70  *
  71  */
  72
  73 static const long double CBRT2 = 1.259921049894873164767210607278228350570251L;
  74 static const long double CBRT4 = 1.587401051968199474751705639272308260391493L;
  75 static const long double CBRT2I = 0.7937005259840997373758528196361541301957467L;
  76 static const long double CBRT4I = 0.6299605249474365823836053036391141752851257L;
  77
  78 long double
  79 cbrtl (long double x)
  80 {
  81   if (isfinite (x) && x != 0.0L)
  82     {
  83       int e, rem, sign;
  84       long double z;
  85
  86       if (x > 0)
  87         sign = 1;
  88       else
  89         {
  90           sign = -1;
  91           x = -x;
  92         }
  93
  94       z = x;
  95       /* extract power of 2, leaving mantissa between 0.5 and 1  */
  96       x = frexpl (x, &e);
  97
  98       /* Approximate cube root of number between .5 and 1,
  99          peak relative error = 1.2e-6  */
 100       x = ((((1.3584464340920900529734e-1L * x
 101               - 6.3986917220457538402318e-1L) * x
 102              + 1.2875551670318751538055e0L) * x
 103             - 1.4897083391357284957891e0L) * x
 104            + 1.3304961236013647092521e0L) * x + 3.7568280825958912391243e-1L;
 105
 106       /* exponent divided by 3 */
 107       if (e >= 0)
 108         {
 109           rem = e;
 110           e /= 3;
 111           rem -= 3 * e;
 112           if (rem == 1)
 113             x *= CBRT2;
 114           else if (rem == 2)
 115             x *= CBRT4;
 116         }
 117       else
 118         {                           /* argument less than 1 */
 119           e = -e;
 120           rem = e;
 121           e /= 3;
 122           rem -= 3 * e;
 123           if (rem == 1)
 124             x *= CBRT2I;
 125           else if (rem == 2)
 126             x *= CBRT4I;
 127           e = -e;
 128         }
 129
 130       /* multiply by power of 2 */
 131       x = ldexpl (x, e);
 132
 133       /* Newton iteration */
 134       x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
 135       x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
 136       x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
 137
 138       if (sign < 0)
 139         x = -x;
 140       return x;
 141     }
 142   else
 143     {
 144 # ifdef __sgi /* so that when x == -0.0L, the result is -0.0L not +0.0L */
 145       return x;
 146 # else
 147       return x + x;
 148 # endif
 149     }
 150 }
 151
 152 #endif