User talk:Lunch/Wilkinson's polynomial

Code to calculate the largest root of Wilkinson's polynomial using just IEEE single precision:

/* for printf */ /* for fabsf */ 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. */ }
 * 1) include 
 * 1) include 

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