Real.README - documentation for real.py. Goes with version 1.11.
Written by Jurjen N.E. Bos (jurjen@q2c.nl)

Introduction
============
real.py is a library that introduces a new class, called Real, of
abitrarily precise numbers, allowing computations with "infinite"
precision.
Each number keeps an approximation of the accuracy of the number.
Routines are available to compute numbers with more and more digits.
real.py allows you to:
- compute with very high precision, e.g. 1000 digits
- compute with very large or small numbers, up to about 10**500000000
- compute without unexpected precision loss: all digits guaranteed
- compute with infinite precision, printing digits on the fly
- printing floating point numbers to their exact accuracy

Just to whet your appetite:
>>> from real import pa,exp,pi,sqrt
>>> pa('exp(pi()*sqrt(163))')
262537412640768743.99999999999925007259719818568887935385633733699086270
753741037821064791011860731295118134618606450419308388794975386404490572
871447719681485232243203911647829148864228272013117831706501045222687801
444841770346969463355707681723887681000923706539519386506362757657888558
223948114276912100830886651107284710623465811298183012459132836100064982
...
and so on.

The function r() in the library constructs Real numbers, and has a
cornucopia of options. Details are shown below.

The library supports the following operators:
+, -, abs
+, -, *, /, **; also both ways with int
<< and >> with int
int(), long()

The library supports the following mathematical functions:
floor, ceil
sqrt, root (n-th root)
exp, log, log10
atan, asin, acos, sin, cos, tan
sinh, cosh, tanh, asinh, acosh, atanh
fact (for all arguments!)
pi and e are functions taking number of digits (with default)

All of these functions are also available as methods. NumPy users
should be happy with this.

The library supports the following constants, implemented as functions
taking a number of digits:
pi
e

Simple examples
---------------
This library allows calculations with arbirary precision numbers:
>>> from real import *
>>> 1/(pi()**3-31)
159.3198876209530096514514069+-3
pi() is a function computing a given number of digits of pi. The default
is real.default (initially 32):
>>> pi(50)
3.141592653589793238462643383279502884197169399375105+-2
Note the precision indication "+-2" at the end: it means that the last
digit may be 2 too high or 2 too low. In fact, the digit is correct:
>>> pi(52)
3.14159265358979323846264338327950288419716939937510582+-2
The precision indicator is "built in" to every number of the Real
class.
In this way, the digits you get are always right (but not neccesarily
plentiful).
For example:
>>> sin(3L**61)
.940+-1
This should be read as "a number between 0.939 and 0.941".

The class Real allows for very large and small numbers:
>>> fact(10000000)
1.2024234005159034561401534879443075698+-1e65657059
>>> exp(exp(exp(3)))
2.05098643605164889510586094208878506609+-5e229520860

The function r() allows to construct numbers with given precision:
>>> r(3L**61,10)
1.27173474827+-5e29
>>> r(3L**61,50)
127173474825648610542883299603.0000000000000000000000+-5
The function can also be used as a "pretty printer" for floating-point
numbers:
>>> import math
>>> math.pi
3.14159265359
>>> r(math.pi)
3.1415926535897931+-5
The precision of the number math.pi is now explicitly shown.

Infinite precision
------------------
The function pa() prints "all" digits of a number. One of the simplest
examples is computing "all" digits of pi:
>>> pa(pi)
3.1415926535897932384626433832795028841971693993751058209749445923078164
062862089986280348253421170679821480865132823066470938446095505822317253
(and so on).
The routine pa() computes a number repeatedly with increasing precision,
each time printing the digits it can guarantee.
If you pass a string to pa(), you can compute arbitrary expressions:
>>> pa('sqrt(pi()**7-e()**8-39)')
.57899977121221107665115128180224906923190232347298078138819681954282919
919476693269700788120335999771483821199512086460383486796180647225032362
...

About precision
---------------
All computations and numbers in the Real class have a very carefully
defined notion of precision. Normally, when writing "precision" I mean
the number of digits in the mantissa of the number. Thus:
>>> r(3)**123456789
1.17852993898070743873540+-4e58903858
has a precision of just under 24 digits. Internally, the precision is
stored in bits, but you don't have to know that.
The default precision is real.default, which is initially set to a
generous 32 digits:
>>> pi()
3.141592653589793238462643383279503+-2
>>> real.default=20
>>> pi()
3.141592653589793238462+-2

The number 0 does not have a precision, so when a precision is needed,
it is taken to be the default precision.

