spigot: an exact real calculator

This manual documents spigot, a command-line exact real calculator.

Chapter 1: What is spigot?

spigot is a calculating program. It supports the usual arithmetic operations, square and cube roots, trigonometric and exponential functions, and a few other special functions such as erf.

spigot differs from the average calculating program in that it is an exact real calculator. This means that it does not suffer from rounding errors; in principle, it can keep generating more and more digits of the number you asked for until it runs out of memory.

In particular, if you ask for a complex expression such as sin(sqrt(pi)), then most calculating systems would compute first pi, then sqrt(pi) and finally sin(sqrt(pi)), accumulating a rounding error at each step, so that the final result had a build-up of error and you would have to do some additional error analysis to decide how much of the output you could trust.

spigot, on the other hand, does not output any digit until it is sure that digit is correct, so if you ask for (say) 100 digits of sin(sqrt(pi)) then you can be sure they are the right 100 digits.

(The downside of doing this is that spigot is slow, compared to the more usual methods of computer arithmetic. You wouldn't want to use even its individual arithmetic operations and functions for any kind of bulk computation. And because it evaluates a complicated expression as a whole, it's especially unsuited to any iterative algorithm which involves looking at the first result you get and then deciding what further arithmetic to do on it: normal arithmetic systems can do the next step of the algorithm with only the cost of the extra operations, but spigot would have to re-evaluate the entire thing from the start every time. So spigot can be useful if you really do need a lot of digits, or if you're doing a computation prone to numerical error, or as a cross-check on other calculating systems, but it wouldn't be usable outside that kind of specialist niche.)

spigot can evaluate expressions starting from numbers it can construct internally (integers, rationals, π, etc), and it can also read a number from a file or a pipe. In the latter mode, it operates ‘on-line’, i.e. it writes output as it goes along, and once it has read enough digits of an input number to know some of the output, it will print it.

Chapter 2: Examples of using spigot

The simplest thing spigot can do is to compute well-known mathematical constants such as π. For example, try running this command:

$ spigot pi
3.1415926535897932384626433832795028841971693993751058209749445923078164
062862089986280348253421170679821480865132823066470938446095505822317253
594081284811174502841027019385211055596446229489549303819644288109756659
(and so on)

In this default mode, spigot generates unbounded output: it will continue writing digits of π to the terminal until you interrupt it (e.g. by pressing Ctrl-C).

As spigot generates more digits, it will slow gradually down and consume more and more memory to keep track of where it's got to, and one of time and memory will ultimately be the limit to how much data it can generate. But in principle, given unlimited time and memory, it could just keep on going.

spigot can also evaluate more complicated mathematical expressions. (Some will be much slower than others.) You could try some of these, for example:

$ spigot pi*e
$ spigot pi/e
$ spigot pi^e
$ spigot -- '-pi^2'
$ spigot 'sin(sqrt(pi))'
$ spigot 'exp(pi*sqrt(163))'

(The above command lines assume that you are invoking spigot from a POSIX shell, so that expressions containing parentheses need to be single-quoted in order for spigot to receive them unmodified without the shell interfering. On other shells or operating systems, quoting conventions may vary.)

Since spigot supports a variety of command-line options beginning with -, if your expression also begins with a - (as one of the examples above does), then you need to precede it with the special argument word --, which instructs spigot to treat further command-line arguments as an expression rather than an option.

See chapter 3 for a full description of the range of expressions you can give to spigot.

2.1 Digit limits and rounding

In its default mode, spigot keeps generating digits until you interrupt it, or until it determines that it has output the exactly correct answer and needs no more digits at all.

You can limit spigot to a finite number of digits by using the -d option. For example:

$ spigot -d40 pi
3.1415926535897932384626433832795028841971

(The argument to -d counts decimal places, rather than significant figures. So here, there are 40 digits after the decimal point.)

By default, spigot does no rounding, i.e. it simply truncates the full decimal expansion of the number you asked for. (This is equivalent to rounding towards zero in all cases.)

If you want to change that, there's a variety of rounding-mode options, detailed in section 4.2. For example, you can ask spigot to round the final digit towards whichever value is nearer to the true result using --rn:

$ spigot --rn -d40 pi
3.1415926535897932384626433832795028841972

(This example rounds up, because the next digit following the 1 is a 6.)

2.2 Generating output in different bases

By default, spigot's output is in base 10. You can use the -b option to output in another base instead:

$ spigot -b2 -d40 pi
11.0010010000111111011010101000100010000101
$ spigot -b3 -d40 pi
10.0102110122220102110021111102212222201112
$ spigot -b16 -d40 pi
3.243f6a8885a308d313198a2e03707344a4093822
$ spigot -b36 -d40 pi
3.53i5ab8p5fsa5jhk72i8asc47wwzlacljj9zn98l

For bases larger than 10, letters of the alphabet are used as additional digits. If you use the -B option in place of -b, those letters will be displayed in upper case:

$ spigot -B16 -d40 pi
3.243F6A8885A308D313198A2E03707344A4093822

2.3 Generating continued fractions of numbers

As well as generating numbers in ordinary base notation, spigot can also generate them in the form of continued fractions. You can generate the continued fraction terms of a number, one per line, using the -c option:

$ spigot -c pi
3
7
15
1
292
1
(and so on)

(For those unfamiliar with continued fractions, this is equivalent to saying that π is equal to 3 + 1 / (7 + 1 / (15 + 1 / ...)).)

This one-per-line format is useful for applications which are consuming the numbers automatically, but for readability, you might prefer to add -l to replace the newlines with commas, so that you can see many more terms at once:

