4 ms·

Computing pi with bc

80 points by lelf 7y ago

jolmg 7y ago
> I imagine bc is using some sort of power series to compute arctan I didn't know or didn't remember that pi could be calculated with a power series. I thought, bc is probably using a constant estimation of pi and this is going to fail, but sure enough, bc/libmath.b:273 in bc's source is the location of the imagined power series in a()'s definition. I found this page explaining the formula: http://blog.recursiveprocess.com/2015/03/12/using-the-arctan-power-series-to-calculate-pi/ http://blog.recursiveprocess.com/2015/03/12/using-the-arctan...
ncmncm 7y ago
It is probably better to use "export BC_ENV_ARGS=-lq" than "alias bc='bc -lq'".
imglorp 7y ago
Fun fact: In the old days, as was unix tradition, bc was a thin wrapper on top of dc (desktop calculator) which is even more spartan than bc. It's quick and tiny, you can pipe things to it. It syntax is postfix, unlike bc which is infix. It's stack based with conditionals, so you can write full programs in dc if you have some free time. It's also got arbitrary input and output radixes. An elegant weapon for a more civilized OS! Eg: echo 16o 5 5 5 ^^p | dc https://linux.die.net/man/1/dc https://linux.die.net/man/1/dc PS: (insane?) people have written some big stuff in dc (because it's there?). I'm drawing a blank if it was Zork or Lisp, hard search to do... PPS: https://esolangs.org/wiki/Dc https://esolangs.org/wiki/Dc
- espadrine 7y ago
  It seems hard to beat bc to compute 10k digits of π. Can you do that with a one-liner?
  - imglorp 7y ago
    Tempting... must sit on hands... But yes, I think you could. The Chudnovsky algorithm is only a half page of Python. dc can be formatted nicelyish with comments but you can also golf it obviously.
  - yesenadam 7y ago
    I'd never used dc before today, or knew anything about it. Here's a first attempt at 1000 digits, using the first algorithm I thought of using[0]: echo 1020k 1sa 1 2v/sb .25st 1sp [lalb+2/sA]s1 [lalb*vsB]s2 [ltlalA-2^lp*-sT]s3 [lp2*sP]s4 [l1xl2xl3xl4xlAsalBsblPsplTst]s7 [l7x7l7xl7xl7x]s8 [l8xl8xl8xl8x]s9 [l9xl9xl9x]sq lqxlqx lAlB+2^lT4*/ p| dc This seems to get about 1017 correct digits. The macros stored with s1, s2, s3, s4 are the 4 equations for a_n, b_n, t_n, p_n. (a_n is put in register a, a_{n+1} in register A etc.) Undoubtedly some other algorithm would be much better. But thanks, I had fun doing this! The man page was unusually easy to read. The program says: Use 1020 digit accuracy. Store 1 in register a. Store 1/sqrt(2) in register b. Store 0.25 in register t. Store 1 in register p. Store the macro "Add register a and b, divide by 2 and store in register A" in register 1. etc. l1x = execute the macro stored in register 1, v = square root, p = print top of stack. [0] https://en.wikipedia.org/wiki/Gauss–Legendre_algorithm https://en.wikipedia.org/wiki/Gauss–Legendre_algorithm
    - yesenadam 7y ago
      Whoops! [l7x7l7xl7xl7x]s8 should be [l7xl7xl7xl7x]s8 EDIT: Gee, that algorithm is way more accurate than I remembered! "rapidly convergent, with only 25 iterations producing 45 million correct digits"! "the number of correct digits doubles with each iteration". This shorter program does 10,000 correct digits from 16 iterations, just takes a few minutes on my (old) machine: echo 10010k 1sa 1 2v/sb .25st 1sp [lalb+2/sA]s1 [lalb*vsB]s2 [ltlalA-2^lp*-sT]s3 [lp2*sP]s4 [l1xl2xl3xl4xlAsalBsblPsplTst]s7 [l7xl7xl7xl7x]s8 l8xl8xl8xl8x lAlB+2^lT4*/ p| dc I checked the result here http://www.eveandersson.com/pi/digits/1000000 http://www.eveandersson.com/pi/digits/1000000 Also - I was curious how old dc is. Old: "dc is the oldest surviving Unix language. When its home Bell Labs received a PDP-11, dc—written in B—was the first language to run on the new computer, even before an assembler. Ken Thompson has opined that dc was the very first program written on the machine." https://en.wikipedia.org/wiki/Dc_(computer_program) https://en.wikipedia.org/wiki/Dc_(computer_program)
      - xelxebar 7y ago
        Hey, I just coded this up and think you might get a kick out of it: [[lipljplkpllplxplmplsp]sD]sT # enable debug tracing [426880 10005v*sc]sC 0si [li1+si]sI 6sk [lk12+sk]sK 13591409sl [ll545140134+sl]sL 1sx [lx_262537412640768000*sx]sX 1sm [lk3^16lk*-li1+3^/lm*sm]sM [lmll*lx/ls+ss]sS [lSxlMxlXxlLxlKxlIx lDx ljli<Z]sZ [ksjlCxlZx lcls/]sP # ( iterations precision -- pi ) It's basically a direct implementation of the Chudnovsky algorithm. Stick the above code in some file, say chud.dc, and you can run it from the command line like this: $ dc -f chud.dc -e '10 100 lPxp' $ dc -f chud.dc -e 'lTx 10 100 lPxp' # outputs tracing info The salient macro is P (for "pi" of course). The comment documents its usage, showing how it transforms the stack. Items on the left of "--" are what it expects; items on the right are its outputs. Unfortunately, this series seems to converge more slowly than Gauss-Legendre.
