← Back to branch summary

~ubuntu-branches/ubuntu/natty/python3.1/natty-security

« back to all changes in this revision

Viewing changes to Python/pymath.c

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2010-07-06 16:52:42 UTC
  • mfrom: (1.2.1 upstream) (2.1.11 sid)
  • Revision ID: james.westby@ubuntu.com-20100706165242-2xv4i019r3et6c0j
Tags: 3.1.2+20100706-1ubuntu1
* Merge with Debian; remaining changes:
  - Regenerate the control file.
  - Add debian/patches/overwrite-semaphore-check for Lucid buildds.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Python/pymath.c
7
7
   thus rounding from extended precision to double precision. */
8
8
double _Py_force_double(double x)
9
9
{
10
 
        volatile double y;
11
 
        y = x;
12
 
        return y;
 
10
    volatile double y;
 
11
    y = x;
 
12
    return y;
13
13
}
14
14
#endif
15
15
 
34
34
#ifndef HAVE_HYPOT
35
35
double hypot(double x, double y)
36
36
{
37
 
        double yx;
 
37
    double yx;
38
38
 
39
 
        x = fabs(x);
40
 
        y = fabs(y);
41
 
        if (x < y) {
42
 
                double temp = x;
43
 
                x = y;
44
 
                y = temp;
45
 
        }
46
 
        if (x == 0.)
47
 
                return 0.;
48
 
        else {
49
 
                yx = y/x;
50
 
                return x*sqrt(1.+yx*yx);
51
 
        }
 
39
    x = fabs(x);
 
40
    y = fabs(y);
 
41
    if (x < y) {
 
42
        double temp = x;
 
43
        x = y;
 
44
        y = temp;
 
45
    }
 
46
    if (x == 0.)
 
47
        return 0.;
 
48
    else {
 
49
        yx = y/x;
 
50
        return x*sqrt(1.+yx*yx);
 
51
    }
52
52
}
53
53
#endif /* HAVE_HYPOT */
54
54
 
56
56
double
57
57
copysign(double x, double y)
58
58
{
59
 
        /* use atan2 to distinguish -0. from 0. */
60
 
        if (y > 0. || (y == 0. && atan2(y, -1.) > 0.)) {
61
 
                return fabs(x);
62
 
        } else {
63
 
                return -fabs(x);
64
 
        }
 
59
    /* use atan2 to distinguish -0. from 0. */
 
60
    if (y > 0. || (y == 0. && atan2(y, -1.) > 0.)) {
 
61
        return fabs(x);
 
62
    } else {
 
63
        return -fabs(x);
 
64
    }
65
65
}
66
66
#endif /* HAVE_COPYSIGN */
67
67
 
73
73
    absx = fabs(x);
74
74
    y = floor(absx);
75
75
    if (absx - y >= 0.5)
76
 
        y += 1.0;
 
76
    y += 1.0;
77
77
    return copysign(y, x);
