LBL "DPI"
; DEEP PI
; calculates the next 5 digits
; of PI at digit X+1.
; ie X=1 gives 141592, X >=1
; press R/S for the next 5 and so on.

STO 22
LBL 25
XEQ 24
1e5
*
INT
STOP
5
ST+ 22
RCL 22
GTO 25
LBL 24
                ; compute digits at X
                ; this is the computation subroutine
STO 17          ; n
20
+
.8847
*
INT
STO 16          ; N
3
*
LN              ; log(3N)
2
LN              ; log(2)
/
INT             ; vmax = ln(3n)/ln(2)
RCL 16
+
RCL 17
-               
STO 12          ; vmax = vmax + N - n
0
STO 13          ; sum=0
RCL 12          ; vmax
X<=0?
GTO 09          ; skip case 2 (l1)

; unrolled case for a=2

2
RCL 12
Y^X             ; 2^vmax
STO 8           ; av=2^vmax
2
RCL 12
1
+
2 
/ 
INT
Y^X             ; 2^(vmax/2)
STO 9           ; B=2^(vmax/2)
0
STO 11          ; s=0
1
STO 10          ; den=1
STO 7           ; num=1
STO 15          ; k=1
RCL 17          
CHS
STO 6           ; v=-n
LBL 10          ; for(k=1; k<=N; ++k)
1
CHS
STO 1           ; dec v
RCL 15          ; k
STO 0           ; t
XEQ 11          ; reduce mod 2
RCL 0           ; t
RCL 15          ; k
2
*
1
-
*               ; t*(2k-1)
RCL 7           ; num
XEQ 18          ; num*t mod av
STO 7           ; num
1
STO 1           ; inc v
RCL 15          ; k
9
*
3
-    
STO 0           ; t=3*(3k-1)
XEQ 11          ; reduce t
RCL 0           ; t
RCL 10          ; den
XEQ 18
STO 10          ; den = den * t mod av
1
STO 1           ; inc v
RCL 15          ; k
3
*
2
-               ; 3k-2
STO 0
XEQ 11
RCL 0
RCL 10
XEQ 18
STO 10          ; den = den *t mov av
RCL 6           ; v
X<=0?
GTO 12
                ; v > 0
RCL 8           ; av
RCL 10          ; den
XEQ 19          ; den^-1 mod av
RCL 7           ; num
XEQ 18          ; num*(inv den) mod av
2
RCL 12          ; vmax
RCL 6           ; v
-
Y^X             ; 2^(vmax-v)
RCL 15          ; k
25
*
3
-     
*               ; 2^(vmax-v) * (25k-3)
XEQ 18          ; (num*(inv den) mod av) * above mod av
RCL 11          ; s
+
RCL 8           ; av
MOD             ; exp + s mod av
STO 11          ; s =
LBL 12
1
ST+ 15          ; k=k+1
RCL 16          ; N
RCL 15          ; k
X<=Y?
GTO 10
RCL 17          ; n
1
-               ; n-1
5
XEQ 21          ; 5^(n-1) mod av
RCL 11          ; s
XEQ 18          ; s*(expr) mod av
RCL 8           ; av
/
STO 13          ; sum=
LBL 09
; cases 3,5,...

3
STO 14          ; a=3
LBL 13          ; for(a=3, a<3N; a=nextprime(a)) ..
RCL 14          ; a
RCL 16          ; N
3
*
LN              ; ln(3N)
RCL 14          ; a
LN        
/             
INT
STO 12          ; vax = ln(3N)/ln(a)
Y^X             ; a^(vmax)
STO 8           ; av=
RCL 17          ; n
2
XEQ 20          ; 2^n mod av
STO 7           ; num=
0
STO 11          ; s=0
STO 6           ; v=0
STO 18          ; kq1 = 0
1
STO 10          ; den=1
STO 15          ; k=1
CHS
STO 19          ; kq2=-1
2
CHS
STO 21          ; kq4=-2
3
CHS
STO 20          ; kq3=-3

LBL 14          ; for (k=1; k<=N; ++k) ..
RCL 15          ; k
2
ST+ 18          ; kq1 += 2
*               ; 2k
STO 0           ; t=2k
18
STO 5           ; use kq1
1
CHS
STO 9           ; dec v
XEQ 15          ; divn(t, a, kq1)
RCL 0
RCL 7           ; num
*
RCL 8           ; av
MOD
STO 7           ; num = num * t mod av
RCL 15          ; k
2
ST+ 19          ; kq2 += 2
*
1
ST+ 5           ; select kq2
-               ; 2k-1
STO 0           ; t=
XEQ 15          ; divn(t, a, kq2)
RCL 0
RCL 7           ; num
*
RCL 8           ; av
MOD
STO 7           ; num = num * t mov av
RCL 15          ; k
9
ST+ 20          ; kq3 += 9
*
3
-
STO 0           ; t=9k-3
1
STO 9           ; inc v signal
ST+ 5           ; select kq3
XEQ 15          ; divn(t, a, kq3)
RCL 0
RCL 10          ; den
*
RCL 8           ; av
MOD
STO 10          ; den = den * t mod av
RCL 15          ; k
3
ST+ 21          ; kq4 += 3
*
2
-
STO 0           ; t=3k-2
1
ST+ 5           ; select kq4
XEQ 15          ; divn(t, a, kq4)
RCL 0
2
*
RCL 10
*
RCL 8
MOD
STO 10          ; den = den*2*t mod av
RCL 6           ; v
X<=0?
GTO 17
RCL 12          ; vmax
RCL 6           ; v
-
RCL 14          ; a
XEQ 20          ; a^(vmax-v) mod av
RCL 15          ; k
25
*
3
-
*               ; *(25k-3)
RCL 8           ; av
MOD
STO 9           ; t1=
RCL 8           ; av
RCL 10          ; den
XEQ 19          ; den^-1 mod av
RCL 7           ; num
*
RCL 8           ; av
MOD
RCL 9           ; t1
*
RCL 8           ; av
MOD
RCL 11          ; s
+
RCL 8
MOD
STO 11          ; s += t mod av
LBL 17
1
ST+ 15          ; k+=1
RCL 16          ; N
RCL 15          ; k
X<=Y?
GTO 14          ; back for for(k=..)

