Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : CBRT.C

 
Output of file : CBRT.C contained in archive : CEPHES22.ZIP
/* cbrt.c
*
* Cube root
*
*
*
* SYNOPSIS:
*
* double x, y, cbrt();
*
* y = cbrt( x );
*
*
*
* DESCRIPTION:
*
* Returns the cube root of the argument, which may be negative.
*
* Range reduction involves determining the power of 2 of
* the argument. A polynomial of degree 2 applied to the
* mantissa, and multiplication by the cube root of 1, 2, or 4
* approximates the root to within about 0.1%. Then Newton's
* iteration is used three times to converge to an accurate
* result.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC 0,8 50000 1.8e-17 5.6e-18
* IEEE 0,1e308 30000 1.5e-16 5.0e-17
*
*/
/* cbrt.c */

/*
Cephes Math Library Release 2.1: December, 1988
Copyright 1984, 1987, 1988 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/


#include "mconf.h"

#ifdef UNK
static double CBRT2 = 1.25992104989487316477;
static double CBRT4 = 1.58740105196819947475;
#endif

#ifdef DEC
static short CBT2[] = {040241,042427,0146153,0112127};
static short CBT4[] = {040313,027765,024753,070744};
#define CBRT2 *(double *)CBT2
#define CBRT4 *(double *)CBT4
#endif

#ifdef IBMPC
static short CBT2[] = {0x728b,0xf98d,0x28a2,0x3ff4};
static short CBT4[] = {0x6e3d,0xa53d,0x65fe,0x3ff9};
#define CBRT2 *(double *)CBT2
#define CBRT4 *(double *)CBT4
#endif

#ifdef MIEEE
static short CBT2[] = {
0x3ff4,0x28a2,0xf98d,0x728b
};
static short CBT4[] = {
0x3ff9,0x65fe,0xa53d,0x6e3d,
};
#define CBRT2 *(double *)CBT2
#define CBRT4 *(double *)CBT4
#endif

double cbrt(x)
double x;
{
int e, rem, sign;
double z;
double frexp(), ldexp();

if( x == 0 )
return( 0.0 );
if( x > 0 )
sign = 1;
else
{
sign = -1;
x = -x;
}

z = x;
/* extract power of 2, leaving
* mantissa between 0.5 and 1
*/
x = frexp( x, &e );

/* Approximate cube root of number between .5 and 1,
* peak relative error = 6.36e-4
*/
x = (-1.9150215751434832257e-1 * x
+ 6.9757045195246484393e-1) * x
+ 4.9329566506409572486e-1;


/* exponent divided by 3 */
if( e >= 0 )
{
rem = e;
e /= 3;
rem -= 3*e;
if( rem == 1 )
x *= CBRT2;
else if( rem == 2 )
x *= CBRT4;
}


/* argument less than 1 */

else
{
e = -e;
rem = e;
e /= 3;
rem -= 3*e;
if( rem == 1 )
x /= CBRT2;
else if( rem == 2 )
x /= CBRT4;
e = -e;
}

/* multiply by power of 2 */
x = ldexp( x, e );

/* Newton iteration */
x -= ( x - (z/(x*x)) )/3.0;
x -= ( x - (z/(x*x)) )/3.0;
x -= ( x - (z/(x*x)) )/3.0;

if( sign < 0 )
x = -x;
return(x);
}


  3 Responses to “Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : CBRT.C

  1. Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

  2. This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

  3. But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/