Creating digits of pi

[brickviking]

Yes, this old chestnut again. Except my guugle-foo has foo'ed and I can't find anything in HP-50g that's relevant. Put simply, I want to have a procedure that generates x digits of PI starting from an offset, i.e. 1:offset, 2:digit_count. An example would be ten digits of PI starting from the 49,238th position.

I don't know of good programs that can generate only a subset of digits, only of programs that can generate x digits of PI as a string (PI50G being one example I came across). Does anyone have a good pointer to what I'd like? (Need is a strong word, perhaps desire just for the sake of it).

(Post 168)

[Joe Horn]

There *is* a way of directly calculating the Nth digit of pi (without generating all the previous digits) but only in hex. Nobody to date has discovered a way of doing it for pi in base 10.

[toml_12953]

(02-06-2018 01:38 AM)Joe Horn Wrote:  There *is* a way of directly calculating the Nth digit of pi (without generating all the previous digits) but only in hex. Nobody to date has discovered a way of doing it for pi in base 10.

Why not just convert the hex to decimal?

[brickviking]

(02-06-2018 09:24 PM)toml_12953 Wrote:  Why not just convert the hex to decimal?

That's not the easiest thing in the world for an inexperienced user like me. I can certainly convert integers, but converting stuff to the right of the hexadecimal point is quite a bit beyond me; though I do understand it uses powers of two, it uses them in a fractional way, i.e. a^½ (SQR()), a^(¼) (fourth root) and so on. The problem actually goes deeper than this though. Converting a fractional power of 16 is only half the issue, and would still result in multiple decimal digits as a result, not the single digit sequence I was after. Let's consider the second hex digit as an example:

\[ \{ 1/256 .. 255/256 \} \notin \{ 1/16 .. 15/16 \} \]

1/256 = 3.90625E-3, 2/256 = 7.8125E-3 …

That's six digits and a decimal point, just for the first hex decimal.

(Post 169)

[TASP]

LOL, this topic makes my brain hurt.

Thanx !!!

[Valentin Albillo]

(02-06-2018 01:38 AM)Joe Horn Wrote:  There *is* a way of directly calculating the Nth digit of pi (without generating all the previous digits) but only in hex.

Nope. It also works for base 2 as well.

Quote: Nobody to date has discovered a way of doing it for pi in base 10.

It can be proved that there's no formula in base 10 for Pi like the one for bases 16 and 2. There are base-10 formulas for many other constants but not for Pi.

See A compendium of BBP-type formulas.

V.
.

[brickviking]

(02-07-2018 10:30 PM)TASP Wrote:  LOL, this topic makes my brain hurt.

Thanx !!!

And I thought I had it tough! I had to google for the set theory symbols for "not in a set", I haven't done this sort of thing since 1984! So far, the only way I've found to partially solve my problem is: generate all N digits of pi with PI50G from 3.14...., keeping the last x digits. I think that's the only feasible way to do it, and is infeasible after a certain number of N, at least on a 50G. It would be larger for calculators with more storage such as the Prime—assuming a PI50G equivalent on that platform, but would still eventually run into problems.

(Post 171)

[toml_12953]

(02-08-2018 02:59 AM)Valentin Albillo Wrote:  

(02-06-2018 01:38 AM)Joe Horn Wrote:  There *is* a way of directly calculating the Nth digit of pi (without generating all the previous digits) but only in hex.

Nope. It also works for base 2 as well.

Quote: Nobody to date has discovered a way of doing it for pi in base 10.

It can be proved that there's no formula in base 10 for Pi like the one for bases 16 and 2. There are base-10 formulas for many other constants but not for Pi.

See A compendium of BBP-type formulas.

V.
.

So what good is it in base 16? What use does it have? Can you calculate areas or volumes with it?

[toml_12953]

(02-08-2018 03:39 PM)Mike (Stgt) Wrote:  

(02-08-2018 01:02 PM)toml_12953 Wrote:  So what good is it in base 16? What use does it have? Can you calculate areas or volumes with it?

The goal behind this is a) hours of fun, b) get workaday tasks done in almost no time.

Fun? I just don't see fun in calculating something useless. Of course I'm an engineer, not a mathematician so anything theoretical has little or no interest for me.

[toml_12953]

(02-08-2018 06:58 PM)Mike (Stgt) Wrote:  Fun is something emotional, and emotions are known to lower the IQ fast and considerable to 50 or so. Strange enough, as for humans the drive to understand something is stronger than sex drive, it is said. So the mainly motivation for many scientists could be just the joy of "now I know", quite possible that an engineer is trained diversely, think and think different nothing but slogans, fun at work something useless. With hindsight worthles 'trouble and sorrow'.

Ciao.....Mike

Now when I see one of my circuits in a spacecraft perform perfectly, I do find that fun so I'm not criticizing how anyone finds fun but just saying that for me, calculating for it's own sake isn't fun. À chacun ses goûts!

[pier4r]

(02-08-2018 03:50 PM)toml_12953 Wrote:  Of course I'm an engineer, not a mathematician so anything theoretical has little or no interest for me.

I have no intention to be polemic, anyway the above sentence sounds to me like "I'm a driver I don't need to know how the engine in my car or drivetrain in my bike works." And while that is completely legit, knowing the tools that one uses a bit better is nothing bad, rather it should improve the results at the end.

[Gerson W. Barbosa]

(02-08-2018 12:26 PM)Mike (Stgt) Wrote:  

(02-08-2018 04:22 AM)brickviking Wrote:  ... on that platform, but would still eventually run into problems.