lelf 7y ago
> I imagine bc is using some sort of power series There is also a cool Haskell package, numbers, that has all arithmetics defined via power series: $ cabal v2-repl -b numbers ⡇ λ> import Data.Number.CReal λ> cos pi + 1 :: CReal 0.0 λ> :set +s λ> showCReal 10000 pi "3.141592653589793238462643383279502884197169399375105820… (2.69 secs, 688,668,864 bytes)
- jandrese 7y ago
  > 688,668,864 bytes That's an eye popping number for what should be a pretty simple calculation.
  - lelf 7y ago
    Welcome to the world of functional programming! This is how much memory the garbage collector has freed during the computation. It does not mean REPL actually used that much at any given moment.
    - giornogiovanna 7y ago
      Yep! Most of the allocations die very young.
DannyB2 7y ago
Use maxima: > maxima (%i1) bfloat(%pi),fpprec:1000; (%o1) 3.1415926535897932384626433832795028841971693993751058209749445923078164\ 062862089986280348253421170679821480865132823066470938446095505822317253594081\ 284811174502841027019385211055596446229489549303819644288109756659334461284756\ 482337867831652712019091456485669234603486104543266482133936072602491412737245\ 870066063155881748815209209628292540917153643678925903600113305305488204665213\ 841469519415116094330572703657595919530921861173819326117931051185480744623799\ 627495673518857527248912279381830119491298336733624406566430860213949463952247\ 371907021798609437027705392171762931767523846748184676694051320005681271452635\ 608277857713427577896091736371787214684409012249534301465495853710507922796892\ 589235420199561121290219608640344181598136297747713099605187072113499999983729\ 780499510597317328160963185950244594553469083026425223082533446850352619311881\ 710100031378387528865875332083814206171776691473035982534904287554687311595628\ 63882353787593751957781857780532171226806613001927876611195909216420199b0 It's done in the blink of an eye. Or use WxMaxima if you need a GUI. Console: (verb) to comfort someone in time of grief because they are forced to use the command line.
- kimixa 7y ago
  That's not /calculating/ pi though, just returning a constant. """ maxima-5.43.0/src/mlisp.lisp: (mdefprop $%pi 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608 $numer) """
saagarjha 7y ago
> I imagine bc is using some sort of power series to compute arctan It uses the Taylor series for the function, as many calculators do.

