
Computing pi with bc - lelf
https://www.johndcook.com/blog/2019/10/29/computing-pi-with-bc/
======
imglorp
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
It seems hard to beat bc to compute 10k digits of π.

Can you do that with a one-liner?

~~~
yesenadam
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
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
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.

~~~
espadrine
Nice!

It is interesting how Chudnovsky is still faster (112sec on my machine, with
parameters 706 10000) than Gauss-Legendre (176s with parameters 15 10005)
despite the slower convergence.

I found it too bad that they still were that much slower than bc’s Machin
(40s) considering how much longer the code was.

I tried it with Unbounded Spigot (Jeremy Gibbons) and got 42s, similar to bc’s
performance, and 2.6× faster than Chudnovsky:

    
    
        # dc -f spigot.dc -e '10000 lPx'
        1sq180sr60st2si
        [27li*12-lq*5lr*+5lt*/sylyn]sY
        [3li*1+3*3li*2+*su]sU
        [5li*2-lq*lr+lylt*-lu*10*sr]sR
        [2li*1-li*lq*10*sq]sQ
        [lult*st]sT [1li+si]sI
        [lYxlUxlRxlQxlTxlIxlj1-sj0lj!<L]sL
        [sjlLx[]p]sP
    

Now it is competitive with bc (on performance, not on code size). However,
compared to Node.js (2.7s), it gets crushed:

    
    
        const pi = s => {
          let y = (s.q*(27n*s.i-12n)+5n*s.r)/(5n*s.t);
          let u = 3n*(3n*s.i+1n)*(3n*s.i+2n);
          s.r = 10n*u*(s.q*(5n*s.i-2n)+s.r-y*s.t);
          s.q = 10n*s.q*s.i*(2n*s.i-1n);
          s.t = s.t*u;
          s.i = s.i+1n;
          return y;
        };
    
        let s = {q:1n, r:180n, t:60n, i:2n};
        for (let i = 0; i < 1e4 + 1; i++) {
          process.stdout.write(''+pi(s));
        }
        console.log();

------
jolmg
> 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...](http://blog.recursiveprocess.com/2015/03/12/using-the-arctan-power-
series-to-calculate-pi/)

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

------
lelf
> _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
> 688,668,864 bytes

That's an eye popping number for what should be a pretty simple calculation.

~~~
lelf
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
Yep! Most of the allocations die _very_ young.

------
DannyB2
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
That's not /calculating/ pi though, just returning a constant.

""" maxima-5.43.0/src/mlisp.lisp:

(mdefprop $%pi
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608
$numer)

"""

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