A little 'look ahead' on another platform (slow laptop using REXX, an interpreter language) showing the last 30 digits, time elapsed in seconds, iterations, and method.
requested 500 digits:
...673518857527248912279381830118358 4.056000 1643 Newton
...673518857527248912279381830119702 1.248000 1034 Madhava-Leibniz

Madhava-Leibniz is interesting, especially because it's from the XIVth century. And yet we can use it on the HP-50g to compute 500 digits, if we are not in a hurry:

500

« PUSH RAD -105 CF -3 CF DUP .653 * 1.74 + IP R->I DUP 2 MOD + DUP 4 * OVER DUPDUP 1 - 1
FOR i i SQ SWAP / PICK3 + ROT SWAP -1
STEP INV NIP UNROT + 1 - 3 0 UNROT
FOR i i INV i 2 - INV - + -4
STEP - 4 * EXPAND FXND DUP SIZE R->I ALOG OVER - PICK3 * SWAP IQUOT + ->STR
DUP HEAD -51 FC? { "." } { "," } IFTE + SWAP TAIL + 1 ROT 2 + SUB POP
»

TEVAL -->

2: "3.
14159265358979323846264338327950288419716939937510582097494459230781640628620899​86280348253421170679
82148086513282306647093844609550582231725359408128481117450284102701938521105559​64462294895493038196
44288109756659334461284756482337867831652712019091456485669234603486104543266482​13393607260249141273
72458700660631558817488152092096282925409171536436789259036001133053054882046652​13841469519415116094
33057270365759591953092186117381932611793105118548074462379962749567351885752724​89122793818301194912"
1:       s:21533.1807


21533 seconds on my HP-50g (almost 6 hours!)

This is much longer than it should be, as I am using no external long-float library (only the built-in Long Int in UserRPL). By using an arbitrary precision library, the running time might be significantly decreased.

Quoting from the Wikipedia article on Madhava of Sangamagrama:

Quote:This text attributes most of the expansions to Madhava, and gives the following infinite series expansion of π, now known as the Madhava-Leibniz series:

\[{\displaystyle {\frac {\pi }{4}}=1-{\frac {1}{3}}+{\frac {1}{5}}-{\frac {1}{7}}+\cdots +{\frac {(-1)^{n}}{2n+1}}+\cdots } {\frac {\pi }{4}}=1-{\frac {1}{3}}+{\frac {1}{5}}-{\frac {1}{7}}+\cdots +{\frac {(-1)^{n}}{2n+1}}+\cdots\]

which he obtained from the power series expansion of the arc-tangent function. However, what is most impressive is that he also gave a correction term, Rn, for the error after computing the sum up to n terms. Madhava gave three expressions for the correction term Rn,[4] namely

Rn = 1/(4n), or
Rn = n/(4n^2+ 1), or
Rn = (n^2 + 1)/(4n^3 + 5n).

where the third correction leads to highly accurate computations of π.

It is not clear how Madhava might have found these correction terms.

Perhaps that Wikipedia article should be updated. The next correction terms are

Rn = (4n^3 + 13n)/(16n^4 + 56n^2 +9),
Rn = (4n^4 + 29n^2 + 16)/(16n^5 + 120n^3 + 89n) and
Rn = (16n^5 + 216n^3 + 389n)/(64n^6 + 880n^4 + 1756n^2 + 225)

Madhava's correction terms are simply the first three approximants to the continued fraction in the following expression, from an old thread:

\[\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+...​​ }}}}}\]

From observation, 2*d/3 terms of the series and the corresponding correction term (or the evaluation of the same number of terms of the continued fraction) are necessary to obtain d correct decimal places of π.

[toml_12953]

(02-09-2018 02:01 AM)Mike (Stgt) Wrote:  

(02-08-2018 11:48 PM)Gerson W. Barbosa Wrote:  Madhava-Leibniz is interesting, [...]

\[{\displaystyle {\frac {\pi }{4}}=1-{\frac {1}{3}}+{\frac {1}{5}}-{\frac {1}{7}}+\cdots +{\frac {(-1)^{n}}{2n+1}}+\cdots } \]

Seems I confused by so many formulas. What you quote is the Madhava-Leibniz formula. In contrast what I use is a more rapidly converging series also from Madhava.

Code:

s = 1
n = 1
DO m = 3 by 2 UNTIL m = m + n
   n = n / -3
   s = s + n / m
end
Pi = 2 * s * Sqrt(3)

This needs about 2 iterations per correct digit of Pi. In my list I marked it falsely with 'Madhava-Leibniz'.

To reduce the loop overhead I merged two steps (what results in a remarkable gain on my system):

Code:

s = 1
n = 1
DO m = 3 by 4 UNTIL m = m + n
   n = n / 9
   s = s + n / (m + 2) - 3 * n / m
end
Pi = 2 * s * Sqrt(3)

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.

Ciao.....Mike

You're still limited to the precision of a variable (Pi), right? How would you go beyond 12 or 13 digits? If you can't, why go to all that trouble to get a number you can get with a single command. That's not a rhetorical question, I'd really like to know.

[toml_12953]

(02-09-2018 12:08 PM)Mike (Stgt) Wrote:  

(02-09-2018 02:42 AM)toml_12953 Wrote:  You're still limited to the precision of a variable (Pi), right? How would you go beyond 12 or 13 digits?

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

Ah! That explains it. Please excuse my obtuseness.

[toml_12953]

(02-09-2018 12:08 PM)Mike (Stgt) Wrote:  

(02-09-2018 02:42 AM)toml_12953 Wrote:  You're still limited to the precision of a variable (Pi), right? How would you go beyond 12 or 13 digits?

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

[Massimo Gnerucci]

(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

[pier4r]

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.

[toml_12953]

(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)

[Gerson W. Barbosa]

(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).

[Gerson W. Barbosa]

(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.