[racket] 80-bit precision in Racket

> We are doing numerical integration of celestial bodies over large periods of time (100 years is a norm).

I'm new to Scheme, so I may be totally wrong about this --- but, isn't a numerical program like this exactly what Scheme is *not* designed for? Isn't Scheme primarily for dynamic-OOP programs? Isn't this mostly a convenience to help people write programs quickly without having to worry about what type each datum is? I mean, isn't this mostly for writing scripts?

I would really write a program like yours in Forth. You don't have dynamic-OOP with all of the data types mushed together so that you don't know what type any particular datum is. You have explicit data types. It is an easy matter to implement another floating-point data stack that is 80-bit and all of the words are distinguished by a prefix --- like DF+ etc. (with the DF prefix meaning double-precision float, as compared to F+ which is just the usual floating point addition that is part of the standard). Also, you could modify your Forth system so that the standard floating-point stack is 80-bit rather than 64-bit, if you can be sure that everything you will be doing will be 80-bit --- that would make for cleaner source-code, and also allow your code to be ANS-Forth compliant.

If I couldn't use Forth, I would use assembly-language. Most modern assemblers have a pretty good macro system (although not as good as Forth's). I wouldn't want to use C, as it doesn't have macros at all (except #define string expansion, which doesn't count). 

If I was going to use Scheme, I would use Gambit or one of the other Scheme systems that compile into C, which does support 80-bit floats. I definitely wouldn't mess with any system that uses a VM and JIT, as that is horribly complicated (and it also tends to be slow). 

For the most part though, I would avoid using any dynamic-OOP language for a numerical program. That seems very clumsy, to have your data getting cast from one type to another without you explicitly casting it --- for numerical programming, you want to know what type (precision) your data is at all times. I actually like Racket as a learning language, and I'm working on a game in Racket --- I'm just saying that a numerical application seems inappropriate. I'm new to Scheme though, so I could be wrong --- maybe what I'm saying is heresy, and Scheme is actually the best language for *every* program that could ever be written. 

I'm certainly not saying that Forth is the best language for *every* program that could ever be written --- over on comp.lang.forth I have argued many times that the Forth community should focus on using Forth for what Forth is good at, and not try to use Forth for what Forth is not good at. I have said that if we try to do this, Forth will become a "jack of all trades and a master of none." Your astronomy program, however, seems to be one of the few cases in which Forth is a good fit.



Message: 1
Date: Sun, 11 Nov 2012 15:11:50 +0400
From: Dmitry Pavlov <dpavlov at ipa.nw.ru>
To: Matthias Felleisen <matthias at ccs.neu.edu>
Cc: users at racket-lang.org
Subject: Re: [racket] 80-bit precision in Racket
Message-ID:
    <CAA-Jc7GMds-Um3p0aaeg-X-Y6FzQB=t1oS3gK-pNVmx5ZgSPpg at mail.gmail.com>
Content-Type: text/plain; charset=UTF-8

Neil, Matthias,

Certainly not #2, and I doubt that #4 is a problem for us
(although it may be, buried somewhere in specific calculations).

#1 and #3 definitely are our case. We are doing numerical
integration of celestial bodies over large periods of time
(100 years is a norm). The forces that act on bodies may be of
quite different orders of magnitude: for example, for near-Earth
objects the acceleration towards the Sun is about 2.9e-4 AU/days^2,
while corrections for peculiar relativistic effects (Lense-Thirring
acceleration etc.) can be as small as 1e-17 AU/days^2 (I do not
have exact numbers now). And those effects must be taken
into account. Yet the numerical integration procedure has its
numerical inaccuracies, too. So when we integrate these
accelerations over ten or more years, it sort of pushes the
limit of double precision representation. Again, I do not
(yet) have a precise mathematical proof of it, but it is
more or less obvious that we need to try 80-bit at least.

In addition to the above: it is known that the DE ephemeris
done at JPL, NASA were integrated with *quadruple* (128-bit)
precision:
http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/a_old_versions/de414.cmt
so it may be a clue that even 80-bit is not enough!
We are going to try 128-bit after we get 80-bit done.

#5 is also our case. Legacy system that we are porting to Racket
uses extended precision numbers.

Thanks for your help!

Best regards,

Dmitry
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.racket-lang.org/users/archive/attachments/20121113/4cf32d59/attachment-0001.html>