amath  1.8.5
Simple command line calculator
Functions | Variables
pow.c File Reference

Expontation function. More...

#include "prim.h"
Include dependency graph for pow.c:

Go to the source code of this file.

Functions

double pow (double x, double y)
 Expontation function. More...
 

Variables

static const double bp []
 
static const double dp_h []
 
static const double dp_l []
 
static const double zero = 0.0
 
static const double one = 1.0
 
static const double two = 2.0
 
static const double two53 = 9007199254740992.0
 
static const double huge = 1.0e300
 
static const double tiny = 1.0e-300
 
static const double L1 = 5.99999999999994648725e-01
 
static const double L2 = 4.28571428578550184252e-01
 
static const double L3 = 3.33333329818377432918e-01
 
static const double L4 = 2.72728123808534006489e-01
 
static const double L5 = 2.30660745775561754067e-01
 
static const double L6 = 2.06975017800338417784e-01
 
static const double P1 = 1.66666666666666019037e-01
 
static const double P2 = -2.77777777770155933842e-03
 
static const double P3 = 6.61375632143793436117e-05
 
static const double P4 = -1.65339022054652515390e-06
 
static const double P5 = 4.13813679705723846039e-08
 
static const double lg2 = 6.93147180559945286227e-01
 
static const double lg2_h = 6.93147182464599609375e-01
 
static const double lg2_l = -1.90465429995776804525e-09
 
static const double ovt = 8.0085662595372944372e-0017
 
static const double cp = 9.61796693925975554329e-01
 
static const double cp_h = 9.61796700954437255859e-01
 
static const double cp_l = -7.02846165095275826516e-09
 
static const double ivln2 = 1.44269504088896338700e+00
 
static const double ivln2_h = 1.44269502162933349609e+00
 
static const double ivln2_l = 1.92596299112661746887e-08
 

Detailed Description

Expontation function.

Definition in file pow.c.

Function Documentation

◆ pow()

double pow ( double  x,
double  y 
)

Expontation function. 

