Code to calculate the largest root of Wilkinson's polynomial using just IEEE single precision:
/* for printf */
#include <stdio.h>
/* for fabsf */
#include <math.h>
int main() {
/* xc is the current iterate
xp is the previous iterate
dx is their difference */
float xc=21.0, xp=20.0, dx=1.0;
/* fcProd is Wilkinson's polynomial at xc
fcMon is the perturbation at xc
fc = fcProd-fcMon
fpProd is Wilkinson's polynomial at xp (inited exactly)
fpMon is the perturbation at xp (inited almost exactly)
NB: 6.25e17 = Plus@@(2^{59,55,53,51,50,46,45,44,41,37,36,30,27,23,20,19,18,17,15})
That is, a 44 bit mantissa is needed (remember the leading
"1" is assumed). IEEE single precision doesn't have the
bits.
xcSq is a temp for holding xc*xc in computing x^19 */
float fcProd, fcMon, fc, fpProd = 0.0, fpMon = 6.25e17, xcSq;
/* zi is the ith root of Wilkinson's polynomial */
float zi;
/* dfProd = fcProd-fpProd
dfMon = fcMon-fpMon
df = dfProd-dfMon */
float dfProd, dfMon, df;
printf("xc = %f dx = %g ",xc,dx);
/* keep going until the relative change is less than two binary ULPs */
while (fabsf(dx/xc)>0x1p-22f) {
/* calculate current function values */
/* fcProd
unrolling this loop is OK. if we assume xc>20, then
this product is calculated by multiplying the largest
factor first down to the smallest. */
for (zi=1.0, fcProd=1.0; zi<20.5; zi+=1.0)
fcProd *= xc-zi;
/* fcMon
calculate x^19 via binary expansion of 19=16+2+1
NB: 0x1p-23 = 1*2^-23 in hexadecimal (for the mantissa;
the exponent is in decimal; go figure) */
xcSq = xc*xc;
fcMon = xcSq*xcSq;
fcMon *= fcMon;
fcMon *= fcMon;
fcMon *= xcSq;
fcMon *= xc;
fcMon *= 0x1p-23f;
/* calculate denominator */
dfProd = fcProd-fpProd;
dfMon = fcMon-fpMon;
df = dfProd-dfMon;
/* calculate fc */
fc = fcProd-fcMon;
printf("fc = %g\n",fc);
/* calculate increment and update previous iterate storage */
dx *= -fc;
dx /= df;
xp = xc;
fpProd = fcProd;
fpMon = fcMon;
/* update current iterate and report it */
xc += dx;
printf("xc = %f dx = %g ",xc,dx);
}
/* all the assignments and stores of intermediate values
are needed to force stores of single precision values
with -ffloat-store. (otherwise many compilers and FPUs
will substitute doubles) if optimization for speed is
desired, most compilers will perform common subexpression
elimination when asked. */
}
The command line used to compile the above:
gcc -o find-largest-root -ffloat-store -O0 find-largest-root.c
And the text the program outputs:
xc = 21.000000 dx = 1 fc = 8.53558e+17 xc = 20.422709 dx = -0.577291 fc = -7.25372e+17 xc = 20.687920 dx = 0.265212 fc = -4.67218e+17 xc = 21.167912 dx = 0.479991 fc = 2.52229e+18 xc = 20.762936 dx = -0.404975 fc = -2.87837e+17 xc = 20.804417 dx = 0.041481 fc = -1.58599e+17 xc = 20.855322 dx = 0.0509049 fc = 3.48639e+16 xc = 20.846148 dx = -0.00917355 fc = -3.09265e+15 xc = 20.846895 dx = 0.00074745 fc = -5.30514e+13 xc = 20.846909 dx = 1.30456e-05 fc = 1.78671e+12 xc = 20.846909 dx = -4.25044e-07