Givaro
examples/Integer/probable_primroot.C

NO DOC

// Copyright(c)'1994-2009 by The Givaro group
// This file is part of Givaro.
// Givaro is governed by the CeCILL-B license under French law
// and abiding by the rules of distribution of free software.
// see the COPYRIGHT file for more details.

#include <iostream>
#include <givaro/givintnumtheo.h>
#include <givaro/givtimer.h>
#include <math.h>

using namespace Givaro;



//  Polynomial-time generation of primitive roots
//  L is number of loops of Pollard partial factorization of n-1
//  10,000,000 gives at least 1-2^{-40} probability of success
//  [Dubrois & Dumas, Industrial-strength primitive roots]
//  Returns the probable primitive root and the probability of error.

int main(int argc, char** argv)
{
    IntNumTheoDom<> IP;
#ifdef __GMP_PLUSPLUS__
    IP.seeding( BaseTimer::seed() );
#endif

    double error;
    IntNumTheoDom<>::Element a,pr,g;
    if (argc > 1) a = IntNumTheoDom<>::Element(argv[1]); else std::cin >> a;
    bool comp ; if ( (comp=(! IP.isprime(a))) )  std::cerr << a << " is not prime, primitive root will have no sense and may loop forever ..." << std::endl;
    double epsilon = argc > 2 ? atof(argv[2]) : 0.0000001;

    Timer tim; tim.clear();
    tim.start();
        //======================================================================
        // Default is partial factorization up to factors of at least 12 digits.
        // with probability of error much less than 2^{-40}
        // IP.probable_prim_root(pr, error, a );


        //======================================================================
        // Choosing L to be O(log^2(p))
        // gives the best probability with O(log^4(p)) complexity
        // Probability of error is approximately O(1/log^4(p))
        // IP.probable_prim_root(pr, error, a, (unsigned long)power(logtwo(a),2));

        //======================================================================
        // Choosing L to be O( \sqrt(epsilon) )
        // gives probability of error at most epsilon
        // Newton-Raphson iteration is used for
        // 1-epsilon = (1+2/(p-1))*(1-1/B)^(ln( (p-1)/2 )/ln(B))
        // So that no factor less than B can be avoided
        // With Pollard's rho factorization, L is chosen to be sqrt(B)
        // Might not be polynomial if epsilon is too big
#define GIVARO_POLLARD
    IP.probable_prim_root(pr, error, a, epsilon );
    tim.stop();

    if (comp) std::cerr << IP.gcd(g,a,pr) << " is a factor of " << a << std::endl;
    IntegerDom().write( std::cout << "Prim root   : ", pr );
    if (error > 0) {
        std::cout << ", correct with probability at least : 1-" << error << std::endl;
        std::cerr << tim << std::endl;

        std::cerr << "Now checking primitivity, this may take some time (complete factorization of n-1) ...";

#define GIVARO_LENSTRA

        Timer verif; verif.clear(); verif.start();
        if ( IP.isorder(a-1, pr, a) ) {
            verif.stop();
            std::cerr << "... Pimitivity checked" << std::endl;
            std::cerr << verif << std::endl;
        }
        else {
            verif.stop();
            std::cerr << "... WARNING : FAILURE" << std::endl;
            std::cerr << verif << std::endl;
        }
    } else {
        std::cout << ", deterministically correct" << std::endl;
        std::cerr << tim << std::endl;
    }

    return 0;
}