bignum / BigFloat

An infinite-precision floating-point calculator and function library for Miranda and type for Haskell

Introduction

There are many arithmetical packages available to perform scientific calculations to frightening precision but, without exception, you have to specify the required precision (usually a fixed number of decimal places) before starting the calculation, with the unfortunate side-effect that the last few places of a result of a complicated calculation are not reliable, depending as they do on the precision of intermediate results. Consequently, such an apparently super-precise result can in fact be misleading, and a rigorous study of the precision lost in intermediate calculations is required to be able to determine how many digits of the result are sure to be correct.

This calculator is different. It evaluates to a depth such that it is always certain that the digits of the results that it prints are correct, with no rounding or intermediate errors. To obtain a result to a certain number of figures it is sufficient to discard all the figures that are not required. In its native form, the calculator outputs an infinite number of digits, limited only by the amount of physical memory available on the host machine, which is a product of the complexity of the calculation in progress.

It works by doing the minimum amount of carry-lookahead required to be certain of the next digit it must print. The syntax in which expressions are entered is currently horrendous, since it consists of calls to functions written in the Miranda functional programming language; it requires a normal infix parser to accept the usual syntax for formulae. That said, it is usable.

News

Newcastle, 11 Sep 2008

Version 4.2 is published, containing a much faster generator for pi;

Newcastle, 14 Nov 2007

Version 4.1 fixes the ubn_add bug mentioned below. It only took me two years to get round to releasing it...

Newcastle, 15 Nov 2005

There is a bug in the implementation of the overly complex ubn_add2 function, which was only detected by the (forthcoming) faster arcsin implementation giving wrong results from the 80th or 99th digit onwards, but it may affect other calculations.

An immediate fix to the source of all versions is to replace the line

ubn_add = ubn_add2

with

ubn_add = ubn_add1

Newcastle, 6 Nov 2005

bignum/BigFloat now has a Home Page and Download Site on SourceForge.

Nijmegen, 5 Oct 2005

BigFloat participated in the "manydigits" competition and came last for speed.

The winning entries obtained inaccurate answers very fast and the competition's rules were changed during the competition to calculate the required number of digits plus ten and lop off the last ten digits.

Professor David Turner said: "I'd have give minus infinity points for getting the wrong answer."

Sittingbourne, 26 Feb 2004

BigFloat has calculated

