Fast Pi Calculator

This board is to show, discuss and archive useful combinator- and logic-creations.
Smart triggering, counters and sensors, useful circuitry, switching as an art :), computers.
Please provide if possible always a blueprint of your creation.
Post Reply
Copy link
User avatar
DaveMcW
Smart Inserter
Smart Inserter
Posts: 3716
Joined: Tue May 13, 2014 11:06 am
Contact:
Contact DaveMcW

Fast Pi Calculator

Post by DaveMcW »

I built a fast calculator for the digits of pi. See it in action here, the video plays at real speed:

Image


Blueprint



Steps to build:

1. Place the calculator
pi1.jpg
pi1.jpg (45.94 KiB) Viewed 8434 times

2. Place the digits, with the 4 constant combinators overlapping. You can keep placing digits until you run out of map space or computing power. *
pi2.jpg
pi2.jpg (87.14 KiB) Viewed 8434 times

3. Turn on the far left combinator labeled "GO".
pi3.jpg
pi3.jpg (39.36 KiB) Viewed 8434 times

You have to place all the digits before starting the calculation, it does not have a reset feature.
Last edited by DaveMcW on Fri Oct 04, 2019 12:25 am, edited 1 time in total.
Top
User avatar
DaveMcW
Smart Inserter
Smart Inserter
Posts: 3716
Joined: Tue May 13, 2014 11:06 am
Contact:
Contact DaveMcW

Re: Fast Pi Calculator

Post by DaveMcW »

Formula

I use Machin's formula, which was used in the world record pi calculation for 250 years. It converges quickly (1.398 digits per term), requires no multiplication, and reuses calculations from previous terms.

machin.png
machin.png (18.12 KiB) Viewed 8424 times


Algorithm

Here is the algorithm I used. Implementing the bigint functions in your favorite language is left as an exercise for the reader.

Code: Select all

pi = bigint(0)
power1 = bigint(16 * 5)
power2 = bigint(4 * 239)
n = 1

for i = 0, DIGITS / 1.398 do
  -- Divide by x^2
  bigint_divide(power1, 5^2)
  bigint_divide(power2, 239^2)

  -- Divide by odd number
  n = n + 2
  term1 = bigint_copy(power1)
  term2 = bigint_copy(power2)
  bigint_divide(term1, n)
  bigint_divide(term2, n)

  -- Adjust sign
  if i % 2 == 1 then
    bigint_negative(term1)
  else
    bigint_negative(term2)
  end

  -- Add to total
  bigint_add(pi, term1)
  bigint_add(pi, term2)
end

Big Integers

I implemented big integers using a word size of 10000. There were several reasons to pick this value.
  • 10000 is a multiple of 10. This makes converting to decimal digits easy.
  • 10000 * 239^2 < 2^31. This means I don't have to worry about int32 overflow until term 107375 of the series.
  • Word size 10000 is 78% faster than word size 1000.
Each group of 4 digits has its own bigint hardware, which means expanding the display automatically expands the bigint size.

pi5.jpg
pi5.jpg (376.28 KiB) Viewed 8424 times


BigInt Division

BigInt division requires 3 combinators per word.
Combinator 1 calculates: word / denominator -> word
Combinator 2 calculates: word % denominator -> remainder
Combinator 3 calculates: remainder * 10000 -> carry
The wires calculate for free: next_word + carry -> next_word

The minimum time per word is 2 ticks, since combinators 2 and 3 must run in series.

With a bigint of 500 words, a division operation takes 1000 ticks. But this can be parallelized! I feed another division calculation into the circuit each tick, resulting in an average runtime of 1 tick per division. To perform the 4 divisions in each term of the series, I use 2 division circuits and 2 ticks.


BigInt Addition

Addition is supposed to be performed right to left, to calculate the carry correctly. But division is performed left to right. I don't want to delay my parallel division scheme, so I simply add the words together and ignore the carry.

The wires calculate for free: pi_word + term_word -> pi_word

This results in words outside the range [0,9999]. I then use 5 combinators to rebalance the words.
Combinator 1 is a timer that only runs the balancer once per 3 ticks.
Combinator 2 calculates: word > 9999 -> carry = 1
Combinator 3 calculates: word > 9999 -> word -= 10000
Combinator 4 calculates: word < 0 -> carry = -1
Combinator 5 calculates: word < 0 -> word += 10000

If you watch carefully, you can spot the corrupted negative digits of pi before the balancer circuit fixes them.


Performance

