simple algorithm for finding primes

[Richard Heathfield]

(e-mail address removed) says:

Well, as you can imagine my considerations of whether or not Richer
Heathfield can find some violation of the ANSI standard in my code is not
exactly high on
my todo list. I've got navel lint that needs sorting.

Oh, I see. Since you posted it in comp.lang.c, I made a not unnatural, but
nevertheless apparently incorrect, assumption that the code was intended to
be portable. If it were not so intended, then of course it's not surprising
that I got compilation errors. What /is/ surprising is that you should post
(a link to) such code here in comp.lang.c.

I am not typically prone to writing non-portable code just for the fun of
it. Using __int64 is not a "performance concern" (it rarely is) -- its a
functionality consideration. Using "long int" which may be implemented by
a compiler to be 32 bits, is not an acceptable compromise -- it would
reduce the correctness of the algorithm to a range of <65536 (the Sieve of
Eratosthenes,
by itself can easily go beyond this.)

Yes, I managed to get a sieve working cheerfully even right up to one
million even under MS-DOS (admittedly I saved lots of space by testing 2
separately and throwing out all other even numbers).

This is the whole "widening
multiply" thing all over again (which I won't regurgitate here, as it
always seems to take 10 posts of explanation before people "get it".)

But since you asked so nicely, I have made a correction which uses
*doubles* in
place of __int64 if compiled with an unrecognized compiler. This allows
it to work for inputs that are up to 26 bits on most typical machines with
a 64 bit IEEE-754 double implementation (it would also work for all 32bits
for typical
80+ bit FP implementations). And I threw in some "guesswork code" about
how it might work on gcc.

Why not just use bignums, which can be written 100% portably and with 100%
accuracy?

So its now more "portable", but that won't help it be any more correct on
platforms without a 64 bit integer (or 80 bit FP.)

Solutions are possible that are correct and portable even to such platforms.

 

[Arthur J. O'Dwyer]

Oh, I see. Since you posted it in comp.lang.c, I made a not unnatural, but
nevertheless apparently incorrect, assumption that the code was intended to
be portable.

You know, Richard, I was *trying* to defuse the situation.
Anyway, the point is that Paul (AFAICT) writes code for Pentiums. It
seems to be "what he does." It is slightly unreasonable to expect
Paul to re-write the Pentium code on his website -- which /is/ *mostly*
portable C code -- every time he wants to give an example of how
something could be done in (*completely* portable) C.
It's not good to post code that isn't portable -- especially if
it's unportable enough to actually refuse to compile on many
implementations -- in comp.lang.c. But Paul certainly added to the
discussion, and it would be fairly trivial to "portabilize" his code
when or if the OP needed it. All that was required, IMHO, was a little
note: "__int64 is not standard C. Replace it with some type that has
64 or more bits of integer precision on your system."


Given a suitably Turing-complete programming language, *everything* is
a performance concern. ;-) In this case, one portable alternative would
be to compute the "widening multiply" in halves (most significant and
least significant 32-bit sections). I don't know how you would do this
in C, but I assure you it can be done.
Richard, in his flamboyant way, suggests bignums.


I think I get it, and I know why you used __int64 in the first
place. I merely think that you should have given a slightly larger
nod to portability before posting the code in a group dedicated to
*portable* code.


<OT> Does gcc have C99's type int64_t (or whatever it really is)? </OT>

-Arthur

 

[Richard Heathfield]

On Wed, 10 Dec 2003, Richard Heathfield wrote:


You know, Richard, I was *trying* to defuse the situation.

So was I, believe it or not. The paragraph you quote above was intended to
explain to Paul that I had jumped to an incorrect conclusion about the
intended scope of the code he posted.

Anyway, the point is that Paul (AFAICT) writes code for Pentiums. It
seems to be "what he does." It is slightly unreasonable to expect
Paul to re-write the Pentium code on his website -- which /is/ *mostly*
portable C code -- every time he wants to give an example of how
something could be done in (*completely* portable) C.

Um, no, it's not unreasonable at all. If he wants to give an example of how
something could be done in *completely* portable C, then nothing is
preventing him from doing that, but I don't see how he can provide a
portable example using /non/-portable code.

