********** * Foxpro ã calculation by spigot algorithm. Purpose here is to show error * bound and to demo infrequent cases where quotient exactly equals radix. * Foxpro LOG() function is natual logarithm, so notes will use same * convention. * * One bit of accuracy is gained for each additional array element. E.g., * an increment of 20 in length of spigot array gains about 6 digits * accuracy since 20*log(2)/log(10) = 20*.30103 = 6.0+. In addition, each * time the array size is doubled, an extra half bit of accuracy is gained. * * Setting last array element to 4 initially (instead of 2) improves the * second order error term to 1.5 bits when array size is doubled, but * that is not demonstrated, since not used in PI386.COM. Results are * similar though (with ûã/2 and ûã*15/16 convergence limits). * * Integer size for array elements must accommodate values up to twice the * array length less one. The dividend and interim quotient values must be * double integer length. The maximum radix used is 100,000,000 since * Foxpro handles integers to 16 digit accuracy. PI386.COM uses 1 billion. * * On output, if quotient is exactly radix, asterisks are printed to mean * all zeroes for that block, but with implied bump by one of digits to left * (already output). This occurs less than one time per every radix output * quotients (assuming ã digits are random and noting that all zeroes * occurs sometimes without bump), so is most noticeable for small radix. * * All observations here are empirical. My starting point was the following * demo C program (by Jim Roth at DEC) downloaded from the Internet---CRH. * * /* 1000 digits of PI */ * /* 'Spigot' algorithm originally due to Stanly Rabinowitz */ * * #include * * main() * { * int d = 4, r = 10000, n = 251, m = 3.322*n*d; * int i, j, k, q; * static int a[3340]; * * for (i = 0; i <= m; i++) a[i] = 2; * a[m] = 4; * * for (i = 1; i <= n; i++) { * q = 0; * for (k = m; k > 0; k--) { * a[k] = a[k]*r+q; * q = a[k]/(2*k+1); * a[k] -= (2*k+1)*q; * q *= k; * } * a[0] = a[0]*r+q; * q = a[0]/r; * a[0] -= q*r; * printf("%04d%s",q, i & 7 ? " " : "\n"); * } * } ********** PRIV nMaxPwr, aErr, cBasePi, cCalcPi, n,; nSqrtPi, nRatio, nDiff1, nProd, nDiff2, nAdjust CLEAR SET TALK OFF nMaxPwr = 11 && Can go higher with FOXPROX.EXE DIMEN aErr[nMaxPwr+1] ? "Compute base ã value for error comparisons that follow..." ? cBasePi = PICALC(2**nMaxPwr+48, 8) ? WAIT "Double array size repeatedly to check error estimate--press key" ? ? " size calculated value of ã with about 12 slack decimals" ? " ---- --------------------------------------------------" FOR n = 0 TO nMaxPwr cCalcPi = PICALC(2**n, 8, .T.) aErr[n+1] = GETERR(n, cCalcPi, cBasePi) NEXT ? WAIT "Ratio of estimate 0.5**(size+n/2) to relative error (ã-picalc)/ã"+; "--press key" ? ? " array error less times less adjusted" ? " n size ratio ûã/2 size ûã*7/16 ratio" ? " -- ---- -------- -------- -------- -------- --------" nSqrtPi = SQRT(2*ACOS(0)) FOR n = 0 TO nMaxPwr nRatio = aErr[n+1] nDiff1 = nRatio-nSqrtPi/2 nProd = nDiff1*2**n nDiff2 = nProd-nSqrtPi*7/16 nAdjust = nRatio*(1-7/8/2**n) ? STR(n, 3), STR(2**n, 5), STR(nRatio, 9, 6), STR(nDiff1, 9, 6),; STR(nProd, 9, 6), STR(nDiff2, 9, 6), STR(nAdjust, 9, 6) NEXT ? WAIT "Slow radix 10, showing some asterisk bumps--press key" ? DO PICALC WITH 1024, 1, .T. ? WAIT "Slow radix 100, showing a double asterisk bump--press key" ? DO PICALC WITH 2048, 2, .T. RETURN ********** * Returns ratio of error estimate 0.5**(2**n+0.5*n) to observed relative * error (ã-picalc)/ã to demo the second order error estimate. Ratio * converges downward to about 0.88623-. The limit seems to be ûã/2. * Third order convergence is also displayed in the output table. * * Absolute error is 0.5**(2**n+0.5*n) * ã/ratio, which is equivalent to * (2**n+0.5*n)*log(2)/log(10) - log(ã/ratio)/log(10) decimal digits. * Since ratio converges downward, using 0.55 > log(2*ûã)/log(10) for * last term gives strict error bound. That is, absolute error ã-picalc * is less than 0.1**((k + 0.5*lg(k))*log(2)/log(10) - 0.55), where k is * array size and lg() denotes logarithm base 2. * * The convergences above were refined by trial and error, beginning with * first order try of 0.5**(2**n), hinted at by factor 3.322 in C program. * The ûã/2 limit was guessed from a scan of mathematical constants in an * appendix to Knuth's "Art of Computer Programming", Vol I. ********** FUNC ERRRATIO PARA n, nDcmls, nFrac && Log base 2 of array size, accurate PRIV nPi && decimals, error fraction nPi = 2*ACOS(0) && Just need a few digits accuracy for factor RETURN EXP(LOG(10)*nDcmls-LOG(2)*(2**n+0.5*n))*nPi/nFrac ********** * Compare base string with shorter calc string to find tail error in calc, * then return error ratio. Calculated approximation is just under ã. ********** FUNC GETERR PARA nTwoPwr, cCalc, cBase PRIV n, nErrLen, cOk, cErr FOR n = LEN(cCalc) TO 0 STEP -1 IF LEFT(cBase, n)=LEFT(cCalc, n) EXIT && Exit guaranteed at n=0 ENDIF NEXT nErrLen = LEN(cCalc)-n && Tail length not in agreement cOk = "."+SUBS(cBase, n+1, nErrLen) cErr = "."+SUBS(cCalc, n+1, nErrLen) RETURN ERRRATIO(nTwoPwr, n-1, VAL(cOk)-VAL(cErr)) ********** * Calculate and display ã approximation, given array/radix sizes. ********** FUNC PICALC PARA nSize, nRdxDgts, lScale && Array size, 1-8 radix digits, scale flag PRIV a, nSlack, n, nIter, nLast, cRet, nDgts, cDgts DIMEN a[nSize] && Foxpro arrays are 1-based STORE 2 TO a && Initialize array to all 2's nSlack = 12 && Slack inaccurate decimals n = (LOG(2)*nSize+0.5*LOG(nSize))/LOG(10)-0.55 n = ROUND(n, 0)+nSlack && Decimal digits to output, with slack nIter = CEILING(n/nRdxDgts) && Iterations (aside from ordinal pass) n = n%nRdxDgts nLast = IIF(n=0, nRdxDgts, n) && Digits to output on last pass cRet = "" && Concatenate digits for return FOR n = 0 TO nIter && Start from zero for ordinal pass * Each array pass returns nRdxDgts digits. If scaling flagged, adjust * second argument downward linearly after a few slop iterations. This * is not quite same as scaling in PI386.COM (which steps down less * frequently, but by larger chunks) but demonstrates the principle. * With enough slop, the error due to scaling is small compared to * the normal algorithm error. The slop here is intentionally large * so that the slack inaccurate digits are same as with no scaling, * which in turn allows accurate error comparisons. PI386.COM walks * a finer line. Scaling cuts run time roughly in half by ignoring * the insignificant portion of array tail. nDgts = n-INT(2*nSlack/nRdxDgts)-1 nDgts = IIF(lScale AND nDgts>0, INT(nSize*nDgts/(nIter+1)), 0) nDgts = GETDGTS(10**nRdxDgts, nSize-nDgts) * Convert result to character digits, with special handling for first * and last iterations, then concatenate to return string. Foxpro * STR() function returns asterisks if format length is exceeded. This * is used to flag output that exactly matches radix. A large radix is * used for error testing so that asterisks do not appear. cDgts = STR(nDgts, nRdxDgts) && Asterisks possible... cDgts = STRTRAN(cDgts, " ", "0") IF n=0 && First pass? cDgts = RIGHT(cDgts, 1) && Grab ordinal digit only ENDIF IF n=nIter && Last pass? cDgts = LEFT(cDgts, nLast) && Shorten output as needed ENDIF cRet = cRet+cDgts && Extend return value DO DISPDGTS WITH cDgts, (n-1)*nRdxDgts, IIF(n=0, nSize, 0) NEXT RETURN cRet ********** * Display next chunk of digits. ********** PROC DISPDGTS PARA cDgts, nCnt, nSize && Digit string, decimals previously output, PRIV n && array size if first pass (else zero) IF nSize > 0 && First pass? Output size and ordinal ? STR(nSize, 5)+" "+cDgts+"." ELSE FOR n = 1 TO LEN(cDgts) ?? SUBS(cDgts, n, 1) && Output next digit nCnt = nCnt+1 IF nCnt%10=0 IF nCnt%60=0 ? SPACE(8) && New line and indent every 60 ELSE ?? " " && Space delimiter every 10 ENDIF ENDIF NEXT ENDIF RETURN ********** * Heart of spigot algorithm. Each call extracts a few more digits of ã * from the array spigot (up to the accuracy limit). Array is updated. * Note that array is 1-based rather than 0-based (as in C demo and in * PI386.COM). ********** FUNC GETDGTS PARA nRadix, nSize && Radix, array size (possibly reduced) PRIV n, nQuot, nDvnd, nDvsr nQuot = 0 FOR n = nSize TO 2 STEP -1 && Work backwards nDvsr = n+n-1 && Divisor--bounds array element size nDvdn = a[n]*nRadix+nQuot && Dividend--double size a[n] = nDvdn%nDvsr && Remainder--single size nQuot = INT(nDvdn/nDvsr) && Quotient--single size nQuot = nQuot*(n-1) && Interim double size, but down to NEXT && single size on loop exit nDvsr = nRadix nDvdn = a[1]*nRadix+nQuot a[1] = nDvdn%nDvsr && Remainder nQuot = INT(nDvdn/nDvsr) && Output quotient--can be radix rarely RETURN nQuot * End program