Post Reply 
Creating digits of pi
02-09-2018, 01:40 PM
Post: #21
RE: Creating digits of pi
(02-09-2018 12:54 PM)toml_12953 Wrote:  
(02-09-2018 12:08 PM)Mike (Stgt) Wrote:  Sorry, I did not mention that you should take it as pseudocode only. Well, it is REXX what I boiled down to the essential loop. So a 'minor' detail went up in smoke:

numeric digits 500

what sets precision of all variables and has its effect on the loop end criterion DO..UNTIL m = m + n what is -- I confess -- an awesome concealment or simple substitution of the compare n/m < epsilon (I quit the loop when it has no influence on the result any more).

How to run REXX? Try Regina, or fast brexx from CERN, or Open Object Rexx (I used it for a HP7470A simulation), or NetRexx, or -- very outstanding -- get a mainframe under Hercules on your machine (which could be a Raspberry Pi - tested, works).

Ciao.....Mike

Here's a version in Decimal BASIC interpreter that returns 1000 digits before you can remove your finger from the Enter key on a 4.2 GHz Core i7:
(PI is a reserved word in Decimal BASIC) Wait a second. That was fun. I didn't just have fun did I? No, not me. I refuse to have fun.

Code:
OPTION ARITHMETIC DECIMAL_HIGH ! 1000 digits precision
LET s = 1
LET n = 1
LET m = 3
DO UNTIL m = m + n
   LET n = n / 9
   LET s = s + n / (m + 2) - 3 * n / m
   LET m = m + 4
LOOP
LET PIE = 2 * s * SQR(3)
PRINT PIE
END

Or the equivalent in bc:

Code:
scale = 1000
s=1
n=1
m=3
while (m != m + n) {
    n = n/9
    s = s+n / (m+2) - 3*n/m
    m = m+4
}
pi = 2*s*sqrt(3)
print pi

Greetings,
    Massimo

-+×÷ ↔ left is right and right is wrong
Visit this user's website Find all posts by this user
Quote this message in a reply
02-09-2018, 02:03 PM
Post: #22
RE: Creating digits of pi
(02-09-2018 12:54 PM)toml_12953 Wrote:  
Code:
OPTION ARITHMETIC DECIMAL_HIGH ! 1000 digits precision
...
And OPTION ARITHMETIC DECIMAL_double_HIGH will get you the next 1000 decimals? With REXX the decimals setting is only limited by available storage.

BTW, forgot to mention an other way you may get in contact with REXX: PC DOS 7. When we like obsolete calculators why not ancient software too.
/M.
Find all posts by this user
Quote this message in a reply
02-09-2018, 02:59 PM (This post was last modified: 02-09-2018 02:59 PM by pier4r.)
Post: #23
RE: Creating digits of pi
wait a moment. REXX for the little that I know should have been designed as a sort of bash (or vbscript & co). How come that it has such a library for arbitrary precision? Or did you add a library yourself?

Edit: nice the equivalent in bc. This let me remember to ask something in another thread.

Wikis are great, Contribute :)
Find all posts by this user
Quote this message in a reply
02-09-2018, 03:35 PM
Post: #24
RE: Creating digits of pi
(02-09-2018 02:03 PM)Mike (Stgt) Wrote:  
(02-09-2018 12:54 PM)toml_12953 Wrote:  
Code:
OPTION ARITHMETIC DECIMAL_HIGH ! 1000 digits precision
...
And OPTION ARITHMETIC DECIMAL_double_HIGH will get you the next 1000 decimals? With REXX the decimals setting is only limited by available storage.

BTW, forgot to mention an other way you may get in contact with REXX: PC DOS 7. When we like obsolete calculators why not ancient software too.
/M.

No, Decimal BASIC only has up to 1000 digits. Not bad considering most versions of BASIC only have single precision (about 7 digits) or double precision (about 16 digits)

Tom L

DM42 SN: 00025 (Beta)
SN: 00221 (Production)
Find all posts by this user
Quote this message in a reply
02-09-2018, 04:26 PM
Post: #25
RE: Creating digits of pi
(02-09-2018 02:59 PM)pier4r Wrote:  wait a moment. REXX for the little that I know should have been designed as a sort of bash (or vbscript & co). How come that it has such a library for arbitrary precision? Or did you add a library yourself?