It's not good to post code that isn't portable -- especially if
it's unportable enough to actually refuse to compile on many
implementations -- in comp.lang.c.
Right.

But Paul certainly added to the
discussion, and it would be fairly trivial to "portabilize" his code
when or if the OP needed it.

I don't dispute either of these points. Nor was I trying to get on Paul's
case, particularly. I was not being disingenuous in my original reply; I
simply made an incorrect assumption (viz. that the posted code was intended
to be portable).

<snip>

 

[Paul Hsieh]

Arthur J. O'Dwyer says...

You know, Richard, I was *trying* to defuse the situation.
Anyway, the point is that Paul (AFAICT) writes code for Pentiums.

As is evidenced by the feedback I get for http://bstring.sf.net/ from people
who use Linux and Mac OS X. I have a bunch of stuff on my website about
Pentiums in my assembly pages. And __int64 is a MSVC/WATCOM-ism, not a
"Pentiumism". Linux on Pentium and Borland C people are likely to have as many
objections with my use of it as anyone else.

The OP was clearly writing DOS or Windows code, and my response was clearly
directed at the OP, not anyone else. I was trying to tell the OP that this is
not a group where algorithms or performance can be sensibly discussed. The
fact that *I* am posting and reading this newsgroup is a complete fluke, and so
I was trying to indicate to the OP should probably broaden his/her horizons
when seeking those kind of answers.

[...] It
seems to be "what he does." It is slightly unreasonable to expect
Paul to re-write the Pentium code on his website -- which /is/ *mostly*
portable C code -- every time he wants to give an example of how
something could be done in (*completely* portable) C.

No, I would prefer to write portable code, because my occasional use of
"Pentium code" is limited to when I can't find a realistic way to do without,
or when I am really just analyzing a Pentium (or one of its clones.) Writing
"Pentium-C" is tantamount to writing assembly -- and I do it sometimes, but its
not because I have some aversion to writing portable code.

It's not good to post code that isn't portable -- especially if it's
unportable enough to actually refuse to compile on many implementations --
in comp.lang.c. But Paul certainly added to the discussion, and it would be
fairly trivial to "portabilize" his code when or if the OP needed it.

No its not, or I would have done so. For example, I've now changed it to code
that will *compile* portably, but will not *function* optimally. On platforms
that don't have a 64 bit integer, the code will produce incorrect results. In
fact, because of the way it was written, it will actually fail using the native
compilers on Sparc, MIPS, PA-RISC and Alpha even though *ALL* those CPUs have
64bit integers with all the right ALUs.

[...] All that was required, IMHO, was a little note: "__int64 is not
standard C. Replace it with some type that has 64 or more bits of integer
precision on your system."

Well, that's what I started with. There still a note there.

Given a suitably Turing-complete programming language, *everything* is a
performance concern. ;-) In this case, one portable alternative would be to
compute the "widening multiply" in halves (most significant and least
significant 32-bit sections).

And how will that help the fact that there is also a shortening divide in
there?

A bignum divide that I have written on top of a bignum multiply is 98 lines
with a big while-loop around most of it. And that's only because I realized
that you could use FP approximations to sneakily perform the estimation step.
Of course, its a modulo divide so there's probably a better approach, but
either way its not going to be very pretty.

That's the alternative to the 1-line divide in the source I showed.

I would have brought this up earlier, but people in this newsgroup have enough
trouble understanding widening multiplies, I didn't think bringing this up was
exactly going to register with anyone.

[...] I don't know how you would do this in C, but I assure you it can be
done.

Its algebra:

((a<<16) + b) * ((c<<16) + d) == ((a*c)<<32) + ((a*d + b*c)<<16) + b*d

I'm the last person that needs to be assured that this can be done.

Richard, in his flamboyant way, suggests bignums.

No -- in his *hyppocritical* way, he suggest bignums. Making bignums portable
is tantamount to reducing it to the Turing machine that you suggested tongue
in cheek. At some point you have to stop and recognize when a programming
language has failed you, and in the case of ANSI C, bignum-like issues show it
for all its inadequacies.

I think I get it, and I know why you used __int64 in the first place.

You do? You seemed to have missed the shortening divide that's in there.

