Edit detail for #167 Infinite Floats Domain revision 1 of 2

changed:
-
We need a way to handle infinite but bounded length floating
point computation. At present we have FLOAT and DFLOAT but
we really need general purpose floating point algorithms
similar in spirit to infinite length integers.

From BillPage Sun Jun 12 14:26:58 -0500 2005
From: Bill Page
Date: Sun, 12 Jun 2005 14:26:58 -0500
Subject: ??? FLOAT
Message-ID: <20050612142658-0500@page.axiom-developer.org>

Tim,

What is the domain 'FLOAT' if not already "infinite but bounded
length"? I thought that FLOAT has internal representation based
on a pair of "infinite but bounder length" INT's.

>From 'float.spad.pamphlet' I read::

<pre>
  )abbrev domain FLOAT Float

  B ==> Boolean
  I ==> Integer
  S ==> String
  PI ==> PositiveInteger
  RN ==> Fraction Integer
  SF ==> DoubleFloat
  N ==> NonNegativeInteger

  ++ Author: Michael Monagan
  ++ Date Created:
  ++   December 1987
  ++ Change History:
  ++   19 Jun 1990
  ++ Basic Operations: outputFloating, outputFixed, outputGeneral, outputSpacing,
  ++   atan, convert, exp1, log2, log10, normalize, rationalApproximation,
  ++   relerror, shift, / , **
  ++ Keywords: float, floating point, number
  ++ Description: \spadtype{Float} implements arbitrary precision floating
  ++ point arithmetic.
  ++ The number of significant digits of each operation can be set
  ++ to an arbitrary value (the default is 20 decimal digits).
  ...
  ++ The underlying representation for floats is binary
  ++ not decimal. The implications of this are described below.
  ++
  ++ The model adopted is that arithmetic operations are rounded to
  ++ to nearest unit in the last place, that is, accurate to within
  ++ \spad{2**(-\spadfunFrom{bits}{FloatingPointSystem})}.
  ++ Also, the elementary functions and constants are
  ++ accurate to one unit in the last place.
  ++ A float is represented as a record of two integers, the mantissa
  ++ and the exponent.  The \spadfunFrom{base}{FloatingPointSystem}
  ++ of the representation is binary, hence
  ++ a \spad{Record(m:mantissa,e:exponent)} represents the number \spad{m * 2 ** e
  ...
  Rep := Record( mantissa:I, exponent:I )
  ...

</pre>
 
<hr>
From Tim Daly Sun, 12 Jun 2005 17:17:54 -0500<br>
From: Tim Daly<br>
Date: Sun, 12 Jun 2005 17:17:54 -0500<br>
Subject: infinite floats domain<br>
Message-ID: <20050612171754-0500@page.axiom-developer.org><br>

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

t

<hr>
From Bill Page, Sun, 12 Jun 2005 21:24:48 -0400<br>
From: "Bill Page" <bill.page1@sympatico.ca><br>
Subject: infinite floats domain<br>
Message-ID: <20050612212448-0400@page.axiom-developer.org>,br>

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.

<hr>
From wyscc Mon, 13 Jun 2005 07:00:00 -4:00

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 <B>is</B> 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 $m \times b^e$, where $m$ is the mantissa (assume finite precision), $b$ is the base (fixed), and $e$ is the exponent (assume also finite precision). In practice, there is an offset when $m$ and $e$ are stored so that $m$ is usually interpreted with a "base-point"  and the most significant digit is normalized to be non-zero, and $e$ has an offset depending on the machine word length and how a word is divided between $m, e$ (and the signs, which for simplicity, we will assume are included in $m$ and $e$). 

Axiom implementation uses the simplified representation of $m \times b^e$ as it mathematically stands ($m$ 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 $b = 10$ 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 $-9$ 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 $-99$ to $-10$, 0, 10 to 99). For simplicity, we will only discuss positive numbers. So one sees that there can only be $19 \times 90$ 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 $b = 2$, 
$m$ 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 $e$ has 11 bits (including sign, ranging from $-1075$ to $971$ 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: 

<B>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)?</B>

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.

\begin{axiom}
bits()$Float
x:Float:=2^(-34)
y:Float:=2^35
x+y
\end{axiom}

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

\begin{axiom}
x:Float:=sqrt(2)
bits(100)$Float
z:Float:=sqrt(2)
z - x
\end{axiom}

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.

\begin{axiom}
x - 1.4142135623730950488
bits(120)$Float
z - 1.4142135623730950488016887242
\end{axiom}

<hr>

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 daly Mon Jun 13 10:34:55 -0500 2005
From: daly
Date: Mon, 13 Jun 2005 10:34:55 -0500
Subject: computable real numbers
Message-ID: <20050613103455-0500@page.axiom-developer.org>

>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




From BillPage Mon Jun 13 11:08:12 -0500 2005
From: Bill Page
Date: Mon, 13 Jun 2005 11:08:12 -0500
Subject: That's cool but ...
Message-ID: <20050613110812-0500@page.axiom-developer.org>

How is this different than what Axiom already does?
I can write:
\begin{axiom}
a:=2*asin(1)
a::Expression Float
digits(100)
a::Expression Float
\end{axiom}

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?

<hr>
From William Sit, Tuesday June 14 20:06:00 -4:00<br>
Subject: Float, Real, RealClosure

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 $2^{(-35)} + 2^{34}$ 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).

