RandomFFTPrime.randomPrime does not respect prime_bound

[cyrilbouvier]

The class RandomFFTPrime has a method randomPrime taking k as parameter and returning a prime such that p-1 is divisible by 2^k. The constructor of the class RandomFFTPrime has a parameter pbound. One could expected that the prime returned by an instance of this class initialized with pbound=2^b returns a random prime less than 2^b but it is not the case. I wrote the simple test below. With the two following sets of arguments, it returns a prime that is larger than 2^b:

• -b 27 -k 16 -s 42 : the test program returns 2^27 + 33554433
• -b 16 -k 15 : the test program returns 2^16+1

The second case is more annoying. Let's assume I want a prime less than 16 bits in order to perform the FFT with Modular<uint16_t, uint32_t>. Calling RandomFFTPrime(1<<16).randomPrime (15) returns 2^16+1 which does not fit in an uint16_t. [ Note that there is no prime p less than 2^16 such that 2^15 | p-1, the function fails to detected this but this is another problem].

After reading the code, it seems that, given pbound as parameter, the actual bound is 2^floor(log_2(pbound)). So with pbound=2^b the actual bound is 2^(b+1).

Another problem is the fact that if k>=b (in the clear case where there is no solution) and that the code is compile with -NDEBUG, the method randomPrime fails with a Floating point exception. It would be better to raise a more explicit exception or to return 0 (instead of a prime).

Finally, if the code is compile with -NDEBUG (this implies that the linbox_check macro is an empty macro), the method randomPrime could return a non prime without any indication that it fails to find a prime.

Code used for testing:

#include "linbox/linbox-config.h"

#include "linbox/randiter/random-fftprime.h"
#include "fflas-ffpack/utils/args-parser.h"

using namespace LinBox;

int main (int argc, char *argv[]) {
    unsigned long bits = 27;
    unsigned long k = 16;
    unsigned long seed = time (NULL);

    Argument args[] = {
        { 'b', "-b nbits", "number of bits of prime.", TYPE_INT, &bits },
        { 'k', "-k k", "bench FFT for n=2^k.", TYPE_INT, &k },
        { 's', "-s seed", "set the seed.", TYPE_INT, &seed },
        END_OF_ARGUMENTS
    };
    FFLAS::parseArguments (argc, argv, args);

    std::cout << "seed = " << seed << std::endl;
    std::cout << "bits = " << bits << std::endl;
    std::cout << "k = " << k << std::endl;

    RandomFFTPrime RandomGen (1<<bits, seed);
    integer p = RandomGen.randomPrime (k);
    std::cout << "p = " << p << std::endl;

    bool b1 = (p >> bits) == 0;
    std::cout << "p < 2^bits ? " << (b1 ? "true" : "false") << std::endl;

    integer mask = 1 << k; mask -= 1;
    bool b2 = ((p-1) & mask) == 0;
    std::cout << "2^k | p-1 ? " << (b2 ? "true" : "false") << std::endl;
    return (b1 & b2) ? 0 : 1;
}

[cyrilbouvier]

Fixed by PR #177