Paul Koning wrote:
>>>>"Jerome" == Jerome H Fine <jhfinexgs2 at compsys.to>
writes:
>>>>
>>>>
Jerome> What I am complaining about is not the internal operations,
Jerome> but the conversion routines. I went to the trouble of
Jerome> printing out the values in OCTAL so that I could figure out
Jerome> the problem. A 56 bit mantissa should allow a value of
Jerome> 36,028,797,018,963,968 when the leading (invisible) bit is
Jerome> one with the remaining bits zero. The displayed value is
Jerome> 36,028,797,018,963,965 instead or less by 3. ...
That would be a bug. Not surprising. Floating point arithmetic
should be, and can be, accurate to the last bit. But very often it is
not, because it's *hard* to do that. DEC had a team of specialists
who worked on the numeric algorithms for VAX. The PDP-11 libraries
didn't get the benefit of that work, so it is not surprising to find
errors in the last few bits.
As I said, if you want exact answers, DO NOT use float.
Jerome Fine replies:
Again, I was not complaining about the internal floating point
arithmetic, but the conversion. As you noted, I also consider
it to be a bug which is actually rather easy to correct. Although
it takes 8 * 16 bit words to convert just the integer portion of
the floating point word, I quickly wrote a routine in assembler
to do just that (well not the whole thing, just when it was able
to fit into a 64 bit unsigned integer variable or only 4 * 16 bit
words) in about 50 assembler instructions. I then moved the
bits back into a separate 64 bit floating point number and used
octal to check that the code was correct. I am sure that just a
few more instructions would have allowed me to do the complete
conversion into a temporary 256 bit value with 128 bits for the
integer and 128 bits for the fraction. That would seem to be the
maximum required for temporary storage while the value was
converted to decimal notation for display with a format of "F80.40".
> 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...
>
Jerome> Any idea where I could download the source code? Given that
Jerome> FORTRAN is not really suitable for shifting and using carry
Jerome> bits, I would port to MACRO-11.
It's right in the GNU software library listing:
http://directory.fsf.org/GNU/gnump.html
Jerome> In general, sieve methods for prime numbers don't ever need
Jerome> division and rarely need multiplication. Addition and
Jerome> comparing are the primary operations. Increment is also
Jerome> needed (with increment by 2 and 4 being a big time saver).
Jerome> So the only other thing needed are accurate conversion
Jerome> routines to allow output in decimal instead of octal. I
Jerome> presume that the latter is best done using subtraction of
Jerome> multiple powers of TEN which are stored in a table unless
Jerome> there are better methods of converting the output? Any
Jerome> suggestions?
I'm sure I/O is part of GnuMP. Decimal? Probably, I haven't looked.
Multiplication is very definitely in there, it has to be, that's a
critical function for crypto. Division I don't know about either.
It certainly should be there if GnuMP was aiming to be a general
arithmetic package. If all else fails, it certainly has modular
arithmetic (again, because of crypto) from which you can synthesize
division if it isn't already there.
I downloaded the source file - THANK YOU for the link. I hope to find
what I may need there.
Oh yes, I'm quite sure it also includes primality
tests --
presumably the Miller-Rabin test.
I am not familiar with the Miller-Rabin test? If it
has anything to do with prime numbers, I would be very
interested.
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.