No additional library, it is by design of REXX to set any precision you like, only limited by the available storage. This is possible as there is only one variable type: text. Text strings of 'unlimited' length. When it comes to math operations it is tested if the arguments may be interpreted as numbers. Slow but smart. Alas only "4-banger-math" plus exponential function, modulo and integer division. For those who need a SIN in a scripting language there are libraries for trig- and log-math, up to 10 digits precision or so. To get more digits I use a Sqrt of my own:
Code:
Sqrt: procedure; arg z
a = RxCalcSqrt(z)       /* get a first approximation from the lib function */
do until b = a
   b = a
   a = .5 * (a + z / a)
end
return a
This is only the essential code, to speed up execution I start with a precision of about 20 digits only and increase it with every iteration up to the desired one. Anyway I am looking for a faster replacement. Then Brent-Salamin formula or a Borwein&Borwein algorithm could outreach the Chudnovsky formula, I assume.

Ciao.....Mike
Find all posts by this user
Quote this message in a reply
02-09-2018, 07:27 PM (This post was last modified: 02-10-2018 01:33 AM by Gerson W. Barbosa.)
Post: #26
RE: Creating digits of pi
(02-09-2018 12:54 PM)toml_12953 Wrote:  Here's a version in Decimal BASIC interpreter that returns 1000 digits before you can remove your finger from the Enter key on a 4.2 GHz Core i7:
(PI is a reserved word in Decimal BASIC) Wait a second. That was fun. I didn't just have fun did I? No, not me. I refuse to have fun.

Code:
OPTION ARITHMETIC DECIMAL_HIGH ! 1000 digits precision
LET s = 1
LET n = 1
LET m = 3
DO UNTIL m = m + n
   LET n = n / 9
   LET s = s + n / (m + 2) - 3 * n / m
   LET m = m + 4
LOOP
LET PIE = 2 * s * SQR(3)
PRINT PIE
END

Perhaps a little more fun here, having written my own code using my own algorithm. Well, actually Madhava's since it was probably known to him, given his exact correction terms. Three digits per iteration, but noticeably slower.
Code:

OPTION ARITHMETIC DECIMAL_HIGH ! 1000 digits precision
LET d=1000                     ! 1000 decimal places
LET n=INT(D/3)+1
LET a=2*n
LET b=8*n  
LET s0=0
LET s1=0
FOR i=n TO 1 STEP -1
   LET s0=s0+2/(16*i*(i-1)+3)
   LET t=4*i
   LET s1=(t*(i-1)+1)/(a+t*i/(b+s1))
NEXT i
LET s1=1/(8*n+s1)
LET p=4*(s0+s1)
PRINT p
END

\[\frac{\pi }{4}= 1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\cdots +\frac{1}{2n-3}-\frac{1}{2n-1}+\frac{1}{4n+\frac{1^{2}}{n+\frac{2^{2}}{4n+\frac{3^{2}}{n+\frac{4^{2}}{4n+...​​​ }}}}}\]

Edited to improve code (removed 2 multiplications out of the loop), but I am not familiar with this BASIC to try other optimizations.

Edited again to include a somewhat more optimized code:

Code:

OPTION ARITHMETIC DECIMAL_HIGH ! 1000 digits precision
LET t1=TIME
LET d=1000                     ! 1000 decimal places
LET n=INT(D/3) + 1
LET a=2*n
LET b=8*n  
LET s0=0
LET s1=0
FOR i=n TO 1 STEP -1
   LET s0=s0 + 2/(16*i*(i - 1) + 3)
   LET t=4*i
   LET s1=((b + s1)*((i - 1)*t + 1))/(a*(b + s1) + i*t)
NEXT i
LET s1=1/(b + s1)
LET p=4*(s0 + s1)
PRINT p
PRINT TIME - t1 
END

334 iterations (but 327 would have sufficed for 1000 decimal places), 342 ms (average of 5 runs). Still slower than Tom's version (243 ms here, Intel E4700 @ 2.6 GHz), but about 130 % faster when compared to my previous version above).
Find all posts by this user
Quote this message in a reply
02-09-2018, 08:03 PM
Post: #27
RE: Creating digits of pi
(02-09-2018 02:01 AM)Mike (Stgt) Wrote:  BTW, instead of tantalizing your real machine for 6 hrs, I suggest a well known emulator of it. Recently I set it to authentic speed for some reason, and I am so glad I can switch it back to full speed. Could help in this case too.