\begin{axiom}
)clear all
bits(68)$Float
x == 2^(-34)::Float
y == 2^(35)::Float
z == y+x
x
y
z
bits(200)$Float
z
bits(68)$Float
2^(-34)+2^35
\end{axiom}

  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.

\begin{axiom}
bits(68)$Float
x == sin(1/sqrt(2)::Float)
x
bits(100)$Float
x
\end{axiom}

> 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, $\sqrt{2}$ and perform exact arithmetic like $1/3 + 1/5$. 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 $k$-precision FPS over all positive $k$; 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.

\begin{axiom}
z == sqrt(3)::Float
bits(10)$Float
mz := mantissa z
ez := exponent z
z10 := mz*2^ez
z10float == z10::Float
mantissa z10float
exponent z10float
mantissa(z-z10float)
\end{axiom}

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.

\begin{axiom}
bits(20)$Float
mantissa z
mantissa z10float
\end{axiom}

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 $K$ is a *real field* if $-1$ is not the sum of squares in $K$.  A field $K$ is *real closed* if $K$ is real, and any algebraic extension of $K$ that is real must be $K$ itself.  Every real closed field has a unique ordering. The *real closure* of a real field $K$ is an algebraic extension of $K$ that is real closed. Every real field $K$ has a real closure $L$, and if $K$ is ordered, $L$ is unique (up to a unique isomorphism). In $L$, every positive element is a square and the algebraic closure of $L$ is $L[\sqrt{-1}]$.

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.

From BillPage Tue Jun 14 23:00:11 -0500 2005
From: Bill Page
Date: Tue, 14 Jun 2005 23:00:11 -0500
Subject: Continued on another page see RealNumbers
Message-ID: <20050614230011-0500@page.axiom-developer.org>

William, thank you very much for your comments on this issue.
Your presentation above is very instructive and useful to me.
I agree completely with your statements about the relevance of

> 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.

Submitted by : (unknown) at: 2007-11-17T22:01:07-08:00 (17 years ago) Name : Axiom Version : default friCAS-20090114 Axiom-20050901 OpenAxiom-20091012 OpenAxiom-20110220 OpenAxiom-Release-141 Category : general Severity : critical serious normal minor wishlist Status : open pending closed Optional subject :   Optional comment :

We need a way to handle infinite but bounded length floating point computation. At present we have FLOAT and DFLOAT but we really need general purpose floating point algorithms similar in spirit to infinite length integers.

??? FLOAT --Bill Page, Sun, 12 Jun 2005 14:26:58 -0500 reply

Tim,

What is the domain FLOAT if not already "infinite but bounded length"? I thought that FLOAT has internal representation based on a pair of "infinite but bounder length" INT's.