$ spigot -c -l pi
3;7,15,1,292,1,1,1,2,1,3,1,14,2,1,1,2,2,2,2,1,84,2,1,1,15,3,13,1,4,2,6,6
,99,1,2,2,6,3,5,1,1,6,8,1,7,1,2,3,7,1,2,1,1,12,1,1,1,3,1,1,8,1,1,2,1,6,1
,1,5,2,2,3,1,2,4,4,16,1,161,45,1,22,1,2,2,1,4,1,2,24,1,2,1,3,1,2,1,1,10,
(and so on)

(The semicolon emphasises that the first number is the integer part.)

Finally, you can use -C to output the continued fraction convergents (that is, the rational numbers obtained by evaluating successive truncations of the continued fraction, which are the best approximations to the target number by a particular metric):

$ spigot -C pi
3/1
22/7
333/106
355/113
103993/33102
104348/33215
(and so on)

In any of these modes, the -d option allows you to limit the number of continued fraction terms or convergents that spigot outputs. Again, the integer part is not counted, so that for example -d3 gives you three terms or convergents as well as the integer one:

$ spigot -c -d3 pi
3
7
15
1
$ spigot -C -d3 pi
3/1
22/7
333/106
355/113

2.4 Reading input numbers from a file or another program

As well as applying its collection of mathematical operators and functions to real numbers that it has derived itself from first principles, spigot can also apply the same operations to numbers received as input from elsewhere, which makes it at least theoretically able to operate on any real number.

To begin with, you can write a number into a file in ordinary decimal notation, and use a special syntax in the spigot expression language to read from that file:

$ echo 1.12358132134558914423337761098715972584418167651094617711 > z
$ spigot -d30 'sin( base10file:z )'
0.901654985409730168388244848164

You might wonder why spigot needs a feature like this: since the file contains only a finite amount of data, you surely could just have pasted the same number on to spigot's command line directly. The difference is that if you write a decimal number on spigot's command line, spigot will assume it is exact, i.e. that the precise number you meant is the rational number with that terminating decimal expansion. But if you read from a file like this, spigot will treat it as a prefix of some longer decimal expansion, and if it reaches the end of the file then it will report a special error saying that it can't compute any more output digits because it ran out of input. So this permits you to compute with imperfectly known inputs, without the risk of generating any incorrect digits of output due to the input imprecision.

You can also tell spigot to read from a file in number bases other than 10, or in continued fraction format. See section 3.5 for full details of all the available options.