• pi to 3,481,598 decimal places (agrees with [pi to 1,000,000 places] (http://3.141592653589793238462643383279502884197169399375105820974944592.com)
• e to 133,705 decimal places (agrees with Gutenberg's [e to 1 million places] (ftp://sailor.gutenberg.org/pub/gutenberg/etext94/ee710.txt)
• The square root of 2 to 1,930,097 decimal places (the first million agree with the [Gutenberg million-digit 2sqrt] (http://www.ibiblio.org/gutenberg/etext93/2sqrt10.txt))
• The golden ratio to 433,174 decimal places.

Bignum's sqrt(2) calculation has exceeded the existing publicly-available value (1,250,000 figures). It managed 1,930,097 decimal places having run from 20 Jan - 26 Feb 2004. Then some idiot managed to turn the power off in the house.

Syntax

Here is the input syntax for the Miranda and KRC versions; the Haskell version just creates a new type BigFloat, to which all the usual operations apply.

Calculated constants

• bn_0: The bignum version of zero
• bn_1: The bignum version of one
• bn_2: The bignum version of two
• bn_e: 2.7182818...
• bn_pi: 3.1415926535...
• bn_pi_2, bn_pi_4 and bn_2_pi: pi/2, pi/4 and twice pi

Functions taking a single parameter

• bn: takes a character string containing a decimal number and converts it to a bignum. The string can have an arbitrary number of decimal places after the decimal point.
• bn: takes an integer and converts it to a bignum (Note: an integer, not a floating point value, nor the result of a floating point calculation. You can try converting those using bn (show x), or more simply with the compound function (bn.show), but entirely at your own risk!)
• showbignum: is the inverse function of bn. It converts a bignum to a printable form. Miranda calls this function automatically whenever it has to display a bignum.
• bn_twice: doubles its parameter (faster than bn_times 2 x)
• bn_half: halves its parameter (faster than bn_quot x 2)
• bn_sqr: returns the square of its argument (faster than bn_mul x x)
• bn_sqrt: returns the square root of its argument
• bn_neg: negates its argument
• bn_rnd: takes a positive integer seed as its parameter and returns a pseudo-random bignum in the range 0 <= n < 1

Functions taking two parameters

• bn_add: takes two bignum parameters and returns their sum.
• bn_sub: subtracts its second bignum parameter from its first, returning the difference.
• bn_times: takes an integer and a bignum and returns their bignum product. The integer is the first of the two parameters, so as to ease partial parameterisazion.
• bn_mul: multiplies two bignums
• bn_quot: divides a bignum by an positive integer. If you want to partially parameterise the integer divisor, use ($bn_quot n).
• bn_div: performs long division of two bignums
• bn_raise: takes a bignum and an integer and returns the bignum raised to that integral power
• bn_pow: takes two bignums and returns the first raised to the power of the second

Trigonometric functions

All these take and return one bignum. All angles are expressed in radians.

• bn_ln: returns the natural logarithm (in base e) of its parameter
• bn_exp: returns the inverse natural logarithm of its parameter
• bn_sin: returns the sine of its parameter. It is fastest for values between -pi/2 and +pi/2
• bn_cos: returns the cosine of its parameter (defined using sin).
• bn_tan: returns the tangent of its parameter (defined using sin/cos).
• bn_asin: returns the arcsine (in the range -pi/2 to pi/2) of its parameter (defined using atan)
• bn_acos: returns the arccosine (in the range -pi/2 to pi/2) of its parameter (defined using atan)
• bn_atan: returns the arctangent (in the range -pi/2 to pi/2) of its parameter.
• bn_sinh, bn_cosh, bn_tanh, bn_asinh, bn_acosh, bn_atanh: return hyperbolic sine, cosine and tangent and their inverses.

Logical functions

For the convenience of Miranda programmers who wish to use the bignum library as part of a larger program, there is also a set of comparison functions of little use in desk-calculator mode:

• bn_is_zero: tells you whether a bignum is zero.
• bn_is_neg: tells you whether a bignum has a negative value.
• bn_cmp: compares two bignums and returns 0 if the two numbers have the same value, 1 if the first number is greater than the second, -1 if the second number is greater than the first.
• bn_eq, bn_ne, bn_gt, bn_lt, bn_ge, bn_le: The six usual comparison functions: equal to, not equal to, greater than, less than, greater than or equal to, less than or equal to. They all take two bignum parameters and return a boolean value.

Examples

bn_sqrt bn_2

The square root of two.

bn_raise (bn_mul (bn "1.5") (bn "0.75")) 6

The interval obtained by tuning up by a perfect fifth and down by a perfect fourth, six times (not quite an octave!)

bn_quot (bn_add bn_1 (bn_sqrt (bn "5"))) 2

The golden ratio, calculated as (1 + sqrt(5)) / 2

Note that, since functions are "greedy" to bind to their parameters, brackets are necessary to make "5" the parameter to bn, rather than bn being taken as the first parameter to bn_sqrt (which would make no sense).

Avoiding bottoming out

Bignum calculator's main drawback is that calculations can "bottom out" which means that the calculation goes on forever without producing a result, or blocks halfway through.

This happens when a calculation involving infinitely-long figures turns out to give a finite-length result. Simple examples are subtracting an infinitely long number from itself:

bn_e $bn_sub bn_e

(a $bn_sub b is the infix version of the regular syntax bn_sub a b) or the classic "sum of three thirds":

bn_times 3 (bn_quot bn_1 3)

because the carry lookahead mechanism has to inspect an infinite list of zeroes to know that the result really is zero, or an infinite list of nines wondering whether there is going to be a carry into it or not.

For a more insidious example of this (and a way of avoiding it), suppose you wanted to check the correctness of Dudeney's solution to the problem of finding two rational numbers, other than 1 and 2, whose cubes add up to 9.

(pn/pd)^3 + (qn/qd)^3 = 9

Dudeney's solution is

415280564497 676702467503
------------ and ------------
348671682660 348671682660

This appeared in his first book, "The Canterbury Puzzles" in 1907, and its denominator is shorter than any other previously published solution. As Martin Gardner observes, he did this without the aid of a modern calculating machine, which makes his achievement almost a miracle!

A naive but useless attempt to verify this would be (as a Miranda script):

%include "bignum"
bn_pn = bn "415280564497"
bn_qn = bn "676702467503"
bn_d = bn "348671682660"
bn_dud = (bn_raise (bn_pn $bn_div bn_d) 3) $bn_add (bn_raise (bn_qn $bn_div bn_d) 3)

but the sum of two infinite-length bignums giving an finite-length result never terminates. It just goes looking forever down an infinite list that starts out 8.99999999999999999999999... for a carry that would resolve its doubts about the first digit.

You can make it work, however, by rearranging the formula

(pn/pd)^3 + (qn/qd)^3 = 9

as

(pn^3 + qn^3) / d^3 = 9

As a Miranda script:

%include "bignum"
bn_pn3 = bn_raise (bn "415280564497") 3
bn_qn3 = bn_raise (bn "676702467503") 3
bn_d3 = bn_raise (bn "348671682660") 3
dudeney = bn_div (bn_add bn_pn3 bn_qn3) bn_d3

which does indeed give 9.

(Though, to be fair, you don't have to use the bignum package to perform this calculation. Miranda's native integers are infinitely precise, so you can use the rearranged formula to get the top and bottom of the division as integers, then check that the integer division is 9 and that the remainder is zero).

Free Bignum calculator in Internet

You can use the bignum calculator at the address

telnet://medialab.freaknet.org:31415

where you can type in a bignum expression and see the result. If the result has an infinite number of digits, you will need to close the telnet session to terminate the calculation.

In reality, this is a free Miranda evaluator preloaded with the bignum package and with the various command- and shell-escapes disabled. For safety, it is limited to 15 minutes of processing time per call and a maximum memory usage of 10 megabytes. Users requiring more processing time should contact martin@freaknet.org.

Here is an example session (from a Unix command prompt):

$ telnet medialab.freaknet.org 31415
Trying 78.134.26.141...
Connected to medialab.freaknet.org.
Escape character is '^]'.
bn_raise (bn_mul (bn "1.5") (bn "0.75")) 6
2.027286529541015625
Connection closed by foreign host.
$ 

If you use telnet from MSWindows, often it doesn't show you what you are typing, the <-- key doesn't work, and the output is sometimes garbled. Try fiddling with the option settings to enable local echo, or telnet to medialab.freaknet.org with login name luther and password luther, and then use:

telnet localhost 31415

from there.

Source code

Bignum is maintained in three versions:

• the Miranda version implements a Bignum type and a set of functions that operate upon them (as described in the Syntax section).
• the Haskell version just implements a BigFloat type, on which all the usual operations can be performed as with any other Floating type. See the comments at the start of the source file.
• the KRC version, a translation of the Miranda version.

Stable releases are on SourceForge

Development is on Codeberg

Other bignum calculation system

At least two other arbitrary-precision calculating systems are included in "standard" Unix: bc and perl both of which make you specify the required precision of the result in advance, and both of which make you wait for the whole calculation to complete before they will tell you any digits of the result.

bc

This classic Unix calculator always works with integers of unbounded size. When invoked with the -l option, you can specify the number of decimal places by setting the magic scale variable.

For example, to calculate 200 decimal places of "pi" as four times the arctangent of one:

$ bc -l
scale=200
4 * a(1)

Authors:

• of the original Unix version: Lorinda L. Cherry & Robert Morris.
• of the GNU reimplementation: Philip A. Nelson, philnelson@acm.org

perl

Perl has a bignum library which replaces all the regular numerical operators with arbitrary-precision versions, and is documented under the Unix manual page bignum(3). It implements square roots and natural logarithms, but fails to provide substitutes for regular perl's sin, cos and atan2 functions.

To calculate the square root of two to 2000 decimal places, you can say:

perl -Mbignum=p,-2000 -le 'print sqrt 2'

An interesting sister package by the same author is bigrat which works in, and gives its results as, vulgar fractions. For example:

perl -Mbigrat -le 'print 1/3 + 1/4'
7/12

This nicely solves the "sum of three thirds" problem for which bc and perl -Mbignum give .999999999... and which just gives bignum/bigfloat a terminal headache.

Author:

• Tels on bloodgate.com in early 2002.

Others

• Genius Math tool and the GEL Language
• GMP: GNU Multiple Precision a extensive library for C programming, "the fastest".
• A list of bignum packages, compiled in 1994, floats around on the internet: search for the string "This is the file BIGNUMS.TXT"
• For the Oct 2005 "manydigits" competition, Milad Niqui compiled another list of bigfloat packages

Author

Martin Guy martinwguy@gmail.com, 1985 - 15 November 2007.