It calculates 42 digits per Factorio second. Converted to real time, this is 42 * m / n digits per second, where n is the number of digits and m is the maximum number of digits your computer can handle before dropping below 60 ups.

* After 150,000 digits the signed integer math will overflow and fail. This can be extended to 1.5 million digits with a small modification (reducing the number of digits per word to 3).
Last edited by DaveMcW on Fri Oct 04, 2019 6:01 pm, edited 1 time in total.
Top
xng
Fast Inserter
Fast Inserter
Posts: 165
Joined: Fri Feb 14, 2014 1:04 pm
Contact:
Contact xng

Re: Fast Pi Calculator

Post by xng »

Incredible!

I don't get one thing though, it says it doesn't need any multiplication..
DaveMcW wrote: ↑Thu Oct 03, 2019 11:54 pm I use Machin's formula, which was used in the world record pi calculation for 250 years. It converges quickly (1.398 digits per term), requires no multiplication, and reuses calculations from previous terms.
But the formula has serious amounts of multiplication going on (both normal multiplication and implicit multiplication from power of base), do you mean you don't have to do any of those?

I am very impressed by this solution! <3
Top
User avatar
darkfrei
Smart Inserter
Smart Inserter
Posts: 2905
Joined: Thu Nov 20, 2014 11:11 pm
Contact:
Contact darkfrei

Re: Fast Pi Calculator

Post by darkfrei »

Why not Bailey–Borwein–Plouffe formula Bellard's formula?
Last edited by darkfrei on Fri Oct 04, 2019 11:29 am, edited 1 time in total.
Top
User avatar
DaveMcW
Smart Inserter
Smart Inserter
Posts: 3716
Joined: Tue May 13, 2014 11:06 am
Contact:
Contact DaveMcW

Re: Fast Pi Calculator

Post by DaveMcW »

All the multiplication is in the denominator, so it is really division. I use constants 25 and 57121 instead of calculating x^2 each time.

BBP Formula does not work in decimal, and is inefficient if you want all the digits. But you can see it here: viewtopic.php?f=193&t=76353
Top
Intangir_V
Inserter
Inserter
Posts: 22
Joined: Wed Oct 02, 2019 3:50 pm
Contact:
Contact Intangir_V

Re: Fast Pi Calculator

Post by Intangir_V »

bravo :shock:
Top
User avatar
DaveMcW
Smart Inserter
Smart Inserter
Posts: 3716
Joined: Tue May 13, 2014 11:06 am
Contact:
Contact DaveMcW

Re: Fast Pi Calculator

Post by DaveMcW »

I spent the past few days studying Machin-like formulas. There are some very fast converging formulas published, but they prefer big numbers and I require all terms to be below 463 to avoid int32 overflow. I finally built a custom program to search for formulas that fit my needs.

I discovered some interesting ones:

Code: Select all

pi/4 = 12*arctan(1/191) + 39*arctan(1/239) + 68*arctan(1/268) + 32*arctan(1/302) + 20*arctan(1/307) + 44*arctan(1/327)
pi/4 = 29*arctan(1/268) + 51*arctan(1/278) + 44*arctan(1/302) + 20*arctan(1/307) + 17*arctan(1/327) + 39*arctan(1/378) + 27*arctan(1/401) + 27*arctan(1/447)
The 8-term formula converges 6.4% faster than the 6-term formula, but requires 7.4% more combinators. So I went with the 6-term formula.
Expanded formula

Here is the updated blueprint using the 6-term formula. It generates 77 digits per second and can handle 480,000 digits before overflowing.

Top
movax20h
Fast Inserter
Fast Inserter
Posts: 164
Joined: Fri Mar 08, 2019 7:07 pm
Contact:
Contact movax20h

Re: Fast Pi Calculator

Post by movax20h »

Word size 10000 is 78% faster than word size 1000.
You mean 1.78 times faster? Or taking 0.78 of the time compared to 1000.
Top
movax20h
Fast Inserter
Fast Inserter
Posts: 164
Joined: Fri Mar 08, 2019 7:07 pm
Contact:
Contact movax20h

Re: Fast Pi Calculator

Post by movax20h »

DaveMcW wrote: ↑Thu Oct 10, 2019 3:32 am I spent the past few days studying Machin-like formulas. There are some very fast converging formulas published, but they prefer big numbers and I require all terms to be below 463 to avoid int32 overflow. I finally built a custom program to search for formulas that fit my needs.

I discovered some interesting ones:

Code: Select all

