← Back to branch summary

~ubuntu-branches/ubuntu/precise/eglibc/precise

« back to all changes in this revision

Viewing changes to sysdeps/ieee754/flt-32/e_lgammaf_r.c

  • Committer: Package Import Robot
  • Author(s): Matthias Klose
  • Date: 2011-10-04 17:48:26 UTC
  • mfrom: (216.1.23 oneiric)
  • Revision ID: package-import@ubuntu.com-20111004174826-2cyb9ewn3ucymlsx
Tags: 2.13-20ubuntu5
https://launchpad.net/bugs/853688
libc6-dev: Don't break the current {gnat,gcj}-4.4-base versons. LP: #853688.

Show diffs side-by-side

added added

removed removed

Lines of Context:
sysdeps/ieee754/flt-32/e_lgammaf_r.c
13
13
 * ====================================================
14
14
 */
15
15
 
 
16
#if defined(LIBM_SCCS) && !defined(lint)
 
17
static char rcsid[] = "$NetBSD: e_lgammaf_r.c,v 1.3 1995/05/10 20:45:47 jtc Exp $";
 
18
#endif
 
19
 
16
20
#include "math.h"
17
21
#include "math_private.h"
18
22
 
 
23
#ifdef __STDC__
19
24
static const float
 
25
#else
 
26
static float
 
27
#endif
20
28
two23=  8.3886080000e+06, /* 0x4b000000 */
21
29
half=  5.0000000000e-01, /* 0x3f000000 */
22
30
one =  1.0000000000e+00, /* 0x3f800000 */
84
92
w5  =  8.3633989561e-04, /* 0x3a5b3dd2 */
85
93
w6  = -1.6309292987e-03; /* 0xbad5c4e8 */
86
94
 
 
95
#ifdef __STDC__
87
96
static const float zero=  0.0000000000e+00;
 
97
#else
 
98
static float zero=  0.0000000000e+00;
 
99
#endif
88
100
 
89
 
static float
90
 
sin_pif(float x)
 
101
#ifdef __STDC__
 
102
        static float sin_pif(float x)
 
103
#else
 
104
        static float sin_pif(x)
 
105
        float x;
 
106
#endif
91
107
{
92
108
        float y,z;
93
109
        int n,ix;
108
124
            y   = (float)2.0*(y - __floorf(y)); /* y = |x| mod 2.0 */
109
125
            n   = (int) (y*(float)4.0);
110
126
        } else {
111
 
            if(ix>=0x4b800000) {
112
 
                y = zero; n = 0;                 /* y must be even */
113
 
            } else {
114
 
                if(ix<0x4b000000) z = y+two23;  /* exact */
 
127
            if(ix>=0x4b800000) {
 
128
                y = zero; n = 0;                 /* y must be even */
 
129
            } else {
 
130
                if(ix<0x4b000000) z = y+two23;  /* exact */
115
131
                GET_FLOAT_WORD(n,z);
116
132
                n &= 1;
117
 
                y  = n;
118
 
                n<<= 2;
119
 
            }
120
 
        }
 
133
                y  = n;
 
134
                n<<= 2;
 
135
            }
 
136
        }
121
137
        switch (n) {
122
138
            case 0:   y =  __kernel_sinf(pi*y,zero,0); break;
123
139
            case 1:
132
148
}
133
149
 
134
150
 
135
 
float
136
 
__ieee754_lgammaf_r(float x, int *signgamp)
 
151
#ifdef __STDC__
 
152
        float __ieee754_lgammaf_r(float x, int *signgamp)
 
153
#else
 
154
        float __ieee754_lgammaf_r(x,signgamp)
 
155
        float x; int *signgamp;
 
156
#endif
137
157
{
138
158
        float t,y,z,nadj,p,p1,p2,p3,q,r,w;
139
159
        int i,hx,ix;
143
163
    /* purge off +-inf, NaN, +-0, and negative arguments */
144
164
        *signgamp = 1;
145
165
        ix = hx&0x7fffffff;
146
 
        if(__builtin_expect(ix>=0x7f800000, 0)) return x*x;
147
 
        if(__builtin_expect(ix==0, 0))
 
166
        if(ix>=0x7f800000) return x*x;
 
167
        if(ix==0)
148
168
          {
149
169
            if (hx < 0)
150
170
              *signgamp = -1;
151
171
            return one/fabsf(x);
152
172
          }
153
 
        if(__builtin_expect(ix<0x1c800000, 0)) {
154
 
            /* |x|<2**-70, return -log(|x|) */
 
173
        if(ix<0x1c800000) {     /* |x|<2**-70, return -log(|x|) */
155
174
            if(hx<0) {
156
 
                *signgamp = -1;
157
 
                return -__ieee754_logf(-x);
 
175
                *signgamp = -1;
 
176
                return -__ieee754_logf(-x);
158
177
            } else return -__ieee754_logf(x);
159
178
        }
160
179
        if(hx<0) {
161
 
            if(ix>=0x4b000000)  /* |x|>=2**23, must be -integer */
 
180
            if(ix>=0x4b000000)  /* |x|>=2**23, must be -integer */
162
181
                return x/zero;
163
182
            t = sin_pif(x);
164
183
            if(t==zero) return one/fabsf(t); /* -integer */
171
190
        if (ix==0x3f800000||ix==0x40000000) r = 0;
172
191
    /* for x < 2.0 */
173
192
        else if(ix<0x40000000) {
174
 
            if(ix<=0x3f666666) {        /* lgamma(x) = lgamma(x+1)-log(x) */
 
193
            if(ix<=0x3f666666) {        /* lgamma(x) = lgamma(x+1)-log(x) */
175
194
                r = -__ieee754_logf(x);
176
195
                if(ix>=0x3f3b4a20) {y = one-x; i= 0;}
177
196
                else if(ix>=0x3e6d3308) {y= x-(tc-one); i=1;}
178
 
                else {y = x; i=2;}
 
197
                else {y = x; i=2;}
179
198
            } else {
180
 
                r = zero;
181
 
                if(ix>=0x3fdda618) {y=(float)2.0-x;i=0;} /* [1.7316,2] */
182
 
                else if(ix>=0x3F9da620) {y=x-tc;i=1;} /* [1.23,1.73] */
 
199
                r = zero;
 
200
                if(ix>=0x3fdda618) {y=(float)2.0-x;i=0;} /* [1.7316,2] */
 
201
                else if(ix>=0x3F9da620) {y=x-tc;i=1;} /* [1.23,1.73] */
183
202
                else {y=x-one;i=2;}
184
203
            }
185
204
            switch(i) {
203
222
                r += (-(float)0.5*y + p1/p2);
204
223
            }
205
224
        }
206
 
        else if(ix<0x41000000) {                        /* x < 8.0 */
 
225
        else if(ix<0x41000000) {                        /* x < 8.0 */
207
226
            i = (int)x;
208
227
            t = zero;
209
228
            y = x-(float)i;
232
251
        if(hx<0) r = nadj - r;
233
252
        return r;
234
253
}
235
 
strong_alias (__ieee754_lgammaf_r, __lgammaf_r_finite)