stylewarning / hypergeometrica Goto Github PK

Livin' like it's 1813 (or 1988).

License: BSD 3-Clause "New" or "Revised" License

Common Lisp 99.75% Makefile 0.25%

hypergeometrica's Introduction

Hypergeometrica

Hypergeometrica aims to democratize the calculation of pi to trillions of digits. As of March 2020, the software used to break world-record computations has remained closed source. This has been a 20+ year trend, and includes famous software such as y-cruncher, TachusPI, PiFast, and SuperPi.

Please watch this introductory video.

CL-USER> (asdf:load-system :hypergeometrica)
CL-USER> (in-package :hypergeometrica)
#<PACKAGE "HYPERGEOMETRICA">
HYPERGEOMETRICA> (compute-pi/ramanujan 100)
31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
HYPERGEOMETRICA> (compute-pi/chudnovsky 100)
31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
HYPERGEOMETRICA> (compute-e 100)
27182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274

Unfortunately, Hypergeometrica cannot yet calculate pi in a completely competent way. What you see above does actually compute pi, but is taking a few very inefficient shortcuts. In order to be efficient, Hypergeometrica needs some additional key routines to eliminate these inefficient shortcuts.

What is the most efficient way to calculate pi with Hypergeometrica?

The call (MPD-PI b) for $b$ bits of precision is the fastest way to compute pi with Hypergeometrica. Here is a call to calculate at least $b=100$ bits of pi.

CL-USER> (in-package :hypergeometrica)
HYPERGEOMETRICA> (mpd-pi 100)

terms = 4
[0 Δ0] split
[0 Δ0] sqrt
[4 Δ4] recip
[8 Δ4] final
#<MPD {10160AD883}>

This output can be interpreted as:

terms = 4 means 4 terms of the Chudnovsky series were calculated.
[x Δy] thing means x milliseconds have elapsed since the start of the computation, and thing took y milliseconds since the last step
#<MPD {...}> is the object returned. Currently there is no base conversion routine to actually show this. One can use mpd-mpfr to convert it into an MPFR object for viewing or checking (as is done in hypergeometrica-tests::test-pi).

What is it?

Hypergeometrica is a Lisp library for performing extremely fast, extremely high-precision arithmetic. At the heart of it all are routines for doing fast multiplication. Hypergeometrica aims to support:

Fast in-core multiplication using a variety of algorithms, from schoolbook to floating-point FFTs.
Fast in-core multiplication for extremely huge numbers using exact convolutions via number-theoretic transforms. This is enabled by extremely fast 64-bit modular arithmetic.
Fast out-of-core multiplication using derivatives of the original Cooley-Tukey algorithm.
Implementation of dyadic rationals for arbitrary precision float-like numbers.
Elementary automatic parallelization when reasonable.

On top of multiplication, one can build checkpointed algorithms for computing a variety of classical constants, like pi.

How is it implemented?

It's a Lisp library that takes advantage of assembly code via SBCL's VOP facilities.

It would probably be easier to get higher performance quicker in C or C++, but there's a lot of non-hot-loop code (such as calculating suitable primes) that are better served without the baggage of a low-level language.

What works and what doesn't?

There's a test suite, I recommend looking at that to see what (should be) working. In any case, a short list of features:

Basic bigint (MPZ) routines.
Basic dyadic rational (MPD) routines.
Code to compute "suitable primes" for number-theoretic transforms.
An in-core number-theoretic transform employing tricks for fast modular arithmetic.
Binary-splitting for the computation of arbitrary hypergeometric series.
Out-of-core/disk-based number representation and automatic upgrading, with specialized algorithms.
In-core computation of pi, with basic asymptotically fast algorithms for division, square root, or inversion.

An implementation of disk-backed bigints exists, but it's not vetted and I'm not sure it's a good architecture.

There's also a broken implementation of out-of-core multiplication called the "matrix Fourier algorithm" following Arndt. Some corner case isn't working, and I'm not even sure this is the best way to do out-of-core multiplication.

Can I contribute?

Please, yes. Even if it's just telling me to document something. File an issue!

I know a lot about {I/O, disks, computer arithmetic, assembly, SBCL, ...} but I'm not really interested in rolling up my sleeves for this project.

Please contact me so I have somebody to ask questions to!

Where can I learn more about arbitrary precision arithmetic?

I'm keeping a list of references.

hypergeometrica's People

Contributors

Stargazers

Watchers

Forkers

salewski svspire hm0880 jwg4 hxzrx paulapatience

hypergeometrica's Issues

implement EQUAL-UPTO and optimize subtraction

Implement multiplication tuning

Implement a benchmarking procedure that tunes multiplication thresholds for different algorithms with different representations.