Method:  Let x =  2   * (1+f)
 1. Compute and return log2(x) in two pieces:
    log2(x) = w1 + w2,
    where w1 has 53-24 = 29 bit trailing zeros.
 2. Perform y*log2(x) = n+y' by simulating muti-precision
    arithmetic, where |y'|<=0.5.
 3. Return x**y = 2**n*exp(y'*log2)
Special cases:
 1.  (anything) ** 0  is 1
 2.  (anything) ** 1  is itself
 3.  (anything) ** NAN is NAN
 4.  NAN ** (anything except 0) is NAN
 5.  +-(|x| > 1) **  +INF is +INF
 6.  +-(|x| > 1) **  -INF is +0
 7.  +-(|x| < 1) **  +INF is +0
 8.  +-(|x| < 1) **  -INF is +INF
 9.  +-1         ** +-INF is NAN
 10. +0 ** (+anything except 0, NAN)               is +0
 11. -0 ** (+anything except 0, NAN, odd integer)  is +0
 12. +0 ** (-anything except 0, NAN)               is +INF
 13. -0 ** (-anything except 0, NAN, odd integer)  is +INF
 14. -0 ** (odd integer) = -( +0 ** (odd integer) )
 15. +INF ** (+anything except 0,NAN) is +INF
 16. +INF ** (-anything except 0,NAN) is +0
 17. -INF ** (anything)  = -0 ** (-anything)
 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
 19. (-anything except 0 and inf) ** (non-integer) is NAN
Accuracy:
 pow(x,y) returns x**y nearly rounded. In particular
        pow(integer,integer)
 always returns the correct integer provided it is
 representable.
Constants :
The hexadecimal values are the intended ones for the following
constants. The decimal values may be used, provided that the
compiler will convert from decimal to binary accurately enough
to produce the hexadecimal values shown.

Definition at line 138 of file pow.c.

References bp, cp, cp_h, cp_l, dp_h, dp_l, fabs(), huge, ivln2, ivln2_h, ivln2_l, L1, L2, L3, L4, L5, L6, lg2, lg2_h, lg2_l, one, ovt, P1, P2, P3, P4, P5, scalbn(), sqrt(), tiny, two, two53, and zero.

Referenced by cpow(), PositionalNumeralSystem::GetText(), PositionalNumeralSystem::Parse(), and RealNumber::Raise().

139 {
140  double z, ax, z_h, z_l, p_h, p_l;
141  double y1, t1, t2, r, s, t, u, v, w;
142  int32_t i, j, k, yisint, n;
143  int32_t hx, hy, ix, iy;
144  uint32_t lx, ly;
145 
146  EXTRACT_WORDS(hx, lx, x);
147  EXTRACT_WORDS(hy, ly, y);
148  ix = hx & 0x7fffffff;
149  iy = hy & 0x7fffffff;
150 
151  /* y==zero: x**0 = 1 */
152  if ((iy | ly) == 0)
153  return one;
154 
155  /* +-NaN return x+y */
156  if (ix > 0x7ff00000 || ((ix == 0x7ff00000) && (lx != 0)) ||
157  iy > 0x7ff00000 || ((iy == 0x7ff00000) && (ly != 0)))
158  return x + y;
159 
160  /* determine if y is an odd int when x < 0
161  * yisint = 0 ... y is not an integer
162  * yisint = 1 ... y is an odd int
163  * yisint = 2 ... y is an even int
164  */
165  yisint = 0;
166  if (hx < 0)
167  {
168  if (iy >= 0x43400000)
169  yisint = 2; /* even integer y */
170  else if (iy >= 0x3ff00000)
171  {
172  k = (iy >> 20) - 0x3ff; /* exponent */
173  if (k > 20)
174  {
175  j = ly >> (52 - k);
176  if ((uint32_t)(j << (52 - k)) == ly)
177  yisint = 2 - (j & 1);
178  }
179  else if (ly == 0)
180  {
181  j = iy >> (20 - k);
182  if ((j << (20 - k)) == iy)
183  yisint = 2 - (j & 1);
184  }
185  }
186  }
187 
188  /* special value of y */
189  if (ly == 0)
190  {
191  if (iy == 0x7ff00000)
192  { /* y is +-inf */
193  if (((ix - 0x3ff00000) | lx) == 0)
194  return y - y; /* inf**+-1 is NaN */
195  else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
196  return (hy >= 0) ? y : zero;
197  else /* (|x|<1)**-,+inf = inf,0 */
198  return (hy < 0) ? -y : zero;
199  }
200  if (iy == 0x3ff00000)
201  { /* y is +-1 */
202  if (hy < 0)
203  return one / x;
204  else
205  return x;
206  }
207  if (hy == 0x40000000)
208  return x * x; /* y is 2 */
209  if (hy == 0x3fe00000)
210  { /* y is 0.5 */
211  if (hx >= 0) /* x >= +0 */
212  return sqrt(x);
213  }
214  }
215 
216  ax = fabs(x);
217  /* special value of x */
218  if (lx == 0)
219  {
220  if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000)
221  {
222  z = ax; /*x is +-0,+-inf,+-1*/
223  if (hy < 0)
224  z = one / z; /* z = (1/|x|) */
225  if (hx < 0)
226  {
227  if (((ix - 0x3ff00000) | yisint) == 0)
228  {
229  z = (z - z) / (z - z); /* (-1)**non-int is NaN */
230  }
231  else if (yisint == 1)
232  z = -z; /* (x<0)**odd = -(|x|**odd) */
233  }
234  return z;
235  }
236  }
237 
238  n = (hx >> 31) + 1;
239 
240  /* (x<0)**(non-int) is NaN */
241  if ((n | yisint) == 0)
242  return (x - x) / (x - x);
243 
244  s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
245  if ((n | (yisint - 1)) == 0)
246  s = -one; /* (-ve)**(odd int) */
247 
248  /* |y| is huge */
249  if (iy > 0x41e00000)
250  { /* if |y| > 2**31 */
251  if (iy > 0x43f00000)
252  { /* if |y| > 2**64, must o/uflow */
253  if (ix <= 0x3fefffff)
254  return (hy < 0) ? huge * huge : tiny * tiny;
255  if (ix >= 0x3ff00000)
256  return (hy > 0) ? huge * huge : tiny * tiny;
257  }
258  /* over/underflow if x is not close to one */
259  if (ix < 0x3fefffff)
260  return (hy < 0) ? s * huge * huge : s * tiny * tiny;
261  if (ix > 0x3ff00000)
262  return (hy > 0) ? s * huge * huge : s * tiny * tiny;
263  /* now |1-x| is tiny <= 2**-20, suffice to compute
264  log(x) by x-x^2/2+x^3/3-x^4/4 */
265  t = ax - one; /* t has 20 trailing zeros */
266  w = (t * t) * (0.5 - t * (0.3333333333333333333333 - t * 0.25));
267  u = ivln2_h * t; /* ivln2_h has 21 sig. bits */
268  v = t * ivln2_l - w * ivln2;
269  t1 = u + v;
270  SET_LOW_WORD(t1, 0);
271  t2 = v - (t1 - u);
272  }
273  else
274  {
275  double ss, s2, s_h, s_l, t_h, t_l;
276  n = 0;
277  /* take care subnormal number */
278  if (ix < 0x00100000)
279  {
280  ax *= two53;
281  n -= 53;
282  GET_HIGH_WORD(ix, ax);
283  }
284  n += ((ix) >> 20) - 0x3ff;
285  j = ix & 0x000fffff;
286  /* determine interval */
287  ix = j | 0x3ff00000; /* normalize ix */
288  if (j <= 0x3988E)
289  k = 0; /* |x|<sqrt(3/2) */
290  else if (j < 0xBB67A)
291  k = 1; /* |x|<sqrt(3) */
292  else
293  {
294  k = 0;
295  n += 1;
296  ix -= 0x00100000;
297  }
298  SET_HIGH_WORD(ax, ix);
299 
300  /* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
301  u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
302  v = one / (ax + bp[k]);
303  ss = u * v;
304  s_h = ss;
305  SET_LOW_WORD(s_h, 0);
306  /* t_h=ax+bp[k] High */
307  t_h = zero;
308  SET_HIGH_WORD(t_h, ((ix >> 1) | 0x20000000) + 0x00080000 + (k << 18));
309  t_l = ax - (t_h - bp[k]);
310  s_l = v * ((u - s_h * t_h) - s_h * t_l);
311  /* compute log(ax) */
312  s2 = ss * ss;
313  r = s2 * s2 * (L1 + s2 * (L2 + s2 * (L3 + s2 * (L4 + s2 * (L5 + s2 * L6)))));
314  r += s_l * (s_h + ss);
315  s2 = s_h * s_h;
316  t_h = 3.0 + s2 + r;
317  SET_LOW_WORD(t_h, 0);
318  t_l = r - ((t_h - 3.0) - s2);
319  /* u+v = ss*(1+...) */
320  u = s_h * t_h;
321  v = s_l * t_h + t_l * ss;
322  /* 2/(3log2)*(ss+...) */
323  p_h = u + v;
324  SET_LOW_WORD(p_h, 0);
325  p_l = v - (p_h - u);
326  z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */
327  z_l = cp_l * p_h + p_l * cp + dp_l[k];
328  /* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */
329  t = (double)n;
330  t1 = (((z_h + z_l) + dp_h[k]) + t);
331  SET_LOW_WORD(t1, 0);
332  t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
333  }
334 
335  /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
336  y1 = y;
337  SET_LOW_WORD(y1, 0);
338  p_l = (y - y1) * t1 + y * t2;
339  p_h = y1 * t1;
340  z = p_l + p_h;
341  EXTRACT_WORDS(j, i, z);
342  if (j >= 0x40900000)
343  { /* z >= 1024 */
344  if (((j - 0x40900000) | i) != 0) /* if z > 1024 */
345  return s * huge * huge; /* overflow */
346  else
347  {
348  if (p_l + ovt > z - p_h)
349  return s * huge * huge; /* overflow */
350  }
351  }
352  else if ((j & 0x7fffffff) >= 0x4090cc00)
353  { /* z <= -1075 */
354  if (((j - 0xc090cc00) | i) != 0) /* z < -1075 */
355  return s * tiny * tiny; /* underflow */
356  else
357  {
358  if (p_l <= z - p_h)
359  return s * tiny * tiny; /* underflow */
360  }
361  }
362  /*
363  * compute 2**(p_h+p_l)
364  */
365  i = j & 0x7fffffff;
366  k = (i >> 20) - 0x3ff;
367  n = 0;
368  if (i > 0x3fe00000)
369  { /* if |z| > 0.5, set n = [z+0.5] */
370  n = j + (0x00100000 >> (k + 1));
371  k = ((n & 0x7fffffff) >> 20) - 0x3ff; /* new k for n */
372  t = zero;
373  SET_HIGH_WORD(t, (n & ~(0x000fffff >> k)));
374  n = ((n & 0x000fffff) | 0x00100000) >> (20 - k);
375  if (j < 0)
376  n = -n;
377  p_h -= t;
378  }
379  t = p_l + p_h;
380  SET_LOW_WORD(t, 0);
381  u = t * lg2_h;
382  v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
383  z = u + v;
384  w = v - (z - u);
385  t = z * z;
386  t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
387  r = (z * t1) / (t1 - two) - (w + z * w);
388  z = one - (r - z);
389  GET_HIGH_WORD(j, z);
390  j += (n << 20);
391  if ((j >> 20) <= 0)
392  z = scalbn(z, n); /* subnormal output */
393  else
394  {
395  uint32_t hz;
396  GET_HIGH_WORD(hz, z);
397  SET_HIGH_WORD(z, hz + (n << 20));
398  }
399  return s * z;
400 }
cp_l
static const double cp_l
Definition: pow.c:87
L6
static const double L6
Definition: pow.c:75
GET_HIGH_WORD
#define GET_HIGH_WORD(i, d)
Get the more significant 32 bit int from a double.
Definition: prim.h:167
lg2_h
static const double lg2_h
Definition: pow.c:82
dp_h
static const double dp_h[]
Definition: pow.c:60
ivln2_h
static const double ivln2_h
Definition: pow.c:89
two53
static const double two53
Definition: pow.c:67
P2
static const double P2
Definition: pow.c:77
P1
static const double P1
Definition: pow.c:76
huge
static const double huge
Definition: pow.c:68
one
static const double one
Definition: pow.c:67
EXTRACT_WORDS
#define EXTRACT_WORDS(ix0, ix1, d)
Get two 32 bit ints from a double.
Definition: prim.h:156
ivln2_l
static const double ivln2_l
Definition: pow.c:90
P3
static const double P3
Definition: pow.c:78
ivln2
static const double ivln2
Definition: pow.c:88
lg2
static const double lg2
Definition: pow.c:81
L3
static const double L3
Definition: pow.c:72
zero
static const double zero
Definition: pow.c:66
scalbn
double scalbn(double x, int n)
ovt
static const double ovt
Definition: pow.c:84
dp_l
static const double dp_l[]
Definition: pow.c:63
P4
static const double P4
Definition: pow.c:79
L1
static const double L1
Definition: pow.c:70
L5
static const double L5
Definition: pow.c:74
cp_h
static const double cp_h
Definition: pow.c:86
lg2_l
static const double lg2_l
Definition: pow.c:83
tiny
static const double tiny
Definition: pow.c:68
two
static const double two
Definition: pow.c:67
SET_LOW_WORD
#define SET_LOW_WORD(d, v)
Set the less significant 32 bits of a double from an int.
Definition: prim.h:209
L4
static const double L4
Definition: pow.c:73
L2
static const double L2
Definition: pow.c:71
cp
static const double cp
Definition: pow.c:85
sqrt
double sqrt(double x)
Square root function.
Definition: sqrt.c:119
SET_HIGH_WORD
#define SET_HIGH_WORD(d, v)
Set the more significant 32 bits of a double from an int.
Definition: prim.h:198
fabs
double fabs(double x)
Returns the absolute value of x.
Definition: fabs.c:51
bp
static const double bp[]
Definition: pow.c:57
P5
static const double P5
Definition: pow.c:80
Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ bp

