28 Jul
2005
28 Jul
'05
1:41 a.m.
Paul Koning wrote:
Jerome> I have been having some fun with Prime
Numbers using FORTRAN
Jerome> 77 under RT-11 on a PDP-11. Since I rapidly run out of the
Jerome> capacity of 32 bit numbers (INTEGER * 4 NUMBER), I have
Jerome> started to use 64 bit values (REAL * 8 NUMBER) but the
Jerome> conversion routines leave a lot to be desired on output.
Jerome> When using F32.0, I can count on only about 15 digits of
Jerome> accuracy even though the values are often accurate to an
Jerome> exact integer.
Of course. Just because it's exactly representable as an integer
doesn't mean it's exactly representable as a 64 bit float. Those have
a 56 bit mantissa, which means you have 2^56 bits of resolution, which
translates to roughly 16 digits. (The VAX architecture handbook
mentions that explicitly, by the way.)
Jerome Fine replies:
What I am complaining about is not the internal operations, but the
conversion routines. I went to the trouble of printing out the values
in OCTAL so that I could figure out the problem. A 56 bit mantissa
should allow a value of 36,028,797,018,963,968 when the leading
(invisible) bit is one with the remaining bits zero. The displayed value
is 36,028,797,018,963,965 instead or less by 3. I checked the
actual 64 bit float at 56000 octal which is an exponent of 270 octal
or 2^55 in excess floating on a PDP-11 (unless my understanding of
exactly 1 which is the 64 bit float value of 40200 which is an
exponent of 201 octal that must be 2^0 = 1 is incorrect).
Jerome> If needed, I can write my own conversion
routines, but I
Jerome> though I might inquire if anyone knows of any libraries which
Jerome> are exact or better still, multi-precision libraries which
Jerome> can handle up to 128 bit integers?
If you want exact arithmetic, run, don't walk, away from floating
point numbers. Use only integers.
There the GNU MP library, and others like it. Long integer arithmetic
libraries are commonly found in crypto library packages. Come to
think of it, you'll find functions that are helpful in stuff relating
to prime numbers -- since prime numbers also figure prominently in
cryptography.
So you could use MP (which is in C) directly, given a suitable
compiler. Or you could port it to Fortran...
Any idea where I could download the source code? Given
that FORTRAN is not really suitable for shifting and
using carry bits, I would port to MACRO-11.
In general, sieve methods for prime numbers don't ever
need division and rarely need multiplication. Addition
and comparing are the primary operations. Increment is
also needed (with increment by 2 and 4 being a big time
saver). So the only other thing needed are accurate
conversion routines to allow output in decimal instead
of octal. I presume that the latter is best done using
subtraction of multiple powers of TEN which are stored
in a table unless there are better methods of converting
the output? Any suggestions?
Sincerely yours,
Jerome Fine
--
If you attempted to send a reply and the original e-mail
address has been discontinued due a high volume of junk
e-mail, then the semi-permanent e-mail address can be
obtained by replacing the four characters preceding the
'at' with the four digits of the current year.