RFC: Arbitrary-precision floats for Guile

RFC: Arbitrary-precision floats for Guile

[Mark H Weaver]

Hello all,

I want to add arbitrary-precision floats to Guile.
I'd use the mpf_t type from GNU MP.

There would be a fluid whose value would determine the minimum precision
to use for inexact operators.  A value of #f (the default) would mean
that normal floats would be used unless one of the operands was a
bigfloat.  Procedures that produce inexact numbers (including read and
exact->inexact) would use the precision of the most precise inexact
operand contributing to the result, or the value of the fluid, whichever
is greater.  There would also be a procedure to reduce the precision of
an inexact number to some specified value.  The syntax for inexact
numbers would allow an optional suffix to specify the precision.

One complication is that we'd have to implement the transcendental
functions.  However, I already have code to do that, as part of a
bigfloat library I wrote a while back.

I'd also like to add another more general representation of complex
numbers whose real and imaginary parts are each of type SCM, like the
way fractions are represented.  The primary motivation would be to
support arbitrary-precision complex numbers, but as an added bonus,
would be able to support exact Gaussian integers.  The existing complex
representation would still be supported in the interests of efficiency.

scm_make_rectangular would automatically choose the right
representation: if the imaginary part was an exact zero, it would simply
return the specified real part.  If both parts were normal floats it
would return the existing fast complex numbers, otherwise it would
return the new more general complex number type.

Now I know we have a release coming soon.  I don't necessarily expect to
get these changes merged before then, although I'm willing and able to
implement all of this in the next week.  However, I was hoping that I
could at least get enough of it into 2.0 so that it could be part of 2.1
without breaking ABI compatibility.  Do you think this is feasible?
If not, would it be possible to delay the release a week or two?

What do you think of this?  Am I crazy? :)

    Best,
     Mark

Re: RFC: Arbitrary-precision floats for Guile

[Hans Aberg]

On 1 Feb 2011, at 14:57, Mark H Weaver wrote:

> I want to add arbitrary-precision floats to Guile.
> I'd use the mpf_t type from GNU MP.
>
> There would be a fluid whose value would determine the minimum  
> precision
> to use for inexact operators.  A value of #f (the default) would mean
> that normal floats would be used unless one of the operands was a
> bigfloat.

Only checking takes a lot of time, so it is best to make a new type.  
Same for integral types. So if speed is an issue, Guile ought to have  
C99 integral types.

The name 'floating' seems free for such a type.

> One complication is that we'd have to implement the transcendental
> functions.  However, I already have code to do that, as part of a
> bigfloat library I wrote a while back.

There is already MPFR, which is a package of top of GMP (see the  
manual of the latter).

> I'd also like to add another more general representation of complex
> numbers whose real and imaginary parts are each of type SCM, like the
> way fractions are represented.

One might have a complexification class with one template argument.  
The current Guile 'complex' type is (in my installation) the  
complexification of 64-bit IEEE floats.

Re: RFC: Arbitrary-precision floats for Guile

[dsmich]

---- Hans Aberg <haberg-1@telia.com> wrote: 
> On 1 Feb 2011, at 14:57, Mark H Weaver wrote:
> 
> > I want to add arbitrary-precision floats to Guile.
> > I'd use the mpf_t type from GNU MP.
> 
> > One complication is that we'd have to implement the transcendental
> > functions.  However, I already have code to do that, as part of a
> > bigfloat library I wrote a while back.
> 
> There is already MPFR, which is a package of top of GMP (see the  
> manual of the latter).

Sounds cool.  Could this be made optional though?  There is some recent interest in cross compiling Guile, and I'd like to try and put it on some embedded device some day.  Something like this probably isn't all that necessary for an embedded system, and even one less dependency helps.

Thanks,
  -Dale

Re: RFC: Arbitrary-precision floats for Guile

[Mark H Weaver]

