gnu / gcc / 4a0fed0c0c7241562eaa5f1a4c916b689429ad86 / . / libgcc / config / arc / ieee-754 / divtab-arc-sf.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; | |

} |