On June 12, 2005 6:18 PM Tim Daly wrote:

you're right, of course. my mind is mush at the moment. i've been driving all week.

Yah, that'll do it to ya ... :) How's CMU?

Seriously, I think that there are some issues that should be dealt with here. There may very well be some situations that could make more efficient and/or effective use of Axiom's "infinite precision" floats. I expect that just as in the case of "infinite" integers, there may be some better ways and some worse ways of doing calculations with this objects.

For that matter, mathematically speaking just what is Axiom's Float type? Is it formally equivalent (though obviously not computationally equivalent) to FRAC INT? How does it relate to standard (IEEE ?) definitions of floating point? How does it differ mathematically from the reals, c.f. RealClosure, etc.

Regards, Bill Page.

Axiom has DoubleFloat (used to be called SmallFloat, hence sf.spad, I think), and Float (float.spad). DoubleFloat is IEEE or "native precision" (according to sf.spad). Any floating point system is only, mathematically speaking, a small subset, and not evenly distributed one for that, of the reals, and for that matter, of the rationals. It is definitely not FRAC INT, which is mathematically equivalent to the field of rational numbers.

In a fixed precision floating point system, there are only finitely many numbers, so it cannot be FRAC INT. In such a system, a floating point number has the form , where is the mantissa (assume finite precision), is the base (fixed), and is the exponent (assume also finite precision). In practice, there is an offset when and are stored so that is usually interpreted with a "base-point" and the most significant digit is normalized to be non-zero, and has an offset depending on the machine word length and how a word is divided between (and the signs, which for simplicity, we will assume are included in and ).

Axiom implementation uses the simplified representation of as it mathematically stands ( is still normalized to give a unique representation, but the base point is to the right of the least significant digit). The intrinsic properties do not depend on the base or the word length, so let's assume for discussion purposes.

Let's say the exponent has precision one (that is, only one digit plus one sign bit). So there are 19 exponents (from to 9). Say the mantissa m has precision 2 (that is, only two digits plus one sign bit). So there are 181 possible mantissas (from to , 0, 10 to 99). For simplicity, we will only discuss positive numbers. So one sees that there can only be positive numbers. But these are not evenly distributed. There are 90 positive numbers for ANY two successive exponents: so 90 positive numbers between 10 and 99, 90 between 100 and 990, and 90 between 1000 and 9900. It is clear that the accuracy of a real number by a floating point number in this system becomes less and less as the exponent goes up. This is the source of round-off errors.

This is basically the situation in DFLOAT (on Intel processors), where , has 53 bits stored in a 52 bit field (not including sign, note that in base 2, the most significant digit normalized must be 1, so no need to store it!) and has 11 bits (including sign, ranging from to which accounts for an offset of 53 because the base point location). This is exactly equivalent to the IEEE-754 standard for 64 bit floating point. The actual arithmetic operations are done via Lisp, which I assume calls the hardware.

In FLOAT, conceptually the infinite precision floating point system, is basically also finite precision floating point system, with the ability to increase precision as requested. However, this promotion is not automatic. In other words, every FLOAT number has its own precision, implicit in the length of its mantissa. In effect, FLOAT is the union of floating point systems with various precisions, but without exact arithmetic. So Tim is not entirely wrong in saying we don't have infinite precision floating point in the spirit of INT. What we have in FLOAT is arbitrary precision, not infinite precision (exact arithmetic). The real question for the developers is then:

Shall we improve FLOAT or shall we implement a new domain of infinite precision float (and if so, how, or rather, what are the specifications of such a domain)?

Example. Addition of two FLOAT with different orders of magnitude will not lead to a higher precision answer. What should be the right answer depends on your point of view: if the precision of the numbers are important, then this example is correctly handled in Axiom. If exact arithmetic with floating point is what you want, this is not correct.

If x and y are FLOAT, but y is much larger than x, say the exponent of y exceeds that of x (when both x and y are normalized) by the current precision of FLOAT, then x+y would be y, not x+y with precision extended (see routine plus in float.spad). Note that 68 bits is approximately 20 decimal digits. The display for y should have 9 trailing zeros after the decimal.

(1) -> bits()$Float

(1)

x:Float:=2^(-34)

(2)

y:Float:=2^35

(3)

x+y

(4)

By the way, in interpreting correctness of displayed answers for floating point computation, there is another source of error besides round off. For example,

x:Float:=sqrt(2)

(5)

bits(100)$Float

(6)

z:Float:=sqrt(2)

(7)

z - x

(8)

Most people would expect the answer of z-x to be 0.16887242 E-20 but this ignores the fact that the display is converted from an internal binary representation to a decimal one. During the conversion there is truncation error (think of 0.1 in base 3 converted to decimal 0.3333... with output at any finite precision). So x is not what it seems, and neither is z. Below, the constants are converted to binary internally before computation, at a higher precision than x (resp. z) to bring out the differences. We can now see that x is larger than z. So Axiom is correct and the expected answer is wrong.

