;  File: PICALC.ASM    System: PI Calc       Version: 1.00  Date: 10-01-94  ;

;-----------------------------------------------------------------------------
; PICALC.COM computes and displays number PI using spigot algorithm (due to S.
; Rabinowitz).  Accuracy here is determined by size of spigot array: 32786
; words => 32768*log(2) = 9864 digits => 2466 iterations.  Output may be
; redirected to file.  For radix 10000 and array size here, each output word
; is luckily in range [0, radix-1].  For smaller radix or larger (dword)
; array, output word (dword) will sometimes match radix exactly, implying
; zeroes and bump by one of previously output digits.
;
; Program takes 2 minutes on a 486/66 MHz machine.  See PI386.COM for a faster
; version algorithm, using double word array elements and capable of finding
; several million digits of pi.  PICALC.COM is 85 bytes.  Goal here was
; smallest possible COM file (with no command-line arguments) capable of
; calculating and displaying nearly 10K decimal digits of pi.
;
; C. Hessel/ER Support/DOE/Germantown MD
;-----------------------------------------------------------------------------

Code_Seg  SEGMENT                  ; Assume 128K free and AH zero from DOS
          ASSUME cs:Code_Seg

          ORG   100h
Begin:    mov   dx,cs              ; Get segment 64K forward, insuring
          add   dh,10h             ;   past stack
          mov   es,dx              ; Save as buffer segment for 64K spigot
          std                      ; Work backwards for duration
          xor   di,di              ; Buffer start
          mov   al,2               ; Store word 2's to buffer, with wrap
          mov   cx,8001h           ;   leaving DI pointing to last word and
          rep   stosw              ;   zeroing CX for main loop entry
          mov   si,10000           ; Fix at largest 10-power fitting word
          mov   bp,2466            ; Maximum iteration value with BX slop
          jmp   SHORT iterloop     ; Best BX is 20000, but any ok--also could
                                   ;   zero BX and make last array word 4
bufloop:  inc   di                 ; 2*k+1
          div   di                 ; temp <- temp/(2*k+1)
          mov   bx,dx              ; a[k] <- remainder
          shr   di,1               ; k
          mul   di                 ; temp <- k*temp
          add   di,di              ; 2*k
          xchg  bx,ax              ; q <- temp and ready a[k] for store
          mov   cx,dx
          stosw                    ; Store a[k], decrementing k
iterloop: mov   ax,es:[di]         ; Fetch next a[k]
          mul   si                 ; temp <- a[k]*10000
          add   ax,bx              ; temp <- temp+q
          adc   dx,cx
          test  di,di              ; Test if k = 0
          jne   bufloop            ; No?  Loop

          div   si                 ; Output word <- temp/10000
          xchg  ax,dx              ; a[0] <- remainder
          stosw                    ; Store a[0], resetting k to buffer end
          xchg  ax,dx              ; Output word back to AX (0-9999)
          mov   cl,4               ; CX was zero on inner loop exit
          mov   bx,10              ; Divisor 10 to extract digits
pushloop: cwd                      ; Zeroes DX since AX < 32768
          div   bx
          add   dl,'0'             ; Make remainder ASCII digit
          push  dx                 ; Save, since must reverse order
          loop  pushloop

          mov   cl,4
poploop:  pop   dx                 ; Pop digits in reverse order
          mov   ah,2
          int   21h                ; DOS display character DL
          loop  poploop

          dec   bp                 ; Should zero CX:BX, but 10 ok
          jne   iterloop           ; Not about 10K digits yet?  Loop

          ret                      ; Exit via int 20h at PSP start

Code_Seg  ENDS
          END  Begin