RCL 17          ; n
1
-
5
XEQ 20          ; 5^(n-1) mod av
RCL 11          ; s
*
RCL 8
MOD
RCL 8
/
RCL 13
+
FRC
STO 13          ; sum

; a<3N... next prime

RCL 14         ; a
XEQ 23         ; next prime
STO 14
RCL 16         ; N
3
*
X<>Y
X<=Y?
GTO 13          ; back to for a=3..

RCL 13          ; sum
RTN             ; display!

LBL 15          ; DIVN(r0, a, kqn)
RCL IND 5       ; kqn
RCL 14          ; a
X>Y?
RTN             ; kqn < a
MOD             ; kqn mod a
STO IND 5
X<>0?
RTN
LBL 16
RCL 14          ; a
ST/ 0           ; r0 /= a
RCL 9
ST+ 6           ; v +/- inc
RCL 0           ; t
RCL 14
MOD             ; t mod a
X=0?
GTO 16          ; while mod = 0
RTN

LBL 11          ; deflate r0 of 2, count in V
RCL 0
2
MOD
X<>0?
RTN
2
ST/ 0
RCL 1
ST+ 6           ; v = v + r1
GTO 11
LBL 18
; Y = b
; X = a
; R9 = B
; R8 = M
; a*b mod M where M=2^k
; USES R0,1
ENTER
ENTER
RCL 9   ; B
/
INT
STO 0   ; a1
RCL 9
*
-
STO 1   ; a0
X<>Y    ; b
ENTER
ENTER
RCL 9   
/
INT     ; b1
ST* 1   ; a0*b1
RCL 9
*
-       ; b0
ST* 0   ; a1*b0
*       ; a0*b0
RCL 0   ; a1*b0
RCL 1
+       ; a1*b0+a0*b1
RCL 9
MOD     ; a1*b0+a0*b1 MOD B
RCL 9
*       ; (a1*b0+a0*b1 MOD B) * B
+       ; (a1*b0+a0*b1 MOD B)*B + a0*b0
RCL 8   ; M
MOD     
RTN
LBL 19
; Y 
; X
; X^-1 mod Y
; USES R012345

STO 0  ; U=X
RDN    ; Y
STO 1  ; V=Y
STO 5  ; Y
1
STO 2   ; C=1
0
STO 3   ; A=0
LBL 00
RCL 3   ; A
RCL 1   ; V
RCL 0   ; U
/   
INT     ; V/U
STO 4   ; Q
RCL 2   ; C
STO 3   ; A=C
*       ; Q*C
-       ; A-Q*C
STO 2   ; c=a-q*c
RCL 0   ; U
RCL 1   ; V
RCL 4   ; Q
RCL 0   ; U
STO 1   ; V=U
*       ; Q*U
-       ; v-qu
STO 0   ; u=v-qu
X<>0?   
GTO 00
RCL 3   ; A
RCL 5   ; Y
MOD
RTN
LBL 20
; R8=M
; Y=b
; X=a
; a^b mod M
; USES R012
STO 0    ; a
X<>Y
STO 1    ; b
1
STO 2   ;r=1
LBL 01
RCL 1   ; b
RCL 1
2
/
INT
STO 1   ;b=b/2
2
*
-       ;b mod 2
X=0?
GTO 02  ;b&1 == 0
RCL 0   ; a
RCL 2  ; r
*       ;r*a
RCL 8
MOD  
STO 2   ; r = r*a mod M
LBL 02
RCL 1   ; b
X=0?
GTO 03
RCL 0
RCL 0
*
RCL 8
MOD
STO 0  ; a=a*a mod M
GTO 01
LBL 03
RCL 2
RTN
LBL 21
; R8=M
; R9=B
; Y=b
; X=a
; a^b mod M, where M=2^k
STO 2    ; a
X<>Y
STO 3    ; b
1
STO 4   ;r=1
LBL 04
RCL 3   ; b
RCL 3
2
/
INT
STO 3   ;b=b/2
2
*
-       ;b mod 2
X=0?
GTO 05  ;b&1 == 0
RCL 2   ; a
RCL 4  ; r
XEQ 18
STO 4  ; r = r*a mod M
LBL 05
RCL 3   ; b
X=0?
GTO 06
RCL 2
RCL 2
XEQ 18
STO 2  ; a=a*a mod M
GTO 04
LBL 06
RCL 4
RTN
LBL 23
; nextprime (X)
; USES R012
STO 2
LBL 08
2
ST+ 2
RCL 2
SQRT
STO 0   ; r = sqrt(x)
3
STO 1   ; i=3
LBL 07
RCL 2   
RCL 1
MOD     ; X mod i
X=0?
GTO 08
RCL 0
RCL 1
2
+
STO 1
X<=Y?
GTO 07
RCL 2
RTN
END