primes returns composites

[b49P23TIvg]

Binary package hint: bsdgames

Yesterday I reported that primes on my hardware yields values beyond the documented range. The man page claims:

DIAGNOSTICS
     Out of range or invalid input results in an appropriate error message
     being written to standard error.

I doubt that it's a security vulnerability since primes isn't on the live CD.

Thus this should result in an error message. Instead primes returns a composite number.

$ primes 7757777773 7757777787
7757777777

In j notation, q: is the factor verb and */ means "insert multiply". The precedence rule for sequences of nouns and verbs is simply right to left evaluation.

   q:7757777777x
65827 117851
   */q:7757777777x
7757777777

(Or you could try Factor[7757777777] at Wolfram Alpha, or use a calculator.)

ProblemType: Bug
DistroRelease: Ubuntu 10.10
Package: bsdgames 2.17-19
ProcVersionSignature: Ubuntu 2.6.35-25.44-generic 2.6.35.10
Uname: Linux 2.6.35-25-generic x86_64
Architecture: amd64
Date: Fri Feb 25 19:45:06 2011
InstallationMedia: Ubuntu 10.10 "Maverick Meerkat" - Release amd64 (20101007)
ProcEnviron:
 PATH=(custom, no user)
 LANG=en_US.utf8
 SHELL=/bin/bash
SourcePackage: bsdgames

[Santtu Lakkala (inz)]

Indeed, the primes program uses Eratosthenes sieve using pre-calculated array of primes up to 65537, leading to possible errors (and, in this case real errors, as 65539 is also a prime) when range hits (65537+2)^2.

There are a few possible fixes, that I can think of:
 a) don't support 64-bit numbers, set upper limit to UINT32_MAX
 b) support 64-bit numbers, but set upper limit to (MAX(primes) + 2)**2 - 1
 c) add some number of pre-calculated primes to array, and repeat b (supporting the whole 64-bit range would take 203280221 primes (which could be stored in a 4-byte uints, leading to ~800 megabytes of memory)
 d) recursively use the sieve to prime itself to reach unlimited range (this adds computation time and algorithm complexity)

[Santtu Lakkala (inz)]

The attachment implements option b) add some to primes table (max 2^17) and limit range to achieve approx. 2^34 bit result range.

[Santtu Lakkala (inz)]

Erm, the previous was obviously solution c), not b) =)

Anyway, here's d), it uses the sieve recursively to init itself. This should allow full 64-bit range to be traversed (and potentially larger, if there is a native integer type).

[Alexander A. Maly (a-a-maly)]

Project https://launchpad.net/primes/ was designed to fix this issue. It contains primes32, that support only 32-bit numbers, and primes64, that supports 64-bit numbers. It is somewhat incompatible with the original primes program, as it, if only 1 argument x given, prints prime numbers from 0 to x, not from x to a kind of infinity.
As well, it has additional keys to dump results in binary or/and difference format to improve overall performance.
Good testing is needed!

[b49P23TIvg (b49p23tivg)]

Thank you for investigating and writing new code. I built the code from
your repository trunk. primes32 might be correct, assuming you added
only a test for the argument size. primes64 didn't work. 7757777775
is, by inspection, divisible by 5.

Thank you again,
Dave Lambert

$ primes32 7757777773 7757777787 # good!
argument '7757777773' incorrect at position 9
$
$ primes64 7757777773 7757777787 # no!
start=7757777773, stop=7757777787, mode=0x0
q_min=3878888886, q_lim=3878888893
sqrt(7757777785)=88078
q_len=7, maxsp=88078
tmul=1, tdiv=255255, sres=8563
p_min=9, p_lim=44039, p_len=44030
sqrt(88078)=296
maxssp=296, si_min=9, si_lim=148
Primes::init _pmax=296
nssevs=55
cnt_d=0, cnt_m=44030
sdevs.size=8542
cnt_dd=0
7757777775
totally 1 primes with sum=7757777775:0
$

[Ubuntu Foundations Team Bug Bot (crichton)]

The attachment "Solution b" of this bug report has been identified as being a patch. The ubuntu-reviewers team has been subscribed to the bug report so that they can review the patch. In the event that this is in fact not a patch you can resolve this situation by removing the tag 'patch' from the bug report and editing the attachment so that it is not flagged as a patch. Additionally, if you are member of the ubuntu-reviewers team please also unsubscribe the team from this bug report.

[This is an automated message performed by a Launchpad user owned by Brian Murray. Please contact him regarding any issues with the action taken in this bug report.]

[Alexander A. Maly (a-a-maly)]

Thank you for assistance!
That seems to be output from an older version (before rev. 9) of the program.
How get you the code? I tried bzr branch lp:primes/ from a clean place, built via make, run it as

./primes64 7757777773 7757777787

and it prints:

start=7757777773, stop=7757777787, mode=0x0, loglevel=1
totally 0 primes, last one is 7757777773
sum is 0 + 2^64 * 0

Or maybe you built it in some unusual environment, like mingw?

[b49P23TIvg (b49p23tivg)]

OH! I had gotten revision 1. ~a-a-maly_primes_trunk-r1.tgz
Now I've got -r9 .

J's implementation of Miller-Rabin test confirms 10,000 ten digit primes.

I originally found the problem when I used primes for a problem at ProjectEuler.net and came up with an incorrect answer. My algorithm for solution looked right, right, right so I finally investigated primes.

[Alexander A. Maly (a-a-maly)]

I wrote this fast implementation for PE problems as well, but mostly for SOLSTRAS challenge at spoj.pl.
It is based on "memoryless" variant of the sieve of Eratosthenes (over priority queue). I heard that sieve of Atkins may be much faster, but cannot investigate it right now.
A good test may be to diff primes64 output with the output of current primes program from bsdgames package, then testing difference by Miller-Rabin.

[b49P23TIvg (b49p23tivg)]

The games package primes, your primes32 and your primes64 agree with start 1, stop 4294967295.

Also, 2892 randomly chosen primes given by

$ ./primes64 4611686018427388039 4611686018428387862

all tested prime by j's Miller-Rabin test.

[Hubert depesz Lubaczewski (depesz)]

I know this is not the most world-shattering bug, but any chance we could see a fix for this?

primes still (in 2016!) shows composites, which can be seen by doing:

=$ primes 35198908481059600 35198908481059602
35198908481059601

=$ factor 35198908481059601
35198908481059601: 181399 405277 478787

=$ bc <<< "181399 * 405277 * 478787 - 35198908481059601"
0