← Back to branch summary

~xnox/ubuntu/trusty/gcc-arm-linux-androideabi/dima

« back to all changes in this revision

Viewing changes to android/bionic/libm/bsdsrc/b_exp.c

  • Committer: Package Import Robot
  • Author(s): Dmitrijs Ledkovs
  • Date: 2013-07-05 10:12:24 UTC
  • Revision ID: package-import@ubuntu.com-20130705101224-6qo3e8jbz8p31aa1
Tags: upstream-0.20130705.1
Import upstream version 0.20130705.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
android/bionic/libm/bsdsrc/b_exp.c
 
1
/*
 
2
 * Copyright (c) 1985, 1993
 
3
 *      The Regents of the University of California.  All rights reserved.
 
4
 *
 
5
 * Redistribution and use in source and binary forms, with or without
 
6
 * modification, are permitted provided that the following conditions
 
7
 * are met:
 
8
 * 1. Redistributions of source code must retain the above copyright
 
9
 *    notice, this list of conditions and the following disclaimer.
 
10
 * 2. Redistributions in binary form must reproduce the above copyright
 
11
 *    notice, this list of conditions and the following disclaimer in the
 
12
 *    documentation and/or other materials provided with the distribution.
 
13
 * 3. All advertising materials mentioning features or use of this software
 
14
 *    must display the following acknowledgement:
 
15
 *      This product includes software developed by the University of
 
16
 *      California, Berkeley and its contributors.
 
17
 * 4. Neither the name of the University nor the names of its contributors
 
18
 *    may be used to endorse or promote products derived from this software
 
19
 *    without specific prior written permission.
 
20
 *
 
21
 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
 
22
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 
23
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 
24
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
 
25
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 
26
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 
27
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 
28
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 
29
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 
30
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 
31
 * SUCH DAMAGE.
 
32
 */
 
33
 
 
34
#ifndef lint
 
35
static char sccsid[] = "@(#)exp.c       8.1 (Berkeley) 6/4/93";
 
36
#endif /* not lint */
 
37
#include <sys/cdefs.h>
 
38
/* __FBSDID("$FreeBSD: src/lib/msun/bsdsrc/b_exp.c,v 1.7 2004/12/16 20:40:37 das Exp $"); */
 
39
 
 
40
 
 
41
/* EXP(X)
 
42
 * RETURN THE EXPONENTIAL OF X
 
43
 * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
 
44
 * CODED IN C BY K.C. NG, 1/19/85;
 
45
 * REVISED BY K.C. NG on 2/6/85, 2/15/85, 3/7/85, 3/24/85, 4/16/85, 6/14/86.
 
46
 *
 
47
 * Required system supported functions:
 
48
 *      scalb(x,n)
 
49
 *      copysign(x,y)
 
50
 *      finite(x)
 
51
 *
 
52
 * Method:
 
53
 *      1. Argument Reduction: given the input x, find r and integer k such
 
54
 *         that
 
55
 *                         x = k*ln2 + r,  |r| <= 0.5*ln2 .
 
56
 *         r will be represented as r := z+c for better accuracy.
 
57
 *
 
58
 *      2. Compute exp(r) by
 
59
 *
 
60
 *              exp(r) = 1 + r + r*R1/(2-R1),
 
61
 *         where
 
62
 *              R1 = x - x^2*(p1+x^2*(p2+x^2*(p3+x^2*(p4+p5*x^2)))).
 
63
 *
 
64
 *      3. exp(x) = 2^k * exp(r) .
 
65
 *
 
66
 * Special cases:
 
67
 *      exp(INF) is INF, exp(NaN) is NaN;
 
68
 *      exp(-INF)=  0;
 
69
 *      for finite argument, only exp(0)=1 is exact.
 
70
 *
 
71
 * Accuracy:
 
72
 *      exp(x) returns the exponential of x nearly rounded. In a test run
 
73
 *      with 1,156,000 random arguments on a VAX, the maximum observed
 
74
 *      error was 0.869 ulps (units in the last place).
 
75
 */
 
76
 
 
77
#include "mathimpl.h"
 
78
 
 
79
static const double p1 = 0x1.555555555553ep-3;
 
80
static const double p2 = -0x1.6c16c16bebd93p-9;
 