A computation may suffer from precision loss. For example,
the above mentioned computation:
>>> real.default=32
>>> sin(3L**61)
.940+-1
The reason this approximation is so inaccurate is that the last digit
of the default approximation of the number 3**61 influences the result
by 0.001. Since the sine function's result cannot be more accurate than
its argument, the result can also not be more accurate than 0.001. So,
the sin function cannot guarantee more digits, if no more digits of
3**61 are known. We can use r() to get more digits:
>>> sin(r(3L**61,50))
.940023765661479919109+-1
Now the result is much more accurate, but still less than the 50 digits
of the argument, of course.
Let's compare this with the built-in floating point numbers:
>>> import math
>>> math.sin(3L**61)
-0.700599464638
What? All digits (and even the sign) are wrong! This is not because
there is a bug in the math library; it is because the floating-point
representation of 3**61 is not accurate to the last digit:
>>> float(3L**61)
1.27173474826e+29
This implies that the sin gives a value based on a number that is
useless. This also shows one of the powers of "real": it does not
give wrong answers.

If a function has very smooth behaviour, it is also possible that it
returns more digits than its argument:
>>> r(10)**10
10000000000.00000000000000000000+-3
>>> atan(r(10)**10)
1.5707963266948966192313216916400847754320+-2
Now we have no less than 41 digits of precision. As a rule, all
functions produce at most twice the precision of their input by rounding
if necessary; this is to prevent excessive computation times.

Mixed arithmetic
----------------
Number in the Real class may be mixed with integers in computations.
Integers and longs are assumed to be _exact_ by the +,_,*,/ operators:
>>> pi()
3.141592653589793238462643383279503+-2
>>> pi()-3
.141592653589793238462643383279503+-2
You can see that no precision loss occurs in the subtraction.
The integer powers -1, 0, 1, 2, 3 are implemented as multiplication or
divisions for speed and accuracy:
>>> 1/pi()
.318309886183790671537767526745029+-1
>>> pi()**-1
.318309886183790671537767526745029+-1
>>> pi()**r(-1)
.3183098861837906715377675267450+-1

Float and strings are converted using the function r() if they are
mixed with Reals: see below.
>>> pi()+math.pi+"3.14159"
9.42474+-4
>>> sqrt(1.5)
1.2247448713915889+-3

The function r(): constructing Reals
====================================
The function r() allows to construct Real numbers from integers, longs,
fractions, floating point numbers and strings.

Constructing Reals from integers and longs
------------------------------------------
r(n) makes n into a real number, with default precision.
r(n,p) makes n into a real number, with precision of p digits. There is
no loss of accuracy in the conversion:
>>> r(1)
1.00000000000000000000000000000000+-1
>>> r(1,10)
1.0000000000+-1

Construction of Reals from fractions
------------------------------------
r(n,d,p) makes n/d into a real number, with precision of p digits. There
is (essentially) no loss of precision in the conversion:
>>> r(1,3,20)
.333333333333333333332+-4
To get the default precision, use real.default or None for p.

Construction of Reals from floating point numbers
-------------------------------------------------
r(f) makes a floating point number into a real with the (at for my
machine) correct amount of digits: 53 bits. This indicates nicely
how accurate the number is:
>>> math.pi
3.14159265359
>>> r(math.pi)
3.1415926535897931+-5
This shows that the default printing routine of floating point numbers
hides some of the significant digits.
It is also allowed to use r(f,p) to get p digits of a floating point
number, but extra digits are "invented" if you increase the precision
beyond that of the floating point number. Be warned!
>>> r(math.pi,30)
3.141592653589793115997963468544+-4
>>> pi(30)
3.1415926535897932384626433832795+-2

Construction of Reals from Reals
--------------------------------
r(x), where x is a real, just returns x.
r(x,p) changes the precision of x to p. The same warning as above
applies: don't increase the precision unless you don't care about the
actual value, or you know what you are doing.
>>> r(pi(30),10)
3.1415926537+-3
>>> r(pi(10),30)
3.141592653584666550159454345703+-4

Construction of Reals from strings
----------------------------------
To allow you to convert the representation format back to a number, you
can construct a Real from a string. The format is:
number = mantissa [error] [exponent]
mantissa = ["+"|"-"] digit* ["." digit*]
error = "+-" digit
exponent = ("E"|"e") ["+"|"-"] digit+

If the period is omitted, it is assumed to be at the end of the
mantissa. If the mantissa is omitted altogether, it is assumed to be
"0".
If the error digits are omitted, it is taken to be +- one half of the
last digit.
If the exponent is omitted, it is taken to be e0.

The string representation of a Real is also in this format. However, due
to roundoff, the error estimation may increase a bit if a number is
converted to string and back.

Also here, r(s,p) may be used to modify the precision of the result
(with the same caveat as above: don't increase the precision if you
want correct accuracy!)
>>> r('1234.567e3')
1234567.0+-5
>>> r('1234.567+-1e3')
1234567.+-1
>>> r('1234.5678e3')
1234567.8+-1
>>> r('1234.5678e3',20)
1234567.81250000000000+-2

