Hacker News new | past | comments | ask | show | jobs | submit login

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!




Guidelines | FAQ | Lists | API | Security | Legal | Apply to YC | Contact

Search: