# Category : C Source Code

Archive : CEPHES22.ZIP

Filename : 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);

}

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

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

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/