Givaro
examples/Polynomial/PolynomialCRT.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 <algorithm>
#include <givaro/givtimer.h>
#include <givaro/givpoly1crt.h>
#include <givaro/givintprime.h>
#include <givaro/montgomery.h>
#include <givaro/extension.h>
#include <givaro/modular.h>
#include <givaro/gfq.h>
#include <givaro/chineseremainder.h> // Chinese Remainder of two elements
#include <givaro/givrns.h> // Chinese Remainder of an array of elements
#include <givaro/givrandom.h>
#include <givaro/givrational.h>
using namespace Givaro;
template <typename Field>
bool tmain(int argc, char ** argv, GivRandom& generator) {
bool pass = true;
typedef Poly1CRT< Field > CRTSystem;
typedef typename CRTSystem::Element Poly;
// typedef typename CRTSystem::Type_t Scal;
// typedef typename CRTSystem::array_E VPoly;
typedef typename CRTSystem::array_T VScal;
Integer a( generator() >>(argc>2?atoi(argv[2]):17) );
Field F(ID.nextprimein( a ));
VScal Primes( argc>1 ? (size_t)atoi(argv[1]):15);
VScal Moduli( Primes.size() );
typename VScal::iterator i = Primes.begin();
typename VScal::iterator e = Moduli.begin();
for(; i != Primes.end(); ++i, ++e) {
do {
F.init(*i, generator());
} while ( (std::find(Primes.begin(), i, *i) != i) || (F.isZero(*i))) ;
F.init(*e, generator());
}
// for(typename VScal::const_iterator it=Primes.begin(); it!=Primes.end();++it)
// F.write(std::cout, *it) << std::endl;
// for(typename VScal::const_iterator it=Moduli.begin(); it!=Moduli.end();++it)
// F.write(std::cout, *it) << std::endl;
CRTSystem CRT( F, Primes, "X" );
Poly res;
Timer tim; tim.clear(); tim.start();
CRT.RnsToRing( res, Moduli );
tim.stop();
F.write( std::cerr << tim << " using ") << std::endl;
if (Primes.size() < 14) {
i = Primes.begin();
e = Moduli.begin();
for( ; i != Primes.end(); ++i, ++e)
if (F.characteristic()>0)
F.write(CRT.getpolydom().write(F.write(std::cout << "subs(X=", *i) << ",", res) << ") mod " << F.characteristic() << " = ", *e) << ';' << std::endl;
else
F.write(CRT.getpolydom().write(F.write(std::cout << "subs(X=", *i) << ",", res) << ") = ", *e) << ';' << std::endl;
}
VScal Verifs( Primes.size() );
CRT.RingToRns( Verifs, res );
typename VScal::const_iterator v = Verifs.begin();
e = Moduli.begin();
for( ; e != Moduli.end(); ++e, ++v)
if (! F.areEqual(*e, *v) ) {
F.write(std::cerr << "incoherency within ") << std::endl;
F.write(std::cerr << "e: ", *e ) << std::endl;
F.write(std::cerr << "v: ", *v ) << std::endl;
pass = false;
break;
}
CRT.getpolydom().random(generator, res, Degree((int64_t)Primes.size()-1));
CRT.RingToRns( Verifs, res );
Poly nres;
tim.clear(); tim.start();
CRT.RnsToRing( nres, Verifs );
tim.stop();
if (! CRT.getpolydom().areEqual(res,nres) ) {
CRT.getpolydom().write(std::cerr << "incoherency within ") << std::endl;
CRT.getpolydom().write(std::cerr << "r: ", res ) << std::endl;
CRT.getpolydom().write(std::cerr << "n: ", nres ) << std::endl;
pass = false;
}
F.write( std::cerr << tim << " using ") << std::endl;
return pass;
}
template <typename Field>
bool tmainext(int argc, char ** argv, GivRandom& generator) {
bool pass = true;
typedef Poly1CRT< Field > CRTSystem;
typedef typename CRTSystem::Element Poly;
// typedef typename CRTSystem::Type_t Scal;
// typedef typename CRTSystem::array_E VPoly;
typedef typename CRTSystem::array_T VScal;
Integer a( generator() >>(argc>2?atoi(argv[2]):17) );
Field F(ID.nextprimein( a ),2);
VScal Primes( argc>1 ? (size_t)atoi(argv[1]):15);
VScal Moduli( Primes.size() );
typename VScal::iterator i = Primes.begin();
typename VScal::iterator e = Moduli.begin();
for(; i != Primes.end(); ++i, ++e) {
F.init(*i); F.init(*e);
do {
F.random(generator,*i);
} while ( (std::find(Primes.begin(), i, *i) != i) || (F.isZero(*i))) ;
F.random(generator,*e);
}
for(typename VScal::const_iterator it=Primes.begin(); it!=Primes.end();++it)
F.write(std::cout, *it) << std::endl;
for(typename VScal::const_iterator it=Moduli.begin(); it!=Moduli.end();++it)
F.write(std::cout, *it) << std::endl;
CRTSystem CRT( F, Primes, "X" );
Poly res;
Timer tim; tim.clear(); tim.start();
CRT.RnsToRing( res, Moduli );
tim.stop();
F.write( std::cerr << tim << " using ") << std::endl;
if (Primes.size() < 14) {
i = Primes.begin();
e = Moduli.begin();
for( ; i != Primes.end(); ++i, ++e)
if (F.characteristic()>0)
F.write(CRT.getpolydom().write(F.write(std::cout << "subs(X=", *i) << ",", res) << ") mod " << F.characteristic() << " = ", *e) << ';' << std::endl;
else
F.write(CRT.getpolydom().write(F.write(std::cout << "subs(X=", *i) << ",", res) << ") = ", *e) << ';' << std::endl;
}
VScal Verifs( Primes.size() );
CRT.RingToRns( Verifs, res );
typename VScal::const_iterator v = Verifs.begin();
e = Moduli.begin();
for( ; e != Moduli.end(); ++e, ++v)
if (! F.areEqual(*e, *v) ) {
F.write(std::cerr << "incoherency within ") << std::endl;
F.write(std::cerr << "e: ", *e ) << std::endl;
F.write(std::cerr << "v: ", *v ) << std::endl;
pass = false;
break;
}
CRT.getpolydom().random(generator, res, Degree((int64_t)Primes.size()-1));
CRT.RingToRns( Verifs, res );
Poly nres;
tim.clear(); tim.start();
CRT.RnsToRing( nres, Verifs );
tim.stop();
if (! CRT.getpolydom().areEqual(res,nres) ) {
CRT.getpolydom().write(std::cerr << "incoherency within ") << std::endl;
CRT.getpolydom().write(std::cerr << "r: ", res ) << std::endl;
CRT.getpolydom().write(std::cerr << "n: ", nres ) << std::endl;
pass = false;
}
F.write( std::cerr << tim << " using ") << std::endl;
return pass;
}
int main(int argc, char ** argv) {
// argv[1] : number of primes
// argv[2] : 2^{32-j} is size of primes
// argv[3] : seed for generator
GivRandom seedor( argc>3 ? (uint64_t)atoi(argv[3]): (uint64_t)BaseTimer::seed() );
uint64_t seed = seedor.seed();
std::cerr << "seed: " << seed << std::endl;
return
tmain<Field1>(argc, argv, *( new GivRandom(seed)))
&& tmain<Field2>(argc, argv, *( new GivRandom(seed)))
&& tmain<Field3>(argc, argv, *( new GivRandom(seed)))
&& tmain<Field4>(argc, argv, *( new GivRandom(seed)))
&& tmain<Field5>(argc, argv, *( new GivRandom(seed)))
&& tmain<Field6>(argc, argv, *( new GivRandom(seed)))
&& tmain<Field7>(argc, argv, *( new GivRandom(seed)))
&& tmain<Field8>(argc, argv, *( new GivRandom(seed)))
&& tmain<Field9>(argc, argv, *( new GivRandom(seed)))
&& tmainext<Field10>(argc, argv, *( new GivRandom(seed)))
;
}