test-ringarith.C
// Copyright(c)'1994-2009,2010 by The Givaro group
// This file is part of Givaro.
// BB: adapted,enhanced from examples/FiniteFields/ff-arith.C
// 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/modular.h>
//#include <givaro/modular-balanced.h>
//#include <givaro/montgomery.h>
#include <givaro/givpoly1.h>
#include <givaro/givinteger.h>
#include <givaro/zring.h>
#include <givaro/gfq.h>
#include <recint/recint.h>
//#include <givaro/modular-extended.h>
//#include <givaro/modular-general.h>
//#include <givaro/modular-integral.h>
//#include <givaro/modular-floating.h>
//#include <givaro/modular-integer.h>
//#include <givaro/modular-ruint.h>
//#include <givaro/modular-log16.h>
//#include <givaro/modular-inttype.h>
#include <givaro/modular.h>
using namespace Givaro;
#define TESTE_EG( a, b ) \
if (!F.areEqual((a),(b))) { \
F.write( F.write(std::cout,a) << "!=",b) \
<< " failed (at line " << __LINE__ << ")" << std::endl; \
return -1; \
}
#define TESTE_T( b ) \
if (!b) { \
F.write(std::cout) \
<< " failed (at line " << __LINE__ << ")" << std::endl; \
return -1; \
}
#define JETESTE( a, s ) \
if (TestRing( (a), (s)) ) { \
std::cout << #a << " failed !" << std::endl; \
return -1; \
}
#define JEPOLTESTE( a, s ) \
if (TestPolRing( (a), (s) ) ) { \
std::cout << #a << " failed !" << std::endl; \
return -1; \
}
#define JEONETESTE( F, x, y ) \
if (TestOneRing(F,x,y)) { \
std::cout << #x << " " << #y << " failed !" << std::endl; \
return -1; \
}
template<class Ring>
int TestOneRing(const Ring& F, const typename Ring::Element& x, const typename Ring::Element& y)
{
#ifdef GIVARO_DEBUG
std::cerr << "Testing " ;
F.write(std::cerr) << " : " << std::endl;
#endif
typename Ring::Element a, b, c, d,a_,b_,c_,d_;
typename Ring::Element e,e_;
F.init(a, 0U);
TESTE_EG(a, F.zero);
TESTE_T(F.isZero(a));
F.init(a, 1U);
// F.write(std::cerr) << std::endl;
// F.write(std::cerr << "a: ", a) << std::endl;
// F.write(std::cerr << "1: ", F.one) << std::endl;
TESTE_EG(a, F.one);
TESTE_T(F.isOne(a));
TESTE_T(F.isUnit(a));
F.init(a, -1);
TESTE_EG(a, F.mOne);
F.assign(a, x);
F.assign(b, y);
F.init(c); // empty constructor
F.init(d); // empty constructor
F.add(c, a, b); // c = a+b
F.assign(c_,c); // c_ <- c
TESTE_EG(c,c_);
F.subin(c_,a);
//F.write(std::cerr) << std::endl;
//F.write(std::cerr << "a: ", a) << std::endl;
//F.write(std::cerr << "b: ", b) << std::endl;
//F.write(std::cerr << "c: ", c) << std::endl;
//F.write(std::cerr << "c_: ", c_) << std::endl;
TESTE_EG(b,c_);
F.axpy(d, a, b, c); // d = a*b + c;
F.init(d_);
F.axmy(d_,a,b,c); // d_ = a*b - c
F.addin(d_,c);
F.subin(d,c);
//F.write(std::cerr) << std::endl;
//F.write(std::cerr << "a: ", a) << std::endl;
//F.write(std::cerr << "b: ", b) << std::endl;
//F.write(std::cerr << "c: ", c) << std::endl;
//F.write(std::cerr << "d: ", d) << std::endl;
//F.write(std::cerr << "d_: ", d_) << std::endl;
TESTE_EG(d_,d);
F.axpy(d, a, b, c); // d = a*b + c;
F.init(d_);
F.assign(d_,c);
F.axpyin(d_,a,b);
//F.write(std::cerr) << std::endl;
//F.write(std::cerr << "a: ", a) << std::endl;
//F.write(std::cerr << "b: ", b) << std::endl;
//F.write(std::cerr << "c: ", c) << std::endl;
//F.write(std::cerr << "d: ", d) << std::endl;
//F.write(std::cerr << "d_: ", d_) << std::endl;
TESTE_EG(d_, d);
F.sub(d,a,b); // d = a -b
F.add(c,a,b); // c = a+b
F.init(e);
F.init(e_);
F.mul(e,d,c); // e = d*c;
F.mul(a_,a,a); // a_ = a*a
F.mul(b_,b,b); // b_ = b*b
F.sub(e_,a_,b_); // e_ = a_ - b_
// F.write(std::cerr) << std::endl;
// F.write(std::cerr << "a: ", a) << std::endl;
// F.write(std::cerr << "a_: ", a_) << std::endl;
// F.write(std::cerr << "b: ", b) << std::endl;
// F.write(std::cerr << "b_: ", b_) << std::endl;
// F.write(std::cerr << "c: ", c) << std::endl;
// F.write(std::cerr << "d: ", d) << std::endl;
// F.write(std::cerr << "e: ", e) << std::endl;
// F.write(std::cerr << "e_: ", e_) << std::endl;
TESTE_EG(e,e_); // a^2 - b^2 = (a-b)(a+b)
F.maxpy(e, a, b, d); // e = d-a*b
F.assign(e_,d);
F.maxpyin(e_, a, b); // e = d - a*b;
TESTE_EG(e,e_);
F.axmy(e, a, b, d); // e = a*b -d;
F.assign(e_,d);
F.maxpyin(e_, a, b); // e = d - a*b;
F.negin(e_);
TESTE_EG(e,e_);
F.maxpy(e, a, b, d); // e = d-a*b;
F.assign(e_,d);
F.axmyin(e_, a, b); // e_ = a*b-e_ = a*b-d
F.negin(e_);
// F.write(std::cerr) << std::endl;
// F.write(std::cerr << "a: ", a) << std::endl;
// F.write(std::cerr << "b: ", b) << std::endl;
// F.write(std::cerr << "c: ", c) << std::endl;
// F.write(std::cerr << "d: ", d) << std::endl;
// F.write(std::cerr << "e: ", e) << std::endl;
// F.write(std::cerr << "e_: ", e_) << std::endl;
TESTE_EG(e,e_);
F.init(a, 3);
F.assign(b , y);
F.mul(c,a,b);
F.subin(c,b);
F.subin(c,b);
F.subin(c,b);
F.init(d, 0U);
TESTE_EG(c, d);
//F.write(std::cerr) << std::endl;
F.init(a, -3); // F.write(std::cerr << "a: ", a) << std::endl;
F.assign(b , y); // F.write(std::cerr << "b: ", b) << std::endl;
F.mul(c,a,b); // F.write(std::cerr << "c: ", c) << std::endl;
F.addin(c,b); // F.write(std::cerr << "c: ", c) << std::endl;
F.addin(c,b); // F.write(std::cerr << "c: ", c) << std::endl;
F.addin(c,b); // F.write(std::cerr << "c: ", c) << std::endl;
F.init(d, 0U);
//F.write(std::cerr << "d: ", d) << std::endl;
TESTE_EG(c, d);
#ifdef GIVARO_DEBUG
F.write(std::cerr );
std::cerr << " done." << std::endl;
#endif
return 0;
}
#ifndef NBITERD
#define NBITER 50
#endif
template<class Ring>
int TestRing(const Ring& F, const uint64_t seed)
{
typename Ring::Element x, y;
typename Ring::RandIter g(F, 0_ui64, seed);
F.init(x, 7U);
F.init(y, -29.0);
JEONETESTE(F,x,y);
F.init(x, Ring::maxCardinality()-1U);
F.init(y, Ring::maxCardinality()-1U);
JEONETESTE(F,x,y);
F.assign(x, F.maxElement());
F.assign(y, F.maxElement());
JEONETESTE(F,x,y);
F.assign(x, F.minElement());
F.assign(y, F.maxElement());
JEONETESTE(F,x,y);
F.assign(x, F.minElement());
F.assign(y, F.minElement());
JEONETESTE(F,x,y);
for (size_t i = 0; i< NBITER; ++i) {
g.random(x); g.random(y);
JEONETESTE(F,x,y);
}
return 0;
}
#ifndef DEGMAX
#define DEGMAX 75
#endif
#ifndef NBITERD
#define NBITERD 10
#endif
template<class Ring>
int TestPolRing(const Ring& F, const uint64_t seed)
{
GivRandom generator(seed);
srand48((int64_t)seed);
for (size_t i = 0; i < NBITERD; ++i) {
int64_t d1 = int64_t (lrand48() % DEGMAX);
int64_t d2 = int64_t (lrand48() % DEGMAX);
typename Ring::Element x, d, z, o;
do { F.random(generator, x, Degree(d1)); } while(F.isZero(x));
do { F.random(generator, d, Degree(d2)); } while(F.isZero(d));
JEONETESTE(F,x,d);
do { F.random(generator, z, Degree(0)); } while(F.isZero(z));
JEONETESTE(F,x,z);
JEONETESTE(F,z,x);
do { F.random(generator, o, Degree(1)); } while(F.isZero(o));
JEONETESTE(F,d,o);
JEONETESTE(F,o,d);
}
return 0;
}
template<class Ring>
int TestInv(const Ring& F, const uint64_t seed)
{
typename Ring::Element a, inva, invinva;
typename Ring::RandIter g(F, 0_ui64, seed);
F.init(a);
F.init(inva);
F.init(invinva);
for (int i=0; i < 100; i++)
{
g.random(a);
while (F.isZero(a)) g.random(a);
F.inv(inva, a);
F.inv(invinva, inva);
// F.write(std::cerr);
// F.write(std::cerr << " => a: ", a);
// F.write(std::cerr << "; inva: ", inva);
// F.write(std::cerr << "; invinva: ", invinva) << std::endl;
TESTE_EG(a, invinva);
}
return 0;
}
int main(int argc, char ** argv)
{
auto seed = static_cast<uint64_t>(argc>1?atoi(argv[1]):BaseTimer::seed());
#ifdef GIVARO_DEBUG
std::cerr << "seed: " << seed << std::endl;
#endif
Integer::seeding(seed);
RecInt::srand(seed);
using ModularCUS = Modular<int8_t, uint16_t>;
using ModularSUZ = Modular<int16_t, uint32_t>;
using ModularZULL = Modular<int32_t, uint64_t>;
using ModularUCUS = Modular<uint8_t, uint16_t>;
using ModularUSUZ = Modular<uint16_t, uint32_t>;
using ModularUZULL = Modular<uint32_t, uint64_t>;
using ModularFD = Modular<float, double>;
#ifdef __GIVARO_HAVE_INT128
using ModularLLULLL = Modular<int64_t, uint128_t>;
using ModularULLULLL = Modular<uint64_t, uint128_t>;
#endif
#define TEST_SPECIFIC(Ring, Name, Modulus...) \
std::cout << "TEST_SPECIFIC: " << #Name << std::endl; \
Ring Name(Modulus); \
JETESTE(Name, seed);
//Modular<float> GF5(5);
//float elt, inv;
//GF5.init(elt, 2);
//GF5.inv(inv, elt);
//GF5.write(std::cout, inv) << " ####\n";
//-------------//
//----- 4 -----//
TEST_SPECIFIC(Modular<int8_t>, C4, 4);
TEST_SPECIFIC(Modular<int16_t>, S4, 4);
TEST_SPECIFIC(Modular<int32_t>, Z4, 4);
TEST_SPECIFIC(Modular<int64_t>, LL4, 4);
TEST_SPECIFIC(Modular<uint8_t>, UC4, 4);
TEST_SPECIFIC(Modular<uint16_t>, US4, 4);
TEST_SPECIFIC(Modular<uint32_t>, UZ4, 4);
TEST_SPECIFIC(Modular<uint64_t>, ULL4, 4);
TEST_SPECIFIC(ModularCUS, CUS4, 4);
TEST_SPECIFIC(ModularSUZ, SUZ4, 4);
TEST_SPECIFIC(ModularZULL, ZULL4, 4);
TEST_SPECIFIC(ModularUCUS, UCUS4, 4);
TEST_SPECIFIC(ModularUSUZ, USUZ4, 4);
TEST_SPECIFIC(ModularUZULL, UZULL4, 4);
TEST_SPECIFIC(ModularFD, FD4, 4);
#ifdef __GIVARO_HAVE_INT128
TEST_SPECIFIC(ModularLLULLL, LLULLL4, 4);
TEST_SPECIFIC(ModularULLULLL, ULLULLL4, 4);
#endif
TEST_SPECIFIC(Modular<float>, F4, 4);
TEST_SPECIFIC(Modular<double>, D4, 4);
TEST_SPECIFIC(Modular<Integer>, I4, 4);
TEST_SPECIFIC(Modular<RecInt::ruint128>, RU4, 4);
TEST_SPECIFIC(Modular<RecInt::rint128>, R4, 4);
//TEST_SPECIFIC(ZRing<Integer>, ZR4, 4);
TEST_SPECIFIC(Modular<Log16>, L5, 5);
TEST_SPECIFIC(Modular<Log16>, L7, 7);
//--------------//
//----- 75 -----//
TEST_SPECIFIC(Modular<Log16>, L79, 79);
TEST_SPECIFIC(Modular<int8_t>, C75, 13);
TEST_SPECIFIC(Modular<int16_t>, S75, 75);
TEST_SPECIFIC(Modular<int32_t>, Z75, 75);
TEST_SPECIFIC(Modular<int64_t>, LL75, 75);
TEST_SPECIFIC(Modular<uint8_t>, UC75, 13);
TEST_SPECIFIC(Modular<uint16_t>, US75, 75);
TEST_SPECIFIC(Modular<uint32_t>, UZ75, 75);
TEST_SPECIFIC(Modular<uint64_t>, ULL75, 75);
TEST_SPECIFIC(ModularCUS, CUS75, 75);
TEST_SPECIFIC(ModularSUZ, SUZ75, 75);
TEST_SPECIFIC(ModularZULL, ZULL75, 75);
TEST_SPECIFIC(ModularUCUS, UCUS75, 75);
TEST_SPECIFIC(ModularUSUZ, USUZ75, 75);
TEST_SPECIFIC(ModularUZULL, UZULL75, 75);
TEST_SPECIFIC(ModularFD, FD75, 75);
#ifdef __GIVARO_HAVE_INT128
TEST_SPECIFIC(ModularLLULLL, LLULLL75, 75);
TEST_SPECIFIC(ModularULLULLL, ULLULLL75, 75);
#endif
TEST_SPECIFIC(Modular<float>, F75, 75);
TEST_SPECIFIC(Modular<double>, D75, 75);
TEST_SPECIFIC(Modular<Integer>, I75, 75);
TEST_SPECIFIC(Modular<RecInt::ruint128>, RU75, 75);
TEST_SPECIFIC(Modular<RecInt::rint128>, R75, 75);
//TEST_SPECIFIC(ModularBalanced<int32_t>, BZ75, 75);
//TEST_SPECIFIC(ModularBalanced<int64_t>, BLL75, 75);
//TEST_SPECIFIC(ModularBalanced<float>, BF75, 75);
//TEST_SPECIFIC(ModularBalanced<double>, BD75, 75);
//TEST_SPECIFIC(Montgomery<int32_t>, MZ75, 75);
//TEST_SPECIFIC(Montgomery<RecInt::ruint128>, MRU75, 75);
//TEST_SPECIFIC(ModularExtended<float>, MEF75, 75);
//TEST_SPECIFIC(ModularExtended<double>, MED75, 75);
#define TEST_POLYNOMIAL(BaseRing, Name, BaseRingName) \
Poly1Dom<BaseRing, Dense> Name(BaseRingName, "X"); \
JEPOLTESTE(Name, seed);
TEST_POLYNOMIAL(Modular<Log16>, PL79, L79);
TEST_POLYNOMIAL(Modular<int8_t>, PC75, C75);
TEST_POLYNOMIAL(Modular<int16_t>, PS75, S75);
TEST_POLYNOMIAL(Modular<int32_t>, PZ75, Z75);
TEST_POLYNOMIAL(Modular<int64_t>, PLL75, LL75);
TEST_POLYNOMIAL(Modular<uint8_t>, PUC75, UC75);
TEST_POLYNOMIAL(Modular<uint16_t>, PUS75, US75);
TEST_POLYNOMIAL(Modular<uint32_t>, PUZ75, UZ75);
TEST_POLYNOMIAL(Modular<uint64_t>, PULL75, ULL75);
TEST_POLYNOMIAL(ModularCUS, PCUS75, CUS75);
TEST_POLYNOMIAL(ModularSUZ, PSUZ75, SUZ75);
TEST_POLYNOMIAL(ModularZULL, PZULL75, ZULL75);
TEST_POLYNOMIAL(ModularUCUS, PUCUS75, UCUS75);
TEST_POLYNOMIAL(ModularUSUZ, PUSUZ75, USUZ75);
TEST_POLYNOMIAL(ModularUZULL, PUZULL75, UZULL75);
TEST_POLYNOMIAL(ModularFD, PFD75, FD75);
TEST_POLYNOMIAL(Modular<float>, PF75, F75);
TEST_POLYNOMIAL(Modular<double>, PD75, D75);
TEST_POLYNOMIAL(Modular<Integer>, PI75, I75);
TEST_POLYNOMIAL(Modular<RecInt::ruint128>, PRU75, RU75);
TEST_POLYNOMIAL(Modular<RecInt::rint128>, PR75, R75);
//TEST_POLYNOMIAL(ModularBalanced<int32_t>, MBZ75, BZ75);
//TEST_POLYNOMIAL(ModularBalanced<int64_t>, MBLL75, BLL75);
//TEST_POLYNOMIAL(ModularBalanced<float>, MBF75, BF75);
//TEST_POLYNOMIAL(ModularBalanced<double>, MBD75, BD75);
//TEST_POLYNOMIAL(Montgomery<int32_t>, PMZ75, MZ75);
// @bug Convert to double inside? //TEST_POLYNOMIAL(Montgomery<RecInt::ruint128>, PMRU75, MRU75);
//TEST_POLYNOMIAL(decltype(PI75), PPI75, PI75);
//TEST_POLYNOMIAL(decltype(PPI75), PPPI75, PPI75);
//--------------------------------//
//----- Modulo maximal prime -----//
#define TEST_LAST(Field, Name) \
std::cout << "TEST_LAST: " << #Name; \
Field Name(Field::maxCardinality()); \
Name.write(std::cout << " (", Field::maxCardinality()) << ")"<< std::endl; \
JETESTE(Name, seed);
TEST_LAST(Modular<Log16>, Lmax);
TEST_LAST(Modular<int8_t>, Cmax);
TEST_LAST(Modular<int16_t>, Smax);
TEST_LAST(Modular<int32_t>, Zmax);
TEST_LAST(Modular<int64_t>, LLmax);
TEST_LAST(Modular<uint8_t>, UCmax);
TEST_LAST(Modular<uint16_t>, USmax);
TEST_LAST(Modular<uint32_t>, UZmax);
TEST_LAST(Modular<uint64_t>, ULLmax);
TEST_LAST(ModularCUS, CUSmax);
TEST_LAST(ModularSUZ, SUZmax);
TEST_LAST(ModularZULL, ZULLmax);
TEST_LAST(ModularUCUS, UCUSmax);
TEST_LAST(ModularUSUZ, USUZmax);
TEST_LAST(ModularUZULL, UZULLmax);
#ifdef __GIVARO_HAVE_INT128
TEST_LAST(ModularLLULLL, LLULLLmax);
TEST_LAST(ModularULLULLL, ULLULLLmax);
#endif
TEST_LAST(Modular<float>, Fmax);
TEST_LAST(Modular<double>, Dmax);
TEST_LAST(ModularFD, FDmax);
TEST_LAST(Modular<RecInt::ruint128>, RUmax);
typedef Modular<RecInt::ruint128,RecInt::ruint256> MyMod;
TEST_LAST(MyMod, RUExtmax);
//TEST_LAST(ModularBalanced<int32_t>, BZmax);
//TEST_LAST(ModularBalanced<int64_t>, BLLmax);
//TEST_LAST(ModularBalanced<float>, BFmax);
//TEST_LAST(ModularBalanced<double>, BDmax);
//TEST_LAST(Montgomery<int32_t>, MZmax);
//TEST_LAST(Montgomery<RecInt::ruint128>, MRUmax);
//TEST_LAST(ModularExtended<float>, MEFmax);
//TEST_LAST(ModularExtended<double>, MEDmax);
// -----------------------
// Test inversions
// -----------------------
#define TEST_INV(Field, Name, Prime) \
std::cout << "TEST_INV: " << #Name; \
Field Name(Prime); \
std::cout << " (" << (Integer)Name.cardinality() << ',' << (Integer)Name.maxCardinality() << ')'<< std::endl; \
if (TestInv( (Name), (seed))) { \
std::cout << #Name << " failed !" << std::endl; \
return -1; \
}
//TEST_INV(Modular<int8_t>, C17, 17); => GF(17) does not fit in Modular<int8_t>
TEST_INV(Modular<int8_t>, C13, 13);
TEST_INV(Modular<int16_t>, S17, 17);
TEST_INV(Modular<int32_t>, Z17, 17);
TEST_INV(Modular<int64_t>, LL17, 17);
TEST_INV(Modular<uint8_t>, UC13, 13);
TEST_INV(Modular<uint16_t>, US17, 17);
TEST_INV(Modular<uint32_t>, UZ17, 17);
TEST_INV(Modular<uint64_t>, ULL17, 17);
TEST_INV(ModularCUS, CUS17, 17);
TEST_INV(ModularSUZ, SUZ17, 17);
TEST_INV(ModularZULL, ZULL17, 17);
TEST_INV(ModularUCUS, UCUS17, 17);
TEST_INV(ModularUSUZ, USUZ17, 17);
TEST_INV(ModularUZULL, UZULL17, 17);
TEST_INV(ModularFD, FD17, 17);
#ifdef __GIVARO_HAVE_INT128
TEST_INV(ModularLLULLL, LLULLL17, 17);
TEST_INV(ModularULLULLL, ULLULLL17, 17);
#endif
TEST_INV(Modular<float>, F17, 17);
TEST_INV(Modular<double>, D17, 17);
TEST_INV(Modular<Integer>, I17, 17);
TEST_INV(Modular<RecInt::ruint128>, RU17, 17);
TEST_INV(Modular<RecInt::rint128>, R17, 17);
//TEST_INV(ZRing<Integer>, ZR17, 17);
TEST_INV(Modular<Log16>, L17, 17);
return 0;
}
/* -*- 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