x - 1.4142135623730950488

(9)

bits(120)$Float

(10)

z - 1.4142135623730950488016887242

(11)

Now that I'm awake the idea is coming back to me. What originally triggered the thought was that we need a way to compute an answer to a variable number of decimal places which could be expanded later. Thus we could compute the coefficients of a polynomial to 3 decimal places for display but expand them to 200 decimal places when we are evaluating the polynomial at a point.

The idea was to have a different representation that stored a closure of the computation, similar to the series mechanism, so we could "continue" to get more digits at a later time.

This raises the same kind of implementation issue that indefinite computation raises except that indefinites are symbolic and infinite precision floats are not.

The issue is that we want to keep the state of a computation in a form that we can expand. This involves storing functions and composition of functions. This would imply that floating point computations are no longer strictly numeric. So the product of two infinite floats (INFL) object would be stored in such a way to keep the original values as well as the result. Ideally we could stop a computation, store it as a closure, and restart it when we need more digits.

I have to give some thought to the more general case of a "closure" based mathematical object that captures the computations. As we push for more "symbolic" answers this issue keeps appearing.

Tim

From Kai Kaminski:

I just read your posts about infinite precision floats on the Axiom list and I recalled that I have seen something like this a while ago. I am not sure if this is what you are looking for but it sounds like it:

http://www.haible.de/bruno/MichaelStoll/reals.html

It is implemented in Common Lisp apparently.

yes, this is exactly the idea.

Tim

a:=2*asin(1)

(12)

a::Expression Float

(13)

digits(100)

(14)

a::Expression Float

(15)

So %pi already has this kind of "closure" built-in. Is it really possible to do this more generally for all possible computations with real numbers?

How are "computable reals" different than actual real numbers?

wyscc wrote:

Any floating point system is only, mathematically speaking, a small subset, and not evenly distributed one for that, of the reals, and for that matter, of the rationals. It is definitely not FRAC INT, which is mathematically equivalent to the field of rational numbers.

But surely there is an isomorphism between the domain of infinite precision floating point numbers and the domain of rationals, no?

Maybe these computable reals are something else? Isn't it related to the RealClosure as already implemented in Axiom?

Tim wrote:

This raises the same kind of implementation issue that indefinite computation raises except that indefinites are symbolic and infinite precision floats are not.

But that is a very big difference. For example, Axiom using FLOAT can simulate arbitrary precision with automatic recomputation right now. All one has to do is instead of storing the value of any identifier in a fixed floating point representation, store it as a function (or macro) to compute that value. This is exactly the idea used by Michael Stoll. In Mathematica, it is using SetDelayed instead of Set and in Axiom, it is using == instead of :=. The defining expression (function) is stored with the identifier, but no composition is actually stored. Instead, when evaluation time comes, the values of the identifiers involved are recursively computed (the actual implementation may be more clever and efficient, but this is discussion on the conceptual level). Evaluation is always possible as long as all the identifiers at the leaves are.

However, in my view, this still is not infinite precision in the sense of that will not compute exactly because the system does not automatically increase precision in FLOAT (it won't, even using macros). But one can compute it by manually increasing the precision (which precision is by no means obvious). Compare this with FRAC INT where there is no need to manually increase the length of integers in numerators or denominators. Floating point systems are designed to compute only significant digits. (More comments on implementing infinite precision below).

)clear all

   All user variables and function definitions have been cleared.
bits(68)$Float

(16)

x == 2^(-34)::Float

y == 2^(35)::Float

z == y+x

x

Compiling body of rule x to compute value of type Float

(17)

y

Compiling body of rule y to compute value of type Float

(18)

z

Compiling body of rule z to compute value of type Float

(19)

bits(200)$Float

(20)

z

(21)

bits(68)$Float

(22)

2^(-34)+2^35

(23)

In contrast, evaluation is not always possible for indefinites (a topic to be discussed separately). It is difficult to evaluate indefinites since that results in the recursive composition of functions (technique involved bears similarity to code optimization by compiler, but unfortunately, we do not just want the code, but a mathematically readable expression). The indefinites are the functions themselves, not their functional values. The output form of an indefinite is either the defining expression (typically useless, since no computation is performed), or recursively composed (typically a mess), or recursively analyzed into cases (proviso-attached expressions, and extremely computationally demanding to simplify provisos). When an indefinite is involved in a loop control construct, it is difficult to "simplify". None of these difficulties are present in infinite precision FLOAT, however it is to be implemented.

Bill Page wrote:

So %pi already has this kind of "closure" built-in. Is it really possible to do this more generally for all possible computations with real numbers?

Yes.

bits(68)$Float

(24)

x == sin(1/sqrt(2)::Float)

   Compiled code for x has been cleared.
   Compiled code for z has been cleared.
   1 old definition(s) deleted for function or rule x

x

