URL https://opencores.org/ocsvn/or1k/or1k/trunk

Subversion Repositories or1k

[/] [or1k/] [trunk/] [newlib/] [newlib/] [libm/] [math/] [s_cbrt.c] - Blame information for rev 40

Go to most recent revision | Details | Compare with Previous | View Log


 
/* @(#)s_cbrt.c 5.1 93/09/24 */
/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 *
 */
 
/*
FUNCTION
        <<cbrt>>, <<cbrtf>>---cube root
 
INDEX
        cbrt
INDEX
        cbrtf
 
ANSI_SYNOPSIS
        #include <math.h>
        double cbrt(double <[x]>);
        float  cbrtf(float <[x]>);
 
TRAD_SYNOPSIS
        #include <math.h>
        double cbrt(<[x]>);
        float  cbrtf(<[x]>);
 
DESCRIPTION
        <<cbrt>> computes the cube root of the argument.
 
RETURNS
        The cube root is returned.
 
PORTABILITY
        <<cbrt>> is in System V release 4.  <<cbrtf>> is an extension.
*/
 
#include "fdlibm.h"
 
#ifndef _DOUBLE_IS_32BITS
 
/* cbrt(x)
 * Return cube root of x
 */
#ifdef __STDC__
static const __uint32_t
#else
static __uint32_t
#endif
        B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
        B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
 
#ifdef __STDC__
static const double
#else
static double
#endif
C =  5.42857142857142815906e-01, /* 19/35     = 0x3FE15F15, 0xF15F15F1 */
D = -7.05306122448979611050e-01, /* -864/1225 = 0xBFE691DE, 0x2532C834 */
E =  1.41428571428571436819e+00, /* 99/70     = 0x3FF6A0EA, 0x0EA0EA0F */
F =  1.60714285714285720630e+00, /* 45/28     = 0x3FF9B6DB, 0x6DB6DB6E */
G =  3.57142857142857150787e-01; /* 5/14      = 0x3FD6DB6D, 0xB6DB6DB7 */
 
#ifdef __STDC__
        double cbrt(double x)
#else
        double cbrt(x)
        double x;
#endif
{
        __int32_t       hx;
        double r,s,t=0.0,w;
        __uint32_t sign;
        __uint32_t high,low;
 
        GET_HIGH_WORD(hx,x);
        sign=hx&0x80000000;             /* sign= sign(x) */
        hx  ^=sign;
        if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */
        GET_LOW_WORD(low,x);
        if((hx|low)==0)
            return(x);          /* cbrt(0) is itself */
 
        SET_HIGH_WORD(x,hx);    /* x <- |x| */
    /* rough cbrt to 5 bits */
        if(hx<0x00100000)               /* subnormal number */
          {SET_HIGH_WORD(t,0x43500000); /* set t= 2**54 */
           t*=x; GET_HIGH_WORD(high,t); SET_HIGH_WORD(t,high/3+B2);
          }
        else
          SET_HIGH_WORD(t,hx/3+B1);
 
 
    /* new cbrt to 23 bits, may be implemented in single precision */
        r=t*t/x;
        s=C+r*t;
        t*=G+F/(s+E+D/s);
 
    /* chopped to 20 bits and make it larger than cbrt(x) */
        GET_HIGH_WORD(high,t);
        INSERT_WORDS(t,high+0x00000001,0);
 
 
    /* one step newton iteration to 53 bits with error less than 0.667 ulps */
        s=t*t;          /* t*t is exact */
        r=x/s;
        w=t+t;
        r=(r-t)/(w+r);  /* r-s is exact */
        t=t+t*r;
 
    /* retore the sign bit */
        GET_HIGH_WORD(high,t);
        SET_HIGH_WORD(t,high|sign);
        return(t);
}
 
#endif /* _DOUBLE_IS_32BITS */