But spigot can also read numeric data from its own standard input (or, if it's compiled on Unix, a numbered file descriptor of your choice), which can point at an unbounded data source such as a pipe. So if you invoke spigot with its standard input pointing at such a pipe, with a program at the other end of the pipe prepared to generate as many digits of the number as requested, then you need only use another special piece of expression syntax to tell spigot to read a real number from it, as far as is necessary.

For example, let's generate the Champernowne constant, which you get by concatenating all the positive integers in order immediately after the decimal point, with no leading zeroes. Here's a piece of Perl which generates that:

$ perl -e 'print"0.";print while++$_'
0.1234567891011121314151617181920212223242526272829303132333435363738394
041424344454647484950515253545556575859606162636465666768697071727374757
677787980818283848586878889909192939495969798991001011021031041051061071
(and so on)

Now we can use spigot to compute with this number (or any other input real), by piping that script's output into spigot, and telling spigot to expect to read a base-10 number from its standard input. As a really simple example, let's ask for spigot to output the same number, but in continued fraction form:

$ perl -e 'print"0.";print while++$_' | spigot -c -l base10stdin
0;8,9,1,149083,1,1,1,4,1,1,1,3,4,1,1,1,15,457540111391031076483646628242
956118599603939710457555000662004393090262659256314937953207747128656313
8641209375503552094607183089984575801469863148833592141783010987,6,1,1,2
1,1,9,1,1,2,3,1,7,2,1,83,1,156,4,58,8,54,4457353800911178833959067671634
293788437292958096324947188556700067877659324583930837874799958333344441
(and so on)

(Apart from it being very easy to generate, another reason I chose this example constant is because it has interestingly large numbers in its continued fraction – the huge number spanning the first three lines occurs because there's an extremely good rational approximation to the constant at the point where numbers change from 2 to 3 digits. The number starting on the fourth line is the corresponding term at the 3/4 digit boundary, and goes on for far longer!)

You can also use the expression base10stdin as part of a larger expression, so that you could compute mathematical functions of an input number, or add it to things (perhaps another input pipe, if you execute spigot in such a way as to assign multiple input fds).

Chapter 3: spigot's expression language

This section gives a full specification of the language spigot uses for mathematical expressions.

spigot's expression language is semantically simpler than most such languages. In deference to the limitations of efficient computing, most expression languages have to carefully distinguish a collection of numeric types representing different subsets of the real numbers with different tradeoffs between range, precision and performance. But spigot, because its whole purpose is to deal in exact real numbers regardless of performance, does not have to do this. Every expression and subexpression in spigot has the same type: real.

So you don't have to remember to avoid writing 1/2 for a half (which in, say, gnuplot, would be interpreted as integer division and yield a rounded-down quotient of zero), or include any explicit syntax to change between formats as appropriate. Just write what looks like the obvious thing, and it will probably mean what you wanted.

3.1 Operators

spigot's expression language provides the following mathematical operators:

The priorities and associativity of these operators, from lowest to highest, are as follows.

3.2 Functions

As well as infix and prefix operators, spigot also supports a range of mathematical functions, all written in the style function(argument) or function(argument1,argument2). The following functions are provided:

Some of these functions are considerably slower than others. The inverses of erf and Phi and the W functions, in particular, are implemented by laborious interval-bisection and are very slow indeed.

3.3 Numeric literals

You can write actual numbers in spigot expressions, of course. By default a string of digits, or a string of digits with a decimal point somewhere in it, will be interpreted in decimal.

spigot supports C-style scientific notation, in which a decimal number (with or without a decimal point) is suffixed with e followed by an optional + or - and then another decimal integer; this denotes the first number times 10 to the power of the second. For example, 1.2 means 1.2, but 1.2e10 (or 1.2e+10) means 12000000000, and 1.2e-10 means 0.00000000012.

spigot also supports hex numbers, by prefixing 0x or 0X to a number, e.g. 0xabc means 2748. Hex numbers can be suffixed with an exponent part similar to the decimal one described above, only using p in place of e as the separator character; this means the hex number should be multiplied by the specified power of two. For example, 0x1.2 means 1.125 (one and an eighth); then 0x1.2p1 (or 0x1.2p+1) means twice that, i.e. 2.25, and 0x1.2p-1 means half of it, i.e. 0.5625.

If you want to input numbers in bases other than 10 and 16, you can prefix a numeric literal with base2:, base3:, …, base36:. For example, you could write base2:11.011, or base36:ZYX.WVU. No exponent suffix is recognised in this mode, because there's no solid convention for what it should look like or even what should be raised to the specified power (since the usual exponent-bases for bases 10 and 16 are 10 and 2 respectively). If you want to specify a number in scientific notation and an arbitrary base, write the exponent as an explicit power of something, e.g. base7:1.234 * 7^13.

A final way you can specify literal numbers to spigot is by specifying them as a hex number interpreted according to IEEE 754. You do this by writing ieee: followed by a number of hex digits, which must be exactly 4, 8, 16 or 32. These lengths correspond to the following formats:

An IEEE hex bit pattern can be followed by a . and further hex digits, in which case those digits are treated as an extension to the mantissa field. For example, ieee:3f800000 represents 1 (in IEEE single precision); ieee:3f800001 is the next representable number, namely 1+2^-23; and spigot will interpret ieee:3f800000.8 as the number half way between those two, i.e. 1+2^-24.

Infinities and NaNs are not permitted as IEEE format input.

3.4 Built-in constants

spigot also supports a few mathematical constants under built-in names.

(All of these numbers can be generated by other methods, such as 4*atan(1) or exp(1), but it's more convenient to have them available as predefined constants.)

3.5 Numbers read from files and file descriptors

spigot can read a number from a file in various formats, and use it as input to its range of mathematical operations and functions (or just output it directly in a different format).

To read a number from a file in ordinary decimal, write base10file: followed by the file name. Any sequence of non-space characters following the : will be assumed to be the file name – even if they're punctuation or delimiters, e.g. if you write sin(base10file:myfile) then spigot will interpret the ) as part of the file name, and (even if a file with that strange name does exist) will complain that the expression is incomplete.

If you need to read a file that does have a space in its name, there's a quoting syntax to permit it. If the first character after the : is either ' or ", then spigot will look for the next occurrence of the same character, and treat all characters in between as the file name. A doubled quote character will be treated as a literal quote. For example:

If your file contains the number in a base other than 10, you can prefix its name with base2file:, base3file:, …, base36file:. Bases larger than 36 are not supported, because after that the letters of the alphabet run out.

When spigot reads a file in base notation, it ignores newlines and white space interspersed between the digits.

You can also prefix a file name with cfracfile:, in which case spigot will expect it to contain a sequence of continued fraction coefficients (written in decimal). The coefficients can be separated by any non-digit characters you like (including newlines, like spigot's -c output mode, or ; and , like spigot's -c -l output mode); a minus sign is permitted before the very first coefficient to indicate a negative number, but ignored thereafter.

As mentioned in section 2.4, spigot will treat numbers read from files using any of the above keywords as if they are not an exact representation of a rational number, but a prefix of an infinitely long expansion of an irrational. So spigot will compute the output value as far as the input precision permits, and if it reaches the point where it can't generate any further output without more input precision, it will stop, print an error message indicating which input file ran out first, and return a special exit status 2.

If you really did want the contents of a file to be interpreted as the exact terminating expansion of a rational number, you can use the keywords base10xfile: (or likewise with bases other than 10) and cfracxfile, where the x stands for ‘exact’. The advantage of doing this rather than just specifying your expansion on spigot's command line directly is performance: if your input file is extremely large, then spigot will only have to read enough of it to generate the amount of output you asked for, whereas if you put the same number directly on the command line then spigot would have to parse all of it before starting.

spigot can also read from its standard input in the same way that it would read from a file, if you write an expression containing the keyword base10stdin (or another base, as above) or cfracstdin. For example, if you write base10stdin, then spigot will expect to read a number in base 10 from its own standard input (so you might pipe the output of another program into it).

If spigot has been compiled with the feature enabled (which it typically will have been on Unix systems; use spigot --version to check), you can also read from a numbered file descriptor of your choice in place of standard input, by writing base10fd:n (or another base, as above) or cfracfd:n. If you do this, spigot will expect its own file descriptor numbered n to already be open and readable, and will expect to read a number from there in the appropriate format. For example, base10fd:0 is another way of writing base10stdin; but, more interestingly, you could invoke spigot with a shell redirection operator such as 3<filename (where filename, in turn, could point at something interesting such as a named pipe, not necessarily a regular file), and then write (say) base10fd:3 to refer to that number.

Numbers read from stdin or numbered file descriptors are always treated as exact if the file descriptor signals EOF. (The expectation is that the whole point of using a file descriptor is so that you can connect it to a program that just keeps generating data.)

(If spigot is told to read from one or more file descriptors which refer to Unix terminal devices, you can use the -T option to set those terminals into raw mode. See section 4.3 for more detail.)

If you want spigot to run in a mode where it will refuse to read from files or file descriptors at all (for example, as a mild safety measure if input expressions are coming from an untrusted source), you can use the --safe option.

3.6 User-defined variables and functions

Sometimes, it's convenient to define your own variables or functions in a spigot expression. For example, suppose you want to evaluate a polynomial at some particular value. You could write something like this:

$ spigot 'sin(1.1)^3 - 2*sin(1.1)^2 + 5*sin(1.1) - 7'

but you probably got annoyed at the repetitiveness just reading that. So instead you could define a variable to have the repeated value:

$ spigot 'let x=sin(1.1) in x^3 - 2*x^2 + 5*x - 7'

or alternatively define a function to represent the polynomial:

$ spigot 'let f(x) = x^3 - 2*x^2 + 5*x - 7 in f(sin(1.1))'

In this example, it's more or less a matter of taste which of those you prefer, but the function syntax becomes more useful if you want to evaluate the function multiple times, e.g.

$ spigot 'let f(x) = x^3 - 2*x^2 + 5*x - 7 in f(f(sin(1.1)))'

The full syntax of the let statement is as follows:

let definitions in expression

where definitions is a comma-separated list of definitions, each in one of the following forms:

variable = expression
function(parameters) = expression

In the latter case, parameters is in turn a comma-separated list of identifiers, and expression is allowed to refer to those identifiers as if they were variables.

Each definition comes into scope as soon as it is complete, but not before. So a sequence of definitions can refer back to each other:

$ spigot 'let x=pi/2, y=x^3, z=exp(y) in sin(z)'
$ spigot 'let f(x)=x+1, g(x)=sin(f(x)) in g(3)'

but a definition cannot refer to itself – in particular, functions cannot be defined recursively. (This is a fundamental limitation of spigot's evaluation system.)

Chapter 4: Full list of spigot's command line options

4.1 Output format options

spigot supports the following options to control the format in which numbers are output:

4.2 Rounding mode options

When spigot is given a digit limit via -d, and is not in continued fraction output mode, it defaults to printing only digits from the real expansion of the number, i.e. it always truncates toward zero and never rounds up.

You can make spigot round the output at the last digit position in various modes, using one of the following options:

If spigot is told to print a limited number of digits via --printf rather than -d, then it defaults to rounding to nearest and tiebreaking to even (i.e. equivalent to --rn), because that's the default behaviour of the C printf function which spigot is imitating. This is done by making the --printf option behave as if --rn had been specified alongside it. So if you want to use a different rounding mode in printf-format output, you will need to specify the rounding-mode option second, or else the implicit one in the --printf will override it. For example, --printf --ru will work, but --ru --printf will revert to the --rn implied by the --printf option.

4.3 Miscellaneous options

Here are some options which didn't fit nicely into the above categories.

Chapter 5: Hazards to computation

So far, this manual has more or less portrayed spigot as a magic device that can compute anything it likes. Perhaps some expressions are evaluated more slowly than others, but it'll get there in the end.

Well … not quite, sorry.

spigot's core algorithm works by finding an interval of rationals bracketing the target number, and gradually narrowing that interval further and further as more information is received or computed. Its various output modes work by waiting until the interval has narrowed to the point that the next digit is uniquely determined, and then writing that digit out.

The problem is: suppose that in order to work out the right output value, spigot has to decide whether a number falls on one side or the other of some boundary value. (This can happen in several different kinds of situation, which I categorise below.) It will do this by narrowing its interval until the boundary value doesn't fall inside it any more, and then it knows which side the number is on.

So the closer to the boundary the number is, the longer spigot will take before its interval narrows enough to work out which side it's on. And here's the problem: if the number it's trying to compute is exactly on the boundary value, then spigot may very well not ever notice – it will just keep narrowing its interval further and further, and the boundary value will still be inside the interval no matter how much it does that, so it will never manage to make a decision.

It's not standard terminology in exact real computation, but I've tended to call this sort of problem an ‘exactness hazard’. The general rule of thumb is that spigot is good at generating difficult output – complicated and fiddly irrational numbers – but can get confused if the answer to any computation is too easy, in particular if it's zero, or an integer, or has a finite number of digits in the output representation. So if you're computing any kind of long and fiddly expression, you need to watch out for too-simple exact numbers cropping up anywhere in it, because if you're unlucky, they can cause problems that prevent the expression as a whole from being computed.

In the following subsections, I list some possible effects of this kind of problem, and go into more detail about what numbers can trigger it.

5.1 Exactness hazards on output

One obvious situation in which spigot has to determine which side of a boundary a number falls is during final output, if the entire result of the computation has a terminating representation in the output base or number system. For example, suppose you run this command:

$ spigot 'sin(asin(0.12345))'

Obviously, we can see that the true answer is 0.12345 exactly. But spigot, with no symbolic algebra system, can't see that. So what can it do? It will manage to print ‘0.1234’ easily enough, but then its interval of possible output values will narrow for ever without allowing it to be certain of the next digit – no matter how much the interval narrows, it will still have one end looking like ‘0.123449999…’ and the other like ‘0.123450000…’, and any further work will just add more 9s to the former and more 0s to the latter, so there will never be a point at which it can be sure of that fifth digit.

spigot mitigates this situation by means of its tentative output feature. (See section 4.3 for options to control this.) What will actually happen if you run the above command, at least on a default Unix terminal, is that spigot will print ‘0.1234’ in the normal way, then follow that by printing the next digit ‘5in red. The red text indicates that the output is tentative, i.e. it might be retracted if more information comes along. So you might see a display along these lines:

$ spigot 'sin(asin(0.12345))'
0.12345 (10^-205)

This should be read as spigot saying: ‘The number definitely starts with 0.1234, and it looks as if it's exactly 0.12345, but I've only looked as far as 205 digits beyond that point, so there might still be surprises as I go further.’ As spigot computes more and more digits without finding a divergence from 0.12345, the power of 10 in the suffix will keep growing. And the red tentative text can be retracted if further information is discovered: if you had instead asked spigot to compute a number that was not exactly equal to 0.12345 but merely very very close to it, then eventually spigot would find that out, delete the red text, and print normal text in its place once it knew what it should be.

So if you see this kind of tentative output from spigot, that's your cue to have a closer look at the expression you asked it to compute, and see if you can see some mathematical reason why the answer is exactly what spigot is suggesting it might be. In this example case above, of course, that reason is pretty obvious – that's what asin means.

One exceptional case in which this kind of thing will work is that spigot makes a special case for numbers it knows from first principles to be rational. So while the above expression with answer 0.12345 ran into trouble, the following much simpler thing does work sensibly, and just prints the right answer immediately:

$ spigot 12345/100000
0.12345

because this time, spigot can remember that the number being output is the same as the known-rational number it got as input, so it can spot that it's an exact digit boundary and output it correctly.

(The special case for rationals is the only reason why section 4.2 needs to have all those different options for methods of tie-breaking in round-to-nearest mode – if no terminating representation could ever be successfully generated, then no tie-breaking options would be required, because spigot could never recognise a tie anyway.)

5.2 Exactness hazards in subexpressions

A more difficult case arises when spigot has the same problem of locating a number on one side or another of a boundary value, but in an intermediate result in a complex expression. For example, suppose you run one of these commands:

$ spigot 'tan(pi/2)'
$ spigot 'floor(sin(pi))'

In each of these situations, spigot will hang completely, and never manage to print any output at all. This is because, in each case, the value of a subexpression (respectively pi/2 and sin(pi)) is exactly on a point of discontinuity of the function it's being fed to next. In order to even start outputting the value of tan(x), spigot first needs to know whether it's positive or negative – and it can't find that out until it knows whether x is on one side or the other of π/2, which in this case it will never manage to decide. Similarly with floor(sin(pi)): to compute floor(x), spigot needs to know whether x < 0 or x >= 0, and again, it can't work that out in the case where x is exactly 0.

In this kind of situation, there really is nothing spigot can do but hang; even tentative output can't mitigate the situation.

As in the previous section, if the number is obviously a rational, then spigot can do better. If you actually ask it for floor(0), then it will handle that fine. But in the example above, the input to ‘floor’ is non-obviously zero, in the sense that in order to know it spigot would have to actually think about maths rather than just computing.

Therefore, if this happens to you, you can sometimes debug it by pulling out subexpressions and trying to evaluate them on their own; when you find one with a simple rational answer (probably signalled in turn by some red tentative output), see if you can prove that to be exactly the number you wanted, and substitute it in. For example:

$ spigot 'floor(sin(pi))'
(spigot hangs)
$ spigot 'sin(pi)'
0 (10^-951)
$ spigot 'floor(0)'
0

Here, the user observes spigot hanging, and suspects an exactness hazard somewhere in their expression. The only interesting subexpression is sin(pi), so the user evaluates that on its own, and gets back tentative output suggesting that the answer is exactly zero. The user thinks about it, realises that of course the answer is exactly zero, and simplifies the original expression by substituting a literal zero in place of the needlessly complicated sin(pi).

(Of course, in this simple example case, once the user had spotted the exact zero it hardly needed to ask spigot for floor of it! But if something more complicated were being done after the difficult point, it might well be easiest to substitute in the simple answer and try again with spigot.)

However, this debugging technique can only work in cases where the troublesome intermediate result is rational, because it's only rationals for which spigot has a special-case handler. There would be no equivalent way to debug the tan(pi/2) example above – spigot actually cannot recognise an input value to ‘tan’ as being an odd multiple of pi/2.

(This is the reason, mentioned in section 3.2, why the tan function can never return an error, even if you try to pass it an invalid input value. Because all the values where tan has no finite answer are irrational, so spigot will always run into this exactness hazard which prevents it attempting to compute an infinity.)

5.3 Internal exactness hazards

A third place where this same problem can occur, potentially, is internal to spigot itself. spigot's internal implementation of its various functions will often need to locate the input value on one or other side of some boundary value, and if spigot does that without sufficient caution, then it might hang forever when passed exactly the boundary value.

This should never happen. Exactness hazards internal to spigot itself are bugs. I've found and fixed all the ones I know of, but in case any more turn up, this section demonstrates how to recognise one, so that you can report it as a bug.

Here's an example of an exactness hazard that used to exist. spigot computes the power function a^b by computing exp(b*log(a)), in most cases. But if a < 0, then the right answer will have magnitude exp(b*log(abs(a))), but could be positive, negative or undefined depending on what kind of number b is. So the implementation of pow, having first ruled out a variety of special cases that don't require taking logs at all, would test the sign of a, and run straight into a hazard if a was non-obviously zero. For example, this could happen:

$ spigot '0^pi'
0
$ spigot 'sin(pi)'
0 (10^-951)
$ spigot 'sin(pi)^pi'
(spigot used to just hang)

If spigot knew that a was zero, then it could easily decide that zero to any power is zero. And it's at least capable of realising that sin(pi) is arbitrarily close to zero – but it wasn't able to produce even tentative output for sin(pi)^pi, because that internal sign test of a hung forever trying to decide which side of zero sin(pi) fell on.

Now the bug is fixed, and the last of those commands produces perfectly good tentative output:

$ spigot 'sin(pi)^pi'
0 (10^-807)

So if you suspect you're in a situation like this, where some spigot computation hangs (without even tentative output) for no reason you can see, try this diagnostic procedure:

For example, you might have reported the above example as follows:

The following spigot command hangs without even tentative output:

$ spigot 'sin(pi)^pi'

I think this is an internal exactness hazard, because spigot can evaluate the two operands to the power function safely (the former produces tentative output and the latter produces definite output):

$ spigot 'sin(pi)'
$ spigot 'pi'

And the power function is mathematically continuous at the point (0,π), so I think spigot should not hang in this case.

That example report contains everything needed to check the facts, reproduce the failure, and decide whether it is indeed a bug.

5.4 Design limitations

One last class of hang in spigot is due not to a bug in the implementation, but a lack of generality in the fundamental design, which couldn't really be fixed without writing a totally different program from scratch.

As mentioned above, spigot works by narrowing an interval of rationals bracketing the number. But sometimes you know that the number is somewhere within one of two intervals of rationals, and can keep narrowing both of those intervals without ever working out which one contains the number. And sometimes, depending on what you want to do with the number next, that ought in principle to be good enough.

For example, it's possible to imagine a system that could cope with this:

$ spigot 'abs(atan2(sin(pi),-1))'

because the more you compute the input values to atan2, the closer you know the output is to one of +π and -π, even though you still don't know which – and since both of those values come out the same once they've been through the subsequent abs function, you could imagine the system being clever enough to cope with the expression as a whole. But spigot's design fundamentally is not.

Chapter 6: Output and return values of spigot

Whatever output format you have chosen (see section 4.1), spigot will write data in the specified format to its standard output channel, i.e. Unix file descriptor 1.

If spigot successfully generates all the precision you asked for, it will print a terminating newline (if the output is any of the one-line types and you did not specify -n; see section 4.3), and then it will terminate with exit status 0 (the usual signal for success).

If spigot cannot generate enough output precision because it was reading a number from an input file and that input file ran out, it will still print its terminating newline (subject to the same conditions as above) on standard output; it will also write an error message to its standard error channel (Unix file descriptor 2) indicating which input file ran out first (in case there was more than one), and terminate with exit status 2.

If spigot encounters any other failure to evaluate the expression, it will print an error message to its standard error, and terminate with exit status 1.

Appendix A: References

spigot's central algorithm is derived from the paper ‘An unbounded spigot algorithm for the digits of π’, by Jeremy Gibbons. (American Mathematical Monthly, 113(4):318-328, 2006.) At the time of writing this, Jeremy Gibbons's web site has a PDF copy of the paper available.

Gibbons's algorithm is readily adapted to produce a continued fraction representation as output in place of base notation, and to accept numbers other than π as input. It forms the core of everything spigot does, and hence it seemed appropriate to name the entire program after it.

(Perhaps ironically, one part of Gibbons's algorithm that spigot does not use any more is his representation of π! It turned out that once spigot had grown the ability to compute arctangents and do arithmetic, a representation based on Machin's formula was much faster.)

My algorithm for basic arithmetic is a sort of hybrid of Gibbons's spigot algorithm and William Gosper's algorithm for doing basic arithmetic on continued fractions (which I formerly used unmodified, but had to change it because using a continued fraction representation introduces exactness hazards). Gosper's algorithm comes originally from ‘HAKMEM’: Memo 239, Artificial Intelligence Laboratory, Massachusetts Institute of Technology, Cambridge, Mass., 1972. HAKMEM can also be downloaded at the time of writing this. Also, I haven't checked with great care, but the algorithm I constructed by combining both ideas may very well be the same as the one used by Imperial College's exact real arithmetic system (see the link below), in which case, they definitely had the same idea before I did.

Other algorithms used in this program came from my own head (though I'm not aware that any of them came from my head first).

spigot is not the only exact real calculator around, and was not really written with the aim of filling any particular gap in the available software. (I mostly wrote it for the fun of playing with Gibbons's elegant algorithm, and only realised some time later that it had become useful enough to be worth publishing.) Hence, if spigot is of interest to you, other implementations might also be of interest. Here are some that I managed to find:

Appendix B: Copyright statement and licence

spigot is copyright 2007-2014 Simon Tatham. All rights reserved.

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

Appendix C: Unix man page for spigot

C.1 NAME

spigot – command-line exact real calculator

C.2 SYNOPSIS

spigot [ options ] expression

C.3 DESCRIPTION

spigot is an exact real calculator: that is, you give it a mathematical expression to evaluate, and it computes it to any desired precision, by default simply printing digits to standard output until it is interrupted.

spigot provides command-line options to control the format of the output, restrict it to a specified number of digits, and apply rounding at the end of those digits. It can produce output in any base between 2 and 36 (after that it runs out of digit characters), or as a continued fraction, and it can read input numbers from files in any of those formats as well.

This man page gives only a brief summary of spigot's functionality. For full detail, you should read the main manual spigot.html; if that is not installed on your system, you can find it on the web at

https://www.chiark.greenend.org.uk/~sgtatham/spigot/spigot.html

C.4 OPTIONS

The following options control spigot's basic output format:

-b base, -B base
Output the number in base base, which must be an integer between 2 and 36 inclusive. Digits above 9 are represented by lower case or upper case letters, for the options -b and -B respectively. The default is -b 10.
-c
Output the number as a list of continued fraction coefficients, as decimal integers, by default one per line.
-C
Output the number's continued fraction convergents, one per line, in the form of two decimal integers with a / between them.
-R
Output the number's value as a rational, in the form of two decimal integers with a / between them, or just one decimal integer if the number is a rational. If spigot does not know the number to be rational immediately, it will start evaluating it to see if it turns out rational later, so if it is not rational then spigot will compute for ever.
-S, -D, -Q, -H
Output the number as a hex representation of an IEEE 754 bit pattern, in 32-bit single precision, 64-bit double, 128-bit quad or 16-bit half precision respectively. If that representation is not exact, a decimal point will be printed followed by further mantissa digits.
--printf format, --printf=format
Format the number in the same way that printf(3) would, given the formatting directive format. format must begin with a % and end with the associated conversion specifier, which must be a floating-point one (one of efgaEFGA).

The following options modify the details of those output formats:

-d limit
Limit the amount of data output. In -b mode, no more than limit digits after the decimal point are printed. In -c or -C mode, no more than limit continued fraction coefficients or convergents are printed, not counting the initial one representing the number's integer part. In the IEEE 754 output modes, no more than limit additional bits of precision are generated after the end of the official mantissa. limit may be negative.
-l
In -c mode, output continued fraction terms all on one line, separated by a ; after the first term and , after each subsequent term.
-w min-int-digits
In -b mode, output at least min-int-digits of the number's integer part, by printing leading zeroes if necessary.
-s exponent-base (or -s b or -sb)
In -b mode, output a prefix which is the largest power of exponent-base less than or equal to the absolute value of the number, and then print the number after dividing by that power. The special string ‘b’ for exponent-base means to use the same base for the exponent as the digit base used for the main output.
--nibble
In --printf mode with the ‘a’ or ‘A’ conversion specifier, choose the output exponent to always be a multiple of 4, instead of the default behaviour of choosing it as large as possible.
-n
In any mode where spigot prints output on a single line, suppress the usual trailing newline if spigot's output terminates.

The following options control rounding, when spigot's output is limited by the -d option. (Rounding does not occur in continued fraction modes.)

--rz
Round towards zero. This is the default.
--ri
Round away from zero.
--ru
Round up (towards positive infinity).
--rd
Round down (towards negative infinity).
--rn, --rne
Round to nearest, breaking ties toward an even last digit.
--rno
Round to nearest, breaking ties toward an odd last digit.
--rnz, --rni, --rnu, --rnd
Round to nearest, breaking ties as if rounding via --rz, --ri, --ru or --rd respectively.

Miscellaneous options:

--tentative=state
Control the printing of ‘tentative output’. Tentative output is printed when spigot does not know for sure what the next digit of the number is because it's starting to look as if it's exactly on a digit boundary. Tentative output is in red, and followed by an indication of about how many digits spigot has examined beyond that point (i.e. how close to exact that digit is known to be); spigot will retract it later if it finds out something definite.

state can be ‘on’, ‘off’ or ‘auto’. ‘auto’ is the default, and means that spigot should only print tentative output if its output is directed to a terminal device.

--safe
Disallow the various features of the expression syntax that tell spigot to read from arbitrary files or from its own file descriptors. Might be useful if an expression is coming from an untrusted source (although, in that situation, you should still beware of other risks such as the expression author forcing spigot into an endless loop).
-T
If instructed to read from a file descriptor which points to a terminal, put the terminal into raw mode (turning off ICANON and ECHO modes) while doing so.

C.5 EXPRESSIONS

spigot's expression language supports the following options, in order of priority from lowest to highest:

+ and -
Addition and subtraction. (Left-associative.)
*, /, %, mod, rem
Multiplication, division and remainder. (Left-associative.) % and mod are synonyms, which both return a remainder between 0 and the denominator; rem returns a remainder of either sign, with absolute value at most half that of the denominator, and ties broken by rounding to even in IEEE 754 style.
Unary - and +
Negation and no-op.
^, **
Power. (Right-associative.)
!
Factorial. (Unary suffix operator.)

You can define variables and functions of your own in subexpressions using the let expression, as follows:

let var=value in expression
Defines the name var to refer to the value of the expression value. The definition is in scope within expression, but not in any other parts of the spigot input.
let fn(params)=defn in expression
Defines the syntax fn(args) to refer to the expression defn with the arguments substituted in for the parameters. params must be a comma-separated list of identifiers; args is a comma-separated list of expressions.

A let expression can contain multiple definitions, separated by commas, e.g. ‘let x=1,y=2 in x+y’. Each definition is in scope for subsequent definitions, so you can write ‘let x=1,y=x+1 in’ expr. But definitions are not in scope for themselves; in particular, functions may not be recursive.

spigot also provides the following built-in functions:

sqrt, cbrt
Square and cube roots.
hypot, atan2 (two arguments, or more for hypot)
Rectangular to polar coordinate conversions: the hypotenuse function (square root of the sum of the squared arguments), and two-variable inverse tangent. hypot can also take a number of arguments other than two.
sin, cos, tan, asin, acos, atan
Trigonometric functions and their inverses.
sind, cosd, tand, asind, acosd, atand, atan2d
Trigonometric functions and their inverses, equivalent to the versions without ‘d’ on the end except that angles are measured in degrees.
sinh, cosh, tanh, asinh, acosh, atanh
Hyperbolic functions and their inverses.
exp, exp2, exp10, log, log2, log10
Exponential and logarithmic functions: raise e, 2 and 10 to a power, or take a log with the same three bases. You can also provide a base of your choice as a second argument to log.
expm1, log1p
Shorthands for exp(x)-1 and log(1+x).
pow (two arguments)
Synonym for the ^ operator.
gamma, tgamma, lgamma
Gamma function (gamma and tgamma are synonyms for this), and the log of the absolute value of the gamma function.
factorial
Synonym for the ! operator.
erf, erfc, Phi, norm
Error-function relatives: the error function itself, 1 minus the error function, and Phi and norm are synonyms for the cumulative normal distribution function.
erfinv, erfcinv, Phiinv, norminv
Inverses of the above error-function relatives.
W, Wn
The Lambert W function, i.e. the inverse of x exp(x). W is the branch with value at least -1, and Wn is the branch with value at most -1.
Ei, En (two arguments), E1, Ein
Exponential integrals, i.e. integrals of things like exp(x)/x. Ei(x) is the indefinite integral of exp(x)/x itself; En(n,x) (for non-negative integer n) is the result of integrating exp(-x)/x n times, flipping the sign each time; E1(x) is shorthand for En(1,x); and Ein(x) is the integral of (1-exp(-x))/x.
Li, li
Logarithmic integrals, i.e. integrals of 1/log(x). Li(x) and li(x) are both the indefinite integral of 1/log(x); only their constants differ, in that Li(2) and li(0) are each defined to be zero.
Si, si, Ci, Cin
Sine and cosine integrals, i.e. integrals of sin(x)/x and cos(x)/x. Si(x) and si(x) are both the indefinite integral of sin(x)/x, differing only in the constant: Si(0)=0, but si(x) has limit 0 as x tends to positive infinity. Ci(x) is the indefinite integral of cos(x)/x, also with limit 0 at positive infinity; Cin(x) is the indefinite integral of (1-cos(x))/x, with Cin(0)=0.
UFresnelS, UFresnelC, FresnelS, FresnelC
Fresnel integrals. UFresnelS and UFresnelC are the indefinite integrals of sin(x^2) and cos(x^2); FresnelS and FresnelC are the ‘normalised’ versions, i.e. integrals of sin(π x^2/2) and cos(π x^2/2). All are zero at the origin.
zeta
The Riemann zeta function (restricted to the real numbers).
abs
Absolute value.
sign
Sign: -1 for a negative input, +1 for a positive input, or 0 for a zero input (provided spigot can determine that it is zero without an exactness hazard).
ceil, floor
Ceiling and floor: smallest integer at least x, and largest integer at most x.
frac
Fractional part, i.e. x - floor(x).
algebraic (variable number of arguments)
Return a root of an arbitrary polynomial with integer coefficients. The first two arguments are the rational bounds of an interval to search, and the rest give the polynomial's coefficients, with constant term first.

spigot supports the following names for built-in constants:

pi, tau
The circle constant π, and the often more useful 2 π.
e
The base of natural logarithms.
phi
The golden ratio, (1+sqrt(5))/2.
eulergamma
The Euler-Mascheroni constant: the limiting difference between the sum and the integral of 1/n.
apery
Apéry's constant: the sum of the reciprocals of the cubes.

Numbers can be input in the following formats:

C.6 RETURN VALUE

spigot returns 0 if its output terminates (because the result is exact, or because it reached the specified -d limit) with no problems.

In case of a parse error, or an invalid operand to a function, or any other kind of fatal error, spigot returns 1.

If spigot is unable to generate output to the desired precision because more precision was needed from a number read from an input file using baseNfile: or cfracfile:, then spigot returns 2, and prints an error message indicating which input file (in case there was more than one) ran out first.

C.7 LIMITATIONS

Due to inherent limitations of its exact real arithmetic strategy, spigot is generally unable to recognise when a number it is computing is exactly equal to a specific boundary value.

One effect of this is that spigot will not behave as you'd like if the output number has a terminating representation in the selected base. For example, asking for sin(asin(0.12345)) will not be able to print 0.12345 and exit. Instead, spigot will get as far as printing ‘0.1234’, and then print tentative output (mentioned above) to indicate that it thinks the next digit might be exactly 5, but it will never reach a point where it's sure of that.

Another effect is that if you ask spigot to evaluate an expression in which an intermediate result is precisely on a point of discontinuity of the function it is passed to, then it may never manage to even start producing output. For example, spigot will hang completely if you ask it for floor(sin(pi)), since sin(pi) = 0 is a point of discontinuity of the floor function, and spigot will never be able to work out that the value of the input to floor is exactly zero, only that it seems to be closer and closer to zero the more it computes.

(An exception is numbers that spigot knows from first principles to be rational. For example, if you ask spigot to evaluate the simpler expressions ‘0.12345’ or ‘floor(0)’, it will print the complete output and terminate successfully, in both cases.)

C.8 LICENCE

spigot is free software, distributed under the MIT licence. Type ‘spigot --licence’ to see the full licence text.


[spigot version 20171116.baea20a]