Compiling body of rule x to compute value of type Float

(25)

bits(100)$Float

(26)

x

(27)

But surely there is an isomorphism between the domain of infinite precision floating point numbers and the domain of rationals, no?

If by such a domain, you meant a floating point representation where the mantissa has infinite precision (not arbitrary precision), then you might as well use infinite sequences of rational numbers in decimal (or binary) expansion, and you could represent, in theory, and perform exact arithmetic like . Since the set of rational numbers is dense in the set of real numbers, and FRAC INT is the field of rational numbers, rational number sequences are better for approximating real numbers, and of course, they would be exact on rationals (unlike FPS). In that case the isomorphism would be with the field of real numbers.

But this is of course not what floating point systems are. In FPS, you may be capturing more and more rational numbers as the precision is increased, but the gaps between two consecutive FLOAT numbers remain large at higher exponents (they amplify the gaps). Also, not every rational number has a finite decimal (or binary) expansion (1/3 would not be representable in an arbitrary precision FPS with base 10). So any rational with a recurring fractional (non-terminating) expansion in the base will not be representable and there cannot be a bijection (not to mention an isomorphism since floating point operations don't respect even the associative law).

Another way of stating the difference is that IPFPS requires the limiting value of mantissa (and associated exponent to locate the fractional point), which does not exist in APFPS (which I view as the infinite union of -precision FPS over all positive ; your view may be different). There are problems with simulating IPFPS by automatically adjusting the precision to yield exact results on floating point computation. The resulting system still cannot perform exact arithmetic on all rational (therefore, real) numbers due to truncation and round off. Whether the error term tends to zero would depend on the particular computation involved (problem may be ill-conditioned). We also need algorithms to automatically adjust the precision, worry about efficiency because of repeated recalculation to increase precision, and perhaps worry about termination. There is also an intrinsic problem in the representation.

Consider the problem of extending the precision of two FP numbers, which have identical FP representation at 10-bit precision. To avoid truncation errors due to binary to decimal conversion, we will work with internal representations of FLOAT only.

z == sqrt(3)::Float

   1 old definition(s) deleted for function or rule z

bits(10)$Float

(28)

mz := mantissa z

Compiling body of rule z to compute value of type Float

(29)

ez := exponent z

(30)

z10 := mz*2^ez

(31)

z10float == z10::Float

mantissa z10float

Compiling body of rule z10float to compute value of type Float

(32)

exponent z10float

(33)

mantissa(z-z10float)

(34)

So here we have two identical 10-bit floating point numbers, both defined as macros. If we extend the precision, they would no longer be the same.

bits(20)$Float

(35)

mantissa z

(36)

mantissa z10float

(37)

This raises some questions. It showed that the current representation of FLOAT is not equipped to handle "infinite precision" floating point computation, whatever that may be, and indeed, not even arbitrary precision, since extending the precision of a number depends not just on the number, but on the functions used to define it. An accompanying question is how to test equality. Similarly, if we are to find the squares of z and z10float, what precisions should be set for the proposed "infinite precision" FPS?

There should be a lot of similarities between IPFP and power series because we should treat IPFP numbers as infinite sequences, not floats. In power series, we start off with the infinite, compute with the infinite (lazy evaluation), but display finitely. In FLOAT, we start with the finite and extend, and there is no unique way to do so. In using FLOAT to simulate IPFP, we start with the infinite when we use macros, compute either finitely (:=) or infinitely (==, the equivalent of lazy evaluation) and display finitely; but the objects are not correctly stored in the data representation. Just as we don't simulate power series using polynomials, we shouldn't simulate IPFP using FPS. This is Axiom, afterall.

Maybe these computable reals are something else? Isn't it related to the RealClosure as already implemented in Axiom?

APFPS is not the same as RealClosure.

For some background material, see Lang: Algebra. Here I recall a few basic definitions. A field is a real field if is not the sum of squares in . A field is real closed if is real, and any algebraic extension of that is real must be itself. Every real closed field has a unique ordering. The real closure of a real field is an algebraic extension of that is real closed. Every real field has a real closure , and if is ordered, is unique (up to a unique isomorphism). In , every positive element is a square and the algebraic closure of is .

RealClosure implements computation in a real closure of a base field. In particular, it does not compute any element that is transcendental over the base field. RealClosure provides two modes of computation, one in terms of minimal polynomials and the other, an interval containing the particular root is used for approximation. Its relation with FLOAT, I believe, would be for solving polynomial equations numerically.

infinite sequences of rational numbers

and also that my presumption of an

isomorphism between the domain of infinite precision floating point numbers and the domain of rationals

was wrong.

In fact I think this subject, on the border between symbolic and numeric computations, is very relevant to Axiom and to computer algebra systems in general. For that reason and since this is not really an issue about Axiom itself, I propose that we continue this elsewhere on another page.

I've create a new page named RealNumbers. Please joint me there. I hope others will also contribute to this subject.