I certainly tried it on Debug4x first (Emu-50g @ 2.6 GHz, 555.7218 seconds). I had estimated it would take 5 or 6 hours on the real 50g, but I wanted to get the exact factor as I was not so sure about that.

Regards,

Gerson.
Find all posts by this user
Quote this message in a reply
02-10-2018, 12:44 PM
Post: #28
RE: Creating digits of pi
Could someone of the interested (Decimal BASIC admirer, etc.) please confirm this Ramanujan–Sato of 1993 works? I fail with it, get Pi = 4.1748739, what is a bit too far off.
I assume a typo in the Wiki.

Ciao.....Mike
Find all posts by this user
Quote this message in a reply
02-10-2018, 03:45 PM (This post was last modified: 02-10-2018 03:48 PM by Gerson W. Barbosa.)
Post: #29
RE: Creating digits of pi
(02-10-2018 12:44 PM)Mike (Stgt) Wrote:  Could someone of the interested (Decimal BASIC admirer, etc.) please confirm this Ramanujan–Sato of 1993 works? I fail with it, get Pi = 4.1748739, what is a bit too far off.
I assume a typo in the Wiki.

Ciao.....Mike


Off by a factor of exactly 400/301:

4.1748739 * 301/400 = 3.1415926

Probably just a coincidence, but I would check this out with more digits.
Find all posts by this user
Quote this message in a reply
02-10-2018, 04:42 PM (This post was last modified: 02-10-2018 05:11 PM by toml_12953.)
Post: #30
RE: Creating digits of pi
(02-10-2018 12:44 PM)Mike (Stgt) Wrote:  Could someone of the interested (Decimal BASIC admirer, etc.) please confirm this Ramanujan–Sato of 1993 works? I fail with it, get Pi = 4.1748739, what is a bit too far off.
I assume a typo in the Wiki.

Ciao.....Mike

The quadratic convergence works:

Code:
OPTION ARITHMETIC DECIMAL_HIGH
LET a=SQR(2)
LET b=0
LET p=2+SQR(2)
LET n=1
DO UNTIL a=1
   LET an=a
   LET bn=b
   LET pn=p
   LET a=(SQR(an)+1/SQR(an))/2
   LET b=((1+bn)*SQR(an))/(an+bn)
   LET p=((1+a)*pn*b)/(1+b)
LOOP
PRINT p
END

Tom L

DM42 SN: 00025 (Beta)
SN: 00221 (Production)
Find all posts by this user
Quote this message in a reply
02-10-2018, 04:45 PM
Post: #31
RE: Creating digits of pi
(02-10-2018 03:45 PM)Gerson W. Barbosa Wrote:  Off by a factor of exactly 400/301:

No, was a very late night typo on my side, one constant 100 times to small. Sad
Tnx anyway.
/M.
Find all posts by this user
Quote this message in a reply
02-10-2018, 04:54 PM
Post: #32
RE: Creating digits of pi
(02-10-2018 04:42 PM)toml_12953 Wrote:  ... Something is off there somewhere.

Yes, have a closer look to the indices.

Try this:
Code:
an = sqrt(2)
bn = 0
pn = 2 + an
do until p = pn
   a = an
   b = bn
   p = pn
   ar = sqrt(a)
   an = .5 * (ar + 1 / ar)
   bn = ar * (1 + b) / (a + b)
   pn = p * bn * (1 + an) / (1 + bn)
end
say p    /* or the last 30 digits if it's too long */
an stands for a(n+1), a stands for a(n), and so on.

Ciao.....Mike
Find all posts by this user
Quote this message in a reply
02-10-2018, 05:19 PM
Post: #33
RE: Creating digits of pi
(02-10-2018 04:54 PM)Mike (Stgt) Wrote:  
(02-10-2018 04:42 PM)toml_12953 Wrote:  ... Something is off there somewhere.

Yes, have a closer look to the indices.

Try this:
Code:
an = sqrt(2)
bn = 0
pn = 2 + an
do until p = pn
   a = an
   b = bn
   p = pn
   ar = sqrt(a)
   an = .5 * (ar + 1 / ar)
   bn = ar * (1 + b) / (a + b)
   pn = p * bn * (1 + an) / (1 + bn)
end
say p    /* or the last 30 digits if it's too long */
an stands for a(n+1), a stands for a(n), and so on.