Line No.	Rev	Author	Line
1	39	lampret
2			`/* @(#)s_cbrt.c 5.1 93/09/24 */`
3			`/*`
4			`* ====================================================`
5			`* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.`
6			`*`
7			`* Developed at SunPro, a Sun Microsystems, Inc. business.`
8			`* Permission to use, copy, modify, and distribute this`
9			`* software is freely granted, provided that this notice`
10			`* is preserved.`
11			`* ====================================================`
12			`*`
13			`*/`
14
15			`/*`
16			`FUNCTION`
17			`<<cbrt>>, <<cbrtf>>---cube root`
18
19			`INDEX`
20			`cbrt`
21			`INDEX`
22			`cbrtf`
23
24			`ANSI_SYNOPSIS`
25			`#include <math.h>`
26			`double cbrt(double <[x]>);`
27			`float cbrtf(float <[x]>);`
28
29			`TRAD_SYNOPSIS`
30			`#include <math.h>`
31			`double cbrt(<[x]>);`
32			`float cbrtf(<[x]>);`
33
34			`DESCRIPTION`
35			`<<cbrt>> computes the cube root of the argument.`
36
37			`RETURNS`
38			`The cube root is returned.`
39
40			`PORTABILITY`
41			`<<cbrt>> is in System V release 4. <<cbrtf>> is an extension.`
42			`*/`
43
44			`#include "fdlibm.h"`
45
46			`#ifndef _DOUBLE_IS_32BITS`
47
48			`/* cbrt(x)`
49			`* Return cube root of x`
50			`*/`
51			`#ifdef __STDC__`
52			`static const __uint32_t`
53			`#else`
54			`static __uint32_t`
55			`#endif`
56			`B1 = 715094163, /* B1 = (682-0.03306235651)220 /`
57			`B2 = 696219795; /* B2 = (664-0.03306235651)220 /`
58
59			`#ifdef __STDC__`
60			`static const double`
61			`#else`
62			`static double`
63			`#endif`
64			`C = 5.42857142857142815906e-01, /* 19/35 = 0x3FE15F15, 0xF15F15F1 */`
65			`D = -7.05306122448979611050e-01, /* -864/1225 = 0xBFE691DE, 0x2532C834 */`
66			`E = 1.41428571428571436819e+00, /* 99/70 = 0x3FF6A0EA, 0x0EA0EA0F */`
67			`F = 1.60714285714285720630e+00, /* 45/28 = 0x3FF9B6DB, 0x6DB6DB6E */`
68			`G = 3.57142857142857150787e-01; /* 5/14 = 0x3FD6DB6D, 0xB6DB6DB7 */`
69
70			`#ifdef __STDC__`
71			`double cbrt(double x)`
72			`#else`
73			`double cbrt(x)`
74			`double x;`
75			`#endif`
76			`{`
77			`__int32_t hx;`
78			`double r,s,t=0.0,w;`
79			`__uint32_t sign;`
80			`__uint32_t high,low;`
81
82			`GET_HIGH_WORD(hx,x);`
83			`sign=hx&0x80000000; /* sign= sign(x) */`
84			`hx ^=sign;`
85			`if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */`
86			`GET_LOW_WORD(low,x);`
87			`if((hx\|low)==0)`
88			`return(x); /* cbrt(0) is itself */`
89
90			`SET_HIGH_WORD(x,hx); /* x <- \|x\| */`
91			`/* rough cbrt to 5 bits */`
92			`if(hx<0x00100000) /* subnormal number */`
93			`{SET_HIGH_WORD(t,0x43500000); /* set t= 2*54 /`
94			`t*=x; GET_HIGH_WORD(high,t); SET_HIGH_WORD(t,high/3+B2);`
95			`}`
96			`else`
97			`SET_HIGH_WORD(t,hx/3+B1);`
98
99
100			`/* new cbrt to 23 bits, may be implemented in single precision */`
101			`r=t*t/x;`
102			`s=C+r*t;`
103			`t*=G+F/(s+E+D/s);`
104
105			`/* chopped to 20 bits and make it larger than cbrt(x) */`
106			`GET_HIGH_WORD(high,t);`
107			`INSERT_WORDS(t,high+0x00000001,0);`
108
109
110			`/* one step newton iteration to 53 bits with error less than 0.667 ulps */`
111			`s=tt; / tt is exact /`
112			`r=x/s;`
113			`w=t+t;`
114			`r=(r-t)/(w+r); /* r-s is exact */`
115			`t=t+t*r;`
116
117			`/* retore the sign bit */`
118			`GET_HIGH_WORD(high,t);`
119			`SET_HIGH_WORD(t,high\|sign);`
120			`return(t);`
121			`}`
122
123			`#endif /* _DOUBLE_IS_32BITS */`

Browse

Tools

Subversion Repositories or1k

[/] [or1k/] [trunk/] [newlib/] [newlib/] [libm/] [math/] [s_cbrt.c] - Blame information for rev 40