From float.spad.pamphlet I read::

  )abbrev domain FLOAT Float
  B ==> Boolean
  I ==> Integer
  S ==> String
  PI ==> PositiveInteger
  RN ==> Fraction Integer
  SF ==> DoubleFloat
  N ==> NonNegativeInteger
  ++ Author: Michael Monagan
  ++ Date Created:
  ++   December 1987
  ++ Change History:
  ++   19 Jun 1990
  ++ Basic Operations: outputFloating, outputFixed, outputGeneral, outputSpacing,
  ++   atan, convert, exp1, log2, log10, normalize, rationalApproximation,
  ++   relerror, shift, / , 
  ++ Keywords: float, floating point, number
  ++ Description: \spadtype{Float} implements arbitrary precision floating
  ++ point arithmetic.
  ++ The number of significant digits of each operation can be set
  ++ to an arbitrary value (the default is 20 decimal digits).
  ...
  ++ The underlying representation for floats is binary
  ++ not decimal. The implications of this are described below.
  ++
  ++ The model adopted is that arithmetic operations are rounded to
  ++ to nearest unit in the last place, that is, accurate to within
  ++ \spad{2(-\spadfunFrom{bits}{FloatingPointSystem})}.
  ++ Also, the elementary functions and constants are
  ++ accurate to one unit in the last place.
  ++ A float is represented as a record of two integers, the mantissa
  ++ and the exponent.  The \spadfunFrom{base}{FloatingPointSystem}
  ++ of the representation is binary, hence
  ++ a \spad{Record(m:mantissa,e:exponent)} represents the number \spad{m  2 * e
  ...
  Rep := Record( mantissa:I, exponent:I )
  ...

---

From Tim Daly Sun, 12 Jun 2005 17:17:54 -0500
From: Tim Daly
Date: Sun, 12 Jun 2005 17:17:54 -0500
Subject: infinite floats domain
Message-ID: <20050612171754-0500@page.axiom-developer.org>

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

t

---

From Bill Page, Sun, 12 Jun 2005 21:24:48 -0400
From: "Bill Page"
Subject: infinite floats domain
Message-ID: <20050612212448-0400@page.axiom-developer.org>,br>

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.

---

From wyscc Mon, 13 Jun 2005 07:00:00 -4:00

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.

axiom

bits()$Float

(1)

Type: PositiveInteger?

axiom

x:Float:=2^(-34)

(2)

Type: Float

axiom

y:Float:=2^35

(3)

Type: Float

axiom

x+y

(4)

Type: Float

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

axiom

x:Float:=sqrt(2)

(5)

Type: Float

axiom

bits(100)$Float

(6)

Type: PositiveInteger?

axiom

z:Float:=sqrt(2)

(7)

Type: Float

axiom

z - x

(8)

Type: Float

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.

axiom

x - 1.4142135623730950488

(9)

Type: Float

axiom

bits(120)$Float

(10)

Type: PositiveInteger?

axiom

z - 1.4142135623730950488016887242

(11)

Type: Float

---

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

computable real numbers --daly, Mon, 13 Jun 2005 10:34:55 -0500 reply

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

That's cool but ... --Bill Page, Mon, 13 Jun 2005 11:08:12 -0500 reply

How is this different than what Axiom already does? I can write:

axiom

a:=2*asin(1)

(12)

Type: Expression(Integer)

axiom

a::Expression Float

(13)

Type: Expression(Float)

axiom

digits(100)

(14)

Type: PositiveInteger?

axiom

a::Expression Float

(15)

Type: Expression(Float)

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?

---

From William Sit, Tuesday June 14 20:06:00 -4:00
Subject: Float, Real, RealClosure?

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).

axiom

)clear all

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

(16)

Type: PositiveInteger?

axiom

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

Type: Void

axiom

y == 2^(35)::Float

Type: Void

axiom

z == y+x

Type: Void

axiom

x

axiom

Compiling body of rule x to compute value of type Float

(17)

Type: Float

axiom

y

axiom

Compiling body of rule y to compute value of type Float

(18)

Type: Float

axiom

z

axiom

Compiling body of rule z to compute value of type Float

(19)

Type: Float

axiom

bits(200)$Float

(20)

Type: PositiveInteger?

axiom

z

(21)

Type: Float

axiom

bits(68)$Float

(22)

Type: PositiveInteger?

axiom

2^(-34)+2^35

(23)

Type: Fraction(Integer)

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.

axiom

bits(68)$Float

(24)

Type: PositiveInteger?

axiom

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

Type: Void

axiom

x

axiom

Compiling body of rule x to compute value of type Float

(25)

Type: Float

axiom

bits(100)$Float

(26)

Type: PositiveInteger?

axiom

x

(27)

Type: Float

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.

axiom

z == sqrt(3)::Float

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

Type: Void

axiom

bits(10)$Float

(28)

Type: PositiveInteger?

axiom

mz := mantissa z

axiom

Compiling body of rule z to compute value of type Float

(29)

Type: PositiveInteger?

axiom

ez := exponent z

(30)

Type: Integer

axiom

z10 := mz*2^ez

(31)

Type: Fraction(Integer)

axiom

z10float == z10::Float

Type: Void

axiom

mantissa z10float

axiom

Compiling body of rule z10float to compute value of type Float

(32)

Type: PositiveInteger?

axiom

exponent z10float

(33)

Type: Integer

axiom

mantissa(z-z10float)

(34)

Type: NonNegativeInteger?

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.

axiom

bits(20)$Float

(35)

Type: PositiveInteger?

axiom

mantissa z

(36)

Type: PositiveInteger?

axiom

mantissa z10float

(37)

Type: PositiveInteger?

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.

Continued on another page see RealNumbers? --Bill Page, Tue, 14 Jun 2005 23:00:11 -0500 reply

William, thank you very much for your comments on this issue. Your presentation above is very instructive and useful to me. I agree completely with your statements about the relevance of

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.