examples/Polynomial/PolynomialCRT.C

NO DOC.

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/qfield.h>
using namespace Givaro;
typedef GFqDom<int64_t> Field1;
typedef Modular<int16_t> Field2;
typedef Modular<Log16> Field3;
typedef Modular<int32_t> Field4;
typedef Modular<int64_t> Field5;
typedef Modular<uint32_t> Field6;
typedef Montgomery<int32_t> Field7;
typedef QField<Rational> Field8;
typedef Modular<Integer> Field9;
typedef Extension<> Field10;
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;
IntPrimeDom ID;
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;
IntPrimeDom ID;
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;
Integer::seeding(seed);
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)))
;
}
/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
Givaro::Degree
Degree type for polynomials.
Definition givdegree.h:28
Givaro::Extension
Extension.
Definition extension.h:79
Givaro::GFqDom
class GFqDom
Definition gfq.h:44
Givaro::GivRandom
GivRandom.
Definition givrandom.h:37
Givaro::IntPrimeDom
Primality tests.
Definition givintprime.h:67
Givaro::Integer
This is the Integer class.
Definition gmp++_int.h:160
Givaro::Integer::seeding
static void seeding(uint64_t s)
Random numbers (no doc).
Definition gmp++_int_rand.inl:44
Givaro::Modular
Forward declaration for Givaro::Modular.
Definition modular-inttype.h:38
Givaro::Poly1CRT
Poly1 CRT.
Definition givpoly1crt.h:26
Givaro::Timer
Timer.
Definition givtimer.h:129
Givaro::Timer::start
void start()
Start timer.
Definition givtimer.C:164
Givaro::Timer::stop
void stop()
Stop timer.
Definition givtimer.C:172
Givaro::Timer::clear
void clear()
Clear timer.
Definition givtimer.C:155
Givaro
Namespace in which the whole Givaro library resides.
Definition all_field.C:23