const double bp[]
static
Initial value:
= {
1.0, 1.5,
}

Definition at line 57 of file pow.c.

Referenced by pow().

◆ cp

const double cp = 9.61796693925975554329e-01
static

Definition at line 85 of file pow.c.

Referenced by pow().

◆ cp_h

const double cp_h = 9.61796700954437255859e-01
static

Definition at line 86 of file pow.c.

Referenced by pow().

◆ cp_l

const double cp_l = -7.02846165095275826516e-09
static

Definition at line 87 of file pow.c.

Referenced by pow().

◆ dp_h

const double dp_h[]
static
Initial value:
= {
0.0, 5.84962487220764160156e-01,
}

Definition at line 60 of file pow.c.

Referenced by pow().

◆ dp_l

const double dp_l[]
static
Initial value:
= {
0.0, 1.35003920212974897128e-08,
}

Definition at line 63 of file pow.c.

Referenced by pow().

◆ huge

const double huge = 1.0e300
static

Definition at line 68 of file pow.c.

Referenced by pow().

◆ ivln2

const double ivln2 = 1.44269504088896338700e+00
static

Definition at line 88 of file pow.c.

Referenced by pow().

◆ ivln2_h

const double ivln2_h = 1.44269502162933349609e+00
static

Definition at line 89 of file pow.c.