80 points by lelf 7y ago

> I imagine bc is using some sort of power series to compute arctan I didn't know or didn't remember that pi could be calculated with a power series. I thought, bc is probably using a constant estimation of pi and this is going to fail, but sure enough, bc/libmath.b:273 in bc's source is the location of the imagined power series in a()'s definition. I found this page explaining the formula: http://blog.recursiveprocess.com/2015/03/12/using-the-arctan-power-series-to-calculate-pi/ http://blog.recursiveprocess.com/2015/03/12/using-the-arctan...

ncmncm 7y ago

It is probably better to use "export BC_ENV_ARGS=-lq" than "alias bc='bc -lq'".

imglorp 7y ago

Fun fact: In the old days, as was unix tradition, bc was a thin wrapper on top of dc (desktop calculator) which is even more spartan than bc. It's quick and tiny, you can pipe things to it. It syntax is postfix, unlike bc which is infix. It's stack based with conditionals, so you can write full programs in dc if you have some free time. It's also got arbitrary input and output radixes. An elegant weapon for a more civilized OS! Eg: echo 16o 5 5 5 ^^p | dc https://linux.die.net/man/1/dc https://linux.die.net/man/1/dc PS: (insane?) people have written some big stuff in dc (because it's there?). I'm drawing a blank if it was Zork or Lisp, hard search to do... PPS: https://esolangs.org/wiki/Dc https://esolangs.org/wiki/Dc

espadrine 7y ago
It seems hard to beat bc to compute 10k digits of π. Can you do that with a one-liner?
- imglorp 7y ago
  Tempting... must sit on hands... But yes, I think you could. The Chudnovsky algorithm is only a half page of Python. dc can be formatted nicelyish with comments but you can also golf it obviously.
- yesenadam 7y ago
  I'd never used dc before today, or knew anything about it. Here's a first attempt at 1000 digits, using the first algorithm I thought of using[0]: echo 1020k 1sa 1 2v/sb .25st 1sp [lalb+2/sA]s1 [lalb*vsB]s2 [ltlalA-2^lp*-sT]s3 [lp2*sP]s4 [l1xl2xl3xl4xlAsalBsblPsplTst]s7 [l7x7l7xl7xl7x]s8 [l8xl8xl8xl8x]s9 [l9xl9xl9x]sq lqxlqx lAlB+2^lT4*/ p| dc This seems to get about 1017 correct digits. The macros stored with s1, s2, s3, s4 are the 4 equations for a_n, b_n, t_n, p_n. (a_n is put in register a, a_{n+1} in register A etc.) Undoubtedly some other algorithm would be much better. But thanks, I had fun doing this! The man page was unusually easy to read. The program says: Use 1020 digit accuracy. Store 1 in register a. Store 1/sqrt(2) in register b. Store 0.25 in register t. Store 1 in register p. Store the macro "Add register a and b, divide by 2 and store in register A" in register 1. etc. l1x = execute the macro stored in register 1, v = square root, p = print top of stack. [0] https://en.wikipedia.org/wiki/Gauss–Legendre_algorithm https://en.wikipedia.org/wiki/Gauss–Legendre_algorithm
  - yesenadam 7y ago
    Whoops! [l7x7l7xl7xl7x]s8 should be [l7xl7xl7xl7x]s8 EDIT: Gee, that algorithm is way more accurate than I remembered! "rapidly convergent, with only 25 iterations producing 45 million correct digits"! "the number of correct digits doubles with each iteration". This shorter program does 10,000 correct digits from 16 iterations, just takes a few minutes on my (old) machine: echo 10010k 1sa 1 2v/sb .25st 1sp [lalb+2/sA]s1 [lalb*vsB]s2 [ltlalA-2^lp*-sT]s3 [lp2*sP]s4 [l1xl2xl3xl4xlAsalBsblPsplTst]s7 [l7xl7xl7xl7x]s8 l8xl8xl8xl8x lAlB+2^lT4*/ p| dc I checked the result here http://www.eveandersson.com/pi/digits/1000000 http://www.eveandersson.com/pi/digits/1000000 Also - I was curious how old dc is. Old: "dc is the oldest surviving Unix language. When its home Bell Labs received a PDP-11, dc—written in B—was the first language to run on the new computer, even before an assembler. Ken Thompson has opined that dc was the very first program written on the machine." https://en.wikipedia.org/wiki/Dc_(computer_program) https://en.wikipedia.org/wiki/Dc_(computer_program)
    - xelxebar 7y ago
      Hey, I just coded this up and think you might get a kick out of it: [[lipljplkpllplxplmplsp]sD]sT # enable debug tracing [426880 10005v*sc]sC 0si [li1+si]sI 6sk [lk12+sk]sK 13591409sl [ll545140134+sl]sL 1sx [lx_262537412640768000*sx]sX 1sm [lk3^16lk*-li1+3^/lm*sm]sM [lmll*lx/ls+ss]sS [lSxlMxlXxlLxlKxlIx lDx ljli<Z]sZ [ksjlCxlZx lcls/]sP # ( iterations precision -- pi ) It's basically a direct implementation of the Chudnovsky algorithm. Stick the above code in some file, say chud.dc, and you can run it from the command line like this: $ dc -f chud.dc -e '10 100 lPxp' $ dc -f chud.dc -e 'lTx 10 100 lPxp' # outputs tracing info The salient macro is P (for "pi" of course). The comment documents its usage, showing how it transforms the stack. Items on the left of "--" are what it expects; items on the right are its outputs. Unfortunately, this series seems to converge more slowly than Gauss-Legendre.