Precision of 0
--------------
As a special rule, if you let r() compute 0 to a given precision, it
will take the second argument as "digits after the point" instead of
"mantissa digits":
>>> r(0,10)
.00000000000+-3

Implicit conversions
--------------------
r() is also called implicitly if a function from real.py is called with
an argument that is not a Real:
>>> fact(10)
3628800.00000000000000000+-1
>>> sin(math.pi)
.000000000000000+-1
>>> cos('0.000010000000')
.999999999951+-2

Using the library
=================
The most useful (and fun) functions of the library are rep() and pa().

The function rep()
------------------
The function rep computes the result of an expression using increasing
precision. rep() accepts a function from number of digits to result, but
also a string which is interpreted as an expression:
>>> rep(e)
2.718281828459045235360287471352663+-2
2.71828182845904523536028747135266249775724709+-2
2.71828182845904523536028747135266249775724709369995957496696+-1
...
>>> rep('sin(1)')
.84147098480789650665250232163030+-2
.8414709848078965066525023216302989996225630+-2
.8414709848078965066525023216302989996225630607983710656728+-2
...
rep() is designed so that the amount of time to compute any result is
not more than twice the time spent since the start of the function. This
means that if you want as much precision as you can get, you don't have
to wait more than twice as long as you would wait computing to only that
precision.
This is very useful if you don't want what precision to choose.

The optional second argument of rep() gives the precision of the first
computation:
>>> rep(pi,3)
3.1416+-2
3.14159+-2
3.1415926+-2
3.141592653+-2
3.141592653590+-2
3.1415926535897932+-2
3.141592653589793238462+-2
3.1415926535897932384626433833+-2
3.1415926535897932384626433832795028842+-2
3.1415926535897932384626433832795028841971693993752+-2
...

The function pa()
-----------------
pa() accepts the same arguments as rep(), but attempts to print the
result only once, adding new digits as soon as they are computed:
>>> pa('sin(1)')
.84147098480789650665250232163029899962256306079837106567275170999191040
439123966894863974354305269585434903790792067429325911892099189888119341
032772921240948079195582676660699990776401197840878273256634748480287029
...
The function pa() cannot print numbers that end with zeroes, since it
can never be sure whether this is not a really long sequence of nines.
In particular, something as simple as pa('1+1') does _not_ give an
answer! For these cases, rep() can be used as a replacement.
>>> pa('r(123456789)')
12345678<dead silence>

To show a nice example why pa() cannot print the next digit, try:
>>> pa('log10(sin(10**-r(5432,100,None)))',8)
-54.32000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000000000000001658185300612831004479789196460
369758485244604552424332996973775238492796390137138858936927630712568532
...
>>> pa('log10(tan(10**-r(5432,100,None)))',8)
-54.31999999999999999999999999999999999999999999999999999999999999999999
999999999999999999999999999999999999999996683629398774337991040421607079
260483029510790895151334006052449523014407219725722282126144738574862934
...
In both cases, after having printed "-54.3", pa() takes some time
computing enough digits to decide if the next digit is a 1 or a 2. The
difference is _very_ small, hence the delay.
Also note the use of r(5432,100,None) to get the constant 54.32 to
default precision.