78
78
}
79
79
#endif /* HAVE_ROUND */
84
84
double
85
85
log1p(double x)
86
86
{
87
 
        /* For x small, we use the following approach.  Let y be the nearest
88
 
           float to 1+x, then
89
 
 
90
 
             1+x = y * (1 - (y-1-x)/y)
91
 
 
92
 
           so log(1+x) = log(y) + log(1-(y-1-x)/y).  Since (y-1-x)/y is tiny,
93
 
           the second term is well approximated by (y-1-x)/y.  If abs(x) >=
94
 
           DBL_EPSILON/2 or the rounding-mode is some form of round-to-nearest
95
 
           then y-1-x will be exactly representable, and is computed exactly
96
 
           by (y-1)-x.
97
 
 
98
 
           If abs(x) < DBL_EPSILON/2 and the rounding mode is not known to be
99
 
           round-to-nearest then this method is slightly dangerous: 1+x could
100
 
           be rounded up to 1+DBL_EPSILON instead of down to 1, and in that
101
 
           case y-1-x will not be exactly representable any more and the
102
 
           result can be off by many ulps.  But this is easily fixed: for a
103
 
           floating-point number |x| < DBL_EPSILON/2., the closest
104
 
           floating-point number to log(1+x) is exactly x.
105
 
        */
106
 
 
107
 
        double y;
108
 
        if (fabs(x) < DBL_EPSILON/2.) {
109
 
                return x;
110
 
        } else if (-0.5 <= x && x <= 1.) {
111
 
                /* WARNING: it's possible than an overeager compiler
112
 
                   will incorrectly optimize the following two lines
113
 
                   to the equivalent of "return log(1.+x)". If this
114
 
                   happens, then results from log1p will be inaccurate
115
 
                   for small x. */
116
 
                y = 1.+x;
117
 
                return log(y)-((y-1.)-x)/y;
118
 
        } else {
119
 
                /* NaNs and infinities should end up here */
120
 
                return log(1.+x);
121
 
        }
 
87
    /* For x small, we use the following approach.  Let y be the nearest
 
88
       float to 1+x, then
 
89
 
 
90
         1+x = y * (1 - (y-1-x)/y)
 
91
 
 
92
       so log(1+x) = log(y) + log(1-(y-1-x)/y).  Since (y-1-x)/y is tiny,
 
93
       the second term is well approximated by (y-1-x)/y.  If abs(x) >=
 
94
       DBL_EPSILON/2 or the rounding-mode is some form of round-to-nearest
 
95
       then y-1-x will be exactly representable, and is computed exactly
 
96
       by (y-1)-x.
 
97
 
 
98
       If abs(x) < DBL_EPSILON/2 and the rounding mode is not known to be
 
99
       round-to-nearest then this method is slightly dangerous: 1+x could
 
100
       be rounded up to 1+DBL_EPSILON instead of down to 1, and in that
 
101
       case y-1-x will not be exactly representable any more and the
 
102
       result can be off by many ulps.  But this is easily fixed: for a
 
103
       floating-point number |x| < DBL_EPSILON/2., the closest
 
104
       floating-point number to log(1+x) is exactly x.
 
105
    */
 
106
 
 
107
    double y;
 
108
    if (fabs(x) < DBL_EPSILON/2.) {
 
109
        return x;
 
110
    } else if (-0.5 <= x && x <= 1.) {
 
111
        /* WARNING: it's possible than an overeager compiler
 
112
           will incorrectly optimize the following two lines
 
113
           to the equivalent of "return log(1.+x)". If this
 
114
           happens, then results from log1p will be inaccurate
 
115
           for small x. */
 
116
        y = 1.+x;
 
117
        return log(y)-((y-1.)-x)/y;
 
118
    } else {
 
119
        /* NaNs and infinities should end up here */
 
120
        return log(1.+x);
 
121
    }
122
122
}
123
123
#endif /* HAVE_LOG1P */
124
124
 
128
128
 *
129
129
 * Developed at SunPro, a Sun Microsystems, Inc. business.
130
130
 * Permission to use, copy, modify, and distribute this
131
 
 * software is freely granted, provided that this notice 
 
131
 * software is freely granted, provided that this notice
132
132
 * is preserved.
133
133
 * ====================================================
134
134
 */
140
140
 
141
141
/* asinh(x)
142
142
 * Method :
143
 
 *      Based on 
144
 
 *              asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
145
 
 *      we have
146
 
 *      asinh(x) := x  if  1+x*x=1,
147
 
 *               := sign(x)*(log(x)+ln2)) for large |x|, else
148
 
 *               := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
149
 
 *               := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))  
 
143
 *      Based on
 
144
 *              asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
 
145
 *      we have
 
146
 *      asinh(x) := x  if  1+x*x=1,
 
147
 *               := sign(x)*(log(x)+ln2)) for large |x|, else
 
148
 *               := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
 
149
 *               := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
150
150
 */
151
151
 