Hans Aberg <haberg-1@telia.com> writes:
>> There would be a fluid whose value would determine the minimum precision
>> to use for inexact operators.  A value of #f (the default) would mean
>> that normal floats would be used unless one of the operands was a
>> bigfloat.
>
> Only checking takes a lot of time, so it is best to make a new
> type. Same for integral types. So if speed is an issue, Guile ought to
> have  C99 integral types.
>
> The name 'floating' seems free for such a type.

I'm not sure I understand.  These type checks would only be added as
additional cases to the generic arithmetic operators, which
unfortunately cannot avoid checking the types of their arguments for
each operation.  Sadly this is already the case.  Currently Guile
supports fixnums (small integers), arbitrary-precision integers,
arbitrary-precision rationals, IEEE 64-bit floats, and complex numbers
composed of IEEE 64-bit floats.  All of these cases must be considered
on each operation.  To make matters worse, all of these types except
small integers must be allocated on on the heap.  Therefore every
floating point operation involves heap allocation, and additionally pays
the amortized cost of garbage collecting the number later on.

The only additional overhead added to existing floating point operations
would be to check the value of the fluid.  Admittedly this is not
entirely trivial, but I think it will be lost in the noise compared with
the overheads already present.

>> One complication is that we'd have to implement the transcendental
>> functions.  However, I already have code to do that, as part of a
>> bigfloat library I wrote a while back.
>
> There is already MPFR, which is a package of top of GMP (see the
> manual of the latter).

Excellent, thanks for the pointer!  MPFR seems to contain all of the
functionality I need, so this will make my job much easier.

>> I'd also like to add another more general representation of complex
>> numbers whose real and imaginary parts are each of type SCM, like the
>> way fractions are represented.
>
> One might have a complexification class with one template
> argument. The current Guile 'complex' type is (in my installation) the
> complexification of 64-bit IEEE floats.

Please forgive me if I misunderstand, but it seems to me that you are
thinking in terms of static type checking, which is wholly different
from how Guile operates internally.  Although I don't see a need for
complex numbers whose real and imaginary parts are of different types,
it would require additional implementation complexity to enforce this
constraint.  I'm not sure it's worth the effort.  My plan is to simply
call the generic arithmetic operators recursively to operate on the real
and imaginary parts.

Have I misunderstood you?  I would very much like to understand.

      Thanks,
        Mark

Re: RFC: Arbitrary-precision floats for Guile

[Hans Aberg]

On 1 Feb 2011, at 18:03, Mark H Weaver wrote:

>>> There would be a fluid whose value would determine the minimum  
>>> precision
>>> to use for inexact operators.  A value of #f (the default) would  
>>> mean
>>> that normal floats would be used unless one of the operands was a
>>> bigfloat.
>>
>> Only checking takes a lot of time, so it is best to make a new
>> type. Same for integral types. So if speed is an issue, Guile ought  
>> to
>> have  C99 integral types.
>>
>> The name 'floating' seems free for such a type.
>
> I'm not sure I understand.  These type checks would only be added as
> additional cases to the generic arithmetic operators, which
> unfortunately cannot avoid checking the types of their arguments for
> each operation.  Sadly this is already the case.  Currently Guile
> supports fixnums (small integers), arbitrary-precision integers,
> arbitrary-precision rationals, IEEE 64-bit floats, and complex numbers
> composed of IEEE 64-bit floats.  All of these cases must be considered
> on each operation.  To make matters worse, all of these types except
> small integers must be allocated on on the heap.  Therefore every
> floating point operation involves heap allocation, and additionally  
> pays
> the amortized cost of garbage collecting the number later on.

Right. I'm told that overflow checks on integral types typically add  
the order of ten cycles per arithmetic operation. So that would  
already be a bottleneck in Guile.

> The only additional overhead added to existing floating point  
> operations
> would be to check the value of the fluid.  Admittedly this is not
> entirely trivial, but I think it will be lost in the noise compared  
> with
> the overheads already present.

All those overheads are really going to bog it down, I suspect, if you  
need something fast.