The factorial
-------------
This will probably not surprise you:
>>> fact(10)
3628800.00000000000000000+-1
In fact, the factorial function does a lot more than multiplying
numbers. The fact() function actually compute the gamma function of
n+1:
>>> fact(10+r(1,10,None))
4593083.5895600181037600919+-2
We can check if the gamma function of 1/2 is indeed sqrt(pi), as is
claimed in the schoolbooks:
>>> fact(r(1,2,50)-1)
1.77245385090551602729816748334114518279754946+-1
>>> sqrt(pi(50))
1.772453850905516027298167483341145182797549456122387+-2
The function can also produce very small numbers:
>>> fact("-10000000.500000000000000")
8.2621381+-5e-65657056
Actually, the function fact() is the most complicated funcion of the
entire library. (And that's just why I wrote it :-)

Precision
---------
Sometimes, strange things happen with precision. Don't let yourself
be fooled. An example:
>>> exp(pi(31)*sqrt(163))
262537412640768744.000000000000+-2
The most plausible conclusion is that exp(pi*sqrt(163)) is a whole
number, but we all know that's pretty unlikely. In fact, this is just
a freak event, as you can see by increasing the precision:
>>> real.default=50
>>> exp(pi()*sqrt(163))
262537412640768743.999999999999250072597198185690+-2

Differential quotients
----------------------
The differential quotient of a function is defined as:
(f(x+h)-f(x))/h
Since math class we know that if h goes to zero, the result of the
quotient goes to the derivative of f in point x: f'(x).
This differential quotient can be used to compute derivatives, if h is
chosen carefully. For numerical reasons too hard to explain here, the
correct value to choose for h is 10**(-default/2). This works for all
functions that have reasonable behaviour (i.e. the second derivative is
smaller than, say, 100).
This trick can be used to print one of the most complicated mathematical
constants, called Euler's constant (denoted by lower case gamma).
We use that gamma = -fact'(0), and fact(0)=1.
(By the way: don't expect 10000 digits from this computation, unless you
have an exceptional computer. This computation is _slow_. I never came
any further than 1331 digits.)
>>> pa('(lambda d:(1-fact(d))/d)(10**-r(real.default/2))')
.57721566490153286060651209008240243104215933593992359880576723488486772
677766467093694706329174674951463144724980708248096050401448654283622417
399764492353625350033374293733773767394279259525824709491600873520394816
567085323315177661152862119950150798479374508570574002992135478614669402
960432542151905877553526733139925401296742051375413954911168510280798423
...

FAQ
===
- I compute pa('1+1'), and nothing happens. Why?
* Computing 1+1 gives 2.0000000000000000..., as we all know. Computing
  this number to any precision cannot say whether it should be printed
  as 1.9999999999999999 or 2.0000000000000000, since the next digit
  may influence the result. Therefore, pa() cannot even print the first
  digit.
- Why does fact(r(1,1000)) take forever?
* It doesn't take forever, but it takes a very long time. This is
  because computing factorials in general is very complicated. As a rule
  of thumb, computing factorials small factorials (<10000) to large
  precisions (more than 1000) digits takes very long.
- How can ceil(x) be defined as floor(x)+1?
* No, I'm not giving it away that easy. Use the hint in the source
  code. If you still don't see it, mail me and I'll tell you.
- What the heck is a Bernoulli number, and what do they have to do
  with factorials?
* Read for example Knuth, volume 1, chapter 1.2.9.
- What does hex() do for Real numbers?
* It gives a representation of the number in "floating point binary",
  represented with a hex mantissa and exponent. There are people who
  actually use this :-)
- Why are there methods for the mathematical functions?
* It allows you to say x.sin() instead of sin(x). It is useful for
  people who use the NumPy extension; see
  <http://www.sls.lcs.mit.edu/jjh/numpy/>

Implementation comments
=======================
I started this project to learn about Python programming. I decided to
do something quite challenging, so I could see if learning the language
distracted me from the real problem. It didn't, and I was quite
impressed.
This ease of use allowed me to write something much more advances than I
expected. I would never have guessed that I would ever write a factorial
routine with arbitrary precision!
To me, this project showed that Python is indeed capable of serious
programming. One thing I did not do, was user interfaces: everything
here is completely machine independent. Therefore, I cannot comment on
the portability aspect.
Of course, I like the syntax of Python, but I especially love its
semantics:
- The cooperation between the operators >>, /, % and & is amazing.
  I do not know of any other language that does this:
    x>>2 = x/4 for all x
    x&3 = x%4 for all x
    x&-2 = x-x%2 for all x
- The "and" and "or" operators are very nice, and much more useful than
  their C counterparts. (But this is not new, of course.)
- The name semantics are so clean that I start to hate C.
- The exceptions.
- Debugger and profiler written in Python.
These things combined made me a total Python convert. You have to force
me to use anything else for programming :-)
BTW, if you like to be amazed, use the profiler to see what happens if
you compute fact(r(1,200)). It gives you an idea of the complexity of
the underlying mathematics.

There were also a few things that I didn't like about python, especially
about its implementation:
- the fact that you get an "outrageous left shift count" error on big
  left shifts, for example 1L<<100000. Why is this? I know the result is
  big, but I really _want_ that. To be precise, there are no less than
  36 places in the program where I want to be able to do that. I had to
  write a function that does it, and use that instead. At one point, I
  even had to inline the function for speed reasons.
- the difficulties computing the number of significant bits of a long
  or int. Check out my function bitsize. It is quite dirty, although it
  will probably never fail, and it will certainly not give wrong
  answers. It took me a long time to write it. There are also four
  places in the program I had to do a weird trick not to have to compute
  it for speed reasons.
- the time cost of a function call. I hope someone will invent a
  modification of the implementation of function call that does not
  change the beautiful semantics.
- O yes, and the long number routines could be faster. But I shouldn't
  complain :-)

Conclusion
==========
Thanks to Tim Peters and Christian Tismer for their interesting and
useful discussions, and their help and advice with making this library
into a mature product.

I hope you enjoy the library. If you like it, send me a postcard:
Jurjen N.E. Bos
Deukelven 48
1963 SR Heemskerk
The Netherlands

I'd like to know what you use it for.
- Jurjen N.E. Bos (jurjen@q2c.nl)