[...] I merely think that you should have given a slightly larger nod to
portability before posting the code in a group dedicated to *portable* code.


<OT> Does gcc have C99's type int64_t (or whatever it really is)? </OT>

Who cares? C99 is *far* less portable than anything I've written.

 

[Richard Heathfield]

Paul Hsieh wrote:

No -- in his *hyppocritical* way, he suggest bignums.

I don't deny being any less susceptible to hypocrisy than the next man (and
it would be hypocritical to claim otherwise), but your claim that a
suggestion to use bignums is somehow hypocritical has completely flummoxed
me.

Making bignums portable is tantamount to reducing it to the Turing
machine that you suggested tongue in cheek.

Actually, it isn't all that difficult to make bignums portable. There is no
particular reason why you shouldn't be able to implement, say, RSA using
code that will work on any conforming hosted implementation (and, indeed,
on freestanding implementations if they have the RAM for it; if you have a
few kilobytes, you can do it).

At some point you have to stop and recognize when a programming
language has failed you, and in the case of ANSI C, bignum-like issues
show it for all its inadequacies.

About the only nuisance I've encountered in this regard is the lack of
operator overloading which, in this particular field, would be very welcome
indeed.

<snip>

 

[osmium]

Paul Hsieh wrote:



I don't deny being any less susceptible to hypocrisy than the next man (and
it would be hypocritical to claim otherwise), but your claim that a
suggestion to use bignums is somehow hypocritical has completely flummoxed
me.


Actually, it isn't all that difficult to make bignums portable. There is no
particular reason why you shouldn't be able to implement, say, RSA using
code that will work on any conforming hosted implementation (and, indeed,
on freestanding implementations if they have the RAM for it; if you have a
few kilobytes, you can do it).

I am just an innocent bystander here. But you, Richard, seem to be
comparing a program that *could* be written with one that *has* been
written. I notice you always refer to bignum in the future tense. If you
are speaking of an actual library it most likely has an actual name and a
parent. My guess is that is where the alleged hypocrisy comes in.

If I were to write a reasonably fast bignum library I would sure like to see
that well hidden carry bit.

 

[Mark McIntyre]

Arthur J. O'Dwyer wrote:

(snip)


It just occurred to me how nice it would be to have a statement,
probably for the preprocessor, that would cause a message to
be printed when it compiled. Such warnings could then be
put in that way.

many compilers offer an extension #warning. I tend to use it to
document places I'd like to do better work on, but don't have time to
address right now.

 

[Jim]

It would seem that preprocessed C could also have messages
that need to be printed.

#ifdef NUM<0
#note Oops, NUM is less than zero!
#else
int x[NUM];
#endif

#if NUM<0
#error Oops, NUM is less than zero!
#else
int x[NUM];
#endif

Not quite what I wanted. That seems to cause a fatal error.

I suppose the example isn't very good, but there are cases when
a warning, or even less. Maybe just a reminder to the person
compiling the program.

-- glen

Some compilers do this with a pragma. Two compilers I'm using allow
this:
#pragma message("remember to initialise this structure!")

Jim

 

[Richard Heathfield]

I am just an innocent bystander here. But you, Richard, seem to be
comparing a program that *could* be written with one that *has* been
written.

Well, if he'd written a portable bignum program already, he'd have said so,
wouldn't he?

I notice you always refer to bignum in the future tense. If you
are speaking of an actual library it most likely has an actual name and a
parent. My guess is that is where the alleged hypocrisy comes in.

Ah, I see. Well, we already have Miracl and GMP, do we not? And yes, I have
my own bignum library, which works (except for Rabin-Miller, which I hope
to fix today), but which is still a bit untidy, which is why I haven't
published it yet. The fastest in the world it ain't (because the emphasis
is on readability and portability), but it's fast enough that it is
unlikely to be the bottleneck in typical usages.

If I were to write a reasonably fast bignum library I would sure like to
see that well hidden carry bit.

Perhaps I'm talking through ignorance here, despite having written one
myself, but why would anyone need to *hide* a carry bit?

 

[Grumble]

glen herrmannsfeldt wrote:
...
It would seem that preprocessed C could also have messages
that need to be printed.

#ifdef NUM<0
#note Oops, NUM is less than zero!
#else
int x[NUM];
#endif