152
152
#ifndef HAVE_ASINH
153
153
double
154
154
asinh(double x)
155
 
{       
156
 
        double w;
157
 
        double absx = fabs(x);
158
 
 
159
 
        if (Py_IS_NAN(x) || Py_IS_INFINITY(x)) {
160
 
                return x+x;
161
 
        }
162
 
        if (absx < two_pow_m28) {       /* |x| < 2**-28 */
163
 
                return x;       /* return x inexact except 0 */
164
 
        } 
165
 
        if (absx > two_pow_p28) {       /* |x| > 2**28 */
166
 
                w = log(absx)+ln2;
167
 
        }
168
 
        else if (absx > 2.0) {          /* 2 < |x| < 2**28 */
169
 
                w = log(2.0*absx + 1.0 / (sqrt(x*x + 1.0) + absx));
170
 
        }
171
 
        else {                          /* 2**-28 <= |x| < 2= */
172
 
                double t = x*x;
173
 
                w = log1p(absx + t / (1.0 + sqrt(1.0 + t)));
174
 
        }
175
 
        return copysign(w, x);
176
 
        
 
155
{
 
156
    double w;
 
157
    double absx = fabs(x);
 
158
 
 
159
    if (Py_IS_NAN(x) || Py_IS_INFINITY(x)) {
 
160
        return x+x;
 
161
    }
 
162
    if (absx < two_pow_m28) {           /* |x| < 2**-28 */
 
163
        return x;               /* return x inexact except 0 */
 
164
    }
 
165
    if (absx > two_pow_p28) {           /* |x| > 2**28 */
 
166
        w = log(absx)+ln2;
 
167
    }
 
168
    else if (absx > 2.0) {              /* 2 < |x| < 2**28 */
 
169
        w = log(2.0*absx + 1.0 / (sqrt(x*x + 1.0) + absx));
 
170
    }
 
171
    else {                              /* 2**-28 <= |x| < 2= */
 
172
        double t = x*x;
 
173
        w = log1p(absx + t / (1.0 + sqrt(1.0 + t)));
 
174
    }
 
175
    return copysign(w, x);
 
176
 
177
177
}
178
178
#endif /* HAVE_ASINH */
179
179
 
180
180
/* acosh(x)
181
181
 * Method :
182
182
 *      Based on
183
 
 *            acosh(x) = log [ x + sqrt(x*x-1) ]
 
183
 *            acosh(x) = log [ x + sqrt(x*x-1) ]
184
184
 *      we have
185
 
 *            acosh(x) := log(x)+ln2, if x is large; else
186
 
 *            acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
187
 
 *            acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
 
185
 *            acosh(x) := log(x)+ln2, if x is large; else
 
186
 *            acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
 
187
 *            acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
188
188
 *
189
189
 * Special cases:
190
190
 *      acosh(x) is NaN with signal if x<1.
195
195
double
196
196
acosh(double x)
197
197
{
198
 
        if (Py_IS_NAN(x)) {
199
 
                return x+x;
200
 
        }
201
 
        if (x < 1.) {                   /* x < 1;  return a signaling NaN */
202
 
                errno = EDOM;
 
198
    if (Py_IS_NAN(x)) {
 
199
        return x+x;
 
200
    }
 
201
    if (x < 1.) {                       /* x < 1;  return a signaling NaN */
 
202
        errno = EDOM;
203
203
#ifdef Py_NAN
204
 
                return Py_NAN;
 
204
        return Py_NAN;
205
205
#else
206
 
                return (x-x)/(x-x);
 
206
        return (x-x)/(x-x);
207
207
#endif
208
 
        }
209
 
        else if (x >= two_pow_p28) {    /* x > 2**28 */
210
 
                if (Py_IS_INFINITY(x)) {
211
 
                        return x+x;
212
 
                } else {
213
 
                        return log(x)+ln2;      /* acosh(huge)=log(2x) */
214
 
                }
215
 
        }
216
 
        else if (x == 1.) {
217
 
                return 0.0;                     /* acosh(1) = 0 */
218
 
        }
219
 
        else if (x > 2.) {                      /* 2 < x < 2**28 */
220
 
                double t = x*x;
221
 
                return log(2.0*x - 1.0 / (x + sqrt(t - 1.0)));
222
 
        }
223
 
        else {                          /* 1 < x <= 2 */
224
 
                double t = x - 1.0;
225
 
                return log1p(t + sqrt(2.0*t + t*t));
226
 
        }
 
208
    }
 
