User talk:Lunch/Wilkinson's polynomial

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