# Category : Assembly Language Source Code

Archive : PICALC.ZIP

Filename : SPIGOT.PRG

* 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

Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/