Implement the BBP hex digit extraction algorithm

Investigate using Garner's algorithm for CRT

Currently, CRT is computed by heavy use of Lisp's bignums. Garner's algorithm was implemented with the garner function. Should we be using that instead of the current CRT method? How do we take multiradix numbers and reconstruct them without consing bignums?

put TEST-SQRT into the tests package

implement efficient to-decimal conversion for MPDs

test general branch of DIVREM algorithm for small sizes

i fear that this branch doesn't actually work in specific cases due to weird tests. For example, 2/2 was resulting in weird quotients.

implement MPD-FRACTION-PART

Implement a function to get the fractional part of an MPD.

Implement efficient to-decimal base conversion for integers

Implement to-decimal base conversion using a sub-quadratic algorithm.

Implement Karatsuba's multiplication algorithm

Implement Karatsuba's multiplication algorithm for the MPZ and MPD data types.

Look into pkhuong's mod code

pkhuong wrote some code to do fast modding, which I understand could be better than Shoup's trick. Look into it.

Remove MPZ type

mpz, in its current form, is somewhat useless. The sign is extra baggage that isn't really needed and ought to not be lugged around. A lot of the arithmetic with mpz ought to just be arithmetic with digit storage. These will ultimately be the foundation for the floating point representation. We don't need both mpz and an mpf with duplicate code between them.

MPZ-SIZE and VEC-* are sometimes incompatible

Just because two MPZ-SIZEs are equal doesn't mean their vector sizes will be equal. This is a bad assumption that was made with MPZ-=. Are there other places where this assumption is being incorrectly made?

SUB128 fails when intrinsics are off

This is because the subtraction leads to a negative number which is asserted to not be possible. Investigate where and why this happens, and whether we need to make SUB128 accept any pair of arguments.

move plethora of NTT test code into tests

All the complex FFT stuff, NTT matrices, etc. don't belong in the main library code

MPZ-DIVREM may suffer from non-terminating expansions

if x/x divides to 0.9999..., it will be truncated to 0 remainder 1. We need to make sure this is handled and fixed.

division with newton sometimes does not converge

sometimes there are runaway cases where it does not converge.

rename DIGIT to LIMB, COEF, or something else

It has become too confusing in the code, in the debug messages, in the error messages, etc. what a "digit" is.

Internally in the code, a "digit" is a number in the radix $base. But when I personally read "digit", I often think of a decimal digit.

GMP uses "limb". These things can also be seen as "coefficients" in a polynomial representation of the number. Maybe there's some other better/short name.

In any case, it needs to change.

Implement sum-of-digits test

Implement sum-of-digits test for multiplication.

MFA should be smarter about how it does NTTs and when to split to RAM

MFA immediately delegates to NTT, which implies a single splitting from disk to RAM. This won't be right in all instances.

Implement non-power-of-2 transforms

Implement non-power-of-2, asymptotically fast NTTs. This can be based on the theory of FFTs.

Confusion over rounding, truncation, significant figures, and precision.

Lots of bugs have popped up due to sloppy conflation of these terms.

Moreover, we've been inconsistent in asking for "bits of precision", "digits of precision", etc.

We need to comb through the code and be very precise on what we mean.

Code for manipulating decimal digits.

Let's move the functions digit-count, digits, carry, and undigits from the test suite into the library itself. Then these functions can themselves be tested. It will facilitate #1 which requires that or similar code.

(TEST-PI 6) is too slow

This (intentionally) uses disk to perform the computation, but intermediate numbers seem to be blowing up to incredible and disproportionate sizes.

intrinsics are broken due to a variety of SBCL changes

SBCL made a variety of changes to assembly, VOPs, etc. which has broken what's there. Intrinsics have been disabled, but they're critical for fast code.

put TEST-RECIP into the tests package

add64 still broken sometimes

ADD64 gives bad sum sometimes

Out-of-core mult doesn't work for 9^9^9

I propose calculating 9^9^9 as the multiplication stress test. It seems to work for in-core multiplies, but not out-of-core.

Implement the Bellard hex digit extraction algorithm

TEST-MPZ-SIZE fails

Running (asdf:test-system :hypergeometrica) results in the following exception:

The function
  #<STANDARD-GENERIC-FUNCTION HYPERGEOMETRICA::INTEGER-MPZ (1)>
requires at least 2 arguments.
   [Condition of type SB-INT:SIMPLE-PROGRAM-ERROR]

