Back to index

glibc  2.9
s_cbrtl.c
Go to the documentation of this file.
00001 /*                                               cbrtl.c
00002  *
00003  *     Cube root, long double precision
00004  *
00005  *
00006  *
00007  * SYNOPSIS:
00008  *
00009  * long double x, y, cbrtl();
00010  *
00011  * y = cbrtl( x );
00012  *
00013  *
00014  *
00015  * DESCRIPTION:
00016  *
00017  * Returns the cube root of the argument, which may be negative.
00018  *
00019  * Range reduction involves determining the power of 2 of
00020  * the argument.  A polynomial of degree 2 applied to the
00021  * mantissa, and multiplication by the cube root of 1, 2, or 4
00022  * approximates the root to within about 0.1%.  Then Newton's
00023  * iteration is used three times to converge to an accurate
00024  * result.
00025  *
00026  *
00027  *
00028  * ACCURACY:
00029  *
00030  *                      Relative error:
00031  * arithmetic   domain     # trials      peak         rms
00032  *    IEEE       -8,8       100000      1.3e-34     3.9e-35
00033  *    IEEE    exp(+-707)    100000      1.3e-34     4.3e-35
00034  *
00035  */
00036 
00037 /*
00038 Cephes Math Library Release 2.2: January, 1991
00039 Copyright 1984, 1991 by Stephen L. Moshier
00040 Adapted for glibc October, 2001.
00041 
00042     This library is free software; you can redistribute it and/or
00043     modify it under the terms of the GNU Lesser General Public
00044     License as published by the Free Software Foundation; either
00045     version 2.1 of the License, or (at your option) any later version.
00046 
00047     This library is distributed in the hope that it will be useful,
00048     but WITHOUT ANY WARRANTY; without even the implied warranty of
00049     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00050     Lesser General Public License for more details.
00051 
00052     You should have received a copy of the GNU Lesser General Public
00053     License along with this library; if not, write to the Free Software
00054     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00055 
00056 
00057 #include "math.h"
00058 #include "math_private.h"
00059 
00060 static const long double CBRT2 = 1.259921049894873164767210607278228350570251L;
00061 static const long double CBRT4 = 1.587401051968199474751705639272308260391493L;
00062 static const long double CBRT2I = 0.7937005259840997373758528196361541301957467L;
00063 static const long double CBRT4I = 0.6299605249474365823836053036391141752851257L;
00064 
00065 
00066 long double
00067 __cbrtl (long double x)
00068 {
00069   int e, rem, sign;
00070   long double z;
00071 
00072   if (!__finitel (x))
00073     return x + x;
00074 
00075   if (x == 0)
00076     return (x);
00077 
00078   if (x > 0)
00079     sign = 1;
00080   else
00081     {
00082       sign = -1;
00083       x = -x;
00084     }
00085 
00086   z = x;
00087  /* extract power of 2, leaving mantissa between 0.5 and 1  */
00088   x = __frexpl (x, &e);
00089 
00090   /* Approximate cube root of number between .5 and 1,
00091      peak relative error = 1.2e-6  */
00092   x = ((((1.3584464340920900529734e-1L * x
00093          - 6.3986917220457538402318e-1L) * x
00094         + 1.2875551670318751538055e0L) * x
00095        - 1.4897083391357284957891e0L) * x
00096        + 1.3304961236013647092521e0L) * x + 3.7568280825958912391243e-1L;
00097 
00098   /* exponent divided by 3 */
00099   if (e >= 0)
00100     {
00101       rem = e;
00102       e /= 3;
00103       rem -= 3 * e;
00104       if (rem == 1)
00105        x *= CBRT2;
00106       else if (rem == 2)
00107        x *= CBRT4;
00108     }
00109   else
00110     {                       /* argument less than 1 */
00111       e = -e;
00112       rem = e;
00113       e /= 3;
00114       rem -= 3 * e;
00115       if (rem == 1)
00116        x *= CBRT2I;
00117       else if (rem == 2)
00118        x *= CBRT4I;
00119       e = -e;
00120     }
00121 
00122   /* multiply by power of 2 */
00123   x = __ldexpl (x, e);
00124 
00125   /* Newton iteration */
00126   x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
00127   x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
00128   x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
00129 
00130   if (sign < 0)
00131     x = -x;
00132   return (x);
00133 }
00134 
00135 weak_alias (__cbrtl, cbrtl)