164 lines
4.9 KiB
C
164 lines
4.9 KiB
C
/*
|
|
* IBM Accurate Mathematical Library
|
|
* written by International Business Machines Corp.
|
|
* Copyright (C) 2001-2017 Free Software Foundation, Inc.
|
|
*
|
|
* This program is free software; you can redistribute it and/or modify
|
|
* it under the terms of the GNU Lesser General Public License as published by
|
|
* the Free Software Foundation; either version 2.1 of the License, or
|
|
* (at your option) any later version.
|
|
*
|
|
* This program is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
* GNU Lesser General Public License for more details.
|
|
*
|
|
* You should have received a copy of the GNU Lesser General Public License
|
|
* along with this program; if not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
/*************************************************************************/
|
|
/* MODULE_NAME:mpexp.c */
|
|
/* */
|
|
/* FUNCTIONS: mpexp */
|
|
/* */
|
|
/* FILES NEEDED: mpa.h endian.h mpexp.h */
|
|
/* mpa.c */
|
|
/* */
|
|
/* Multi-Precision exponential function subroutine */
|
|
/* ( for p >= 4, 2**(-55) <= abs(x) <= 1024 ). */
|
|
/*************************************************************************/
|
|
|
|
#include "endian.h"
|
|
#include "mpa.h"
|
|
#include <assert.h>
|
|
|
|
#ifndef SECTION
|
|
# define SECTION
|
|
#endif
|
|
|
|
/* Multi-Precision exponential function subroutine (for p >= 4,
|
|
2**(-55) <= abs(x) <= 1024). */
|
|
void
|
|
SECTION
|
|
__mpexp (mp_no *x, mp_no *y, int p)
|
|
{
|
|
int i, j, k, m, m1, m2, n;
|
|
mantissa_t b;
|
|
static const int np[33] =
|
|
{
|
|
0, 0, 0, 0, 3, 3, 4, 4, 5, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6,
|
|
6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 8
|
|
};
|
|
|
|
static const int m1p[33] =
|
|
{
|
|
0, 0, 0, 0,
|
|
17, 23, 23, 28,
|
|
27, 38, 42, 39,
|
|
43, 47, 43, 47,
|
|
50, 54, 57, 60,
|
|
64, 67, 71, 74,
|
|
68, 71, 74, 77,
|
|
70, 73, 76, 78,
|
|
81
|
|
};
|
|
static const int m1np[7][18] =
|
|
{
|
|
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
|
{0, 0, 0, 0, 36, 48, 60, 72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
|
|
{0, 0, 0, 0, 24, 32, 40, 48, 56, 64, 72, 0, 0, 0, 0, 0, 0, 0},
|
|
{0, 0, 0, 0, 17, 23, 29, 35, 41, 47, 53, 59, 65, 0, 0, 0, 0, 0},
|
|
{0, 0, 0, 0, 0, 0, 23, 28, 33, 38, 42, 47, 52, 57, 62, 66, 0, 0},
|
|
{0, 0, 0, 0, 0, 0, 0, 0, 27, 0, 0, 39, 43, 47, 51, 55, 59, 63},
|
|
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 43, 47, 50, 54}
|
|
};
|
|
mp_no mps, mpk, mpt1, mpt2;
|
|
|
|
/* Choose m,n and compute a=2**(-m). */
|
|
n = np[p];
|
|
m1 = m1p[p];
|
|
b = X[1];
|
|
m2 = 24 * EX;
|
|
for (; b < HALFRAD; m2--)
|
|
b *= 2;
|
|
if (b == HALFRAD)
|
|
{
|
|
for (i = 2; i <= p; i++)
|
|
{
|
|
if (X[i] != 0)
|
|
break;
|
|
}
|
|
if (i == p + 1)
|
|
m2--;
|
|
}
|
|
|
|
m = m1 + m2;
|
|
if (__glibc_unlikely (m <= 0))
|
|
{
|
|
/* The m1np array which is used to determine if we can reduce the
|
|
polynomial expansion iterations, has only 18 elements. Besides,
|
|
numbers smaller than those required by p >= 18 should not come here
|
|
at all since the fast phase of exp returns 1.0 for anything less
|
|
than 2^-55. */
|
|
assert (p < 18);
|
|
m = 0;
|
|
for (i = n - 1; i > 0; i--, n--)
|
|
if (m1np[i][p] + m2 > 0)
|
|
break;
|
|
}
|
|
|
|
/* Compute s=x*2**(-m). Put result in mps. This is the range-reduced input
|
|
that we will use to compute e^s. For the final result, simply raise it
|
|
to 2^m. */
|
|
__pow_mp (-m, &mpt1, p);
|
|
__mul (x, &mpt1, &mps, p);
|
|
|
|
/* Compute the Taylor series for e^s:
|
|
|
|
1 + x/1! + x^2/2! + x^3/3! ...
|
|
|
|
for N iterations. We compute this as:
|
|
|
|
e^x = 1 + (x * n!/1! + x^2 * n!/2! + x^3 * n!/3!) / n!
|
|
= 1 + (x * (n!/1! + x * (n!/2! + x * (n!/3! + x ...)))) / n!
|
|
|
|
k! is computed on the fly as KF and at the end of the polynomial loop, KF
|
|
is n!, which can be used directly. */
|
|
__cpy (&mps, &mpt2, p);
|
|
|
|
double kf = 1.0;
|
|
|
|
/* Evaluate the rest. The result will be in mpt2. */
|
|
for (k = n - 1; k > 0; k--)
|
|
{
|
|
/* n! / k! = n * (n - 1) ... * (n - k + 1) */
|
|
kf *= k + 1;
|
|
|
|
__dbl_mp (kf, &mpk, p);
|
|
__add (&mpt2, &mpk, &mpt1, p);
|
|
__mul (&mps, &mpt1, &mpt2, p);
|
|
}
|
|
__dbl_mp (kf, &mpk, p);
|
|
__dvd (&mpt2, &mpk, &mpt1, p);
|
|
__add (&__mpone, &mpt1, &mpt2, p);
|
|
|
|
/* Raise polynomial value to the power of 2**m. Put result in y. */
|
|
for (k = 0, j = 0; k < m;)
|
|
{
|
|
__sqr (&mpt2, &mpt1, p);
|
|
k++;
|
|
if (k == m)
|
|
{
|
|
j = 1;
|
|
break;
|
|
}
|
|
__sqr (&mpt1, &mpt2, p);
|
|
k++;
|
|
}
|
|
if (j)
|
|
__cpy (&mpt1, y, p);
|
|
else
|
|
__cpy (&mpt2, y, p);
|
|
return;
|
|
}
|