>>> One complication is that we'd have to implement the transcendental
>>> functions.  However, I already have code to do that, as part of a
>>> bigfloat library I wrote a while back.
>>
>> There is already MPFR, which is a package of top of GMP (see the
>> manual of the latter).
>
> Excellent, thanks for the pointer!  MPFR seems to contain all of the
> functionality I need, so this will make my job much easier.

Why reinventing the wheel? There is a group of maintainers of this  
package, hard focusing on getting the numerics right.

>>> I'd also like to add another more general representation of complex
>>> numbers whose real and imaginary parts are each of type SCM, like  
>>> the
>>> way fractions are represented.
>>
>> One might have a complexification class with one template
>> argument. The current Guile 'complex' type is (in my installation)  
>> the
>> complexification of 64-bit IEEE floats.
>
> Please forgive me if I misunderstand, but it seems to me that you are
> thinking in terms of static type checking, which is wholly different
> from how Guile operates internally.

I just it as a language to describe it. It would be function of types  
with one parameter. One can implement template using lambda calculus  
as well, check Template Haskell.

> Although I don't see a need for
> complex numbers whose real and imaginary parts are of different types,
> it would require additional implementation complexity to enforce this
> constraint.  I'm not sure it's worth the effort.  My plan is to simply
> call the generic arithmetic operators recursively to operate on the  
> real
> and imaginary parts.
>
> Have I misunderstood you?  I would very much like to understand.

The Complexification C(R) of a ring R is a mathematical construction,  
the set R x R of pairs, written as x + iy where i is a formal symbol,  
with component-wise addition, and multiplication given by using the  
rule i^ = -1.

So instead of making a new complex type for every numeric type, the  
idea is to do a one parameter complexification class. I'm not sure how  
that works out Scheme, but Haskell has such a module.

Re: RFC: Arbitrary-precision floats for Guile

[Andy Wingo]

Hi Mark,

On Tue 01 Feb 2011 14:57, Mark H Weaver <mhw@netris.org> writes:

> I want to add arbitrary-precision floats to Guile.
> I'd use the mpf_t type from GNU MP.
>
> There would be a fluid whose value would determine the minimum
> precision to use for inexact operators.  A value of #f (the default)
> would mean that normal floats would be used unless one of the operands
> was a bigfloat.  Procedures that produce inexact numbers (including
> read and exact->inexact) would use the precision of the most precise
> inexact operand contributing to the result, or the value of the fluid,
> whichever is greater.

What other procedures would use this fluid?  Can you enumerate them?

In (* inum flonum bigflonum), with what precision would the first
multiplication be performed?  Note that currently the compiler compiles
it as (* (* inum flonum) bigflownum).

I would like to avoid adding more core numeric types during the 2.0.  I
don't think we have time for a full implementation of this before 2.0,
either.  Is there a minimal change that you can make that adds this
type?

The same consideration applies to exact complex numbers.  I wonder
though: could we keep things simple and just have complex components as
SCM values?  In the common flonum case, we can allocate the SCM values
inline, at the end of the complex SCM object.  That's probably enough of
a performance win that we can claim victory there, and avoid two kinds
of complex numbers.

Also, I would like for the mpfr dependency to be optional, if possible.

> What do you think of this?  Am I crazy? :)

Possibly :)  If it makes it in, great.  ~Small changes can go into 2.0,
at this point.  But if it doesn't, that will just be a good reason to
get a 2.2 out soonish.

Cheers,

Andy
-- 
http://wingolog.org/

Re: RFC: Arbitrary-precision floats for Guile

[Nelson H. F. Beebe]

Before rushing off and implementing arbitrary-precision floats in
Guile, please study this recent book carefully first, particularly
Chapter 3 that discusses IEEE Standards requirements for
floating-point arithmetic:

@String{pub-BIRKHAUSER-BOSTON   = "Birkh{\"a}user Boston Inc."}
@String{pub-BIRKHAUSER-BOSTON:adr = "Cambridge, MA, USA"}

