Computing pi with bc | Hacker News Reader

Computing pi with bc

johndcook.com

80 pointslelf17 comments

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

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

It seems hard to beat bc to compute 10k digits of π.

Can you do that with a one-liner?

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.

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

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

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)

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.

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

And here's my version of Gauss-Legendre:

    [[lmp lnp lap lbp ltp lqp lOxp]sD]sT # enable debug tracing

    [0sn [ln 1+ sn]sN
     1sa [la dsc lb + 2/ sa]sA
     1 2v/sb [lc lb *v sb]sB
     1 4/st [lt lc la- 2^ lq* - st]sT
     1sq [lq 2* sq]sQ
    ]sI
    [lDx lNx lAx lBx lTx lQx lnlm>Z]sZ
    [la lb+ 2^ lt 4* /]sO

    [Ksk k sm lIx lZx lOx lkk]sP # ( iterations precision -- pi )

Usage is the same as with the Chudnovsky version from my other comment. Assuming the above script is located in gl.dc:

    $ dc -f gl.dc -e '15 10010 lPxp'

That gives about 10k digits accurately.

It's pretty much just a the implementation one would expect. The only tricky part is ordering the variable updates, since those have interdependencies. Both 'a' and 'b' depend ond each other, so we're force to create an intermediate variable, that's what 'c' is in the above.

Anyway!

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

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

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

> 688,668,864 bytes

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

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.

Yep! Most of the allocations die very young.

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.

That's not /calculating/ pi though, just returning a constant.

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

(mdefprop $%pi 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608 $numer)

"""

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

Reply on news.ycombinator.com