81
static const double p3 = 0x1.1566aaf25de2cp-14;
 
82
static const double p4 = -0x1.bbd41c5d26bf1p-20;
 
83
static const double p5 = 0x1.6376972bea4d0p-25;
 
84
static const double ln2hi = 0x1.62e42fee00000p-1;
 
85
static const double ln2lo = 0x1.a39ef35793c76p-33;
 
86
static const double lnhuge = 0x1.6602b15b7ecf2p9;
 
87
static const double lntiny = -0x1.77af8ebeae354p9;
 
88
static const double invln2 = 0x1.71547652b82fep0;
 
89
 
 
90
#if 0
 
91
double exp(x)
 
92
double x;
 
93
{
 
94
        double  z,hi,lo,c;
 
95
        int k;
 
96
 
 
97
#if !defined(vax)&&!defined(tahoe)
 
98
        if(x!=x) return(x);     /* x is NaN */
 
99
#endif  /* !defined(vax)&&!defined(tahoe) */
 
100
        if( x <= lnhuge ) {
 
101
                if( x >= lntiny ) {
 
102
 
 
103
                    /* argument reduction : x --> x - k*ln2 */
 
104
 
 
105
                        k=invln2*x+copysign(0.5,x);     /* k=NINT(x/ln2) */
 
106
 
 
107
                    /* express x-k*ln2 as hi-lo and let x=hi-lo rounded */
 
108
 
 
109
                        hi=x-k*ln2hi;
 
110
                        x=hi-(lo=k*ln2lo);
 
111
 
 
112
                    /* return 2^k*[1+x+x*c/(2+c)]  */
 
113
                        z=x*x;
 
114
                        c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5))));
 
115
                        return  scalb(1.0+(hi-(lo-(x*c)/(2.0-c))),k);
 
116
 
 
117
                }
 
118
                /* end of x > lntiny */
 
119
 
 
120
                else
 
121
                     /* exp(-big#) underflows to zero */
 
122
                     if(finite(x))  return(scalb(1.0,-5000));
 
123
 
 
124
                     /* exp(-INF) is zero */
 
125
                     else return(0.0);
 
126
        }
 
127
        /* end of x < lnhuge */
 
128
 
 
129
        else
 
130
        /* exp(INF) is INF, exp(+big#) overflows to INF */
 
131
            return( finite(x) ?  scalb(1.0,5000)  : x);
 
132
}
 
133
#endif
 
134
 
 
135
/* returns exp(r = x + c) for |c| < |x| with no overlap.  */
 
136
 
 
137
double __exp__D(x, c)
 
138
double x, c;
 
139
{
 
140
        double  z,hi,lo;
 
141
        int k;
 
142
 
 
143
        if (x != x)     /* x is NaN */
 
144
                return(x);
 
145
        if ( x <= lnhuge ) {
 
146
                if ( x >= lntiny ) {
 
147
 
 
148
                    /* argument reduction : x --> x - k*ln2 */
 
149
                        z = invln2*x;
 
150
                        k = z + copysign(.5, x);
 
151
 
 
152
                    /* express (x+c)-k*ln2 as hi-lo and let x=hi-lo rounded */
 
153
 
 
154
                        hi=(x-k*ln2hi);                 /* Exact. */
 
155
                        x= hi - (lo = k*ln2lo-c);
 
156
                    /* return 2^k*[1+x+x*c/(2+c)]  */
 
157
                        z=x*x;
 
158
                        c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5))));
 
159
                        c = (x*c)/(2.0-c);
 
160
 
 
161
                        return  scalb(1.+(hi-(lo - c)), k);
 
162
                }
 
163
                /* end of x > lntiny */
 
164
 
 
165
                else
 
166
                     /* exp(-big#) underflows to zero */
 
167
                     if(finite(x))  return(scalb(1.0,-5000));
 
168
 
 
169
                     /* exp(-INF) is zero */
 
170
                     else return(0.0);
 
171
        }
 
172
        /* end of x < lnhuge */
 
173
 
 
174
        else
 
175
        /* exp(INF) is INF, exp(+big#) overflows to INF */
 
176
            return( finite(x) ?  scalb(1.0,5000)  : x);
 
177
}