OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-old/] [newlib-1.17.0/] [newlib/] [libm/] [common/] [s_cbrt.c] - Blame information for rev 853

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

Line No. Rev Author Line
1 148 jeremybenn
 
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)*2**20 */
57
        B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
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=t*t;          /* t*t 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 */

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.