pi/4 = 12*arctan(1/191) + 39*arctan(1/239) + 68*arctan(1/268) + 32*arctan(1/302) + 20*arctan(1/307) + 44*arctan(1/327)
pi/4 = 29*arctan(1/268) + 51*arctan(1/278) + 44*arctan(1/302) + 20*arctan(1/307) + 17*arctan(1/327) + 39*arctan(1/378) + 27*arctan(1/401) + 27*arctan(1/447)
The 8-term formula converges 6.4% faster than the 6-term formula, but requires 7.4% more combinators. So I went with the 6-term formula.
Expanded formula

Here is the updated blueprint using the 6-term formula. It generates 77 digits per second and can handle 480,000 digits before overflowing.

FYI. In python, you can generate permutations and combinations using a bit simpler approach that fixed nested loops:

Code: Select all

# Generate all 5-term combinations of primes
for x1 in range(0, len(primes)):
    for x2 in range(x1 + 1, len(primes)):
        for x3 in range(x2 + 1, len(primes)):
            for x4 in range(x3 + 1, len(primes)):
                for x5 in range(x4 + 1, len(primes)):
                    primes_subset = {2: True, primes[x1]: True, primes[x2]: True,
                                     primes[x3]: True, primes[x4]: True, primes[x5]: True}
can be replaced by:

Code: Select all

for prime in itertools.combinations(primes, 5):
    primes_subset = set(prime)
    primes_subset[2] = True
    ...
It is not optimal, because it will return same primes sometimes just in different order, but this can be fixed easily with a bit smarter algorithm I show a sketch at the bottom.

Similarly:

Code: Select all

                    # Generate all 6-term combinations of squares
                    for y1 in range(0, len(squares)):
                        for y2 in range(y1 + 1, len(squares)):
                            for y3 in range(y2 + 1, len(squares)):
                                for y4 in range(y3 + 1, len(squares)):
                                    for y5 in range(y4 + 1, len(squares)):
                                        for y6 in range(y5 + 1, len(squares)):
                                            squares_subset = [squares[y1], squares[y2], squares[y3], squares[y4],
                                                              squares[y5], squares[y6]]
can be replaced by:

Code: Select all

for squares_subset in itertools.combinations(squares, 6):
  ...

You can also implement similar concepts using recursive yield generators using yield from, for example, this is a code I was using to generate all permutations for assembler inputs in my program to optimize layout of factories and belts:

Code: Select all

def generate_all_multi_permutations(possibilities, current_prefix=[]):
  """Argument possibilities is a list of lists, i.e. [[0, 1, 2], [0, 1], [0, 1, 2, 3]].
  This function will return a similar list of lists, with the set of results being
  a Cartesian product of all possible permutation of each sublist.
  Example output:
  [[0, 1, 2], [0, 1], [0, 1, 2, 3]].
  [[0, 1, 2], [0, 1], [0, 1, 3, 2]].
  [[0, 1, 2], [0, 1], [0, 2, 3, 1]].
  [[0, 1, 2], [0, 1], [0, 2, 1, 3]].
  ...
  [[2, 0, 1], [1, 0], [3, 1, 0, 2]].
  ...
  [[2, 1, 0], [1, 0], [3, 2, 1, 0]].
  
  The number of returned (iterated) elements is product of factorials.
  Example from above: 3! * 2! * 4! = 288
  
  For long lists (6+ elements), or where each sublists is long (4+ elements),
  this is will grow fast!
  """
  if not possibilities:
    yield current_prefix
  else:
    for p in permutations(possibilities[0]):
      yield from generate_all_multi_permutations(possibilities[1:], current_prefix=current_prefix + [p])
Your code could be something like this:

Code: Select all

# Generate all unique monotonic sequences (combinations) from L.
def generate_all_sequences(L, length, current_prefix=[]):
  if length == 0:
    yield current_prefix
  else:
    for i, x in enumerate(L):
      yield from generate_all_sequences(L[i+1:], length=length-1, current_prefix=current_prefix + [x])
Example output:

Code: Select all

>>> for seq in generate_all_sequences([2, 3, 5, 7, 11, 13], length=4): print(seq)
... 
[2, 3, 5, 7]
[2, 3, 5, 11]
[2, 3, 5, 13]
[2, 3, 7, 11]
[2, 3, 7, 13]
[2, 3, 11, 13]
[2, 5, 7, 11]
[2, 5, 7, 13]
[2, 5, 11, 13]
[2, 7, 11, 13]
[3, 5, 7, 11]
[3, 5, 7, 13]
[3, 5, 11, 13]
[3, 7, 11, 13]
[5, 7, 11, 13]
Top
Post Reply

Return to “Combinator Creations”