Category : Assembly Language Source Code
Archive   : PICALC.ZIP
Filename : SPIGOT.PRG

 
Output of file : SPIGOT.PRG contained in archive : PICALC.ZIP
**********
* 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


  3 Responses to “Category : Assembly Language Source Code
Archive   : PICALC.ZIP
Filename : SPIGOT.PRG

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

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

  3. 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/