Referenced by pow().

◆ ivln2_l

const double ivln2_l = 1.92596299112661746887e-08
static

Definition at line 90 of file pow.c.

Referenced by pow().

◆ L1

const double L1 = 5.99999999999994648725e-01
static

Definition at line 70 of file pow.c.

Referenced by pow().

◆ L2

const double L2 = 4.28571428578550184252e-01
static

Definition at line 71 of file pow.c.

Referenced by pow().

◆ L3

const double L3 = 3.33333329818377432918e-01
static

Definition at line 72 of file pow.c.

Referenced by pow().

◆ L4

const double L4 = 2.72728123808534006489e-01
static

Definition at line 73 of file pow.c.

Referenced by pow().

◆ L5

const double L5 = 2.30660745775561754067e-01
static

Definition at line 74 of file pow.c.

Referenced by pow().

◆ L6

const double L6 = 2.06975017800338417784e-01
static

Definition at line 75 of file pow.c.

Referenced by pow().

◆ lg2

const double lg2 = 6.93147180559945286227e-01
static

Definition at line 81 of file pow.c.

Referenced by pow().

◆ lg2_h

const double lg2_h = 6.93147182464599609375e-01
static

Definition at line 82 of file pow.c.

Referenced by pow().

◆ lg2_l

const double lg2_l = -1.90465429995776804525e-09
static