#if NUM<0
#error Oops, NUM is less than zero!
#else
int x[NUM];
#endif

Not quite what I wanted. That seems to cause a fatal error.

I suppose the example isn't very good, but there are cases when
a warning, or even less. Maybe just a reminder to the person
compiling the program.

Some compilers do this with a pragma.

Some have a #warning preprocessor directive.
http://gcc.gnu.org/onlinedocs/cpp/Diagnostics.html

 

[osmium]

Well, if he'd written a portable bignum program already, he'd have said so,
wouldn't he?


Ah, I see. Well, we already have Miracl and GMP, do we not? And yes, I have
my own bignum library, which works (except for Rabin-Miller, which I hope
to fix today), but which is still a bit untidy, which is why I haven't
published it yet. The fastest in the world it ain't (because the emphasis
is on readability and portability), but it's fast enough that it is
unlikely to be the bottleneck in typical usages.


Perhaps I'm talking through ignorance here, despite having written one
myself, but why would anyone need to *hide* a carry bit?

From the Miracl page:

"MIRACL is compact, fast and efficient and its
now easier than ever to get the same near-optimal
performance from any processor. Although
essentially a portable library, inline assembly
and special techniques can be invoked for
blistering speed"

It would seem to me that the designers of Miracl also wanted to see that
carry out bit which the designers of the C language hid. And the inline
assembly destroys the portability. It is not clear to me that a sensible,
fast, binary, arithmetic big numbers library is possible. I say binary to
exclude BCD arithmetic which is a different kettle of fish. I read
"blistering speed" as meaning reasonably fast.

I will let someone else research GMP. My algorithm is first in, first out.

If I were designing C, I would have considered adding the carry out (and
perhaps other things) as a standard signal. But even so, the resulting big
numbers library would still be slower than necessary.

 

[Wolfgang Riedel]

Arthur J. O'Dwyer says...

No -- in his *hyppocritical* way, he suggest bignums.

Paul,

maybe, he has every right, to be critical of hippos?!

Btw., I think, you do yourself no good with those hypocritical
ad hominem attacks.

Let's stay civilized.

Wolfgang

 

[Chris Torek]

If I were designing C, I would have considered adding the carry out (and
perhaps other things) as a standard signal.

Carry-out of unsigned addition of two operands is trivial to pick
up:

unsigned long a, b, result;
int carry;

result = a + b;
carry = result < a; /* or result < b */

In assembly languages for CPUs without carry flags (e.g., MIPS),
this C code is often a direct translation of the required
assembly code:

addu rRESULT, rA, RB # set register RESULT = A + B
slt rCARRY, rRESULT, rA # set register CARRY = RESULT < A

and on CPUs with carry flags, a compiler should be able to optimize
the "carry = " operation appropriately. If the carry is then fed
into a subsequent op:

r2 = a2 + b2 + carry;

it should even be able to emit an "ADC" instruction, if one exists.
Of course, obtaining carry-out of 3-or-more operand add instructions
is more difficult to achieve, but quite a few modern CPUs have no
three-operand-input operations (including no "add with carry"
instructions; indeed, MIPS, as mentioned above, has no carry flag
at all). Still:

carry = r2 < a2 || (r2 == a2 && carry);
/* or equivalently, carry = (r2 < a2) | ((r2 == a2) & carry); */

(since "carry" is either 0 or 1 initially).

(Whether your compiler does such optimizations is another question
entirely. If those optimizations are important to you, I suggest
that you could put some time, money, and/or effort into seeking
out -- or creating, e.g., by modifying gcc -- compilers that do
them.)

(V9 SPARC points out one reason MIPS does not have a carry flag.
V8 and earlier had only the "condition codes" register, %cc, with
the usual four condition codes: N V C and Z => negative, overflow,
carry, and zero. V9 has %icc, the 32-bit condition codes, and
%xcc, the 64-bit condition codes. A 64-bit register holding a
result of 0x0000000087654321 is negative in 32-bits, but positive
in 64-bits. The sum of 0x00000000ffffffff and 0x00000000ffffffff
is 0x00000001fffffffe, which has no carry in 64-bits but a carry
in 32-bits.)