amath  1.8.5
Simple command line calculator
log1p.c File Reference
#include "prim.h"
Include dependency graph for log1p.c:

Go to the source code of this file.

Functions

double log1p (double x)
 

Variables

static const double ln2_hi = 6.93147180369123816490e-01
 
static const double ln2_lo = 1.90821492927058770002e-10
 
static const double two54 = 1.80143985094819840000e+16
 
static const double Lp1 = 6.666666666666735130e-01
 
static const double Lp2 = 3.999999999940941908e-01
 
static const double Lp3 = 2.857142874366239149e-01
 
static const double Lp4 = 2.222219843214978396e-01
 
static const double Lp5 = 1.818357216161805012e-01
 
static const double Lp6 = 1.531383769920937332e-01
 
static const double Lp7 = 1.479819860511658591e-01
 
static double zero = 0.0
 

Function Documentation

◆ log1p()

double log1p ( double  x)

Definition at line 122 of file log1p.c.

References ln2_hi, ln2_lo, Lp1, Lp2, Lp3, Lp4, Lp5, Lp6, Lp7, two54, and zero.

Referenced by acosh(), asinh(), and atanh().

123 {
124  double hfsq,f,c,s,z,R,u;
125  int32_t k,hx,hu,ax;
126 
127  f = 0.0;
128  c = 0.0;
129  hu = 0;
130 
131  GET_HIGH_WORD(hx,x); /* high word of x */
132  ax = hx&0x7fffffff;
133 
134  k = 1;
135  if (hx < 0x3FDA827A) { /* x < 0.41422 */
136  if(ax>=0x3ff00000) { /* x <= -1.0 */
137  if(x==-1.0) return -two54/zero; /* log1p(-1)=+inf */
138  else return (x-x)/(x-x); /* log1p(x<-1)=NaN */
139  }
140  if(ax<0x3e200000) { /* |x| < 2**-29 */
141  if(two54+x>zero /* raise inexact */
142  &&ax<0x3c900000) /* |x| < 2**-54 */
143  return x;
144  else
145  return x - x*x*0.5;
146  }
147  if(hx>0||hx<=((int)0xbfd2bec3)) {
148  k=0;
149  f=x;
150  hu=1;
151  } /* -0.2929<x<0.41422 */
152  }
153  if (hx >= 0x7ff00000) return x+x;
154  if(k!=0) {
155  if(hx<0x43400000) {
156  u = 1.0+x;
157  GET_HIGH_WORD(hu,u); /* high word of u */
158  k = (hu>>20)-1023;
159  c = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */
160  c /= u;
161  } else {
162  u = x;
163  GET_HIGH_WORD(hu,u); /* high word of u */
164  k = (hu>>20)-1023;
165  c = 0;
166  }
167  hu &= 0x000fffff;
168  if(hu<0x6a09e) {
169  SET_HIGH_WORD(u, hu|0x3ff00000); /* normalize u */
170  } else {
171  k += 1;
172  SET_HIGH_WORD(u, hu|0x3fe00000); /* normalize u/2 */
173  hu = (0x00100000-hu)>>2;
174  }
175  f = u-1.0;
176  }
177  hfsq=0.5*f*f;
178  if(hu==0) { /* |f| < 2**-20 */
179  if(f==zero) {
180  if(k==0) return zero;
181  else {
182  c += k*ln2_lo;
183  return k*ln2_hi+c;
184  }
185  }
186  R = hfsq*(1.0-0.66666666666666666*f);
187  if(k==0) return f-R;
188  else
189  return k*ln2_hi-((R-(k*ln2_lo+c))-f);
190  }
191  s = f/(2.0+f);
192  z = s*s;
193  R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))));
194  if(k==0) return f-(hfsq-s*(hfsq+R));
195  else
196  return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
197 }
static const double Lp6
Definition: log1p.c:117
static double zero
Definition: log1p.c:120
#define GET_HIGH_WORD(i, d)
Get the more significant 32 bit int from a double.
Definition: prim.h:167
static const double Lp4
Definition: log1p.c:115
static const double ln2_lo
Definition: log1p.c:110
static const double Lp3
Definition: log1p.c:114
static const double Lp5
Definition: log1p.c:116
static const double two54
Definition: log1p.c:111
static const double Lp2
Definition: log1p.c:113
static const double ln2_hi
Definition: log1p.c:109
static const double Lp7
Definition: log1p.c:118
static const double Lp1
Definition: log1p.c:112
#define SET_HIGH_WORD(d, v)
Set the more significant 32 bits of a double from an int.
Definition: prim.h:198
Here is the caller graph for this function:

Variable Documentation

◆ ln2_hi

const double ln2_hi = 6.93147180369123816490e-01
static

Definition at line 109 of file log1p.c.

Referenced by log1p().

◆ ln2_lo

const double ln2_lo = 1.90821492927058770002e-10
static

Definition at line 110 of file log1p.c.

Referenced by log1p().

◆ Lp1

const double Lp1 = 6.666666666666735130e-01
static

Definition at line 112 of file log1p.c.

Referenced by log1p().

◆ Lp2

const double Lp2 = 3.999999999940941908e-01
static

Definition at line 113 of file log1p.c.

Referenced by log1p().

◆ Lp3

const double Lp3 = 2.857142874366239149e-01
static

Definition at line 114 of file log1p.c.

Referenced by log1p().

◆ Lp4

const double Lp4 = 2.222219843214978396e-01
static

Definition at line 115 of file log1p.c.

Referenced by log1p().

◆ Lp5

const double Lp5 = 1.818357216161805012e-01
static

Definition at line 116 of file log1p.c.

Referenced by log1p().

◆ Lp6

const double Lp6 = 1.531383769920937332e-01
static

Definition at line 117 of file log1p.c.

Referenced by log1p().

◆ Lp7

const double Lp7 = 1.479819860511658591e-01
static

Definition at line 118 of file log1p.c.

Referenced by log1p().

◆ two54

const double two54 = 1.80143985094819840000e+16
static

Definition at line 111 of file log1p.c.

Referenced by log1p().

◆ zero

double zero = 0.0
static

Definition at line 120 of file log1p.c.

Referenced by log1p().