lelf 7y ago

> I imagine bc is using some sort of power series There is also a cool Haskell package, numbers, that has all arithmetics defined via power series: $ cabal v2-repl -b numbers ⡇ λ> import Data.Number.CReal λ> cos pi + 1 :: CReal 0.0 λ> :set +s λ> showCReal 10000 pi "3.141592653589793238462643383279502884197169399375105820… (2.69 secs, 688,668,864 bytes)

jandrese 7y ago
> 688,668,864 bytes That's an eye popping number for what should be a pretty simple calculation.
- lelf 7y ago
  Welcome to the world of functional programming! This is how much memory the garbage collector has freed during the computation. It does not mean REPL actually used that much at any given moment.
  - giornogiovanna 7y ago
    Yep! Most of the allocations die very young.

DannyB2 7y ago

Use maxima: > maxima (%i1) bfloat(%pi),fpprec:1000; (%o1) 3.1415926535897932384626433832795028841971693993751058209749445923078164\ 062862089986280348253421170679821480865132823066470938446095505822317253594081\ 284811174502841027019385211055596446229489549303819644288109756659334461284756\ 482337867831652712019091456485669234603486104543266482133936072602491412737245\ 870066063155881748815209209628292540917153643678925903600113305305488204665213\ 841469519415116094330572703657595919530921861173819326117931051185480744623799\ 627495673518857527248912279381830119491298336733624406566430860213949463952247\ 371907021798609437027705392171762931767523846748184676694051320005681271452635\ 608277857713427577896091736371787214684409012249534301465495853710507922796892\ 589235420199561121290219608640344181598136297747713099605187072113499999983729\ 780499510597317328160963185950244594553469083026425223082533446850352619311881\ 710100031378387528865875332083814206171776691473035982534904287554687311595628\ 63882353787593751957781857780532171226806613001927876611195909216420199b0 It's done in the blink of an eye. Or use WxMaxima if you need a GUI. Console: (verb) to comfort someone in time of grief because they are forced to use the command line.

kimixa 7y ago
That's not /calculating/ pi though, just returning a constant. """ maxima-5.43.0/src/mlisp.lisp: (mdefprop $%pi 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608 $numer) """

saagarjha 7y ago

> I imagine bc is using some sort of power series to compute arctan It uses the Taylor series for the function, as many calculators do.