Omnimaga
Calculator Community => Other Calc-Related Projects and Ideas => TI Z80 => Topic started by: JWinslow23 on May 04, 2015, 01:27:55 pm
-
...in vanilla TI-BASIC. To a precision of more than 15 places.
Well, with the help of Jeremy Gibbons (http://www.cs.ox.ac.uk/jeremy.gibbons/publications/spigot.pdf) and myself, now there's hope!
EDIT: Now there's MORE hope...click here (https://www.omnimaga.org/ti-z80-calculator-projects/how-i-wish-i-could-calculate-pi/msg400899/#msg400899)
This program is as fast as I could make it (though not too fast...this ain't Asm), and I welcome any improvements on size and speed. It uses simulated arbitrary precision arithmetic with lists. I have verified it for 100 digits, and the algorithm itself has been verified for 1,000 digits in Haskell.
Download link coming soon attached.
Code is as follows:
PICALC.8xp{1→ʟQ
{180→ʟR
{60→ʟT
2→I
ClrHome
501→D
"CHANGE THE ABOVE NUMBER TO CALCULATE MORE OR LESS DIGITS
"?→Str1
For(J,1,D
Output(1,1,Str1
Output(8,1,J
"ABOVE FOR TESTING
3(3I+1)(3I+2→U
5ʟT→L1
prgmAADJUST
L1→L3
5ʟR→L1
prgmAADJUST
L1→L2
(27I-12)ʟQ→L1
prgmAADJUST
L1+L2→L1
prgmAADJUST
L3→L2
prgmADIVIDE
θ→Y
Str1+sub("0123456789",Y+1,1→Str1
(5I-2)ʟQ→L1
prgmAADJUST
L1→L2
YʟT→L1
prgmAADJUST
L2-L1+ʟR→L1
prgmAADJUST
10UL1→L1
prgmAADJUST
L1→ʟR
10I(2I-1)ʟQ→L1
prgmAADJUST
L1→ʟQ
UʟT→L1
prgmAADJUST
L1→ʟT
I+1→I
EndAADJUST.8xpDelVar BDelVar C1→A
Repeat B or A>dim(L1
C+L1(A→L1(A
If not(Ans
A=dim(L1)+1-sum(int(1/(1+cumSum(abs(seq(L1(Z),Z,dim(L1),1,-1→B
int(10^(7)⁻¹Ans→C
If 0>L1(A
-10^(7)Ans+L1(A→L1(A
10^(7)fPart(10^(7)⁻¹L1(A→L1(A
A+1→A
If C and Ans>dim(L1
Then
0→L1(A
1+dim(ʟQ→dim(ʟQ
Ans→dim(ʟR
Ans→dim(ʟT
Ans→dim(L2
Ans→dim(L3
End
EndACOMPARE.8xpDelVar Gdim(L1→A
Repeat G or A<1
(L1(A)>L2(A))-(L1(A)<L2(A→G
A-1→A
EndADIVIDE.8xpDelVar θprgmACOMPARE
While G≥0
Output(6,1,"
Output(6,1,θ
"ABOVE FOR TESTING
L1-L2→L1
prgmAADJUST
θ+1→θ
prgmACOMPARE
End
-
This sound awesome!
Would you mind explaining some logic behind it / posting the Haskel program?
-
This sound awesome!
Would you mind explaining some logic behind it / posting the Haskel program?
All the logic is explained in the paper linked at the top. The program I am referring to is the last program listed in the paper:
> piG3 = g(1,180,60,2) where
> g(q,r,t,i) = let (u,y)=(3*(3*i+1)*(3*i+2),div(q*(27*i-12)+5*r)(5*t))
> in y : g(10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)I'm not a math wizard or anything :P , but I think it's based off of this series for pi:
(http://latex.codecogs.com/gif.latex?%5Cpi%20%3D%203+%5Cfrac%7B1%5Ctimes1%7D%7B3%5Ctimes4%5Ctimes5%7D%5Ctimes%288+%5Cfrac%7B2%5Ctimes3%7D%7B3%5Ctimes7%5Ctimes8%7D%5Ctimes%28%5Ccdots5i-2+%5Cfrac%7Bi%282i-1%29%7D%7B3%283i+1%29%283i+2%29%7D%5Ctimes%5Ccdots%29%29)
And it views that as an infinite composition of linear fractional transforms. The program also has an optimization based off a conjecture that has not yet been fully proven, but it has been verified to hold for calculation of the first 1,000 digits.
Hope I helped! Check the paper for more: this guy knows more about the subject than I ever would hope to tell you :P .
Of course, if you don't mind a previous commitment to the number of digits calculated, I have one that fits entirely within the precision of the calculator's numbers... that one fails before the first 0 and I don't know why
-
I have new code! With this one, I tested that the maximum amount of digits it can supply on a fully RAM cleared calculator is 377, and this is done so in 3 hours 45 minutes. (don't worry, it doesn't go at the speed of 1.675 minutes per digit...it just gets slower with digits calculated, as expected)
Download links attached. You will be amazed or else :P
Code is as follows:
PICALC.8xpSetUpEditor Q,T,R
{1->|LQ
{60->|LT
3Ans->|LR
ClrHome
377->D
"?->Str1
For(J,2,D+1
Output(1,1,Str1
Output(6,1,"
Output(8,1,J-1
5|LT->L1
prgmAADJUST
L1->L2
(27J-12)|LQ+5|LR->L1
prgmAADJUST
prgmADIVIDE
Str1+sub("0123456789",theta+1,1->Str1
30((5J-2)|LQ-theta|LT+|LR->L1
prgmAADJUST
(3J+1)L1->L1
prgmAADJUST
(3J+2)L1->L1
prgmAADJUST
L1->|LR
10J(2J-1)|LQ->L1
prgmAADJUST
L1->|LQ
3(3J+1)|LT->L1
prgmAADJUST
(3J+2)L1->L1
prgmAADJUST
L1->|LT
EndAADJUST.8xpDelVar C
For([recursiven],1,dim(L1
C+L1([recursiven]->|N
int(Ans|E~7->C
|N-Ans|E7->L1([recursiven]
End
If C
C->L1([recursiven]
dim(L1->dim(L2
Ans->dim(|LQ
Ans->dim(|LR
Ans->dim(|LTACOMPARE.8xp1+dim(L1->A
Repeat Ans or A<2
A-1->A
(L1(A)>L2(A))-(L1(A)<L2(A
EndADIVIDE.8xpDelVar thetaprgmACOMPARE
While Ans>=0
Output(6,1,"
Output(6,1,theta
L1-L2->L1
prgmAADJUST
theta+1->theta
prgmACOMPARE
End
-
Bump.
Even newer code. Now, no divide program. Also a bit faster.
Download links attached, here is the code to type in:
PICALC.8xpSetUpEditor Q,T,R
{1->|LQ
{60->|LT
3Ans->|LR
ClrHome
Output(7,1,"CLEAR TO EXIT. # NEXT: 0...
"?->Str1
2->I%
Repeat getKey=45
If 1<length(Str1
Then
"
If not(fPart((length(Str1)-2)/16
Output(6,1,Ans+Ans+Ans+Ans
2+16int((length(Str1)-82)/16)(98<=length(Str1
Output(1,1,sub(Str1,Ans,length(Str1)+1-Ans
End
Output(8,2,I%-2
Output(8,12,0
5|LT->L1
prgmAADJUST
L1->L2
(9I%-4)|LQ->L1
prgmAADJUST
3L1+5|LR->L1
prgmAADJUST
0->PV
prgmACOMPARE
While Ans>=0
PV+1->PV
Output(8,12,PV
L1-L2->L1
prgmAADJUST
prgmACOMPARE
End
sub("0123456789",PV+1,1
If I%=2
Ans+".
Str1+Ans->Str1
(5I%-2)|LQ+|LR-PV|LT->L1
prgmAADJUST
30L1->L1
prgmAADJUST
If I%<33
Then
9I%^^2+9I%+2
Else
(3I%+1)L1->L1
prgmAADJUST
3I%+2
End
AnsL1->L1
prgmAADJUST
L1->|LR
10I%|LQ->L1
prgmAADJUST
(2I%-1)L1->L1
prgmAADJUST
L1->|LQ
3|LT->L1
prgmAADJUST
If I%<33
Then
9I%^^2+9I%+2
Else
(3I%+1)L1->L1
prgmAADJUST
3I%+2
End
AnsL1->L1
prgmAADJUST
L1->|LT
I%+1->I%
End
ClrHome
sub(Str1,2,length(Str1)-1->Str1
Disp "# OF DIGITS:","DIGITS STORED IN","Str1
Output(1,14,I%-2AADJUST.8xpIf min(L1)>0 and |E10>=max(L1
Return
0->PMT
For([recursiven],1,dim(L1
PMT+L1([recursiven]->|N
int(|N|E~10->PMT
|N-PMT|E10->L1([recursiven]
End
If not(PMT
Return
[recursiven]->dim(L2
Ans->dim(|LQ
Ans->dim(|LR
Ans->dim(|LT
PMT->L1(AnsACOMPARE.8xp1+dim(L1->|N
Repeat Ans or |N=1
|N-1->|N
L1(|N)-L2(|N
End