209
    else if (x >= two_pow_p28) {        /* x > 2**28 */
 
210
        if (Py_IS_INFINITY(x)) {
 
211
            return x+x;
 
212
        } else {
 
213
            return log(x)+ln2;                  /* acosh(huge)=log(2x) */
 
214
        }
 
215
    }
 
216
    else if (x == 1.) {
 
217
        return 0.0;                             /* acosh(1) = 0 */
 
218
    }
 
219
    else if (x > 2.) {                          /* 2 < x < 2**28 */
 
220
        double t = x*x;
 
221
        return log(2.0*x - 1.0 / (x + sqrt(t - 1.0)));
 
222
    }
 
223
    else {                              /* 1 < x <= 2 */
 
224
        double t = x - 1.0;
 
225
        return log1p(t + sqrt(2.0*t + t*t));
 
226
    }
227
227
}
228
228
#endif /* HAVE_ACOSH */
229
229
 
231
231
 * Method :
232
232
 *    1.Reduced x to positive by atanh(-x) = -atanh(x)
233
233
 *    2.For x>=0.5
234
 
 *                1           2x                          x
 
234
 *                1           2x                          x
235
235
 *      atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
236
 
 *                2          1 - x                    1 - x
 
236
 *                2          1 - x                    1 - x
237
237
 *
238
238
 *      For x<0.5
239
239
 *      atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
248
248
double
249
249
atanh(double x)
250
250
{
251
 
        double absx;
252
 
        double t;
 
251
    double absx;
 
252
    double t;
253
253
 
254
 
        if (Py_IS_NAN(x)) {
255
 
                return x+x;
256
 
        }
257
 
        absx = fabs(x);
258
 
        if (absx >= 1.) {               /* |x| >= 1 */
259
 
                errno = EDOM;
 
254
    if (Py_IS_NAN(x)) {
 
255
        return x+x;
 
256
    }
 
257
    absx = fabs(x);
 
258
    if (absx >= 1.) {                   /* |x| >= 1 */
 
259
        errno = EDOM;
260
260
#ifdef Py_NAN
261
 
                return Py_NAN;
 
261
        return Py_NAN;
262
262
#else
263
 
                return x/zero;
 
263
        return x/zero;
264
264
#endif
265
 
        }
266
 
        if (absx < two_pow_m28) {       /* |x| < 2**-28 */
267
 
                return x;
268
 
        }
269
 
        if (absx < 0.5) {               /* |x| < 0.5 */
270
 
                t = absx+absx;
271
 
                t = 0.5 * log1p(t + t*absx / (1.0 - absx));
272
 
        } 
273
 
        else {                          /* 0.5 <= |x| <= 1.0 */
274
 
                t = 0.5 * log1p((absx + absx) / (1.0 - absx));
275
 
        }
276
 
        return copysign(t, x);
 
265
    }
 
266
    if (absx < two_pow_m28) {           /* |x| < 2**-28 */
 
267
        return x;
 
268
    }
 
269
    if (absx < 0.5) {                   /* |x| < 0.5 */
 
270
        t = absx+absx;
 
271
        t = 0.5 * log1p(t + t*absx / (1.0 - absx));
 
272
    }
 
273
    else {                              /* 0.5 <= |x| <= 1.0 */
 
274
        t = 0.5 * log1p((absx + absx) / (1.0 - absx));
 
275
    }
 
276
    return copysign(t, x);
277
277
}
278
278
#endif /* HAVE_ATANH */