@Book{Muller:2010:HFP,
  author =       "Jean-Michel Muller and Nicolas Brisebarre and Florent
                 de Dinechin and Claude-Pierre Jeannerod and Vincent
                 Lef{\`e}vre and Guillaume Melquiond and Nathalie Revol
                 and Damien Stehl{\'e} and Serge Torres",
  title =        "Handbook of Floating-Point Arithmetic",
  publisher =    pub-BIRKHAUSER-BOSTON,
  address =      pub-BIRKHAUSER-BOSTON:adr,
  pages =        "xxiii + 572",
  year =         "2010",
  DOI =          "http://dx.doi.org/10.1007/978-0-8176-4704-9",
  ISBN =         "0-8176-4704-X",
  ISBN-13 =      "978-0-8176-4704-9",
  LCCN =         "QA76.9.C62 H36 2010",
  bibdate =      "Thu Jan 27 16:18:58 2011",
  price =        "US\$90 (est.)",
  acknowledgement = ack-nhfb,
}

The authors of that book are a group of leading French researchers who
have done a lot of important work in the area of hardware and software
for floating-point arithmetic, and some are co-authors of the MPFR
package that builds on top of the facilities provided by GMP.  MPFR
and its companion, MPC (from

	http://mpc.multiprecision.org/

) for complex arithmetic, represents the state of the art in the area
of correctly-rounded arithmetic.

Another recent book addresses arithmetic and the computation of
elementary functions in software arbitrary-precision arithmetic:

@String{pub-CAMBRIDGE           = "Cambridge University Press"}
@String{pub-CAMBRIDGE:adr       = "Cambridge, UK"}

@Book{Brent:2011:MCA,
  author =       "Richard P. Brent and Paul Zimmermann",
  title =        "Modern computer arithmetic",
  volume =       "18",
  publisher =    pub-CAMBRIDGE,
  address =      pub-CAMBRIDGE:adr,
  pages =        "xvi + 221",
  year =         "2011",
  ISBN =         "0-521-19469-5 (hardcover)",
  ISBN-13 =      "978-0-521-19469-3 (hardcover)",
  LCCN =         "QA76.9.C62 BRE 2011",
  bibdate =      "Sat Jan 15 12:25:22 MST 2011",
  bibsource =    "library.ox.ac.uk:210/ADVANCE",
  series =       "Cambridge monographs on applied and computational
                 mathematics",
  acknowledgement = ack-nhfb,
  subject =      "Computer arithmetic",
}

Recent releases of the gcc-4.x compiler family require GMP, MPFR, and
MPC so that correctly-rounded compile-time conversions can be
guaranteed; thus, recent GNU/Linux systems are likely to have all
three libraries already installed.  All are extremely portable to
systems with IEEE 754 arithmetic; don't expect them to run on older
architectures (VAX, PDP-10, etc.).

The programming-language compiler field is littered with far too many
code optimizations that are completely wrong for floating-point
arithmetic; please do not add to that mess!  For example,
floating-point arithmetic is not associative, so (+ (a (+ b c)) is not
the same as (+ (b (+ a c)), and (* (a (* b c))) may differ (possibly
wildly) from (* (b (* a c))) etc.

Reciprocation is rarely exact, so (/ x y) is NOT the same as
(* x (/ 1 y)), unless y is known in advance to be a power of the base.

Importantly, (equal x x) MUST return #f if x is a NaN (quiet or
signaling), and cannot be optimized away to #t.

Comparisons of floating-point values may return UNORDERED, so (< a b)
is NOT the same as (not (>= b a)).

It is imperative that the floating-point programmer have control over
whether higher precision is used in intermediate evaluations.  The
C/C++/Java/C# family, and Fortran 20xx, provide a volatile type
qualifier that forces intermediate results into memory in their normal
storage size.  This can be an issue in the Lisp family, where type
declarations may be omitted, or variable types may change at run time.

If interval arithmetic is to eventually be supported in guile, then
the floating-point system MUST provide access to rounding-mode
control, or else provide additional arithmetic operators that produce
upper and lower bounds for basic arithmetic operations (+, -, /, *, sqrt).
Interval arithmetic is already available in C++ and Fortran 90 in the
freely-available Sun (now Oracle) compilers, and they should be used
as a test bed for any development of interval arithmetic in guile.

Careful study of the C99 Standard and subsequent documents is also
desirable:

@String{pub-ISO                 = "International Organization for
                                  Standardization"}
@String{pub-ISO:adr             = "Geneva, Switzerland"}

@Book{ISO:1999:IIP,
  author =       "{ISO}",
  title =        "{ISO\slash IEC 9899:1999}: Programming Languages ---
                 {C}",
  publisher =    pub-ISO,
  address =      pub-ISO:adr,
  pages =        "538",
  day =          "16",
  month =        dec,
  year =         "1999",
  ISBN =         "????",
  ISBN-13 =      "????",
  LCCN =         "????",
  bibdate =      "Tue Dec 12 06:46:19 2000",
  note =         "Available in electronic form for online purchase at
                 \path=http://webstore.ansi.org/= and
                 \path=http://www.cssinfo.com/=.",
  price =        "US\$18 (electronic), US\$225 (print)",
  URL =          "http://www.iso.ch/cate/d29237.html;
                 http://anubis.dkuug.dk/JTC1/SC22/open/n2620/n2620.pdf;
                 http://anubis.dkuug.dk/JTC1/SC22/WG14/www/docs/n897.pdf;
                 http://webstore.ansi.org/ansidocstore/product.asp?sku=ISO%2FIEC+9899%3A1999",
  acknowledgement = ack-nhfb,
}

There are at least three critical subsequent Technical Corrigenda to
that document.

There are also recommendations for support of decimal arithmetic, and
additional library functions, although in my opinion, more work is
required in both cases (my forthcoming book discusses these issues in
detail):

@Misc{ISO:2006:IIJa,
  key =          "C",
  title =        "{ISO\slash IEC JTC1 SC22 WG14 N1154}: Extension for
                 the programming language {C} to support decimal
                 floating-point arithmetic",
  howpublished = "World-Wide Web document",
  day =          "27",
  month =        feb,
  year =         "2006",
  bibdate =      "Tue Mar 07 18:23:45 2006",
  URL =          "http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1154.pdf",
  acknowledgement = ack-nhfb,
  keywords =     "decimal floating-point arithmetic",
}

@Misc{ISO:2006:IIJb,
  key =          "C",
  title =        "{ISO\slash IEC JTC1 SC22 WG14 N1161}: Rationale for
                 {TR 24732}: Extension to the programming language {C}:
                 Decimal Floating-Point Arithmetic",
  howpublished = "World-Wide Web document",
  day =          "27",
  month =        feb,
  year =         "2006",
  bibdate =      "Tue Mar 07 18:23:45 2006",
  URL =          "http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1161.pdf",
  acknowledgement = ack-nhfb,
  keywords =     "decimal floating-point arithmetic",
}

@Misc{ISO:2006:IIJc,
  key =          "C",
  title =        "{ISO\slash IEC JTC1 SC22 WG14 N1176}: Extension for
                 the programming language {C} to support decimal
                 floating-point arithmetic",
  howpublished = "World-Wide Web document",
  pages =        "iii + 33",
  day =          "24",
  month =        may,
  year =         "2006",
  bibdate =      "Sat Feb 24 20:00:36 2007",
  URL =          "http://open-std.org/jtc1/sc22/wg14/www/docs/n1176.pdf",
  acknowledgement = ack-nhfb,
  keywords =     "decimal floating-point arithmetic",
}

@Book{ISO:2009:IIT,
  key =          "C",
  title =        "{ISO/IEC TR 24732:2009} Information technology ---
                 Programming languages, their environments and system
                 software interfaces --- Extension for the programming
                 language {C} to support decimal floating-point
                 arithmetic",
  publisher =    pub-ISO,
  address =      pub-ISO:adr,
  year =         "2009",
  bibdate =      "Thu Nov 25 08:56:44 2010",
  bibsource =    "http://www.iso.org/iso/search.htm",
  series =       "Technical report",
  URL =          "http://www.iso.org/iso/iso_catalogue/catalogue_tc/catalogue_detail.htm?csnumber=38842",
  acknowledgement = ack-nhfb,
  subject =      "programming languages (electronic computers)",
}

@Book{ISO:2010:IIIa,
  key =          "C++",
  title =        "{ISO/IEC 29124:2010}: Information technology ---
                 Programming languages, their environments and system
                 software interfaces --- Extensions to the {C++ Library}
                 to support mathematical special functions",
  publisher =    pub-ISO,
  address =      pub-ISO:adr,
  year =         "2010",
  bibdate =      "Thu Nov 25 08:56:44 2010",
  bibsource =    "http://www.iso.org/iso/search.htm",
  series =       "Technical report",
  URL =          "http://www.iso.org/iso/iso_catalogue/catalogue_tc/catalogue_detail.htm?csnumber=50511",
  acknowledgement = ack-nhfb,
  subject =      "programming languages (electronic computers)",
}

Recently a draft of the developing ISO C 201x standard that will
eventually replace C99 became available:

	http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1539.pdf

-------------------------------------------------------------------------------
- Nelson H. F. Beebe                    Tel: +1 801 581 5254                  -
- University of Utah                    FAX: +1 801 581 4148                  -
- Department of Mathematics, 110 LCB    Internet e-mail: beebe@math.utah.edu  -
- 155 S 1400 E RM 233                       beebe@acm.org  beebe@computer.org -
- Salt Lake City, UT 84112-0090, USA    URL: http://www.math.utah.edu/~beebe/ -
-------------------------------------------------------------------------------

Re: RFC: Arbitrary-precision floats for Guile

[Hans Aberg]

On 1 Feb 2011, at 21:37, Andy Wingo wrote:

> In (* inum flonum bigflonum), with what precision would the first
> multiplication be performed?  Note that currently the compiler  
> compiles
> it as (* (* inum flonum) bigflownum).

An idea that comes to my mind is to set a minimum float precision,  
which might be IEEE 64-bit or corresponding. Then, for multiprecision  
floats, I think a user-friendly solution would be the output precision  
of an operation being a function of the operation and the precision of  
the values put to it based on an error analysis with a margin, as to  
not introduce too big round-off errors. Experts might want to have  
something else, but that would be less user-friendly.

Re: RFC: Arbitrary-precision floats for Guile

[Pierpaolo Bernardi]

On Wed, Feb 2, 2011 at 02:53, Nelson H. F. Beebe <beebe@math.utah.edu> wrote:

> Another recent book addresses arithmetic and the computation of
> elementary functions in software arbitrary-precision arithmetic:
>
> @String{pub-CAMBRIDGE           = "Cambridge University Press"}
> @String{pub-CAMBRIDGE:adr       = "Cambridge, UK"}
>
> @Book{Brent:2011:MCA,
>  author =       "Richard P. Brent and Paul Zimmermann",
>  title =        "Modern computer arithmetic",
>  volume =       "18",
>  publisher =    pub-CAMBRIDGE,
>  address =      pub-CAMBRIDGE:adr,
>  pages =        "xvi + 221",
>  year =         "2011",
>  ISBN =         "0-521-19469-5 (hardcover)",
>  ISBN-13 =      "978-0-521-19469-3 (hardcover)",
>  LCCN =         "QA76.9.C62 BRE 2011",
>  bibdate =      "Sat Jan 15 12:25:22 MST 2011",
>  bibsource =    "library.ox.ac.uk:210/ADVANCE",
>  series =       "Cambridge monographs on applied and computational
>                 mathematics",
>  acknowledgement = ack-nhfb,
>  subject =      "Computer arithmetic",
> }

In case this is not universally known: this book is also available for
free from Zimmermann home page.

http://www.loria.fr/~zimmerma/mca/pub226.html

Cheers
P.