Restarts:
 0: [CONTINUE] Skip the rest of the test TEST-MPZ-SIZE and continue by returning (values)
 1: [RETEST] Rerun the test TEST-MPZ-SIZE
 2: [CONTINUE] Skip the rest of the test FIASCO-SUITES::HYPERGEOMETRICA-TESTS and continue by returning (values)
 3: [RETEST] Rerun the test FIASCO-SUITES::HYPERGEOMETRICA-TESTS
 4: [CONTINUE-WITHOUT-DEBUGGING] Turn off debugging for this test session and invoke the first CONTINUE restart
 5: [CONTINUE-WITHOUT-DEBUGGING-ERRORS] Do not stop at unexpected errors for the rest of this test session and continue by invoking the first CONTINUE restart
 --more--

Backtrace:
0: (SB-INT:%PROGRAM-ERROR "~@<The function ~2I~_~S ~I~_requires at least ~W argument~:P.~:>" #<STANDARD-GENERIC-FUNCTION HYPERGEOMETRICA::INTEGER-MPZ (1)> 2)
1: (SB-PCL::ERROR-NEED-AT-LEAST-N-ARGS #<STANDARD-GENERIC-FUNCTION HYPERGEOMETRICA::INTEGER-MPZ (1)> 2)
2: (SB-PCL::INITIAL-DFUN #<STANDARD-GENERIC-FUNCTION HYPERGEOMETRICA::INTEGER-MPZ (1)> (0))
3: ((LABELS TEST-MPZ-SIZE :IN TEST-MPZ-SIZE))
4: ((LABELS FIASCO::RUN-TEST-BODY :IN FIASCO::RUN-TEST-BODY-IN-HANDLERS))
5: (FIASCO::CALL-WITH-TEST-HANDLERS #<FUNCTION (LAMBDA NIL :IN FIASCO::RUN-TEST-BODY-IN-HANDLERS) {10077C93CB}>)
6: (FIASCO::PRETTY-RUN-TEST #<test TEST-MPZ-SIZE> #<FUNCTION (LABELS TEST-MPZ-SIZE :IN TEST-MPZ-SIZE) {538CB92B}>)
7: ((LABELS #:BODY-SYM0 :IN TEST-MPZ-SIZE))
8: (TEST-MPZ-SIZE)
9: ((LABELS FIASCO-SUITES::HYPERGEOMETRICA-TESTS :IN FIASCO-SUITES::HYPERGEOMETRICA-TESTS))
10: ((LABELS FIASCO::RUN-TEST-BODY :IN FIASCO::RUN-TEST-BODY-IN-HANDLERS))
11: (FIASCO::CALL-WITH-TEST-HANDLERS #<FUNCTION (LAMBDA NIL :IN FIASCO::RUN-TEST-BODY-IN-HANDLERS) {1003EC993B}>)
12: (FIASCO::PRETTY-RUN-TEST #<test FIASCO-SUITES::HYPERGEOMETRICA-TESTS :tests 22> #<FUNCTION (LABELS FIASCO-SUITES::HYPERGEOMETRICA-TESTS :IN FIASCO-SUITES::HYPERGEOMETRICA-TESTS) {538B267B}>)
13: ((LABELS #:BODY-SYM1 :IN FIASCO-SUITES::HYPERGEOMETRICA-TESTS))
14: (FIASCO-SUITES::HYPERGEOMETRICA-TESTS)
15: (RUN-SUITE-TESTS #<test FIASCO-SUITES::HYPERGEOMETRICA-TESTS :tests 22> :VERBOSE NIL :STREAM #<SB-SYS:FD-STREAM for "socket 127.0.0.1:36697, peer: 127.0.0.1:48318" {1001358393}> :INTERACTIVE T)
16: (RUN-TESTS :HYPERGEOMETRICA-TESTS :DESCRIBE-FAILURES T :VERBOSE NIL :STREAM #<SB-SYS:FD-STREAM for "socket 127.0.0.1:36697, peer: 127.0.0.1:48318" {1001358393}> :INTERACTIVE T)
17: (TEST-HYPERGEOMETRICA)
18: ((SB-PCL::EMF ASDF/ACTION:PERFORM) #<unused argument> #<unused argument> #<ASDF/LISP-ACTION:TEST-OP > #<ASDF/SYSTEM:SYSTEM "hypergeometrica/tests">)
19: ((LAMBDA NIL :IN ASDF/ACTION:CALL-WHILE-VISITING-ACTION))
--more--

I tried adding 'mpz/ram as the second argument in the INTEGER-MPZ calls in TEST-MPZ-SIZE, to no avail:

There is no applicable method for the generic function
  #<STANDARD-GENERIC-FUNCTION HYPERGEOMETRICA::INTEGER-MPZ (1)>
when called with arguments
  (0 MPZ/RAM).
   [Condition of type SB-PCL::NO-APPLICABLE-METHOD-ERROR]

Restarts:
 0: [RETRY] Retry calling the generic function.
 1: [CONTINUE] Skip the rest of the test TEST-MPZ-SIZE and continue by returning (values)
 2: [RETEST] Rerun the test TEST-MPZ-SIZE
 3: [CONTINUE] Skip the rest of the test FIASCO-SUITES::HYPERGEOMETRICA-TESTS and continue by returning (values)
 4: [RETEST] Rerun the test FIASCO-SUITES::HYPERGEOMETRICA-TESTS
 5: [CONTINUE-WITHOUT-DEBUGGING] Turn off debugging for this test session and invoke the first CONTINUE restart
 --more--

Backtrace:
0: ((:METHOD NO-APPLICABLE-METHOD (T)) #<STANDARD-GENERIC-FUNCTION HYPERGEOMETRICA::INTEGER-MPZ (1)> 0 MPZ/RAM) [fast-method]
1: (SB-PCL::CALL-NO-APPLICABLE-METHOD #<STANDARD-GENERIC-FUNCTION HYPERGEOMETRICA::INTEGER-MPZ (1)> (0 MPZ/RAM))
2: ((LABELS TEST-MPZ-SIZE :IN TEST-MPZ-SIZE))
3: ((LABELS FIASCO::RUN-TEST-BODY :IN FIASCO::RUN-TEST-BODY-IN-HANDLERS))
4: (FIASCO::CALL-WITH-TEST-HANDLERS #<FUNCTION (LAMBDA NIL :IN FIASCO::RUN-TEST-BODY-IN-HANDLERS) {100B00560B}>)
5: (FIASCO::PRETTY-RUN-TEST #<test TEST-MPZ-SIZE> #<FUNCTION (LABELS TEST-MPZ-SIZE :IN TEST-MPZ-SIZE) {5390762B}>)
6: ((LABELS #:BODY-SYM0 :IN TEST-MPZ-SIZE))
7: (TEST-MPZ-SIZE)
8: ((LABELS FIASCO-SUITES::HYPERGEOMETRICA-TESTS :IN FIASCO-SUITES::HYPERGEOMETRICA-TESTS))
9: ((LABELS FIASCO::RUN-TEST-BODY :IN FIASCO::RUN-TEST-BODY-IN-HANDLERS))
10: (FIASCO::CALL-WITH-TEST-HANDLERS #<FUNCTION (LAMBDA NIL :IN FIASCO::RUN-TEST-BODY-IN-HANDLERS) {10076F2F4B}>)
11: (FIASCO::PRETTY-RUN-TEST #<test FIASCO-SUITES::HYPERGEOMETRICA-TESTS :tests 22> #<FUNCTION (LABELS FIASCO-SUITES::HYPERGEOMETRICA-TESTS :IN FIASCO-SUITES::HYPERGEOMETRICA-TESTS) {538B267B}>)
12: ((LABELS #:BODY-SYM1 :IN FIASCO-SUITES::HYPERGEOMETRICA-TESTS))
13: (FIASCO-SUITES::HYPERGEOMETRICA-TESTS)
14: (RUN-SUITE-TESTS #<test FIASCO-SUITES::HYPERGEOMETRICA-TESTS :tests 22> :VERBOSE NIL :STREAM #<SB-SYS:FD-STREAM for "socket 127.0.0.1:36697, peer: 127.0.0.1:48318" {1001358393}> :INTERACTIVE T)
15: (RUN-TESTS :HYPERGEOMETRICA-TESTS :DESCRIBE-FAILURES T :VERBOSE NIL :STREAM #<SB-SYS:FD-STREAM for "socket 127.0.0.1:36697, peer: 127.0.0.1:48318" {1001358393}> :INTERACTIVE T)
16: (TEST-HYPERGEOMETRICA)
17: ((SB-PCL::EMF ASDF/ACTION:PERFORM) #<unused argument> #<unused argument> #<ASDF/LISP-ACTION:TEST-OP > #<ASDF/SYSTEM:SYSTEM "hypergeometrica/tests">)
18: ((LAMBDA NIL :IN ASDF/ACTION:CALL-WHILE-VISITING-ACTION))
19: ((:METHOD ASDF/ACTION:PERFORM-WITH-RESTARTS :AROUND (T T)) #<ASDF/LISP-ACTION:TEST-OP > #<ASDF/SYSTEM:SYSTEM "hypergeometrica/tests">) [fast-method]
--more--

I know basically nothing about this library, as I just discovered it, so perhaps the error is trivial, and my attempted fix ignorant. I may investigate later. I'm not sure if this error is related to #23.

Test MPD against MPFR

Create a robust set of unit tests that checks for the correctness of MPDs against MPFR.