Definition at line 83 of file pow.c.

Referenced by pow().

◆ one

const double one = 1.0
static

Definition at line 67 of file pow.c.

Referenced by pow().

◆ ovt

const double ovt = 8.0085662595372944372e-0017
static

Definition at line 84 of file pow.c.

Referenced by pow().

◆ P1

const double P1 = 1.66666666666666019037e-01
static

Definition at line 76 of file pow.c.

Referenced by pow().

◆ P2

const double P2 = -2.77777777770155933842e-03
static

Definition at line 77 of file pow.c.

Referenced by pow().

◆ P3

const double P3 = 6.61375632143793436117e-05
static

Definition at line 78 of file pow.c.

Referenced by pow().

◆ P4

const double P4 = -1.65339022054652515390e-06
static

Definition at line 79 of file pow.c.

Referenced by pow().

◆ P5

const double P5 = 4.13813679705723846039e-08
static

Definition at line 80 of file pow.c.

Referenced by pow().

◆ tiny

const double tiny = 1.0e-300
static

Definition at line 68 of file pow.c.

Referenced by pow().

◆ two

const double two = 2.0
static

Definition at line 67 of file pow.c.

Referenced by pow().

◆ two53

const double two53 = 9007199254740992.0
static

Definition at line 67 of file pow.c.

Referenced by pow().

◆ zero

const double zero = 0.0
static

Definition at line 66 of file pow.c.

Referenced by pow().