Float comparison

[Keith Thompson]

I see nothing in 5.2.4.2.2 that supports your claim that a stored
floating-point value represents a range of mathematical real
values rather than one specific mathematical real value. In
particular, see 5.2.4.2.2p2:

A floating-point number (x) is defined by the following model:

x = s b^e summation(k=1, p) f(sub k) b^-k, emin <= e <= emax

That's a clumsy rendition of the mathematical notation used in
the standard; you'll need to use a PDF (or hardcopy) version to
see it.

The model defines a single exact value, not a range of values,
for any floating-point number.

Repeat. In 5.2.4.2.2:

[#10] The values given in the following list shall be
replaced by implementation-defined constant expressions with
(positive) values that are less than or equal to those
shown:

-- the difference between 1 and the least value greater
than 1 that is representable in the given floating
point type, b1-p

FLT_EPSILON 1E-5
DBL_EPSILON 1E-9
LDBL_EPSILON 1E-9

Combine that with the method of expressing non-1.0 values.

Yes, and???

Did you even look at 5.2.4.2.2p2, which defines C's model for
floating-point values? (That's a serious question.)

Yes, I know what FLT_EPSILON, DBL_EPSILON, and LDBL_EPSILON are. In
type double, 1.0 is exactly representable, 1.0+DBL_EPSILON is exactly
representable, and the infinitely many real values between them are
not.

This does not support your claim that a floating-point value
represents a range of real values. A floating-point value, according
to C's model for floating-point numbers, represents a single exact
real value.

 

[Keith Thompson]

BartC wrote: [...]

I think he means: if the range of any floating point value is say 0 to 10%
(eg. 10.0 represents real values 10 to 11), then
multiplying 10 x 10 gives 100, which with a range of 10% gives 100 to 110.

However (10...11) x (10...11) has results in the range 100 to 121. 121 is
not in the range 100 to 110.

This is a chimera. The errors are small, compared with the
values. The sort of effect you are worrying about is of the form
error_ratio squared, and thus negligible.

I'm truly at a loss to understand how that's relevant. What's being
discussed is what a floating-point value means, by definition. What
does the magnitude of an error have to do with that?

 

[Phil Carmody]

Phil Carmody wrote:

[Absolutely none of the quoted matter]

I wouldn't recommend using floating point math (which is what IEEE 754
is about) to do a search for prime numbers.

Do you have any particular insights into the field that Yves Gallot,
Jean Penne, George Woltman, Ernst Mayer, Guillermo Ballester Valor,
Richard Crandall, Dan Bernstein and myself don't have?

Just append " unless you know what you're doing" and I'll have no
problems agreeing, but the people who write the code that helps
people find prime numbers (i.e. the above listed people) do know
what they're doing, and exclusively use floating point.

Phil

 

[Phil Carmody]

Yep laughing. Because now you are saying it has error

You don't define what you mean by 'it' here. There are several things it
could refer to.

and before it was said it didn't.

I'm pretty sure you're making that up. But as you're not capable of
identifying clearly to what you're refering, it's hard to know.

Got an exact quote of me saying what you're claiming I said?

"Excat" is the word that was used.

Please engage brain before posting to usenet.

[SNIP - jibber-jabber by someone who knows nothing non-superficial about
the field]

Phil

 

[Keith Thompson]

Phil Carmody wrote:

[Absolutely none of the quoted matter]

I wouldn't recommend using floating point math (which is what IEEE 754
is about) to do a search for prime numbers.

Do you have any particular insights into the field that Yves Gallot,
Jean Penne, George Woltman, Ernst Mayer, Guillermo Ballester Valor,
Richard Crandall, Dan Bernstein and myself don't have?

Just append " unless you know what you're doing" and I'll have no
problems agreeing, but the people who write the code that helps
people find prime numbers (i.e. the above listed people) do know
what they're doing, and exclusively use floating point.

I find that very surprising; I would have thought that searching for
large prime numbers would call for the use of very large integers, not
floating-point (I assume we're talking about 64-bit floating-point).
To oversimplify tremendously, if you want to know whether
89587810616373014464315520632597 is prime, I don't see how being able
to represent 89587810616373014464315520632597.5 is helful. I presume
I'm missing something. (That's just a random number; it's probably
not prime.)

Can you cite any references for this, particularly anything that
doesn't require higher math to understand?

 

[Phil Carmody]

Phil Carmody wrote:

[Absolutely none of the quoted matter]

I wouldn't recommend using floating point math (which is what IEEE 754
is about) to do a search for prime numbers.

Do you have any particular insights into the field that Yves Gallot,
Jean Penne, George Woltman, Ernst Mayer, Guillermo Ballester Valor,
Richard Crandall, Dan Bernstein and myself don't have?

Just append " unless you know what you're doing" and I'll have no
problems agreeing, but the people who write the code that helps
people find prime numbers (i.e. the above listed people) do know
what they're doing, and exclusively use floating point.

