amath  1.6.2
Simple command line calculator
Functions | Variables
cbrt.c File Reference
#include "prim.h"
Include dependency graph for cbrt.c:

Go to the source code of this file.

Functions

double cbrt (double x)
 Cube root function. More...
 

Variables

static const unsigned B1 = 715094163
 
static const unsigned B2 = 696219795
 
static const double C = 5.42857142857142815906e-01
 
static const double D = -7.05306122448979611050e-01
 
static const double E = 1.41428571428571436819e+00
 
static const double F = 1.60714285714285720630e+00
 
static const double G = 3.57142857142857150787e-01
 

Function Documentation

double cbrt ( double  x)

Cube root function.

Version
1.3
Date
95/01/18
Copyright
Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. Developed at SunSoft, a Sun Microsystems, Inc. business. Permission to use, copy, modify, and distribute this software is freely granted, provided that this notice is preserved.

Definition at line 67 of file cbrt.c.

References B1, B2, C, D, E, F, and G.

68 {
69  sword hx, lx, ht;
70  double r,s,t=0.0,w;
71  uword sign;
72 
73  GET_HIGH_WORD(hx,x); /* high word of x */
74  sign=hx&0x80000000; /* sign= sign(x) */
75  hx ^=sign;
76  if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */
77  GET_LOW_WORD(lx, x);
78  if((hx|lx)==0)
79  return(x); /* cbrt(0) is itself */
80 
81  SET_HIGH_WORD(x,hx); /* x <- |x| */
82  /* rough cbrt to 5 bits */
83  if(hx<0x00100000) /* subnormal number */
84  {
85  SET_HIGH_WORD(t,0x43500000); /* set t= 2**54 */
86  t*=x;
87  GET_HIGH_WORD(ht,t);
88  SET_HIGH_WORD(t,ht/3+B2);
89  }
90  else
91  SET_HIGH_WORD(t,hx/3+B1);
92 
93  /* new cbrt to 23 bits, may be implemented in single precision */
94  r=t*t/x;
95  s=C+r*t;
96  t*=G+F/(s+E+D/s);
97 
98  /* chopped to 20 bits and make it larger than cbrt(x) */
99  SET_LOW_WORD(t,0);
100  GET_HIGH_WORD(ht,t);
101  SET_HIGH_WORD(t,ht + 0x00000001);
102 
103  /* one step newton iteration to 53 bits with error less than 0.667 ulps */
104  s=t*t; /* t*t is exact */
105  r=x/s;
106  w=t+t;
107  r=(r-t)/(w+r); /* r-s is exact */
108  t=t+t*r;
109 
110  /* retore the sign bit */
111  GET_HIGH_WORD(ht,t);
112  SET_HIGH_WORD(t,ht|sign);
113 
114  return(t);
115 }
D
static const double D
Definition: cbrt.c:52
E
static const double E
Definition: cbrt.c:53
GET_HIGH_WORD
#define GET_HIGH_WORD(i, d)
Get the more significant 32 bit int from a double.
Definition: prim.h:165
GET_LOW_WORD
#define GET_LOW_WORD(i, d)
Get the less significant 32 bit int from a double.
Definition: prim.h:176
F
static const double F
Definition: cbrt.c:54
sword
signed int sword
32 bit signed integer.
Definition: prim.h:107
uword
unsigned int uword
32 bit unsigned integer.
Definition: prim.h:101
B2
static const unsigned B2
Definition: cbrt.c:48
G
static const double G
Definition: cbrt.c:55
B1
static const unsigned B1
Definition: cbrt.c:47
C
static const double C
Definition: cbrt.c:51
SET_LOW_WORD
#define SET_LOW_WORD(d, v)
Set the less significant 32 bits of a double from an int.
Definition: prim.h:211
SET_HIGH_WORD
#define SET_HIGH_WORD(d, v)
Set the more significant 32 bits of a double from an int.
Definition: prim.h:199

Variable Documentation

const unsigned B1 = 715094163
static

Definition at line 47 of file cbrt.c.

Referenced by cbrt().

const unsigned B2 = 696219795
static

Definition at line 48 of file cbrt.c.

Referenced by cbrt().

const double C = 5.42857142857142815906e-01
static

Definition at line 51 of file cbrt.c.

Referenced by cbrt().

const double D = -7.05306122448979611050e-01
static

Definition at line 52 of file cbrt.c.

Referenced by cbrt().

const double E = 1.41428571428571436819e+00
static

Definition at line 53 of file cbrt.c.

Referenced by cbrt().

const double F = 1.60714285714285720630e+00
static

Definition at line 54 of file cbrt.c.

Referenced by cbrt().

const double G = 3.57142857142857150787e-01
static

Definition at line 55 of file cbrt.c.

Referenced by cbrt().