Ciao.....Mike
I corrected my original code. It works now. Thanks!

Tom L

DM42 SN: 00025 (Beta)
SN: 00221 (Production)
Find all posts by this user
Quote this message in a reply
02-10-2018, 05:25 PM
Post: #34
RE: Creating digits of pi
(02-10-2018 12:44 PM)Mike (Stgt) Wrote:  Could someone of the interested (Decimal BASIC admirer, etc.) please confirm this Ramanujan–Sato of 1993 works? I fail with it, get Pi = 4.1748739, what is a bit too far off.
I assume a typo in the Wiki.

Nope, the Wiki is correct. I used the constants featured there and this is the result I got:

list

10 word -520:point -160
20 A=63365028312971999585426220+28337702140800842046825600*sqrt(5)+384*sqrt
(5)*sqrt(10891728551171178200467436212395209160385656017+48709290865788102250773​
38534541688721351255040*sqrt(5))
30 B=7849910453496627210289749000+3510586678260932028965606400*sqrt(5)+2515
968*sqrt(3110)*sqrt(6260208323789001636993322654444020882161+2799650273060444296​
577206890718825190235*sqrt(5))
40 C=-214772995063512240-96049403338648032*sqrt(5)-1296*sqrt(5)*sqrt(109852
34579463550323713318473+4912746253692362754607395912*sqrt(5))
50 clr time : S=0:N=0
60 T=S+!(6*N)/!(3*N)/(!(N)^3)*(A+N*B)/C^(3*N):if T<>S then S=T:inc N:goto 60
70 print sqrt((-C)^3)/S:print time1000


run

3.14159265358979323846264338327950288419716939937510582097494459230781640628620
89986280348253421170679821480865132823066470938446095505822317253594081284811174​
50284102701938521105559644622948954930381964428810975665933446128475648233786783​
16527120190914564856692346034861045432664821339360726024914127372458700660631558​
81748815209209628292540917153643678925903600113305305488204665213841469519415116​
09433057270365759591953092186117381932611793105118548074462379962749567351885752​
72489122793818301194912983367336244065664308602139494639522473719070217986094370​
27705392171762931767523846748184676694051320005681271452635608277857713427577896​
09173637178721468440901224953430146549585371050792279689258923542019956112129021​
96086403441815981362977477130996051870721134999999837

which gives 770 correct decimals (all featured) in 15-16 iterations, 8 milliseconds.

Have a nice weekend.
V.
.
Find all posts by this user
Quote this message in a reply
02-10-2018, 06:03 PM
Post: #35
RE: Creating digits of pi
(02-10-2018 05:25 PM)Valentin Albillo Wrote:  Nope, the Wiki is correct.
Thank you, I've found my error.

Quote:Have a nice weekend.
Same to you.
/M.
Find all posts by this user
Quote this message in a reply
02-10-2018, 09:14 PM (This post was last modified: 02-10-2018 09:20 PM by emece67.)
Post: #36
RE: Creating digits of pi
I've tested, using Python + the mpmath package, the nonic convergence method. It gives
  • 1 correct digit at the 1st iteration
  • 22 at the 2nd
  • 217 at the 3rd
  • 1984 at the 4th
  • 17898 at the 5th
  • 161123 at the 6th
  • 1450163 at the 7th
Really awesome!

Regards.

César - Information must flow.
Find all posts by this user
Quote this message in a reply
02-11-2018, 12:12 AM
Post: #37
RE: Creating digits of pi
Yeah, not helpful at all, but for some random digit of Pi in base 2, just guess.

50/50 chance of being right is good enough for me and I don't need to put the batteries back in my HP41 to do it.

Wink

2speed HP41CX,int2XMEM+ZEN, HPIL+DEVEL, HPIL+X/IO, I/R, 82143, 82163, 82162 -25,35,45,55,65,67,70,80
Find all posts by this user
Quote this message in a reply
02-11-2018, 11:21 PM (This post was last modified: 02-11-2018 11:35 PM by Gerson W. Barbosa.)
Post: #38
RE: Creating digits of pi
(02-09-2018 07:27 PM)Gerson W. Barbosa Wrote:  
Code:

OPTION ARITHMETIC DECIMAL_HIGH ! 1000 digits precision
LET t1=TIME
LET d=1000                     ! 1000 decimal places
LET n=INT(D/3) + 1
LET a=2*n
LET b=8*n  
LET s0=0
LET s1=0
FOR i=n TO 1 STEP -1
   LET s0=s0 + 2/(16*i*(i - 1) + 3)
   LET t=4*i
   LET s1=((b + s1)*((i - 1)*t + 1))/(a*(b + s1) + i*t)
NEXT i
LET s1=1/(b + s1)
LET p=4*(s0 + s1)
PRINT p
PRINT TIME - t1 
END

334 iterations (but 327 would have sufficed for 1000 decimal places), 342 ms (average of 5 runs). Still slower than Tom's version (243 ms here, Intel E4700 @ 2.6 GHz), but about 130 % faster when compared to my previous version above).

Now I get running times as lows as 0.2 seconds, but sometimes as high as 0.35 seconds.

Code:

OPTION ARITHMETIC DECIMAL_HIGH ! 1000 digits precision
LET t = TIME
LET d = 1000                   ! 1000 decimal places
LET n = INT(D/3+1/2)
LET a = 2*n
LET b = 8*n 
LET c = 4*n*n
LET d = 4*n - 1 
LET s0 = 0
LET s1 = 0
FOR i = 1 TO n
   LET s0 = s0 + 2/((d - 2)*d)
   LET s1 = (b + s1)*(c - d)/(a*(b + s1) + c)
   LET c = c - d - d + 2   
   LET d = d - 4
NEXT i
LET s1 = 1/(b + s1)
LET p=4*(s0 + s1)
PRINT p
PRINT
PRINT TIME - t; "seconds"  
END

3.141592653589793238462643383279502884197169399375105820974944592307816406286208​99862803482534211706798214808651328230664709384460955058223172535940812848111745​02841027019385211055596446229489549303819644288109756659334461284756482337867831​65271201909145648566923460348610454326648213393607260249141273724587006606315588​17488152092096282925409171536436789259036001133053054882046652138414695194151160​94330572703657595919530921861173819326117931051185480744623799627495673518857527​24891227938183011949129833673362440656643086021394946395224737190702179860943702​77053921717629317675238467481846766940513200056812714526356082778577134275778960​91736371787214684409012249534301465495853710507922796892589235420199561121290219​60864034418159813629774771309960518707211349999998372978049951059731732816096318​59502445945534690830264252230825334468503526193118817101000313783875288658753320​83814206171776691473035982534904287554687311595628638823537875937519577818577805​32171226806613001927876611195909216420199

.2 seconds


These are the more straightforward expressions for s0 and s1, but they're slower:

Code:

LET s0 = s0 - 1/d + 1/(d - 2)
LET s1 = (c - d)/(a + c/(b + s1))
Find all posts by this user
Quote this message in a reply
02-12-2018, 01:38 AM
Post: #39
RE: Creating digits of pi
(02-09-2018 07:27 PM)Gerson W. Barbosa Wrote:  334 iterations (but 327 would have sufficed for 1000 decimal places), 342 ms

Try this
Code:
if arg() then numeric digits arg(1)      /* no argument -> default */
numeric fuzz RxCalcLog10(digits()) % 1   /* ignore some digits at compare */
e = 10 ** -(digits() - fuzz() - 1)       /* epsilon */
a = 1
b = 1 / 2
r = 1
s = b
x = a
DO UNTIL k < e
   b = Sqrt(x * b)
   x = a
   a = (a + b) / 2
   c = (a - b) ** 2
   r = r + r
   k = c * r
   s = s - k
end
pi = 2 * a * a / s
1o interations for 1000 digits.

BTW, 1000 decimals in 10 minutes on an HP-49, that gives everyone the chance to keep count of it.

Ciao.....Mike
Find all posts by this user
Quote this message in a reply
02-12-2018, 03:40 AM
Post: #40
RE: Creating digits of pi
(02-12-2018 01:38 AM)Mike (Stgt) Wrote:  BTW, 1000 decimals in 10 minutes on an HP-49, that gives everyone the chance to keep count of it.

Really fast (and compact), 5 seconds on the emulator!

For real fun, I would try a 2500-year-old algorithm (with a small adjustment for speed) on my 50g. No problem waiting three hours for the glorious thousand digits :-)

http://www.hpmuseum.org/cgi-sys/cgiwrap/...443#189221

There is also a 4-minute alternative (certainly more than 10 minutes on the 49G). Both programs require the LongFloat library, though.
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)