1*05a0b428SJohn Marino /* s_log1pf.c -- float version of s_log1p.c.
2*05a0b428SJohn Marino * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
3*05a0b428SJohn Marino */
4*05a0b428SJohn Marino
5*05a0b428SJohn Marino /*
6*05a0b428SJohn Marino * ====================================================
7*05a0b428SJohn Marino * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8*05a0b428SJohn Marino *
9*05a0b428SJohn Marino * Developed at SunPro, a Sun Microsystems, Inc. business.
10*05a0b428SJohn Marino * Permission to use, copy, modify, and distribute this
11*05a0b428SJohn Marino * software is freely granted, provided that this notice
12*05a0b428SJohn Marino * is preserved.
13*05a0b428SJohn Marino * ====================================================
14*05a0b428SJohn Marino */
15*05a0b428SJohn Marino
16*05a0b428SJohn Marino #include "math.h"
17*05a0b428SJohn Marino #include "math_private.h"
18*05a0b428SJohn Marino
19*05a0b428SJohn Marino static const float
20*05a0b428SJohn Marino ln2_hi = 6.9313812256e-01, /* 0x3f317180 */
21*05a0b428SJohn Marino ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */
22*05a0b428SJohn Marino two25 = 3.355443200e+07, /* 0x4c000000 */
23*05a0b428SJohn Marino Lp1 = 6.6666668653e-01, /* 3F2AAAAB */
24*05a0b428SJohn Marino Lp2 = 4.0000000596e-01, /* 3ECCCCCD */
25*05a0b428SJohn Marino Lp3 = 2.8571429849e-01, /* 3E924925 */
26*05a0b428SJohn Marino Lp4 = 2.2222198546e-01, /* 3E638E29 */
27*05a0b428SJohn Marino Lp5 = 1.8183572590e-01, /* 3E3A3325 */
28*05a0b428SJohn Marino Lp6 = 1.5313838422e-01, /* 3E1CD04F */
29*05a0b428SJohn Marino Lp7 = 1.4798198640e-01; /* 3E178897 */
30*05a0b428SJohn Marino
31*05a0b428SJohn Marino static const float zero = 0.0;
32*05a0b428SJohn Marino
33*05a0b428SJohn Marino float
log1pf(float x)34*05a0b428SJohn Marino log1pf(float x)
35*05a0b428SJohn Marino {
36*05a0b428SJohn Marino float hfsq,f,c,s,z,R,u;
37*05a0b428SJohn Marino int32_t k,hx,hu,ax;
38*05a0b428SJohn Marino
39*05a0b428SJohn Marino GET_FLOAT_WORD(hx,x);
40*05a0b428SJohn Marino ax = hx&0x7fffffff;
41*05a0b428SJohn Marino
42*05a0b428SJohn Marino k = 1;
43*05a0b428SJohn Marino if (hx < 0x3ed413d7) { /* x < 0.41422 */
44*05a0b428SJohn Marino if(ax>=0x3f800000) { /* x <= -1.0 */
45*05a0b428SJohn Marino if(x==(float)-1.0) return -two25/zero; /* log1p(-1)=+inf */
46*05a0b428SJohn Marino else return (x-x)/(x-x); /* log1p(x<-1)=NaN */
47*05a0b428SJohn Marino }
48*05a0b428SJohn Marino if(ax<0x31000000) { /* |x| < 2**-29 */
49*05a0b428SJohn Marino if(two25+x>zero /* raise inexact */
50*05a0b428SJohn Marino &&ax<0x24800000) /* |x| < 2**-54 */
51*05a0b428SJohn Marino return x;
52*05a0b428SJohn Marino else
53*05a0b428SJohn Marino return x - x*x*(float)0.5;
54*05a0b428SJohn Marino }
55*05a0b428SJohn Marino if(hx>0||hx<=((int32_t)0xbe95f61f)) {
56*05a0b428SJohn Marino k=0;f=x;hu=1;} /* -0.2929<x<0.41422 */
57*05a0b428SJohn Marino }
58*05a0b428SJohn Marino if (hx >= 0x7f800000) return x+x;
59*05a0b428SJohn Marino if(k!=0) {
60*05a0b428SJohn Marino if(hx<0x5a000000) {
61*05a0b428SJohn Marino u = (float)1.0+x;
62*05a0b428SJohn Marino GET_FLOAT_WORD(hu,u);
63*05a0b428SJohn Marino k = (hu>>23)-127;
64*05a0b428SJohn Marino /* correction term */
65*05a0b428SJohn Marino c = (k>0)? (float)1.0-(u-x):x-(u-(float)1.0);
66*05a0b428SJohn Marino c /= u;
67*05a0b428SJohn Marino } else {
68*05a0b428SJohn Marino u = x;
69*05a0b428SJohn Marino GET_FLOAT_WORD(hu,u);
70*05a0b428SJohn Marino k = (hu>>23)-127;
71*05a0b428SJohn Marino c = 0;
72*05a0b428SJohn Marino }
73*05a0b428SJohn Marino hu &= 0x007fffff;
74*05a0b428SJohn Marino if(hu<0x3504f7) {
75*05a0b428SJohn Marino SET_FLOAT_WORD(u,hu|0x3f800000);/* normalize u */
76*05a0b428SJohn Marino } else {
77*05a0b428SJohn Marino k += 1;
78*05a0b428SJohn Marino SET_FLOAT_WORD(u,hu|0x3f000000); /* normalize u/2 */
79*05a0b428SJohn Marino hu = (0x00800000-hu)>>2;
80*05a0b428SJohn Marino }
81*05a0b428SJohn Marino f = u-(float)1.0;
82*05a0b428SJohn Marino }
83*05a0b428SJohn Marino hfsq=(float)0.5*f*f;
84*05a0b428SJohn Marino if(hu==0) { /* |f| < 2**-20 */
85*05a0b428SJohn Marino if(f==zero) if(k==0) return zero;
86*05a0b428SJohn Marino else {c += k*ln2_lo; return k*ln2_hi+c;}
87*05a0b428SJohn Marino R = hfsq*((float)1.0-(float)0.66666666666666666*f);
88*05a0b428SJohn Marino if(k==0) return f-R; else
89*05a0b428SJohn Marino return k*ln2_hi-((R-(k*ln2_lo+c))-f);
90*05a0b428SJohn Marino }
91*05a0b428SJohn Marino s = f/((float)2.0+f);
92*05a0b428SJohn Marino z = s*s;
93*05a0b428SJohn Marino R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))));
94*05a0b428SJohn Marino if(k==0) return f-(hfsq-s*(hfsq+R)); else
95*05a0b428SJohn Marino return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
96*05a0b428SJohn Marino }
97