I find that very surprising; I would have thought that searching for
large prime numbers would call for the use of very large integers, not
floating-point (I assume we're talking about 64-bit floating-point).
To oversimplify tremendously, if you want to know whether
89587810616373014464315520632597 is prime, I don't see how being able
to represent 89587810616373014464315520632597.5 is helful. I presume
I'm missing something. (That's just a random number; it's probably
not prime.)

You certainly want the numbers that you start with, and end with,
in a multi-precision multiplication to be exact integers, split over
many limbs.

However, in order to keep the Big-Oh of multiplication down (and
prime hunters typically deal with large enough numbers that the
Big-Oh dominates constants) you want to be using a FFT-based
scheme. And as you probably know, that involves hitting floating
point numbers. (There are things called NTTs which don't, but
nobody's got them to actually be as quick as FP methods yet, due
to CPU designers putting all their effort into FPUs recently.)

Square(a) = InvFFT(AutoConvolution(FFT(a)))
Mult(a,b) = InvFFT(Convolution(FFT(a),FFT(b)))

Ignoring loglog terms, FFTs are O(n.log(n)), convolution's O(n).
So you've got an 'almost linear' O(n^(1+eps)) multiplication.

For comparison, high-school multiplication is O(n^2), and
Karatsuba and Toom Cook are O(n^(1+significant))

So it's FFTs for the win every time (if your numbers are big enough).

The problem is that every operation in the above adds what's
usually a small amount of random noise to the values. Log(n)
steps in the FFT, an extra square, and log(N) steps in the
InvFFT, and you can see that the error's not going to grow too
quickly on average. However, it only takes one outlier with an
error over 0.5, and you will round to the wrong value (these
can sometimes be detected). If you're accumulating several
thousand products into 1 limb, each original limb can't be
over about 20 bits, so it initially looks like the FFTs are
wasteful on space, but they're worth it in the long run.

Can you cite any references for this, particularly anything that
doesn't require higher math to understand?

V. shallow - just implies that FFTs are useful:
http://en.wikipedia.org/wiki/Multiplication_algorithm#Fourier_transform_methods

With comments about rounding errors, but in unreadable MS-specific HTML:
http://numbers.computation.free.fr/Constants/Algorithms/fft.html
(ISTR it's vaguely fixable with either John Walker's or my demoroniser,
http://fatphil.org/perl/)

These have everything, including dozens more references.
Crandall & Pomerance's /Prime Numbers, a Computational Perspective/

Dan Bernstein - /Multidigit Multiplication for Mathematicians/
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.20.8270


Phil

 

[Richard Tobin]

It isn't. The sub-500 prime factors are 13, 41, 59, and 197. This
leaves us with 14461171494310710271575383, which I can't be
bothered to factor so I'll just announce that it's prime and see
what happens.

I'm rather disappointed that googling for "factor 14461171494310710271575383"
doesn't give the answer.

-- Richard

 

[Rich Webb]

Keith Thompson said:



It isn't. The sub-500 prime factors are 13, 41, 59, and 197. This
leaves us with 14461171494310710271575383, which I can't be
bothered to factor so I'll just announce that it's prime and see
what happens.

1201*1553*2605544659*35117*84737

 

[Ben Bacarisse]

Keith Thompson said:



It isn't. The sub-500 prime factors are 13, 41, 59, and 197. This
leaves us with 14461171494310710271575383, which I can't be
bothered to factor so I'll just announce that it's prime and see
what happens.

14461171494310710271575383 = 1201 * 1553 * 35117 * 84737 * 2605544659

in case anyone cares.

 

[Rich Webb]

Rich Webb said:


Proof by laziness. Gotta love it.

A wired world, indeed. Correlates to the famous Linus quote "Only wimps
use tape backup: _real_ men just upload their important stuff on ftp,
and let the rest of the world mirror it."

 

[James Kuyper]

Phil Carmody wrote:

[Absolutely none of the quoted matter]

I wouldn't recommend using floating point math (which is what IEEE 754
is about) to do a search for prime numbers.

Do you have any particular insights into the field that Yves Gallot,
Jean Penne, George Woltman, Ernst Mayer, Guillermo Ballester Valor,
Richard Crandall, Dan Bernstein and myself don't have?

Just append " unless you know what you're doing" and I'll have no
problems agreeing, but the people who write the code that helps
people find prime numbers (i.e. the above listed people) do know
what they're doing, and exclusively use floating point.

I claim no special expertise in the field; if you have indeed have it,
the choice of floating point arithemetic may indeed make sense; but it's
limited range and precision, and tendency to roundoff errors, makes it
seem a poor fit in the search for large integers possessing an
inherently integral property.

 

[CBFalconer]

.... snip ...

Repeat. In 5.2.4.2.2:

[#10] The values given in the following list shall be
replaced by implementation-defined constant expressions with
(positive) values that are less than or equal to those
shown:

-- the difference between 1 and the least value greater
than 1 that is representable in the given floating
point type, b1-p

FLT_EPSILON 1E-5
DBL_EPSILON 1E-9
LDBL_EPSILON 1E-9

Combine that with the method of expressing non-1.0 values.

Yes, and???

Did you even look at 5.2.4.2.2p2, which defines C's model for
floating-point values? (That's a serious question.)

Yes, I know what FLT_EPSILON, DBL_EPSILON, and LDBL_EPSILON are.
In type double, 1.0 is exactly representable, 1.0+DBL_EPSILON is
exactly representable, and the infinitely many real values
between them are not.

Also below, down to at least 1.0 - DBL_EPSILON/2.0.

You are agreeing with my claim. Those values 'between them that
are not' are what form the range of a FP value. Remember that I am
always talking about what you can get from the FP value, not from
the program as a whole.

[OT]
In mathematics you have various classes of numbers, including
positive integers, signed integers, rationals and reals, and
transcendentals. Many operations can be used on all, but some
cannot. In particular, some classes cannot express values in other
classes. It is at least 50 years since I took courses in all
these, but some ideas stick. In particular rationals, reals,
integers can all be put in one to one correspondence with the set
of integers. They form one form of infinity. Transcendentals
cannot, and are thus infinitely more 'dense' than the other forms.
PI, e, root(2) are all transcendental. Integers are not closed
sets under division. Reals are. FP values are another class.
[/OT]

5.2.4.2.2

[#4] The accuracy of the floating-point operations (+, -, *,
/) and of the library functions in <math.h> and <complex.h>
that return floating-point results is implementation
defined. The implementation may state that the accuracy is
unknown.

This handles all operations on FP values. EPSILONS provide an
accuracy limit.

 

[Keith Thompson]

... snip ...

Repeat. In 5.2.4.2.2:

[#10] The values given in the following list shall be
replaced by implementation-defined constant expressions with
(positive) values that are less than or equal to those
shown:

-- the difference between 1 and the least value greater
than 1 that is representable in the given floating
point type, b1-p

FLT_EPSILON 1E-5
DBL_EPSILON 1E-9
LDBL_EPSILON 1E-9

Combine that with the method of expressing non-1.0 values.

Yes, and???

Did you even look at 5.2.4.2.2p2, which defines C's model for
floating-point values? (That's a serious question.)

Which you haven't answered.

Also below, down to at least 1.0 - DBL_EPSILON/2.0.
Ok.

You are agreeing with my claim.

I certainly am not.

Those values 'between them that
are not' are what form the range of a FP value.

You keep asserting that. You have yet to support your assertion with
reference to the standard. The things you've cited from the standard
do not support your claims.

[snip]

 

[Guest]

... snip ...

Repeat. In 5.2.4.2.2:
  [#10] The values  given  in  the  following  list  shall  be
  replaced by implementation-defined constant expressions with
  (positive) values that are  less  than  or  equal  to  those
  shown:
    -- the  difference  between  1 and the least value greater
       than 1 that is  representable  in  the  given  floating
       point type, b1-p
       FLT_EPSILON                   1E-5
       DBL_EPSILON                   1E-9
       LDBL_EPSILON                  1E-9
Combine that with the method of expressing non-1.0 values.
Yes, and???
Did you even look at 5.2.4.2.2p2, which defines C's model for
floating-point values?  (That's a serious question.)

Which you haven't answered.

Also below, down to at least 1.0 - DBL_EPSILON/2.0.
Ok.

You are agreeing with my claim.

I certainly am not.

                                 Those values 'between them that
are not' are what form the range of a FP value.

You keep asserting that.  You have yet to support your assertion with
reference to the standard.  The things you've cited from the standard
do not support your claims.

But it's *obvious*, therefore it doesn't need to said in the standard.

 

[TonyMc]

PI, e, root(2) are all transcendental.

No - root(2) is irrational but not transcendental.


Tony

 

[Richard Tobin]

No - root(2) is irrational but not transcendental.

In fact it must be the most obvious non-transcendental irrational.

-- Richard

 

[Dik T. Winter]

>
> Huh? What do you mean. Seems to me that if you multiply two
> 'ranges' together you will get a result with a third 'range'. Do
> the algebra.

But the new range is *not* the range you get by multiplying the result
by x_EPSILON. Do the algebra.

 

[Dik T. Winter]

> Yes, but that is not necessarily what was offered to be stored. For
> example, use double to compute:
>
> double x = 1.0 + (FLT_EPSILON/2.0)
>
> which is calculable, and is not equal to the double 1.0.
> the result in a float.
>
> float y = x;
>
> Now try to extract something other than 1.0 from y.
>
> I haven't read it for many years, but there is a whole section in
> Knuth about floating point errors, and keeping track.

There is a whole book by J. H. Wilkinson about this subject. Keeping
track will *not* give you reasonable results in many cases. You should
analyse.

 

[Dik T. Winter]

>
> This is a chimera. The errors are small, compared with the
> values. The sort of effect you are worrying about is of the form
> error_ratio squared, and thus negligible.

You are completely wrong. Do some study about error-analysis,
interval-arithmetic and so on.

 

[Carl]

In fact it must be the most obvious non-transcendental irrational.

-- Richard

If only we used a base-pi number system...