128 lines
4.0 KiB
C
128 lines
4.0 KiB
C
/* Copyright (C) 2008-2021 Free Software Foundation, Inc.
|
|
Contributor: Joern Rennecke <joern.rennecke@embecosm.com>
|
|
on behalf of Synopsys Inc.
|
|
|
|
This file is part of GCC.
|
|
|
|
GCC is free software; you can redistribute it and/or modify it under
|
|
the terms of the GNU General Public License as published by the Free
|
|
Software Foundation; either version 3, or (at your option) any later
|
|
version.
|
|
|
|
GCC 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 General Public License
|
|
for more details.
|
|
|
|
Under Section 7 of GPL version 3, you are granted additional
|
|
permissions described in the GCC Runtime Library Exception, version
|
|
3.1, as published by the Free Software Foundation.
|
|
|
|
You should have received a copy of the GNU General Public License and
|
|
a copy of the GCC Runtime Library Exception along with this program;
|
|
see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
|
|
<http://www.gnu.org/licenses/>. */
|
|
|
|
/* We use a polynom similar to a Tchebycheff polynom to get an initial
|
|
seed, and then use a newton-raphson iteration step to get an
|
|
approximate result
|
|
If this result can't be rounded to the exact result with confidence, we
|
|
round to the value between the two closest representable values, and
|
|
test if the correctly rounded value is above or below this value.
|
|
|
|
Because of the Newton-raphson iteration step, an error in the seed at X
|
|
is amplified by X. Therefore, we don't want a Tchebycheff polynom
|
|
or a polynom that is close to optimal according to the maximum norm
|
|
on the errro of the seed value; we want one that is close to optimal
|
|
according to the maximum norm on the error of the result, i.e. we
|
|
want the maxima of the polynom to increase linearily.
|
|
Given an interval [X0,X2) over which to approximate,
|
|
with X1 := (X0+X2)/2, D := X1-X0, F := 1/D, and S := D/X1 we have,
|
|
like for Tchebycheff polynoms:
|
|
P(0) := 1
|
|
but then we have:
|
|
P(1) := X + S*D
|
|
P(2) := 2 * X^2 + S*D * X - D^2
|
|
Then again:
|
|
P(n+1) := 2 * X * P(n) - D^2 * P (n-1)
|
|
*/
|
|
|
|
int
|
|
main (void)
|
|
{
|
|
long double T[5]; /* Taylor polynom */
|
|
long double P[5][5];
|
|
int i, j;
|
|
long double X0, X1, X2, S;
|
|
long double inc = 1./64;
|
|
long double D = inc*0.5;
|
|
long i0, i1, i2;
|
|
|
|
memset (P, 0, sizeof (P));
|
|
P[0][0] = 1.;
|
|
for (i = 1; i < 5; i++)
|
|
P[i][i] = 1 << i-1;
|
|
P[2][0] = -D*D;
|
|
for (X0 = 1.; X0 < 2.; X0 += inc)
|
|
{
|
|
X1 = X0 + inc * 0.5;
|
|
X2 = X1 + inc;
|
|
S = D / X1;
|
|
T[0] = 1./X1;
|
|
for (i = 1; i < 5; i++)
|
|
T[i] = T[i-1] * -T[0];
|
|
#if 0
|
|
printf ("T %1.8f %f %f %f %f\n", (double)T[0], (double)T[1], (double)T[2],
|
|
(double)T[3], (double)T[4]);
|
|
#endif
|
|
P[1][0] = S*D;
|
|
P[2][1] = S*D;
|
|
for (i = 3; i < 5; i++)
|
|
{
|
|
P[i][0] = -D*D*P[i-2][0];
|
|
for (j = 1; j < i; j++)
|
|
P[i][j] = 2*P[i-1][j-1]-D*D*P[i-2][j];
|
|
}
|
|
#if 0
|
|
printf ("P3 %1.8f %f %f %f %f\n", (double)P[3][0], (double)P[3][1], (double)P[3][2],
|
|
(double)P[3][3], (double)P[3][4]);
|
|
printf ("P4 %1.8f %f %f %f %f\n", (double)P[4][0], (double)P[4][1], (double)P[4][2],
|
|
(double)P[4][3], (double)P[4][4]);
|
|
#endif
|
|
for (i = 4; i > 1; i--)
|
|
{
|
|
long double a = T[i]/P[i][i];
|
|
|
|
for (j = 0; j < i; j++)
|
|
T[j] -= a * P[i][j];
|
|
}
|
|
#if 0
|
|
printf ("A %1.8f %f %f\n", (double)T[0], (double)T[1], (double)T[2]);
|
|
#endif
|
|
#if 0
|
|
i2 = T[2]*512;
|
|
long double a = (T[2]-i/512.)/P[2][2];
|
|
for (j = 0; j < 2; j++)
|
|
T[j] -= a * P[2][j];
|
|
#else
|
|
i2 = 0;
|
|
#endif
|
|
for (i = 0, i0 = 0; i < 4; i++)
|
|
{
|
|
long double T0, Ti1;
|
|
|
|
i1 = T[1]*8192. + i0 / (long double)(1 << 19) - 0.5;
|
|
i1 = - (-i1 & 0x1fff);
|
|
Ti1 = ((unsigned)(-i1 << 19) | i0) /-(long double)(1LL<<32LL);
|
|
T0 = T[0] - (T[1]-Ti1)/P[1][1] * P[1][0] - (X1 - 1) * Ti1;
|
|
i0 = T0 * 512 * 1024 + 0.5;
|
|
i0 &= 0x7ffff;
|
|
}
|
|
#if 0
|
|
printf ("A %1.8f %f %f\n", (double)T[0], (double)T[1], (double)T[2]);
|
|
#endif
|
|
printf ("\t.long 0x%x\n", (-i1 << 19) | i0);
|
|
}
|
|
return 0;
|
|
}
|