[givaro] 01/07: Imported Upstream version 4.0.2
Doug Torrance
dtorrance-guest at moszumanska.debian.org
Wed Aug 10 05:09:37 UTC 2016
This is an automated email from the git hooks/post-receive script.
dtorrance-guest pushed a commit to branch master
in repository givaro.
commit 6ffd038226519a8675bef35561e22c8ec1978b08
Author: Doug Torrance <dtorrance at piedmont.edu>
Date: Tue Aug 9 19:10:50 2016 -0400
Imported Upstream version 4.0.2
---
.gitignore | 3 +-
.travis.yml | 52 +++
AUTHORS | 11 +-
ChangeLog | 1 -
Makefile.am | 2 +-
README.md | 23 +-
autogen.sh | 6 +-
benchmarks/Makefile.am | 2 +-
benchmarks/Makefile.tests | 78 ++++
benchmarks/benchmark-recint_exp.C | 42 ++-
benchmarks/benchmark-recint_inadd.C | 65 ++++
benchmarks/benchmark-recint_inmul.C | 65 ++++
benchmarks/benchmark-recint_inv_arazi.C | 25 +-
benchmarks/benchmark-recint_mul.C | 26 +-
benchmarks/benchmark-recint_mulm.C | 74 ++++
benchmarks/benchmark-recint_tmul.C | 63 ++++
benchmarks/generic.gnuplot | 24 ++
benchmarks/perfpublisher.sh | 26 +-
configure.ac | 34 +-
examples/Integer/primitiveroot.C | 54 +--
examples/Polynomial/Makefile.am | 3 +-
examples/Polynomial/isgenerator.C | 53 +++
examples/RecInt/Makefile.am | 34 +-
examples/RecInt/extended-int-types.C | 62 ++++
examples/RecInt/iterator.C | 29 --
.../rint.h => examples/RecInt/recint-iterator.C | 70 ++--
givaro.pc.in | 4 +-
macros/CodeChunk/Makefile.am | 18 +
macros/CodeChunk/avx.C | 11 +
macros/CodeChunk/sse.C | 12 +
macros/ax_check_x86_features.m4 | 77 ++++
macros/ax_gcc_x86_cpu_supports.m4 | 104 ++++++
macros/simd-check.m4 | 126 +++++++
src/Makefile.am | 2 +-
src/kernel/field/extension.h | 2 +-
src/kernel/field/gf2.h | 4 +-
src/kernel/field/gf2.inl | 0
src/kernel/field/gfq.h | 10 +-
src/kernel/field/gfq.inl | 71 ++++
src/kernel/recint/recdefine.h | 1 +
src/kernel/recint/rint.h | 2 +
src/kernel/recint/rmgmodule.h | 10 +-
src/kernel/recint/rmint.h | 2 +
src/kernel/recint/rmintmg.h | 2 +
src/kernel/recint/ruconvert.h | 51 ++-
src/kernel/recint/ruint.h | 2 +
src/kernel/recint/rumul.h | 6 +-
src/kernel/recint/ruruint.h | 11 +-
src/kernel/ring/Makefile.am | 5 +-
src/kernel/ring/modular-double.h | 5 +-
src/kernel/ring/modular-extended.h | 305 +++++++++++++++
src/kernel/ring/modular-extended.inl | 285 ++++++++++++++
src/kernel/ring/modular-float.h | 8 +-
src/kernel/ring/modular-int16.h | 412 +++++++++++----------
src/kernel/ring/modular-int32.h | 330 +++++++++--------
src/kernel/ring/modular-int32.inl | 11 +-
src/kernel/ring/modular-int64.h | 391 ++++++++++---------
src/kernel/ring/modular-int64.inl | 30 +-
src/kernel/ring/modular-int8.h | 407 ++++++++++----------
src/kernel/ring/modular-integer.h | 6 +-
src/kernel/ring/modular-inttype.h | 8 +-
src/kernel/ring/modular-mulprecomp.inl | 95 +++++
src/kernel/ring/modular-ruint.h | 4 +-
src/kernel/ring/modular-uint16.h | 340 +++++++++--------
src/kernel/ring/modular-uint32.h | 335 +++++++++--------
src/kernel/ring/modular-uint64.h | 327 ++++++++--------
src/kernel/ring/modular-uint64.inl | 31 +-
src/kernel/ring/modular-uint8.h | 327 ++++++++--------
src/kernel/ring/modular-uint8.inl | 5 +-
src/kernel/ring/modular.h | 4 +-
src/kernel/ring/montgomery-int32.h | 3 +-
src/kernel/ring/montgomery-ruint.h | 8 +-
src/kernel/system/givconfig.h | 12 +-
tests/Makefile.am | 4 +-
tests/jenkins-maker.sh | 85 +++++
tests/perfpublisher.sh | 25 +-
tests/test-ffarith.C | 27 +-
tests/test-modularmulprecomp.C | 211 +++++++++++
tests/test-ringarith.C | 94 +++--
tests/test-ruint_operators.C | 6 +-
80 files changed, 3993 insertions(+), 1608 deletions(-)
diff --git a/.gitignore b/.gitignore
index cd46ce1..071fcc5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -68,11 +68,12 @@ examples/Polynomial/pol_factor
examples/Polynomial/trunc_arith
examples/Rational/iratrecon
examples/Rational/polydouble
-examples/RecInt/iterator
+examples/RecInt/recint-iterator
examples/RecInt/rsa
givaro-config
givaro-config.h
givaro-makefile
+givaro.pc
libtool
macros/libtool.m4
macros/ltoptions.m4
diff --git a/.travis.yml b/.travis.yml
new file mode 100644
index 0000000..9a84868
--- /dev/null
+++ b/.travis.yml
@@ -0,0 +1,52 @@
+sudo: required
+dist: trusty
+language: cpp
+
+matrix:
+ include:
+ - compiler: gcc
+ env: COMPILER=g++-4.9
+ addons:
+ apt:
+ sources:
+ - ubuntu-toolchain-r-test
+ packages:
+ - g++-4.9
+ - libgmp-dev
+ - libgmpxx4ldbl
+
+ - compiler: gcc
+ env: COMPILER=g++-5
+ addons:
+ apt:
+ sources:
+ - ubuntu-toolchain-r-test
+ packages:
+ - g++-5
+ - libgmp-dev
+ - libgmpxx4ldbl
+
+ - compiler: clang
+ env: COMPILER=clang++
+ addons:
+ apt:
+ packages:
+ - libgmp-dev
+ - libgmpxx4ldbl
+
+install:
+ export CXX="$COMPILER"
+
+before_script:
+ - ./autogen.sh
+
+script:
+ - make
+ - make check
+ - make benchmarks
+ - make examples
+ - make dist
+
+notifications:
+ on_success: change
+ on_failure: always
\ No newline at end of file
diff --git a/AUTHORS b/AUTHORS
index afd5081..d97e99f 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -1,9 +1,12 @@
Authors of the Givaro library
-Thierry Gautier
-Jean-Louis Roch
-Gilles Villard
+Brice Boyer
+Alexis Breust
Jean-Guillaume Dumas
+Thierry Gautier
Pascal Giorgi
+Romain Lebreton
Clément Pernet
-Brice Boyer
+Jean-Louis Roch
+Bastien Vialla
+Gilles Villard
diff --git a/ChangeLog b/ChangeLog
index 776de0f..3ade94e 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,4 +1,3 @@
2003-12-11 Pascal Giorgi <pascal.giorgi at ens-lyon.fr>
* /src/kernel/zpz/givzpz.h add struct Std64 with some definition and redefine in the ggod way int64 & uint64
-Hi
diff --git a/Makefile.am b/Makefile.am
index abc9f3d..2b1eb10 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -66,7 +66,7 @@ docs_dev:docs/givaro-dev-html/index.html
docs/givaro-dev-html/index.html:
(cd docs; ${MAKE} docs_dev)
-VERSION=4.0.1
+VERSION=4.0.2
git:
git commit -a; git pull; git push
diff --git a/README.md b/README.md
index 4aa0492..1229eb6 100644
--- a/README.md
+++ b/README.md
@@ -1,10 +1,12 @@
Givaro
======
+[](https://ci.inria.fr/linbox/job/Givaro/)
+
Download and install
--------------------
-For lastest releases, please check out [this website](https://forge.imag.fr/frs/?group_id=187).
+For lastest releases, please check out [this website](http://github.com/linbox-team/givaro); older releases can be found on [that website](https://forge.imag.fr/frs/?group_id=187).
Then, you can install doing:
```
@@ -23,10 +25,23 @@ Then, you can install doing:
Compile your own files
----------------------
-An optional compilation help file is provided: just add the following line to your Makefile. Then a simple call will compile your C and C++ files.
+Givaro uses pkgconfig to expose the compilation flags it requires.
+You will get the compilation flags by calling
```
-include ##GIVAROROOT##/bin/givaro-makefile
+pkg-config --cflags givaro
+```
+and the linking flags by calling
+```
+pkg-config --libs givaro
```
-However, if you want to do it without this tool, you should add `-I##GIVAROROOT##/include` to the CXX compilation flags, and `-L##GIVAROROOT##/lib -lgivaro` to the LD link flags, along those for GMP.
\ No newline at end of file
+If you have installed givaro in a non-standard directory (such as `/usr/local`), make sure to have added the path where to find givaro's .pc file to the PKG_CONFIG_PATH environment variable.
+```
+PKG_CONFIG_PATH=${PKG_CONFIG_PATH}:<path to your givaro install>/lib/pkgconfig
+```
+
+An alternative option is to just add the following line to your Makefile. Then a simple call will compile your C and C++ files.
+```
+include ##GIVAROROOT##/bin/givaro-makefile
+```
diff --git a/autogen.sh b/autogen.sh
index 1f1a9e3..82e031a 100755
--- a/autogen.sh
+++ b/autogen.sh
@@ -52,8 +52,10 @@ LIBTOOLIZE=libtoolize
(uname -a|grep -v Darwin) < /dev/null > /dev/null 2>&1 ||
{
echo "....Adding fix for OSX"
-LIBTOOL=glibtool
-LIBTOOLIZE=glibtoolize
+if command -v "glibtoolize" >/dev/null; then
+ LIBTOOL=glibtool
+ LIBTOOLIZE=glibtoolize
+fi
}
diff --git a/benchmarks/Makefile.am b/benchmarks/Makefile.am
index 6a7eb7f..8e44711 100644
--- a/benchmarks/Makefile.am
+++ b/benchmarks/Makefile.am
@@ -15,7 +15,7 @@ OPTFLAGS=
OPTLINKS=
AM_CXXFLAGS = @DEFAULT_CFLAGS@
-AM_CPPFLAGS += $(OPTFLAGS) $(GMP_CFLAGS) -I$(top_srcdir)/src/kernel/system -I$(top_srcdir)/src/kernel/recint
+AM_CPPFLAGS += $(OPTFLAGS) $(GMP_CFLAGS) -I$(top_srcdir)/src/kernel/system -I$(top_srcdir)/src/kernel/recint -I$(top_srcdir)/src/kernel/integer -I$(top_srcdir)/src/kernel/gmp++ -I$(top_srcdir)/src/kernel/ring
LDADD = $(OPTLINKS) -L${top_srcdir}/src -L${top_srcdir}/src/.libs -lgivaro $(GMP_LIBS) $(LDFLAGS)
AM_LDFLAGS=-static
diff --git a/benchmarks/Makefile.tests b/benchmarks/Makefile.tests
new file mode 100644
index 0000000..c2335ee
--- /dev/null
+++ b/benchmarks/Makefile.tests
@@ -0,0 +1,78 @@
+ifneq ($(origin RMOD), undefined)
+MODME=${RMOD}
+else
+MODME=mul mulm # exp
+endif
+
+ifneq ($(origin RDIR), undefined)
+DIRME=${RDIR}
+else
+DIRME=inmul tmul inadd
+endif
+
+ifneq ($(origin RTES), undefined)
+ORTES=${RTES}
+else
+ORTES=${MODME} ${DIRME}
+endif
+
+THREADS=4
+
+ifneq ($(origin LOOPS), undefined)
+ OLOOPS=${LOOPS}
+else
+ OLOOPS=22
+endif
+
+MGTES=$(filter-out ${DIRME},${ORTES})
+DETES=$(filter ${DIRME},${ORTES})
+
+SUFFIXES=mul mulm inmul tmul inadd exp inv_arazi
+MEANINGS=Modular_multiplication In_place_modular_multiplication In_place_multiplication Multiplication In_place_addition Modular_exponentiation Base_inverse
+
+EXEC=${ORTES:%=benchmark-recint_%}
+WSRC=${EXEC:%=-W %.C}
+OUTP=output.rint
+MODEL=$(shell cat /proc/cpuinfo | grep "model name" | head -1|cut -d':' -f2| tr -s ' '|sed 's/^ //')
+
+
+SHELL := /bin/bash
+index = $(words $(shell a="$(2)";echo $${a/$(1)*/$(1)} ))
+swap = $(word $(call index,$(1),${SUFFIXES}),${MEANINGS})
+remun = $(shell a="$(1)";echo $${a}|sed 's/_/ /g')
+
+SIZES=6 7 8 9 10 11 12 13
+MONTG=MG_INACTIVE
+PREFI=
+
+all: outclean ${SIZES:%=all_%} outcleanmg ${SIZES:%=allmg_%} plot plotmg
+
+mkplot = $(foreach fil, $(1), sed 's/FUNCTION/${PREFI}${fil}/g;s/MODEL/${MODEL}/g;s/MEANING/$(call remun,$(call swap,${fil}) $(2))/' generic.gnuplot > ${PREFI}${fil}.gnuplot; gnuplot ${PREFI}${fil}.gnuplot;)
+
+mkruns = make -j ${THREADS} "OPTFLAGS=-DSTD_RECINT_SIZE=$(1) -DMG_DEFAULT=${MONTG}" $(addprefix ./benchmark-recint_,$(2)) ${WSRC}; $(foreach fil, $(2), ./benchmark-recint_${fil} `echo '10^((${OLOOPS}-$(1))/2)'|bc` >> ${OUTP}.${PREFI}${fil};)
+
+mkocl = $(foreach fil, $(1), - rm ${OUTP}.${PREFI}${fil};)
+
+all_%:
+ $(call mkruns,$*, ${ORTES})
+
+plot:
+ $(call mkplot, ${DETES})
+ $(call mkplot, ${MGTES}, (${MONTG}))
+
+outclean:
+ $(call mkocl, ${ORTES})
+
+allmg_%: MONTG=MG_ACTIVE
+allmg_%: PREFI=mg_
+allmg_%:
+ $(call mkruns,$*, ${MGTES})
+
+plotmg: MONTG=MG_ACTIVE
+plotmg: PREFI=mg_
+plotmg:
+ $(call mkplot, ${MGTES}, (${MONTG}))
+
+outcleanmg: PREFI=mg_
+outcleanmg:
+ $(call mkocl, ${MGTES})
\ No newline at end of file
diff --git a/benchmarks/benchmark-recint_exp.C b/benchmarks/benchmark-recint_exp.C
index 866687a..02bb326 100644
--- a/benchmarks/benchmark-recint_exp.C
+++ b/benchmarks/benchmark-recint_exp.C
@@ -4,44 +4,68 @@
#include <recint/recint.h>
#include <givaro/givtimer.h>
-#define STD_RECINT_SIZE 9
-#define LOOPS 500
+#if not defined(STD_RECINT_SIZE)
+#define STD_RECINT_SIZE 8
+#endif
+
+#if not defined(LOOPS)
+#define LOOPS 100000
+#endif
#define ALEA_MAX 64
#define ALEA_MASK 63
using namespace RecInt;
-int main(void)
+int main(int argc, char ** argv)
{
+ size_t nbloops = static_cast<size_t>((argc > 1)? atoi(argv[1]) : LOOPS);
+
rmint<STD_RECINT_SIZE> m[ALEA_MAX];
ruint<STD_RECINT_SIZE> u[ALEA_MAX], module;
- Givaro::Timer tim;
+ mpz_class b[ALEA_MAX], c[ALEA_MAX], gmod;
+ Givaro::Timer tim, gmp;
// For montgomery algorithm, the module must be odd
RecInt::srand(42);
rand(module);
if (module % 2 == 0) module++;
rmint<STD_RECINT_SIZE>::init_module(module);
+ ruint_to_mpz(gmod, module);
// Randomness
for (unsigned int i = 0; i < ALEA_MAX; i++) {
- rand(m[i]);
- rand(u[i]);
+ rand(m[i]); ruint_to_mpz(b[i],m[i].Value);
+ rand(u[i]); ruint_to_mpz(c[i],u[i]);
}
// Main loop
tim.clear(); tim.start();
- for (unsigned int l = 0; l < LOOPS; l++) {
+ for (unsigned int l = 0; l < nbloops; l++) {
exp(m[l & ALEA_MASK], m[(l+2) & ALEA_MASK], u[l & ALEA_MASK]);
exp(m[(l+2) & ALEA_MASK], m[(l+1) & ALEA_MASK], u[l & ALEA_MASK]);
}
tim.stop();
+
+ gmp.clear(); gmp.start();
+ for (unsigned int l = 0; l < nbloops; l++) {
+ mpz_powm(b[l & ALEA_MASK].get_mpz_t(),b[(l+2) & ALEA_MASK].get_mpz_t(),c[l & ALEA_MASK].get_mpz_t(), gmod.get_mpz_t());
+ mpz_powm(b[(l+2) & ALEA_MASK].get_mpz_t(),b[(l+1) & ALEA_MASK].get_mpz_t(),c[l & ALEA_MASK].get_mpz_t(), gmod.get_mpz_t());
+ }
+ gmp.stop();
+
+
+ rand(module);
// -----------
// Standard output for benchmark - Alexis Breust 2014/12/11
- std::cout << "Time: " << tim.usertime()
- << " Gflops: " << "Irrelevant" << std::endl;
+ std::cout
+ << "Time: " << tim.usertime()
+ << " Mflops: " << std::scientific << (double(2*nbloops))/tim.usertime()/1000.0/1000.0 << ' ' << (double(2*nbloops))/gmp.usertime()/1000.0/1000.0
+ << " SIZE: " << STD_RECINT_SIZE
+ << " GMP time: " << gmp.usertime()
+ << ' ' << m[(int)(module )& ALEA_MASK] << ' ' << b[(int)(module)& ALEA_MASK] << std::endl ;
+
return 0;
}
diff --git a/benchmarks/benchmark-recint_inadd.C b/benchmarks/benchmark-recint_inadd.C
new file mode 100644
index 0000000..babb974
--- /dev/null
+++ b/benchmarks/benchmark-recint_inadd.C
@@ -0,0 +1,65 @@
+#include <iostream>
+#include <cstdlib>
+
+#include <recint/recint.h>
+#include <givaro/givtimer.h>
+#include <givaro/givinteger.h>
+
+#if not defined(STD_RECINT_SIZE)
+#define STD_RECINT_SIZE 8
+#endif
+
+#if not defined(LOOPS)
+#define LOOPS 1000000
+#endif
+
+#define ALEA_MAX 64
+#define ALEA_MASK 63
+
+using namespace RecInt;
+
+int main(int argc, char ** argv)
+{
+ size_t nbloops = static_cast<size_t>((argc > 1)? atoi(argv[1]) : LOOPS);
+// std::cerr << "nbloops: " << nbloops << std::endl;
+
+ ruint<STD_RECINT_SIZE> a[ALEA_MAX], d[ALEA_MAX];
+ mpz_class b[ALEA_MAX], c[ALEA_MAX];
+ Givaro::Timer tim,gmp;
+
+
+ // Randomness
+ for (unsigned int i = 0; i < ALEA_MAX; i++) {
+ rand(a[i]);
+ ruint_to_mpz(b[i],a[i]);
+ }
+
+ // Main loop
+ tim.clear(); tim.start();
+ for (unsigned int l = 0; l < nbloops; l++) {
+ d[l & ALEA_MASK] = a[l & ALEA_MASK];
+ d[l & ALEA_MASK] += a[l & ALEA_MASK];
+ }
+ tim.stop();
+
+ // Main loop
+ gmp.clear(); gmp.start();
+ for (unsigned int l = 0; l < nbloops; l++) {
+ c[l & ALEA_MASK] = b[l & ALEA_MASK];
+ c[l & ALEA_MASK] += b[l & ALEA_MASK];
+ }
+ gmp.stop();
+
+ ruint<STD_RECINT_SIZE> module; rand(module);
+
+ // -----------
+ // Standard output for benchmark - Alexis Breust 2014/12/11
+ std::cout
+ << "SIZE: " << STD_RECINT_SIZE
+ << " Time: " << tim.usertime() << ' ' << gmp.usertime()
+ << " Mflops: " << std::scientific << (double(nbloops))/tim.usertime()/1000.0/1000.0 << ' ' << (double(nbloops))/gmp.usertime()/1000.0/1000.0
+ << ' ' << a[(int)( module )& ALEA_MASK] << ' ' << b[(int)(module)& ALEA_MASK] << std::endl ;
+
+ return 0;
+}
+
diff --git a/benchmarks/benchmark-recint_inmul.C b/benchmarks/benchmark-recint_inmul.C
new file mode 100644
index 0000000..4d2e173
--- /dev/null
+++ b/benchmarks/benchmark-recint_inmul.C
@@ -0,0 +1,65 @@
+#include <iostream>
+#include <cstdlib>
+
+#include <recint/recint.h>
+#include <givaro/givtimer.h>
+#include <givaro/givinteger.h>
+
+#if not defined(STD_RECINT_SIZE)
+#define STD_RECINT_SIZE 8
+#endif
+
+#if not defined(LOOPS)
+#define LOOPS 1000000
+#endif
+
+#define ALEA_MAX 64
+#define ALEA_MASK 63
+
+using namespace RecInt;
+
+int main(int argc, char ** argv)
+{
+ size_t nbloops = static_cast<size_t>((argc > 1)? atoi(argv[1]) : LOOPS);
+// std::cerr << "nbloops: " << nbloops << std::endl;
+
+ ruint<STD_RECINT_SIZE> a[ALEA_MAX], d[ALEA_MAX];
+ mpz_class b[ALEA_MAX], c[ALEA_MAX];
+ Givaro::Timer tim,gmp;
+
+
+ // Randomness
+ for (unsigned int i = 0; i < ALEA_MAX; i++) {
+ rand(a[i]);
+ ruint_to_mpz(b[i],a[i]);
+ }
+
+ // Main loop
+ tim.clear(); tim.start();
+ for (unsigned int l = 0; l < nbloops; l++) {
+ d[l & ALEA_MASK] = a[l & ALEA_MASK];
+ d[l & ALEA_MASK] *= a[l & ALEA_MASK];
+ }
+ tim.stop();
+
+ // Main loop
+ gmp.clear(); gmp.start();
+ for (unsigned int l = 0; l < nbloops; l++) {
+ c[l & ALEA_MASK] = b[l & ALEA_MASK];
+ c[l & ALEA_MASK] *= b[l & ALEA_MASK];
+ }
+ gmp.stop();
+
+ ruint<STD_RECINT_SIZE> module; rand(module);
+
+ // -----------
+ // Standard output for benchmark - Alexis Breust 2014/12/11
+ std::cout
+ << "SIZE: " << STD_RECINT_SIZE
+ << " Time: " << tim.usertime() << ' ' << gmp.usertime()
+ << " Mflops: " << std::scientific << (double(nbloops))/tim.usertime()/1000.0/1000.0 << ' ' << (double(nbloops))/gmp.usertime()/1000.0/1000.0
+ << ' ' << a[(int)(module)& ALEA_MASK] << ' ' << b[(int)(module)& ALEA_MASK] << std::endl ;
+
+ return 0;
+}
+
diff --git a/benchmarks/benchmark-recint_inv_arazi.C b/benchmarks/benchmark-recint_inv_arazi.C
index 902d7da..da083bf 100644
--- a/benchmarks/benchmark-recint_inv_arazi.C
+++ b/benchmarks/benchmark-recint_inv_arazi.C
@@ -29,14 +29,19 @@ using namespace RecInt;
int main(int argc, char ** argv)
{
size_t nbloops = static_cast<size_t>((argc > 1)? atoi(argv[1]) : LOOPS);
- Givaro::Timer tim;
+ Givaro::Timer tim, gmp;
ruint<STD_RECINT_SIZE> a[ALEA_MAX];
ruint<STD_RECINT_SIZE> pinv[ALEA_MAX];
+ mpz_class b[ALEA_MAX],c[ALEA_MAX], gmod(1);
+ gmod <<= (1<<STD_RECINT_SIZE);
+// std::cerr << "gmod: " << gmod << std::endl;
+
// Randomness
for (unsigned int i = 0; i < ALEA_MAX; i++) {
rand(a[i]);
if (a[i] % 2 == 0) ++a[i];
+ ruint_to_mpz(b[i],a[i]);
}
// Random
@@ -48,11 +53,23 @@ int main(int argc, char ** argv)
}
tim.stop();
+ gmp.clear(); gmp.start();
+ for (UDItype l = 0; l < nbloops; l++) {
+ mpz_invert(c[l & ALEA_MASK].get_mpz_t(),b[l & ALEA_MASK].get_mpz_t(), gmod.get_mpz_t());
+ }
+ gmp.stop();
+
+ ruint<STD_RECINT_SIZE> module; rand(module);
+
// -----------
// Standard output for benchmark - Alexis Breust 2014/12/11
- std::cout << "Time: " << tim.usertime()
- << " Gflops: " << std::scientific << (double(nbloops))/tim.usertime()/1000.0/1000.0/1000.0 << ' ' << a[(int)(rand(a[0]))& ALEA_MASK] << std::endl ;
-
+ std::cout
+ << "Time: " << tim.usertime()
+ << " Mflops: " << std::scientific << (double(nbloops))/tim.usertime()/1000.0/1000.0 << ' ' << (double(nbloops))/gmp.usertime()/1000.0/1000.0
+ << " SIZE: " << STD_RECINT_SIZE
+ << " GMP time" << gmp.usertime()
+ << ' ' << a[(int)(module)& ALEA_MASK] << ' ' << b[(int)(module)& ALEA_MASK] << std::endl ;
+
return 0;
}
diff --git a/benchmarks/benchmark-recint_mul.C b/benchmarks/benchmark-recint_mul.C
index 0ca5173..10fb28c 100755
--- a/benchmarks/benchmark-recint_mul.C
+++ b/benchmarks/benchmark-recint_mul.C
@@ -3,9 +3,10 @@
#include <recint/recint.h>
#include <givaro/givtimer.h>
+#include <givaro/givinteger.h>
#if not defined(STD_RECINT_SIZE)
-#define STD_RECINT_SIZE 8
+#define STD_RECINT_SIZE 7
#endif
#if not defined(LOOPS)
@@ -23,8 +24,9 @@ int main(int argc, char ** argv)
// std::cerr << "nbloops: " << nbloops << std::endl;
rmint<STD_RECINT_SIZE> a[ALEA_MAX];
+ mpz_class b[ALEA_MAX], gmod;
ruint<STD_RECINT_SIZE> module;
- Givaro::Timer tim;
+ Givaro::Timer tim,gmp;
// For montgomery algorithm, the module must be odd
RecInt::srand(42);
@@ -32,10 +34,13 @@ int main(int argc, char ** argv)
if (module % 2 == 0) module++;
rmint<STD_RECINT_SIZE>::init_module(module);
// std::cerr << "module: " << module << std::endl;
+ ruint_to_mpz(gmod, module);
// Randomness
- for (unsigned int i = 0; i < ALEA_MAX; i++)
+ for (unsigned int i = 0; i < ALEA_MAX; i++) {
rand(a[i]);
+ ruint_to_mpz(b[i],a[i].Value);
+ }
// Main loop
tim.clear(); tim.start();
@@ -44,10 +49,21 @@ int main(int argc, char ** argv)
}
tim.stop();
+ // Main loop
+ gmp.clear(); gmp.start();
+ for (unsigned int l = 0; l < nbloops; l++) {
+ b[(l+1) & ALEA_MASK] = (b[l & ALEA_MASK] * b[l & ALEA_MASK]) % gmod;
+ }
+ gmp.stop();
+
// -----------
// Standard output for benchmark - Alexis Breust 2014/12/11
- std::cout << "Time: " << tim.usertime()
- << " Gflops: " << std::scientific << (double(nbloops))/tim.usertime()/1000.0/1000.0/1000.0 << ' ' << a[(int)(rand(module))& ALEA_MASK] << std::endl ;
+ std::cout
+ << "Time: " << tim.usertime()
+ << " Mflops: " << std::scientific << (double(nbloops))/tim.usertime()/1000.0/1000.0 << ' ' << (double(nbloops))/gmp.usertime()/1000.0/1000.0
+ << " SIZE: " << STD_RECINT_SIZE
+ << " GMP time: " << gmp.usertime()
+ << ' ' << a[(int)(rand(module))& ALEA_MASK] << ' ' << b[(int)(rand(module))& ALEA_MASK] << std::endl ;
return 0;
}
diff --git a/benchmarks/benchmark-recint_mulm.C b/benchmarks/benchmark-recint_mulm.C
new file mode 100644
index 0000000..88d732f
--- /dev/null
+++ b/benchmarks/benchmark-recint_mulm.C
@@ -0,0 +1,74 @@
+#include <iostream>
+#include <cstdlib>
+#include <recint/recint.h>
+#include <givaro/givtimer.h>
+#include <givaro/givinteger.h>
+
+#if not defined(STD_RECINT_SIZE)
+#define STD_RECINT_SIZE 8
+#endif
+
+#if not defined(LOOPS)
+#define LOOPS 1000000
+#endif
+
+#define ALEA_MAX 64
+#define ALEA_MASK 63
+
+using namespace RecInt;
+
+int main(int argc, char ** argv)
+{
+ size_t nbloops = static_cast<size_t>((argc > 1)? atoi(argv[1]) : LOOPS);
+// std::cerr << "nbloops: " << nbloops << std::endl;
+
+ rmint<STD_RECINT_SIZE> a[ALEA_MAX],d[ALEA_MAX];
+ mpz_class b[ALEA_MAX],c[ALEA_MAX], gmod;
+ ruint<STD_RECINT_SIZE> module;
+ Givaro::Timer tim,gmp;
+
+ // For montgomery algorithm, the module must be odd
+ RecInt::srand(42);
+ rand(module);
+ if (module % 2 == 0) module++;
+ rmint<STD_RECINT_SIZE>::init_module(module);
+// std::cerr << "module: " << module << std::endl;
+ ruint_to_mpz(gmod, module);
+
+ // Randomness
+ for (unsigned int i = 0; i < ALEA_MAX; i++) {
+ rand(a[i]);
+ ruint_to_mpz(b[i],a[i].Value);
+ }
+
+ // Main loop
+ tim.clear(); tim.start();
+ for (unsigned int l = 0; l < nbloops; l++) {
+ //mul(a[l & ALEA_MASK], a[l & ALEA_MASK], a[(l+1) & ALEA_MASK]);
+ d[l & ALEA_MASK] = a[l & ALEA_MASK];
+ d[l & ALEA_MASK] *= a[l & ALEA_MASK];
+ }
+ tim.stop();
+
+ // Main loop
+ gmp.clear(); gmp.start();
+ for (unsigned int l = 0; l < nbloops; l++) {
+ c[l & ALEA_MASK] = b[l & ALEA_MASK];
+ c[l & ALEA_MASK] *= b[l & ALEA_MASK];
+ c[l & ALEA_MASK] %= gmod;
+ }
+ gmp.stop();
+
+ rand(module);
+
+ // -----------
+ // Standard output for benchmark - Alexis Breust 2014/12/11
+ std::cout
+ << "SIZE: " << STD_RECINT_SIZE
+ << " Time: " << tim.usertime() << ' ' << gmp.usertime()
+ << " Mflops: " << std::scientific << (double(nbloops))/tim.usertime()/1000.0/1000.0 << ' ' << (double(nbloops))/gmp.usertime()/1000.0/1000.0
+ << ' ' << a[(int)(module)& ALEA_MASK] << ' ' << b[(int)(module)& ALEA_MASK] << std::endl ;
+
+ return 0;
+}
+
diff --git a/benchmarks/benchmark-recint_tmul.C b/benchmarks/benchmark-recint_tmul.C
new file mode 100644
index 0000000..27dfdfa
--- /dev/null
+++ b/benchmarks/benchmark-recint_tmul.C
@@ -0,0 +1,63 @@
+#include <iostream>
+#include <cstdlib>
+
+#include <recint/recint.h>
+#include <givaro/givtimer.h>
+#include <givaro/givinteger.h>
+
+#if not defined(STD_RECINT_SIZE)
+#define STD_RECINT_SIZE 8
+#endif
+
+#if not defined(LOOPS)
+#define LOOPS 1000000
+#endif
+
+#define ALEA_MAX 64
+#define ALEA_MASK 63
+
+using namespace RecInt;
+
+int main(int argc, char ** argv)
+{
+ size_t nbloops = static_cast<size_t>((argc > 1)? atoi(argv[1]) : LOOPS);
+// std::cerr << "nbloops: " << nbloops << std::endl;
+
+ ruint<STD_RECINT_SIZE> a[ALEA_MAX], d[ALEA_MAX];
+ mpz_class b[ALEA_MAX],c[ALEA_MAX];
+ Givaro::Timer tim,gmp;
+
+
+ // Randomness
+ for (unsigned int i = 0; i < ALEA_MAX; i++) {
+ rand(a[i]);
+ ruint_to_mpz(b[i],a[i]);
+ }
+
+ // Main loop
+ tim.clear(); tim.start();
+ for (unsigned int l = 0; l < nbloops; l++) {
+ d[l & ALEA_MASK] = a[l & ALEA_MASK] * a[l & ALEA_MASK];
+ }
+ tim.stop();
+
+ // Main loop
+ gmp.clear(); gmp.start();
+ for (unsigned int l = 0; l < nbloops; l++) {
+ c[l & ALEA_MASK] = b[l & ALEA_MASK] * b[l & ALEA_MASK];
+ }
+ gmp.stop();
+
+ ruint<STD_RECINT_SIZE> module; rand(module);
+
+ // -----------
+ // Standard output for benchmark - Alexis Breust 2014/12/11
+ std::cout
+ << "SIZE: " << STD_RECINT_SIZE
+ << " Time: " << tim.usertime() << ' ' << gmp.usertime()
+ << " Mflops: " << std::scientific << (double(nbloops))/tim.usertime()/1000.0/1000.0 << ' ' << (double(nbloops))/gmp.usertime()/1000.0/1000.0
+ << ' ' << d[(int)(module)& ALEA_MASK] << ' ' << c[(int)(module)& ALEA_MASK] << std::endl ;
+
+ return 0;
+}
+
diff --git a/benchmarks/generic.gnuplot b/benchmarks/generic.gnuplot
new file mode 100644
index 0000000..1875012
--- /dev/null
+++ b/benchmarks/generic.gnuplot
@@ -0,0 +1,24 @@
+set xlabel "Integer bit size"
+set ylabel "Speed (M. arith. op./s)"
+set title "MEANING on an MODEL" # tc lt 2
+#set key below
+set logscale y 10
+#set logscale x 2
+#set ytics (10,1000)
+set xtics ("64" 6,"128" 7,"256" 8,"512" 9,"1024" 10,"2048" 11,"4096" 12, "8192" 13)
+#set xtics tc lt 2
+#set ytics tc lt 2
+#set grid noxtics ytics lt 2
+set grid noxtics ytics
+#set border 4095 lt 7
+#set style line 5 pt 2
+
+
+
+set terminal pdf enhanced color solid lw 2 size 6,4
+set output "rint_FUNCTION.pdf"
+plot [6:13] "output.rint.FUNCTION" using 7:($5) title "GMP-6" with linespoint lt 2 lc 1
+
+set terminal pdf enhanced color solid lw 2 size 6,4
+set output "rint_FUNCTION.pdf"
+replot "output.rint.FUNCTION" using 7:($4) title "RecInt" with fsteps lt 3
diff --git a/benchmarks/perfpublisher.sh b/benchmarks/perfpublisher.sh
index 3f872d4..b010b39 100755
--- a/benchmarks/perfpublisher.sh
+++ b/benchmarks/perfpublisher.sh
@@ -8,12 +8,24 @@ XMLFILE=$1
benchmarks=$2
COMPILER=$3
+# choose gdate on OS X
+if command -v "gdate" >/dev/null; then
+ DATE=gdate
+else
+ DATE=date
+fi
#=================#
# Plateform infos #
#=================#
COMPILERVERSION=$($COMPILER --version 2>&1 | head -1)
-CPUFREQ=$(cat /proc/cpuinfo | grep "MHz" | head -1 |rev | cut -f1 -d' ' | rev)
+
+if command -v "lscpu" >/dev/null; then
+ CPUFREQ=$(lscpu | grep "MHz" | rev | cut -f1 -d' ' | rev)
+else
+ CPUFREQ=$((`sysctl -n hw.cpufrequency`/1000000))
+fi
+
ARCH=$(uname -m)
OSNAME=$(uname -s)
OSVERSION=$(uname -r)
@@ -45,8 +57,8 @@ echo '<report name="benchmarks-report" categ="benchmarks">' >> $XMLFILE
#=======#
echo '<start>' >> $XMLFILE
-echo '<date format="YYYYMMDD" val="'$(date +%Y%m%d)'" />' >> $XMLFILE
-echo '<time format="HHMMSS" val="'$(date +%H%M%S)'" />' >> $XMLFILE
+echo '<date format="YYYYMMDD" val="'$($DATE +%Y%m%d)'" />' >> $XMLFILE
+echo '<time format="HHMMSS" val="'$($DATE +%H%M%S)'" />' >> $XMLFILE
echo '</start>' >> $XMLFILE
#============#
@@ -59,9 +71,9 @@ do
then
#File does not exist: compile it
echo '[Compiling]' $benchmark
- COMPILESTART=$(date +%s%3N)
+ COMPILESTART=$($DATE +%s%3N)
COMPILELOG=$(make $benchmark 2>&1; echo 'Returned state: '$?)
- COMPILEEND=$(date +%s%3N)
+ COMPILEEND=$($DATE +%s%3N)
COMPILETIME=$(($COMPILEEND - $COMPILESTART))
COMPILECHECK=$(echo $COMPILELOG | grep -o '[^ ]*$')
COMPILETIMERELEVANT='true'
@@ -96,7 +108,7 @@ do
EXECUTED='yes'
EXECUTIONLOG=$(./$benchmark 2>&1)
- if [[ ${EXECUTIONLOG,,} != "time:"* ]]
+ if [[ ${EXECUTIONLOG} != "Time:"* ]]
then
#Execution failure
PASSED='no'
@@ -114,7 +126,7 @@ do
EXECUTIONTIME=$(echo $EXECUTIONLOG | cut -d' ' -f2)
PERFORMANCEFLOPS=$(echo $EXECUTIONLOG | cut -d' ' -f4)
EXECUTIONTIMERELEVANT='true'
- if [[ ${PERFORMANCEFLOPS,,} != "irrelevant" ]]
+ if [[ ${PERFORMANCEFLOPS} != "Irrelevant" ]]
then
PERFORMANCEFLOPSRELEVANT='true'
else
diff --git a/configure.ac b/configure.ac
index 6b50b6d..05d18b8 100644
--- a/configure.ac
+++ b/configure.ac
@@ -7,13 +7,16 @@
AC_PREREQ([2.61])
-AC_INIT([Givaro],[4.0.1],[Jean-Guillaume.Dumas at imag.fr],[givaro],
+AC_INIT([Givaro],[4.0.2],[Jean-Guillaume.Dumas at imag.fr],[givaro],
[http://ljk.imag.fr/CASYS/LOGICIELS/givaro])
AC_CONFIG_MACRO_DIR([macros])
AC_CONFIG_AUX_DIR([build-aux])
-AM_INIT_AUTOMAKE([1.8 gnu no-dependencies -Wall -Wno-portability foreign])
AC_CONFIG_HEADERS([config.h])
+
+AC_CANONICAL_TARGET
+
+AM_INIT_AUTOMAKE([1.8 gnu no-dependencies -Wall -Wno-portability foreign])
AX_PREFIX_CONFIG_H(givaro-config.h, __GIVARO)
AC_PATH_PROG(RM, rm, $FALSE)
RM="$RM -f"
@@ -158,6 +161,23 @@ echo "-----------------------------------------------"
echo "-----------------------------------------------"
+# checkes which SIMD instructions are available and defines HAVE_{SSE_4_1,AVX,AVX2}_INSTRUCTIONS and compiler flags
+CUSTOM_SIMD="no"
+FF_CHECK_SIMD
+arch=`echo $target | cut -d"-" -f1`
+if [[ "x$CUSTOM_SIMD" = "xno" ]] ; then
+ AX_CHECK_X86_FEATURES([][])
+else
+ CXXFLAGS="${CXXFLAGS} ${SSEFLAGS} ${AVXFLAGS}"
+fi
+
+dnl forcing -mpfmath=sse flag if SSE is available: Veltkamp split would return a different result on the x87 fpu
+AS_IF([ test "x$SSEFLAGS" != "x" -o "x$X86_FEATURE_CFLAGS" != "x" ],[CXXFLAGS="${CXXFLAGS} -mfpmath=sse"])
+
+dnl With GCC's default ABI version, a __m128 or __m256 are the same types and therefore we cannot
+dnl have overloads for both types without linking error.
+AS_IF([test "x$CCNAM" = "xgcc48"],[CXXFLAGS="${CXXFLAGS} -fabi-version=6"],[])
+
# Machine characteristics
AC_CHECK_SIZEOF(char, 8)
@@ -165,7 +185,13 @@ AC_CHECK_SIZEOF(short, 16)
AC_CHECK_SIZEOF(int, 32)
AC_CHECK_SIZEOF(long, 32)
AC_CHECK_SIZEOF(long long, 64)
-AC_CHECK_SIZEOF(__int64, 64)
+AC_CHECK_SIZEOF(__int64_t, 64)
+
+AC_LANG_CPLUSPLUS
+
+AC_CHECK_TYPE([__int128_t], [AC_TRY_COMPILE([#include <type_traits>], [std::make_unsigned<__int128_t>::type y;],[AC_DEFINE(HAVE_INT128, 1, [Define that compiler allows int128_t types])])])
+
+
# Checks for header files.
AC_HEADER_STDC
@@ -203,8 +229,6 @@ AC_DEFINE_UNQUOTED(INT64, $GIVARO_INT64, Canonical 64-bit data type)
echo "-----------------------------------------------"
# Feature checks
-AC_LANG_CPLUSPLUS
-
echo "-----------------------------------------------"
diff --git a/examples/Integer/primitiveroot.C b/examples/Integer/primitiveroot.C
index 2d2475a..cebf2ad 100644
--- a/examples/Integer/primitiveroot.C
+++ b/examples/Integer/primitiveroot.C
@@ -32,33 +32,33 @@ int main(int argc, char** argv)
#endif
IntNumTheoDom<>::Element a,pr;
if (argc > 1) a = IntNumTheoDom<>::Element(argv[1]); else std::cin >> a;
-
- uint64_t runs;
- Timer tim; tim.clear();
- if (IP.isprime(a)) {
- Integer phin; IP.sub(phin,a,IP.one);
- std::vector<Integer> Lf;
- IP.write(std::cout << "Totient : ", Lf,phin) << std::endl;
- tim.start();
- for(uint64_t i = 0; i < TIMING; ++i)
- IP.prim_root_of_prime(pr, a);
- tim.stop();
- IP.write( std::cout << "Deterministic : ", pr ) << std::endl;
- std::cerr << tim << std::endl;
- }
- tim.start();
-for(uint64_t i = 0; i < TIMING; ++i)
- IP.prim_root(pr, runs, a);
- tim.stop();
- IP.write( std::cout << "Random : ", pr ) << std::endl;
- std::cerr << tim << " (" << runs << " runs)" << std::endl;
- tim.start();
-for(uint64_t i = 0; i < TIMING; ++i)
- IP.lowest_prim_root(pr, a);
- tim.stop();
- IP.write( std::cout << "Lowest : ", pr ) << std::endl;
- std::cerr << tim << std::endl;
-
+
+ uint64_t runs;
+ Timer tim; tim.clear();
+ if (IP.isprime(a)) {
+ Integer phin; IP.sub(phin,a,IP.one);
+ std::vector<Integer> Lf;
+ IP.write(std::cout << "Totient : ", Lf,phin) << std::endl;
+ tim.start();
+ for(uint64_t i = 0; i < TIMING; ++i)
+ IP.prim_root_of_prime(pr, a);
+ tim.stop();
+ IP.write( std::cout << "Deterministic : ", pr ) << std::endl;
+ std::cerr << tim << std::endl;
+ }
+ tim.start();
+ for(uint64_t i = 0; i < TIMING; ++i)
+ IP.prim_root(pr, runs, a);
+ tim.stop();
+ IP.write( std::cout << "Random : ", pr ) << std::endl;
+ std::cerr << tim << " (" << runs << " runs)" << std::endl;
+ tim.start();
+ for(uint64_t i = 0; i < TIMING; ++i)
+ IP.lowest_prim_root(pr, a);
+ tim.stop();
+ IP.write( std::cout << "Lowest : ", pr ) << std::endl;
+ std::cerr << tim << std::endl;
+
return 0;
}
diff --git a/examples/Polynomial/Makefile.am b/examples/Polynomial/Makefile.am
index 06ab40d..3bfc986 100755
--- a/examples/Polynomial/Makefile.am
+++ b/examples/Polynomial/Makefile.am
@@ -30,12 +30,13 @@ LDADD = -L${top_srcdir}/src -L${top_srcdir}/src/.libs -lgivaro $(GMP_LIBS) $(LDF
AM_LDFLAGS=-static
-EXTRA_PROGRAMS=isirred isprimitive trunc_arith pol_arith pol_eval pol_factor interpolate PolynomialCRT highorder bivariate
+EXTRA_PROGRAMS=isirred isprimitive trunc_arith pol_arith pol_eval pol_factor interpolate PolynomialCRT highorder bivariate isgenerator
CLEANFILES=$(EXTRA_PROGRAMS)
interpolate_SOURCES = interpolate.C
isirred_SOURCES = isirred.C
+isgenerator_SOURCES = isgenerator.C
isprimitive_SOURCES = isprimitive.C
pol_arith_SOURCES = pol_arith.C
trunc_arith_SOURCES = trunc_arith.C
diff --git a/examples/Polynomial/isgenerator.C b/examples/Polynomial/isgenerator.C
new file mode 100755
index 0000000..4376349
--- /dev/null
+++ b/examples/Polynomial/isgenerator.C
@@ -0,0 +1,53 @@
+// 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.
+/*! @file examples/Polynomial/isprimitive.C
+ * @ingroup examples
+ * @ingroup polynomials
+ * @example examples/Polynomial/isprimitive.C
+ * @brief NO DOC
+ */
+
+#include <iostream>
+#include <stdlib.h>
+#include <givaro/gfq.h>
+#include <givaro/givpoly1factor.h>
+#include <givaro/givtimer.h>
+
+// using namespace std;
+
+using namespace Givaro;
+
+
+int main(int argc, char** argv)
+{
+ typedef GFqDom<int64_t>::Residu_t UT;
+ UT MOD;
+ if (argc > 1)
+ MOD = (UT)atoi(argv[1]);
+ else
+ std::cin >> MOD;
+ uint64_t expo = 1;
+ if (argc > 2) expo = (uint64_t)atoi(argv[2]);
+
+ GFqDom<int64_t> F(MOD, expo);
+
+ Poly1FactorDom<GFqDom<int64_t>, Dense> FD(F,Indeter("X"));
+ Poly1FactorDom<GFqDom<int64_t>, Dense>::Element P, G;
+ FD.read( std::cin, P );
+ FD.read( std::cin, G );
+
+ Timer tim; tim.clear(); tim.start();
+ bool f = FD.is_prim_root(G, P );
+ tim.stop();
+
+ FD.write(F.write( FD.write( std::cout, G ) << " is " << (f?"":"not ") << " a generator in " ) << " defined with ", P) << std::endl;
+
+ // std::cout << f << std::endl;
+ std::cerr << tim << std::endl;
+
+ return 0;
+}
+
diff --git a/examples/RecInt/Makefile.am b/examples/RecInt/Makefile.am
index 18cce52..594ce09 100755
--- a/examples/RecInt/Makefile.am
+++ b/examples/RecInt/Makefile.am
@@ -4,41 +4,23 @@
# and abiding by the rules of distribution of free software.
# see the COPYRIGHT file for more details.
examples: $(EXTRA_PROGRAMS)
-
AM_CPPFLAGS=-I$(top_srcdir)
-
-OPTFLAGS=
-OPTLINKS=
-
-# for gcc ...
-#OPTFLAGS+= -O7 -funroll-all-loops -felide-constructors -fstrict-aliasing
-#OPTFLAGS+= -frerun-loop-opt -fexpensive-optimizations
-#OPTFLAGS+= -fomit-frame-pointer
-#OPTFLAGS+= -fprefetch-loop-arrays -floop-optimize
-#OPTFLAGS+= -malign-double
-#OPTFLAGS+= -falign-loops -falign-jumps -falign-functions -falign-labels
-#OPTFLAGS+= -fschedule-insns2
-#OPTFLAGS+= -fforce-addr -fforce-mem -fstrength-reduce
-#OPTFLAGS+= -ffast-math
-
-# for icc ...
-#OPTFLAGS+= -fast -Ob2 -ipo_obj -unroll
-#OPTFLAGS+= -parallel -par_report1
-#OPTLINKS+= -parallel
-# icc for itanium2
-#OPTFLAGS+= -tpp2 -mcpu=itanium2
-
AM_CXXFLAGS = @DEFAULT_CFLAGS@
-AM_CPPFLAGS += $(OPTFLAGS) $(GMP_CFLAGS) -I$(top_srcdir)/src/kernel/system -I$(top_srcdir)/src/kernel/recint
+AM_CPPFLAGS += $(OPTFLAGS) $(GMP_CFLAGS)
+AM_CPPFLAGS += -I$(top_srcdir)/src/kernel
+
+#OPTFLAGS=
+#OPTLINKS=
LDADD = $(OPTLINKS) -L${top_srcdir}/src -L${top_srcdir}/src/.libs -lgivaro $(GMP_LIBS) $(LDFLAGS)
AM_LDFLAGS=-static
-EXTRA_PROGRAMS=iterator rsa
+EXTRA_PROGRAMS=recint-iterator rsa extended-int-types
CLEANFILES=$(EXTRA_PROGRAMS)
-iterator_SOURCES = iterator.C
+recint_iterator_SOURCES = recint-iterator.C
+extended_int_types_SOURCES = extended-int-types.C
rsa_SOURCES = rsa.C
# for compilation of new examples
diff --git a/examples/RecInt/extended-int-types.C b/examples/RecInt/extended-int-types.C
new file mode 100644
index 0000000..9355e31
--- /dev/null
+++ b/examples/RecInt/extended-int-types.C
@@ -0,0 +1,62 @@
+/* Copyright Université Grenoble Alpes
+Contributors :
+ Jean-Guillaume DUMAS (Jean-Guillaume.Dumas at imag.fr)
+
+Time-stamp: <27 Jul 16 15:23:32 Jean-Guillaume.Dumas at imag.fr>
+
+This software is a computer program whose purpose is to provide a
+fixed precision arithmetic library.
+
+This software is governed by the CeCILL-B license under French law and
+abiding by the rules of distribution of free software. You can use,
+modify and/ or redistribute the software under the terms of the CeCILL-B
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info".
+
+As a counterpart to the access to the source code and rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty and the software's author, the holder of the
+economic rights, and the successive licensors have only limited
+liability.
+
+In this respect, the user's attention is drawn to the risks associated
+with loading, using, modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean that it is complicated to manipulate, and that also
+therefore means that it is reserved for developers and experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or
+data to be ensured and, more generally, to use and operate it in the
+same conditions as regards security.
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL-B license and that you accept its terms.
+*/
+#include <iostream>
+#include <recint/ruint.h>
+#include <recint/rint.h>
+
+int main(void)
+{
+ RecInt::ruint64 a(1234567890);
+ RecInt::ruint128 b("12345678901234567890123456789012345678901234567890123456789012345678901234567890");
+ RecInt::ruint256 c("12345678901234567890123456789012345678901234567890123456789012345678901234567890");
+ RecInt::ruint512 d("12345678901234567890123456789012345678901234567890123456789012345678901234567890");
+ RecInt::rint64 r(1234567890);
+ RecInt::rint128 s("12345678901234567890123456789012345678901234567890123456789012345678901234567890");
+ RecInt::rint256 t("12345678901234567890123456789012345678901234567890123456789012345678901234567890");
+ RecInt::rint512 u("12345678901234567890123456789012345678901234567890123456789012345678901234567890");
+
+ std::cerr << " 64 bits; " << sizeof(a) << " octets: " << a << std::endl;
+ std::cerr << "128 bits;" << sizeof(b) << " octets: " << b << std::endl;
+ std::cerr << "256 bits;" << sizeof(c) << " octets: " << c << std::endl;
+ std::cerr << "512 bits;" << sizeof(d) << " octets: " << d << std::endl;
+ std::cerr << " 64 bits; " << sizeof(r) << " octets: " << r << std::endl;
+ std::cerr << "128 bits;" << sizeof(s) << " octets: " << s << std::endl;
+ std::cerr << "256 bits;" << sizeof(t) << " octets: " << t << std::endl;
+ std::cerr << "512 bits;" << sizeof(u) << " octets: " << u << std::endl;
+
+ return 0;
+}
+
diff --git a/examples/RecInt/iterator.C b/examples/RecInt/iterator.C
deleted file mode 100644
index a2e8365..0000000
--- a/examples/RecInt/iterator.C
+++ /dev/null
@@ -1,29 +0,0 @@
-#include <iostream>
-#include <cstdlib>
-
-#include <recint/ruint.h>
-
-#if not defined(STD_RECINT_SIZE)
-#define STD_RECINT_SIZE 9
-#endif
-
-using namespace RecInt;
-
-int main(void)
-{
- ruint<STD_RECINT_SIZE> n;
- ruint<STD_RECINT_SIZE>::cr_iterator itr;
-
- RecInt::srand(time(NULL));
- rand(n);
-
- std::cout << "Our " << NBBITS<STD_RECINT_SIZE>::value << "-bit-wide number is: " << std::endl << std::hex << n << std::endl;
-
- const ruint<STD_RECINT_SIZE> n2(n);
- std::cout << std::endl << "Reverse:" << std::endl;
- for (itr = n2.rbegin(); itr != n2.rend(); itr++)
- std::cout << std::hex << *itr << std::endl;
-
- return 0;
-}
-
diff --git a/src/kernel/recint/rint.h b/examples/RecInt/recint-iterator.C
similarity index 60%
copy from src/kernel/recint/rint.h
copy to examples/RecInt/recint-iterator.C
index d107d1e..3267acb 100644
--- a/src/kernel/recint/rint.h
+++ b/examples/RecInt/recint-iterator.C
@@ -1,14 +1,11 @@
-/* recint.h - RecInt library for fixed precision arithmetic
-
-Copyright Université Joseph Fourier - Grenoble
+/* Copyright Université Grenoble Alpes
Contributors :
- Alexis BREUST (alexis.breust at gmail.com 2014)
- Christophe CHABOT (christophechabotcc at gmail.com 2011)
- Jean-Guillaume Dumas
+ Jean-Guillaume DUMAS (Jean-Guillaume.Dumas at imag.fr)
-Time-stamp: <20 Jun 12 10:31:24 Jean-Guillaume.Dumas at imag.fr>
+Time-stamp: <27 Jul 16 12:06:56 Jean-Guillaume.Dumas at imag.fr>
-This software is a computer program whose purpose is to provide an fixed precision arithmetic library.
+This software is a computer program whose purpose is to provide a
+fixed precision arithmetic library.
This software is governed by the CeCILL-B license under French law and
abiding by the rules of distribution of free software. You can use,
@@ -36,47 +33,32 @@ same conditions as regards security.
The fact that you are presently reading this means that you have had
knowledge of the CeCILL-B license and that you accept its terms.
*/
+#include <iostream>
+#include <cstdlib>
+#include <recint/ruint.h>
-#ifndef RINT_H
-#define RINT_H
-
-/* Class definition */
-#include "rrint.h"
-
-/* Arithmetics in rint */
-#include "radd.h"
-#include "rsub.h"
-#include "rmul.h"
-#include "rdiv.h"
-//#include "rexp.h"
-//#include "rgcd.h"
-//#include "rinvmod.h"
-//#include "raddmul.h"
-
-/* Comparisons of rint */
-#include "rcmp.h"
-
-/* Bit manipulation of rint */
-#include "rfiddling.h"
-
-/* Limb manipulation in rint */
-//#include "rmanip.h"
-
-/* Random in rint */
-#include "rrandom.h"
+#if not defined(STD_RECINT_SIZE)
+#define STD_RECINT_SIZE 9
+#endif
-/* Shift for rint */
-//#include "rshift.h"
+using namespace RecInt;
-/* Display for rint */
-#include "rdisplay.h"
+int main(void)
+{
+ ruint<STD_RECINT_SIZE> n;
+ ruint<STD_RECINT_SIZE>::cr_iterator itr;
+
+ RecInt::srand(time(NULL));
+ rand(n);
-/* Internal use for rint */
-//#include "rtools.h"
+ std::cout << "Our " << NBBITS<STD_RECINT_SIZE>::value << "-bit-wide number is: " << std::endl << std::hex << n << std::endl;
-/* GMP conversion system */
-#include "rconvert.h"
+ const ruint<STD_RECINT_SIZE> n2(n);
+ std::cout << std::endl << "Reverse:" << std::endl;
+ for (itr = n2.rbegin(); itr != n2.rend(); itr++)
+ std::cout << std::hex << *itr << std::endl;
-#endif
+ return 0;
+}
diff --git a/givaro.pc.in b/givaro.pc.in
index ed2c60f..274d87e 100644
--- a/givaro.pc.in
+++ b/givaro.pc.in
@@ -1,6 +1,6 @@
/------------------ givaro.pc ------------------------
prefix=@prefix@
-exec_prefix=@prefix@/bin
+exec_prefix=@prefix@
libdir=@prefix@/lib
includedir=@prefix@/include
@@ -10,5 +10,5 @@ URL: http://givaro.forge.imag.fr
Version: @VERSION@
Requires:
Libs: -L at libdir@ -lgivaro @LIBS@
-Cflags: @DEFAULT_CFLAGS@ @CXXFLAGS@
+Cflags: -I at includedir@ @CXXFLAGS@
\-------------------------------------------------------
\ No newline at end of file
diff --git a/macros/CodeChunk/Makefile.am b/macros/CodeChunk/Makefile.am
new file mode 100644
index 0000000..7e799d6
--- /dev/null
+++ b/macros/CodeChunk/Makefile.am
@@ -0,0 +1,18 @@
+# Copyright (c) 2013
+# written by Brice Boyer (briceboyer) <boyer.brice at gmail.com>
+# copied from fflas-ffpack by Clement Pernet
+# adapted from LinBox configuration
+#
+# ========LICENCE========
+# Copyright(c)'1994-2016 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.
+# ========LICENCE========
+#/
+
+
+EXTRA_DIST= \
+ sse.C \
+ avx.C \
diff --git a/macros/CodeChunk/avx.C b/macros/CodeChunk/avx.C
new file mode 100644
index 0000000..2c8bb49
--- /dev/null
+++ b/macros/CodeChunk/avx.C
@@ -0,0 +1,11 @@
+#include <immintrin.h>
+int main() {
+ __m256d P ;
+ double p = 0;
+ P = _mm256_set1_pd(p);
+ P = _mm256_add_pd(P,P);
+#ifdef __try_avx2
+ P = _mm256_fnmadd_pd(P,P,P);
+#endif
+ return 0;
+}
diff --git a/macros/CodeChunk/sse.C b/macros/CodeChunk/sse.C
new file mode 100644
index 0000000..0ace724
--- /dev/null
+++ b/macros/CodeChunk/sse.C
@@ -0,0 +1,12 @@
+#include <immintrin.h>
+
+int main() {
+ // SSE 2
+ __m128d P ;
+ double p = 0;
+ P = _mm_set1_pd(p);
+ P = _mm_add_pd(P,P);
+ // SSE 4.1
+ P = _mm_floor_pd(P);
+ return 0;
+}
diff --git a/macros/ax_check_x86_features.m4 b/macros/ax_check_x86_features.m4
new file mode 100644
index 0000000..22e030b
--- /dev/null
+++ b/macros/ax_check_x86_features.m4
@@ -0,0 +1,77 @@
+# ===========================================================================
+# http://www.gnu.org/software/autoconf-archive/ax_check_x86_features.html
+# ===========================================================================
+#
+# SYNOPSIS
+#
+# AX_CHECK_X86_FEATURES([ACTION-IF-FOUND],[ACTION-IF-NOT-FOUND])
+#
+# DESCRIPTION
+#
+# Checks if the host cpu supports various x86 instruction set, the
+# instructions that will get tested are "mmx, popcnt, sse, sse2, sse3,
+# sse4.1, sse4.2, sse4a, avx, avx2, avx512f, fma, fma4, bmi, bmi2". If the
+# instruction set is supported by the host cpu, the C preprocessor macro
+# HAVE_XXX_INSTRUCTIONS is set to 1. The XXX is up-cased instruction case
+# with dot replaced by underscore. For example, the test for "sse4.2"
+# would export HAVE_SSE4_2_INSTRUCTIONS=1. Also the compiler flag
+# "-msse4.2" would be added to X86_FEATURE_CFLAGS variable, that can be
+# obtained in Makefile.am using @X86_FEATURE_CFLAGS at .
+#
+# If any of the test for the instruction set were succeeded, the configure
+# script would run ACTION-IF-FOUND if it is specified, or append
+# X86_FEATURE_CFLAGS to CFLAGS. If none of the instruction were found,
+# ACTION-IF-NOT-FOUND hook is triggered.
+#
+# This macro requires gcc extended builtin function "__builtin_cpu_init"
+# and "__builtin_cpu_supports" to detect the cpu features. It will error
+# out if the compiler doesn't has these builtins.
+#
+# See also AX_GCC_X86_CPU_SUPPORTS, which is the actual macro that perform
+# the checks for the instruction sets.
+#
+# LICENSE
+#
+# Copyright (c) 2016 Felix Chern <idryman at gmail.com>
+#
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+# Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+# As a special exception, the respective Autoconf Macro's copyright owner
+# gives unlimited permission to copy, distribute and modify the configure
+# scripts that are the output of Autoconf when processing the Macro. You
+# need not follow the terms of the GNU General Public License when using
+# or distributing such scripts, even though portions of the text of the
+# Macro appear in them. The GNU General Public License (GPL) does govern
+# all other use of the material that constitutes the Autoconf Macro.
+#
+# This special exception to the GPL applies to versions of the Autoconf
+# Macro released by the Autoconf Archive. When you make and distribute a
+# modified version of the Autoconf Macro, you may extend this special
+# exception to the GPL to apply to your modified version as well.
+
+#serial 1
+
+AC_DEFUN([AX_CHECK_X86_FEATURES],
+ [m4_foreach_w(
+ [ax_x86_feature],
+ [mmx popcnt sse sse2 sse3 sse4.1 sse4.2 sse4a avx avx2 avx512f fma fma4 bmi bmi2],
+ [AX_GCC_X86_CPU_SUPPORTS(ax_x86_feature,
+ [X86_FEATURE_CFLAGS="$X86_FEATURE_CFLAGS -m[]ax_x86_feature"],
+ [])
+ ])
+ AC_SUBST([X86_FEATURE_CFLAGS])
+ m4_ifval([$1],[$1],
+ [CXXFLAGS="$CXXFLAGS $X86_FEATURE_CFLAGS"])
+ $2
+])
diff --git a/macros/ax_gcc_x86_cpu_supports.m4 b/macros/ax_gcc_x86_cpu_supports.m4
new file mode 100644
index 0000000..a61a14a
--- /dev/null
+++ b/macros/ax_gcc_x86_cpu_supports.m4
@@ -0,0 +1,104 @@
+# ===========================================================================
+# http://www.gnu.org/software/autoconf-archive/ax_gcc_x86_cpu_supports.html
+# ===========================================================================
+#
+# SYNOPSIS
+#
+# AX_GCC_X86_CPU_SUPPORTS(X86-INSTRUCTION-SET,
+# [ACTION-IF-FOUND],[ACTION-IF-NOT-FOUND])
+#
+# DESCRIPTION
+#
+# Checks if the host cpu supports X86-INSTRUCTION-SET. The instruction set
+# that can be tested are "mmx, popcnt, sse, sse2, sse3, sse4.1, sse4.2,
+# sse4a, avx, avx2, avx512f, fma, fma4, bmi, bmi2". If the instruction set
+# is supported by the host cpu, the C preprocessor macro
+# HAVE_XXX_INSTRUCTIONS is set to 1. The XXX is up-cased instruction case
+# with dot replaced by underscore. For example, the test for "sse4.2"
+# would export HAVE_SSE4_2_INSTRUCTIONS=1. This macro requires gcc
+# extended builtin function "__builtin_cpu_init" and
+# "__builtin_cpu_supports" to detect the cpu features. It will error out
+# if the compiler doesn't has these builtins.
+#
+# If the test for the instruction set succeeded, the hook ACTION-IF-FOUND
+# would run. Otherwise the hook ACTION-IF-NOT-FOUND would run if
+# specified.
+#
+# See also AX_CHECK_X86_FEATURES, which checks all the possible
+# instruction set and export the corresponding CFLAGS.
+#
+# LICENSE
+#
+# Copyright (c) 2016 Felix Chern <idryman at gmail.com>
+#
+# This program is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+# Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+# As a special exception, the respective Autoconf Macro's copyright owner
+# gives unlimited permission to copy, distribute and modify the configure
+# scripts that are the output of Autoconf when processing the Macro. You
+# need not follow the terms of the GNU General Public License when using
+# or distributing such scripts, even though portions of the text of the
+# Macro appear in them. The GNU General Public License (GPL) does govern
+# all other use of the material that constitutes the Autoconf Macro.
+#
+# This special exception to the GPL applies to versions of the Autoconf
+# Macro released by the Autoconf Archive. When you make and distribute a
+# modified version of the Autoconf Macro, you may extend this special
+# exception to the GPL to apply to your modified version as well.
+
+#serial 1
+
+AC_DEFUN_ONCE([_AX_GCC_X86_CPU_INIT],
+ [AC_LANG_PUSH([C])
+ AC_CACHE_CHECK([for gcc __builtin_cpu_init function],
+ [ax_cv_gcc_check_x86_cpu_init],
+ [AC_RUN_IFELSE(
+ [AC_LANG_PROGRAM([#include <stdlib.h>],
+ [__builtin_cpu_init ();])
+ ],
+ [ax_cv_gcc_check_x86_cpu_init=yes],
+ [ax_cv_gcc_check_x86_cpu_init=no])])
+ AS_IF([test "X$ax_cv_gcc_check_x86_cpu_init" = "Xno"],
+ [AC_MSG_ERROR([Need GCC to support X86 CPU features tests])])
+])
+
+AC_DEFUN([AX_GCC_X86_CPU_SUPPORTS],
+ [AC_REQUIRE([AC_PROG_CC])
+ AC_REQUIRE([_AX_GCC_X86_CPU_INIT])
+ AC_LANG_PUSH([C])
+ AS_VAR_PUSHDEF([gcc_x86_feature], [AS_TR_SH([ax_cv_gcc_x86_cpu_supports_$1])])
+ AC_CACHE_CHECK([for x86 $1 instruction support],
+ [gcc_x86_feature],
+ [AC_RUN_IFELSE(
+ [AC_LANG_PROGRAM( [#include <stdlib.h> ],
+ [ __builtin_cpu_init ();
+ if (__builtin_cpu_supports("$1"))
+ return 0;
+ return 1;
+ ])],
+ [gcc_x86_feature=yes],
+ [gcc_x86_feature=no]
+ )]
+ )
+ AC_LANG_POP([C])
+ AS_VAR_IF([gcc_x86_feature],[yes],
+ [AC_DEFINE(
+ AS_TR_CPP([HAVE_$1_INSTRUCTIONS]),
+ [1],
+ [Define if $1 instructions are supported])
+ $2],
+ [$3]
+ )
+ AS_VAR_POPDEF([gcc_x86_feature])
+])
diff --git a/macros/simd-check.m4 b/macros/simd-check.m4
new file mode 100644
index 0000000..864eeb1
--- /dev/null
+++ b/macros/simd-check.m4
@@ -0,0 +1,126 @@
+dnl Check for SIMD
+dnl Created by BB, 2014-03-25
+dnl modified by CP, 2016-07-11
+dnl copied from fflas-ffpack by CP
+dnl ========LICENCE========
+dnl Copyright(c)'1994-2016 by The Givaro group
+dnl This file is part of Givaro.
+dnl Givaro is governed by the CeCILL-B license under French law
+dnl and abiding by the rules of distribution of free software.
+dnl see the COPYRIGHT file for more details.
+dnl ========LICENCE========
+dnl
+
+dnl FF_CHECK_SIMD
+dnl
+dnl turn on SSE4.1 AVX, AVX2 extensions if available
+
+AC_DEFUN([FF_CHECK_SIMD],
+[
+ AC_ARG_ENABLE(simd,[AC_HELP_STRING([--disable-simd], [ Disable vectorized instructions: SSE4.1, AVX, AVX2])])
+ AS_IF([ test "x$enable_simd" != "xno" ],
+ [
+ AS_ECHO("SIMD enabled")
+ arch=`echo $target | cut -d"-" -f1`
+ # if we are on a x86 (32 or 64 bits) with gcc>=4.8 then run the AX_CHECK_X86_FEATURES macro
+ AS_IF([test "x$arch" = "xx86_64" -o "x$arch" = "xi686"],
+ [archx86="yes"],
+ [archx86="no"]
+ )
+ AS_IF([ test "x$CCNAM" != "xgcc48" -o "x$archx86" = "xno" ],
+ [
+ CUSTOM_SIMD="yes"
+ echo "Compiling with $CCNAM for a $arch target: running custom checks for SSE4.1 and AVX1,2"
+ AC_MSG_CHECKING(for SSE 4.1)
+ BACKUP_CXXFLAGS=${CXXFLAGS}
+ SSEFLAGS="-msse4.1"
+ CXXFLAGS="${BACKUP_CXXFLAGS} ${SSEFLAGS}"
+ CODE_SSE=`cat macros/CodeChunk/sse.C`
+ AC_TRY_RUN([ ${CODE_SSE} ],
+ [ sse_found="yes" ],
+ [ sse_found="no" ],
+ [
+ echo "cross compiling...disabling"
+ sse_found="no"
+ ])
+ AS_IF([ test "x$sse_found" = "xyes" ],
+ [
+ AC_DEFINE(HAVE_SSE4_1_INSTRUCTIONS,1,[Define if SSE is available])
+ AC_MSG_RESULT(yes)
+ ],
+ [
+ SSEFLAGS=""
+ AC_MSG_RESULT(no)
+ ])
+ CXXFLAGS=${BACKUP_CXXFLAGS}
+
+ dnl Check for AVX
+ AC_MSG_CHECKING(for AVX)
+ CODE_AVX=`cat macros/CodeChunk/avx.C`
+ dnl Intel compilers usually do not require option to enable avx
+ dnl Thus, we test with no option on
+ for switch_avxflags in "" "-mavx"; do
+ CXXFLAGS="${BACKUP_CXXFLAGS} -O0 ${switch_avxflags}"
+ AC_TRY_RUN([ ${CODE_AVX} ],
+ [
+ avx_found="yes"
+ AVXFLAGS=${switch_avxflags}
+ break
+ ],
+ [ avx_found="no" ],
+ [
+ echo "cross compiling...disabling"
+ avx_found="no"
+ break
+ ])
+ done
+
+ dnl Is AVX found?
+ AS_IF([ test "x$avx_found" = "xyes" ],
+ [
+ AC_MSG_RESULT(yes)
+ AC_DEFINE(HAVE_AVX_INSTRUCTIONS,1,[Define if AVX is available])
+
+ dnl Check for AVX2
+ AC_MSG_CHECKING(for AVX2)
+ for switch_avx2flags in "" "-mfma -mavx2"; do
+ CXXFLAGS="${BACKUP_CXXFLAGS} -O0 ${switch_avx2flags}"
+ AC_TRY_RUN(
+ [
+ #define __try_avx2
+ ${CODE_AVX}
+ ],
+ [
+ avx2_found="yes"
+ AVX2FLAGS="${switch_avx2flags}"
+ break
+ ],
+ [ avx2_found="no" ],
+ [
+ echo "cross compiling...disabling"
+ avx2_found = "no"
+ break
+ ])
+ done
+
+ dnl Is AVX2 found?
+ AS_IF([ test "x$avx2_found" = "xyes" ],
+ [
+ AC_MSG_RESULT(yes)
+ AC_DEFINE(HAVE_AVX2_INSTRUCTIONS,1,[Define if AVX2 is available])
+ AVXFLAGS=${AVX2FLAGS}
+ ],
+ [ AC_MSG_RESULT(no) ]
+ )
+ ],
+ [
+ dnl No AVX
+ AC_MSG_RESULT(no)
+ ])
+
+ CXXFLAGS=${BACKUP_CXXFLAGS}
+ ],
+ [ ])
+ ],[ AS_ECHO("SIMD disabled")
+ CUSTOM_SIMD="yes" ])
+])
diff --git a/src/Makefile.am b/src/Makefile.am
index 78a5f39..80ba848 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -20,6 +20,6 @@ libgivaro_la_LIBADD= \
kernel/gmp++/libgmppp.la kernel/bstruct/libgivbstruct.la kernel/integer/libgivinteger.la kernel/memory/libgivmemory.la kernel/rational/libgivrational.la kernel/system/libgivsystem.la kernel/field/libgivfield.la library/tools/libgivtools.la library/poly1/libgivpoly1.la
# soname of libgivaro
-libgivaro_la_LDFLAGS = -version-info 8:1:0
+libgivaro_la_LDFLAGS = -version-info 9:0:0
libgivaro_la_LDFLAGS+= -no-undefined
libgivaro_la_LDFLAGS+=$(LDFLAGS)
diff --git a/src/kernel/field/extension.h b/src/kernel/field/extension.h
index cb8068c..5b7aca9 100755
--- a/src/kernel/field/extension.h
+++ b/src/kernel/field/extension.h
@@ -557,7 +557,7 @@ namespace Givaro {
Element& random(Element& elt) const
{
// Create new random Elements
- elt.resize( (size_t)(_field.order()));
+ elt.resize( (uint64_t)(_field.order()));
for(typename Element::iterator it = elt.begin(); it != elt.end() ; ++ it) {
int64_t tmp = static_cast<int64_t>((double (_givrand()) / double(_GIVRAN_MODULO_)) * double(_size));
(_field.base_field()).init(*it , tmp);
diff --git a/src/kernel/field/gf2.h b/src/kernel/field/gf2.h
old mode 100644
new mode 100755
index e5e8a02..3e7eee5
--- a/src/kernel/field/gf2.h
+++ b/src/kernel/field/gf2.h
@@ -61,8 +61,8 @@ namespace Givaro
GF2 () {}
GF2 (int p, int exp = 1)
{
- assert(p != 2);
- assert(exp != 1);
+ assert(p == 2);
+ assert(exp == 1);
}
/** Copy constructor.
diff --git a/src/kernel/field/gf2.inl b/src/kernel/field/gf2.inl
old mode 100644
new mode 100755
diff --git a/src/kernel/field/gfq.h b/src/kernel/field/gfq.h
index baa2d09..04ff2ca 100755
--- a/src/kernel/field/gfq.h
+++ b/src/kernel/field/gfq.h
@@ -5,7 +5,7 @@
// and abiding by the rules of distribution of free software.
// see the COPYRIGHT file for more details.
// file: gfq.h
-// Time-stamp: <28 Oct 15 20:48:41 Jean-Guillaume.Dumas at imag.fr>
+// Time-stamp: <17 Jul 16 16:12:52 Jean-Guillaume.Dumas at imag.fr>
// date: 1999
// version:
// author: Jean-Guillaume.Dumas
@@ -98,6 +98,14 @@ public:
template<typename Vector>
GFqDom(const UTT P, const UTT e, const Vector& modPoly);
+ // Construction with prescribed irreducible polynomial
+ // and with prescribed generator polynomial
+ // coefficients of the vector should be integers-like
+ // there will be a call to this->init to build the
+ // representation of both polynomials
+ template<typename Vector>
+ GFqDom( const UTT P, const UTT e, const Vector& modPoly, const Vector& genPoly);
+
GFqDom( const GFqDom<TT>& F)
: zero(F.zero),
one(F.one),
diff --git a/src/kernel/field/gfq.inl b/src/kernel/field/gfq.inl
index 1b462a7..0e6dcdd 100755
--- a/src/kernel/field/gfq.inl
+++ b/src/kernel/field/gfq.inl
@@ -1115,6 +1115,77 @@ namespace Givaro {
_plus1[size_t(mOne)] = 0;
}
+ // Construction with prescribed irreducible polynomial
+ // and with prescribed generator polynomial
+ // coefficients of the vector should be integers-like
+ // there will be a call to this->init to build the
+ // representation of both polynomials
+ template<typename TT>
+ template<typename Vector>
+ inline GFqDom<TT>::GFqDom(const UTT P, const UTT e,
+ const Vector& modPoly, const Vector& genPoly):
+ zero(0)
+ , one ((TT) power(P,e) - 1 )
+ , mOne( (P==2)? (one) : ( one >> 1) ) // 1 == -1 in GF(2^k)
+ , _characteristic(P)
+ , _exponent(e)
+ , _q( (UTT) one + 1 )
+ , _qm1 ( (UTT)one )
+ , _log2pol( (UT)_q )
+ , _pol2log( (UT)_q )
+ , _plus1( (UT)_q )
+ , _dcharacteristic( (double)P )
+ {
+
+ // 1 is represented by q-1, zero by 0
+ _log2pol[0] = (UTT)zero;
+
+ GFqDom<TT> Zp(P,1);
+ typedef Poly1FactorDom< GFqDom<TT>, Dense > PolDom;
+ PolDom Pdom( Zp );
+ typename PolDom::Element Ft, F(e+1), G(genPoly.size()), H;
+
+ for( size_t i = 0; i < F.size(); ++i )
+ Zp.init( F[i], modPoly[i]);
+
+ for( size_t i = 0; i < G.size(); ++i )
+ Zp.init( G[i], genPoly[i]);
+
+ Pdom.assign(H,G);
+
+ typedef Poly1PadicDom< GFqDom<TT>, Dense > PadicDom;
+ PadicDom PAD(Pdom);
+
+ PAD.eval(_log2pol[1],H);
+ PAD.eval(_irred, F);
+
+ for (UTT i = 2; i < _qm1; ++i) {
+ Pdom.mulin(H, G);
+ Pdom.modin(H, F);
+ PAD.eval(_log2pol[i], H);
+ }
+ _log2pol[_qm1] = 1;
+
+ _log2pol[0] = 0;
+
+ for (UTT i = 0; i < _q; ++i)
+ _pol2log[ _log2pol[i] ] = i;
+
+ _plus1[0] = 0;
+
+ UTT a,b,r;
+ for (UTT i = 1; i < _q; ++i) {
+ a = _log2pol[i];
+ r = a % P;
+ if (r == (P - 1))
+ b = a - r;
+ else
+ b = a + 1;
+ _plus1[i] = (TT)_pol2log[b] - (TT)_qm1;
+ }
+
+ _plus1[size_t(mOne)] = 0;
+ }
template<typename TT> inline void GFqDom<TT>::Init() {}
diff --git a/src/kernel/recint/recdefine.h b/src/kernel/recint/recdefine.h
index 72addea..126df73 100644
--- a/src/kernel/recint/recdefine.h
+++ b/src/kernel/recint/recdefine.h
@@ -47,6 +47,7 @@ knowledge of the CeCILL-B license and that you accept its terms.
// FIXME Get info at configure-time - A.B.
// Check for anonymous unions + anonymous structs + __uint128_t
// + constructor variables in aggregate
+// NOTE : __GIVARO_HAVE_INT128 is now available to test for __uint128_t
// #define __RECINT_USE_FAST_128
namespace RecInt
diff --git a/src/kernel/recint/rint.h b/src/kernel/recint/rint.h
index d107d1e..c9d0d75 100644
--- a/src/kernel/recint/rint.h
+++ b/src/kernel/recint/rint.h
@@ -41,6 +41,8 @@ knowledge of the CeCILL-B license and that you accept its terms.
#ifndef RINT_H
#define RINT_H
+#include "givaro-config.h"
+
/* Class definition */
#include "rrint.h"
diff --git a/src/kernel/recint/rmgmodule.h b/src/kernel/recint/rmgmodule.h
index cd1f46b..d45e62e 100644
--- a/src/kernel/recint/rmgmodule.h
+++ b/src/kernel/recint/rmgmodule.h
@@ -59,11 +59,11 @@ namespace RecInt
}
-// --------------------------------------------------------------
-// ---------- Arazi & Qi with fast basecase ---------------------
-// -- See [JG Dumas, On Newton-Raphson iteration for ---------
-// -- multiplicative inverses modulo prime powers. ---------
-// -- IEEE Transactions on Computers, 2013]. ---------
+// ---------------------------------------------------------------------
+// ----------------- Arazi & Qi with fast basecase ---------------------
+// -- See [JG Dumas, On Newton-Raphson iteration for --
+// -- multiplicative inverses modulo prime powers. --
+// -- IEEE Transactions on Computers, 63(8), pp 2106-2109, 2014]. --
namespace RecInt
{
diff --git a/src/kernel/recint/rmint.h b/src/kernel/recint/rmint.h
index dd8e7e8..70b0f56 100644
--- a/src/kernel/recint/rmint.h
+++ b/src/kernel/recint/rmint.h
@@ -41,6 +41,8 @@ knowledge of the CeCILL-B license and that you accept its terms.
#ifndef RMINT_H
#define RMINT_H
+#include "givaro-config.h"
+
/* Does not use Montgomery */
#include "rmbrmint.h"
#include "rmbmodule.h"
diff --git a/src/kernel/recint/rmintmg.h b/src/kernel/recint/rmintmg.h
index 7367a81..3201a75 100644
--- a/src/kernel/recint/rmintmg.h
+++ b/src/kernel/recint/rmintmg.h
@@ -41,6 +41,8 @@ knowledge of the CeCILL-B license and that you accept its terms.
#ifndef RMINT_MG_H
#define RMINT_MG_H
+#include "givaro-config.h"
+
/* Use Montgomery */
#include "rmgrmint.h"
#include "rmgmodule.h"
diff --git a/src/kernel/recint/ruconvert.h b/src/kernel/recint/ruconvert.h
index 727ab07..b2b55ac 100644
--- a/src/kernel/recint/ruconvert.h
+++ b/src/kernel/recint/ruconvert.h
@@ -40,7 +40,6 @@ knowledge of the CeCILL-B license and that you accept its terms.
#define RUINT_CONVERT_H
#include <gmpxx.h>
-
#include "rumanip.h" /* reset() */
// --------------------------------------------------------------
@@ -62,8 +61,23 @@ namespace RecInt
// Convert a GMP integer into a ruint
template <size_t K>
inline ruint<K>& mpz_to_ruint(ruint<K>& a, const mpz_class& b) {
- unsigned int i;
- mpz_class c(b);
+ /*
+ mpz_t m0; mpz_init(m0);
+ mpz_init(m0);
+ mpz_set(m0, b.get_mpz_t());
+ size_t bitsize = std::min(m0->_mp_alloc*GMP_LIMB_BITS, 1<<K);
+ size_t n= bitsize/GMP_LIMB_BITS + (bitsize%GMP_LIMB_BITS?1:0);
+ mp_limb_t *target = reinterpret_cast<mp_limb_t*> (&a);
+ reset(a);
+ for (size_t i=0;i<n;i++){
+ target[i]=m0->_mp_d[i];
+ }
+
+ std::cout<<"converting to ruint "<<b<<" to "<<a<<std::endl;
+ */
+
+ unsigned int i;
+ mpz_class c(b);
reset(a);
for (i = 0; i < NBLIMB<K>::value; i++) {
@@ -76,14 +90,30 @@ namespace RecInt
set_limb(a, l, i); c >>= 64;
#endif
}
-
+
+
return a;
}
// Convert a ruint into a GMP integer
template <size_t K>
inline mpz_class& ruint_to_mpz(mpz_class& a, const ruint<K>& b) {
- a = 0;
+ /*
+ mpz_t m;
+ mpz_init2(m,1<<K);
+ m->_mp_size=m->_mp_alloc;
+ size_t n= 1<<(K-6);
+ const limb *src = reinterpret_cast<const limb*> (&b);
+ limb *target = reinterpret_cast<limb*> (m->_mp_d);
+ for (size_t i=n-1;i<(size_t)-1;i--){
+ target[i]=src[i];
+ if (m->_mp_size == (int)i+1 &&src[i]==0 ) m->_mp_size--;
+ }
+
+ a=mpz_class(m);
+ */
+
+ a = 0;
for (auto it(b.rbegin()); it != b.rend(); ++it) {
#if GMP_LIMB_BITS != 64
// GMP does not handle uint64_t, need to break it
@@ -100,7 +130,7 @@ namespace RecInt
#endif
#endif
}
-
+
return a;
}
@@ -119,8 +149,15 @@ namespace RecInt
inline ruint<K>& mpz_t_to_ruint(ruint<K>& a, mpz_srcptr b) {
// TODO Optimize...
mpz_class r(b);
- return RecInt::mpz_to_ruint(a, r);
+ return mpz_to_ruint(a, r);
}
+
+ template <size_t K>
+ inline ruint<K>::ruint(const char* b) {
+ mpz_class m(b);
+ mpz_to_ruint(*this, m);
+ }
+
}
#endif
diff --git a/src/kernel/recint/ruint.h b/src/kernel/recint/ruint.h
index 628b789..c7b62ae 100644
--- a/src/kernel/recint/ruint.h
+++ b/src/kernel/recint/ruint.h
@@ -41,6 +41,8 @@ knowledge of the CeCILL-B license and that you accept its terms.
#ifndef RUINT_H
#define RUINT_H
+#include "givaro-config.h"
+
/* Class definition */
#include "ruruint.h"
diff --git a/src/kernel/recint/rumul.h b/src/kernel/recint/rumul.h
index b561c32..b16754b 100644
--- a/src/kernel/recint/rumul.h
+++ b/src/kernel/recint/rumul.h
@@ -148,11 +148,12 @@ namespace RecInt
// a = ahal = b*c with naive method
// Note: this function is safe, ah|al is correctly computed
// even if b, c are really ah or al
- template <size_t K>
+
+ template <size_t K>
inline void lmul_naive(ruint<K>& ah, ruint<K>& al, const ruint<K>& b, const ruint<K>& c) {
bool rmid, rlow;
ruint<K> bcmid, blcl;
-
+
// Low part
lmul_naive(blcl, b.Low, c.Low);
@@ -170,6 +171,7 @@ namespace RecInt
if (rlow) add_1(ah);
if (rmid) add_1(rmid, ah.High);
}
+
template<>
inline void lmul_naive(ruint<__RECINT_LIMB_SIZE>& ah, ruint<__RECINT_LIMB_SIZE>& al, const ruint<__RECINT_LIMB_SIZE>& b, const ruint<__RECINT_LIMB_SIZE>& c) {
umul_ppmm(ah.Value, al.Value, b.Value, c.Value);
diff --git a/src/kernel/recint/ruruint.h b/src/kernel/recint/ruruint.h
index cb90931..8430806 100644
--- a/src/kernel/recint/ruruint.h
+++ b/src/kernel/recint/ruruint.h
@@ -51,11 +51,12 @@ namespace RecInt
template <size_t K> class ruint {
public:
// this = High|Low
- ruint<K-1> High, Low;
-
+ //ruint<K-1> High, Low;
+ ruint<K-1> Low, High;
+
// Constructors
ruint() {}
- ruint(const ruint<K>& r) : High(r.High), Low(r.Low) {}
+ ruint(const ruint<K>& r) : Low(r.Low), High(r.High) {}
ruint(const ruint<K-1>& rl) : Low(rl) {}
ruint(const double b) : Low((b < 0)? -b : b) { if (b < 0) *this = -*this; }
template <typename T, __RECINT_IS_UNSIGNED(T, int) = 0> ruint(const T b) : Low(b) {}
@@ -63,6 +64,8 @@ namespace RecInt
{ if (b < 0) *this = -*this; }
template <typename T, __RECINT_IS_NOT_FUNDAMENTAL(T, int) = 0> ruint(const T& b)
{ *this = b.operator ruint<K>(); } // Fix for Givaro::Integer
+
+ ruint(const char* s);
// Cast
// Note: Templated operators and specialization make compilers clang + icpc
@@ -127,6 +130,7 @@ namespace RecInt
template <typename T, __RECINT_IS_SIGNED(T, int) = 0> ruint(const T b) : Value(limb(b)) {}
template <typename T, __RECINT_IS_NOT_FUNDAMENTAL(T, int) = 0> ruint(const T& b)
{ *this = b.operator ruint<__RECINT_LIMB_SIZE>(); } // Fix for Givaro::Integer
+ ruint(const char* s);
// Cast
// Brutal too, but icc is kind of peaky - AB 2015/02/11
@@ -189,6 +193,7 @@ namespace RecInt
template <typename T, __RECINT_IS_SIGNED(T, int) = 0> ruint(const T b) : Value(__uint128_t(b)) {}
template <typename T, __RECINT_IS_NOT_FUNDAMENTAL(T, int) = 0> ruint(const T& b)
{ *this = b.operator ruint<__RECINT_LIMB_SIZE+1>(); } // Fix for Givaro::Integer
+ ruint(const char* s);
// Cast
operator float() const { return (float)(Value); }
diff --git a/src/kernel/ring/Makefile.am b/src/kernel/ring/Makefile.am
index f970dab..68d8d8f 100644
--- a/src/kernel/ring/Makefile.am
+++ b/src/kernel/ring/Makefile.am
@@ -26,6 +26,7 @@ pkginclude_HEADERS= \
modular-integer.h \
modular-inttype.h \
modular-log16.h \
+ modular-mulprecomp.inl \
modular-ruint.h \
modular-balanced.h \
modular-balanced-int32.h \
@@ -58,6 +59,8 @@ pkginclude_HEADERS= \
modular-balanced-double.inl \
montgomery-int32.inl\
montgomery-ruint.inl\
- ring-interface.h
+ ring-interface.h \
+ modular-extended.h \
+ modular-extended.inl
EXTRA_DIST=ring.doxy
diff --git a/src/kernel/ring/modular-double.h b/src/kernel/ring/modular-double.h
index 72a70c2..2025b9d 100644
--- a/src/kernel/ring/modular-double.h
+++ b/src/kernel/ring/modular-double.h
@@ -1,3 +1,5 @@
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
// ==========================================================================
// Copyright(c)'1994-2015 by The Givaro group
// This file is part of Givaro.
@@ -30,6 +32,7 @@ public:
// ----- Exported Types and constantes
typedef Modular<double> Self_t;
typedef uint64_t Residu_t;
+ using Compute_t = double;
enum { size_rep = sizeof(Residu_t) };
// ----- Constantes
@@ -63,7 +66,7 @@ public:
template<class T> inline T& characteristic(T& p) const { return p = _lp; }
inline Residu_t cardinality() const { return _lp; }
template<class T> inline T& cardinality(T& p) const { return p = _lp; }
- static inline Residu_t maxCardinality() { return 67108864; }
+ static inline Residu_t maxCardinality() { return 94906266; } // Biggest int n s.t. (n-1)^2 < 2^53
static inline Residu_t minCardinality() { return 2; }
// ----- Checkers
diff --git a/src/kernel/ring/modular-extended.h b/src/kernel/ring/modular-extended.h
new file mode 100644
index 0000000..f8fc99d
--- /dev/null
+++ b/src/kernel/ring/modular-extended.h
@@ -0,0 +1,305 @@
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
+// ==========================================================================
+// Copyright(c)'1994-2016 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.
+// Authors: Bastien Vialla <bastien.vialla at lirmm.fr>
+// ==========================================================================
+
+#ifndef __GIVARO_MODULAR_EXTENDED_H
+#define __GIVARO_MODULAR_EXTENDED_H
+
+#include "givaro/givconfig.h"
+
+#include "givaro/givranditer.h"
+#include "givaro/ring-interface.h"
+#include "givaro/modular-general.h"
+
+namespace Givaro{
+ /*
+ *
+ * Modular double/float allowing big moduli
+ * !!: RandIter does not works, use your own random
+ *
+ */
+ template<class _Element>
+ class ModularExtended : public virtual FiniteFieldInterface<_Element>
+ {
+ public:
+
+ typedef _Element Element;
+ typedef Element* Element_ptr ;
+ typedef const Element ConstElement;
+ typedef const Element* ConstElement_ptr;
+ // ----- Exported Types and constantes
+ typedef ModularExtended<_Element> Self_t;
+ using Compute_t = _Element;
+ typedef uint64_t Residu_t;
+ enum { size_rep = sizeof(Residu_t) };
+
+ private:
+ // Verkampt Split
+ inline void split(const Element x, Element &x_h, Element &x_l) const {
+ Element c;
+ if(std::is_same<Element, double>::value){
+ c = (Element)((1 << 27)+1);
+ }else if(std::is_same<Element, float>::value){
+ c = (Element)((1 << 13)+1);
+ }
+ c = c*x;
+ x_h = c-(c-x);
+ x_l = x - x_h;
+ }
+
+ // Dekker mult, a * b = s + t
+ inline void mult_dekker(const Element a, const Element b, Element &s, Element &t) const{
+ s = a*b;
+ Element ah, al, bh, bl;
+ split(a, ah, al);
+ split(b, bh, bl);
+ t = (al*bl-(((s-ah*bh)-al*bh)-ah*bl));
+ }
+
+ public:
+ // ----- Constantes
+ const Element zero = 0.0;
+ const Element one = 1.0;
+ const Element mOne = -1.0;
+
+ // ----- Constructors
+ ModularExtended() = default;
+
+ template<class XXX> ModularExtended(const XXX& p)
+ : zero(0.0), one(1.0), mOne((Element)p - 1.0), _p((Element)p), _invp((Element)1/(Element)_p), _negp(-_p), _lp((Residu_t)p)
+ {
+ assert(_p >= minCardinality());
+ assert(_p <= maxCardinality());
+ }
+
+ // ----- Accessors
+ inline Element minElement() const override { return zero; }
+ inline Element maxElement() const override { return mOne; }
+
+ // ----- Access to the modulus
+ inline Residu_t residu() const { return _lp; }
+ inline Residu_t size() const { return _lp; }
+ inline Residu_t characteristic() const { return _lp; }
+ template<class T> inline T& characteristic(T& p) const { return p = _lp; }
+ inline Residu_t cardinality() const { return _lp; }
+ template<class T> inline T& cardinality(T& p) const { return p = _lp; }
+ static inline Residu_t maxCardinality() {
+ if(std::is_same<Element, double>::value)
+ return 1125899906842623; // 2^(52-2) - 1
+ else if(std::is_same<Element, float>::value)
+ return 2097151; // 2^(23-2) - 1
+ }
+ static inline Residu_t minCardinality() { return 2; }
+
+ // ----- Checkers
+ inline bool isZero(const Element& a) const override { return a == zero; }
+ inline bool isOne (const Element& a) const override { return a == one; }
+ inline bool isMOne(const Element& a) const override { return a == mOne; }
+ inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
+ inline size_t length(const Element a) const { return size_rep; }
+
+ // ----- Ring-wise operators
+ inline bool operator==(const Self_t& F) const { return _p == F._p; }
+ inline bool operator!=(const Self_t& F) const { return _p != F._p; }
+ inline Self_t& operator=(const Self_t& F)
+ {
+ F.assign(const_cast<Element&>(one), F.one);
+ F.assign(const_cast<Element&>(zero), F.zero);
+ F.assign(const_cast<Element&>(mOne), F.mOne);
+ _p = F._p;
+ _negp = F._negp;
+ _invp = F._invp;
+ _lp= F._lp;
+ return *this;
+ }
+
+
+
+ // ----- Initialisation
+ inline Element& init (Element& x) const
+ {
+ return x = zero;
+ }
+
+ template<typename T>
+ Element& init (Element& r, T a) const{
+ r = Caster<Element>(a);
+ return reduce(r);
+ }
+
+ Element &assign (Element &x, const Element &y) const{
+ return x = y;
+ }
+
+ // ----- Convert and reduce
+ template<typename T> T& convert(T& r, const Element& a) const
+ { return r = static_cast<T>(a); }
+
+ Element& reduce (Element& x, const Element& y) const{
+ x=y;
+ return reduce(x);
+ }
+
+ Element& reduce (Element& x) const ;
+
+ // ----- Classic arithmetic
+ Element& mul(Element& r, const Element& a, const Element& b) const override;
+
+ Element& div(Element& r, const Element& a, const Element& b) const override{
+ return mulin(inv(r, a), b);
+ }
+ Element& add(Element& r, const Element& a, const Element& b) const override {
+ r = a + b;
+ if(r >= _p)
+ r += _negp;
+ return r;
+ }
+ Element& sub(Element& r, const Element& a, const Element& b) const override {
+ r = a - b;
+ if(r < 0)
+ r += _p;
+ return r;
+ }
+ Element& neg(Element& r, const Element& a) const override {
+ r = -a;
+ if(r < 0)
+ r += _p;
+ return r;
+ }
+ Element& inv(Element& x, const Element& y) const override{
+ if(std::is_same<Element, double>::value){
+ int64_t x_int, y_int, tx, ty;
+ x_int = int64_t(_lp);
+ y_int = int64_t(y);
+ tx = 0;
+ ty = 1;
+
+ while (y_int != 0) {
+ // always: gcd (modulus,residue) = gcd (x_int,y_int)
+ // sx*modulus + tx*residue = x_int
+ // sy*modulus + ty*residue = y_int
+ int64_t q = x_int / y_int; // integer quotient
+ int64_t temp = y_int; y_int = x_int - q * y_int;
+ x_int = temp;
+ temp = ty; ty = tx - q * ty;
+ tx = temp;
+ }
+
+ if (tx < 0) tx += int64_t(_p);
+
+ // now x_int = gcd (modulus,residue)
+ return x = Element(tx);
+ }else if(std::is_same<Element, float>::value){
+ int32_t x_int, y_int, tx, ty;
+ x_int = int32_t(_lp);
+ y_int = int32_t(y);
+ tx = 0;
+ ty = 1;
+
+ while (y_int != 0) {
+ // always: gcd (modulus,residue) = gcd (x_int,y_int)
+ // sx*modulus + tx*residue = x_int
+ // sy*modulus + ty*residue = y_int
+ int32_t q = x_int / y_int; // integer quotient
+ int32_t temp = y_int; y_int = x_int - q * y_int;
+ x_int = temp;
+ temp = ty; ty = tx - q * ty;
+ tx = temp;
+ }
+
+ if (tx < 0) tx += int32_t(_p);
+
+ // now x_int = gcd (modulus,residue)
+ return x = Element(tx);
+ }
+ }
+
+ Element& mulin(Element& r, const Element& a) const override {
+ return mul(r, r, a);
+ }
+ Element& divin(Element& r, const Element& y) const override{
+ Element iy;
+ return mulin(r, inv(iy, y));
+ }
+ Element& addin(Element& r, const Element& a) const override {
+ return add(r, r, a);
+ }
+ Element& subin(Element& r, const Element& a) const override {
+ return sub(r, r, a);
+ }
+ Element& negin(Element& r) const override {
+ return neg(r, r);
+ }
+ Element& invin(Element& r) const override {
+ return inv(r, r);
+ }
+
+ // -- axpy: r <- a * x + y
+ // -- axpyin: r <- a * x + r
+ Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override {
+ Element tmp;
+ mul(tmp, a, x);
+ return add(r, tmp, y);
+ }
+ Element& axpyin(Element& r, const Element& a, const Element& x) const override {
+ Element tmp(r);
+ return axpy(r, a, x, tmp);
+ }
+
+ // -- axmy: r <- a * x - y
+ // -- axmyin: r <- a * x - r
+ Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override {
+ Element tmp;
+ mul(tmp, a, x);
+ return sub(r, tmp, y);
+ }
+ Element& axmyin(Element& r, const Element& a, const Element& x) const override {
+ return axmy(r, a, x, r);
+ }
+
+ // -- maxpy: r <- y - a * x
+ // -- maxpyin: r <- r - a * x
+ Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override {
+ Element tmp;
+ mul(tmp, a, x);
+ return sub(r, y, tmp);
+ }
+ Element& maxpyin(Element& r, const Element& a, const Element& x) const override {
+ return maxpy(r, a, x, r);
+ }
+
+ // ----- Random generators
+ typedef ModularRandIter<Self_t> RandIter;
+ typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
+ template< class Random > Element& random(const Random& g, Element& r) const { return init(r, g()); }
+ template< class Random > Element& nonzerorandom(const Random& g, Element& a) const {
+ while (isZero(init(a, g())));
+ return a;
+ }
+
+ // --- IO methods
+ std::istream& read (std::istream& s);
+ std::ostream& write(std::ostream& s) const;
+ std::istream& read (std::istream& s, Element& a) const;
+ std::ostream& write(std::ostream& s, const Element& a) const;
+
+ protected:
+ Element _p = 0;
+ Element _invp = 0;
+ Element _negp = 0;
+ Residu_t _lp = 0;
+
+ };
+
+} // Givaro
+
+#include "givaro/modular-extended.inl"
+
+#endif // __GIVARO_MODULAR_EXTENDED_H
diff --git a/src/kernel/ring/modular-extended.inl b/src/kernel/ring/modular-extended.inl
new file mode 100644
index 0000000..dac373e
--- /dev/null
+++ b/src/kernel/ring/modular-extended.inl
@@ -0,0 +1,285 @@
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
+// ==========================================================================
+// Copyright(c)'1994-2016 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.
+// Authors: Bastien Vialla <bastien.vialla at lirmm.fr>
+// ==========================================================================
+
+#ifndef __GIVARO_MODULAR_EXTENDED_INL
+#define __GIVARO_MODULAR_EXTENDED_INL
+
+namespace Givaro {
+
+ // ----------------
+ // ----- IO methods
+
+ template<>
+ inline
+ std::ostream &ModularExtended<float>::write (std::ostream &os) const
+ {
+ return os << "ModularExtended<float> mod " << _p;
+ }
+
+ template<>
+ inline
+ std::ostream &ModularExtended<double>::write (std::ostream &os) const
+ {
+ return os << "ModularExtended<double> mod " << _p;
+ }
+
+ template<typename _Element>
+ inline
+ std::istream &ModularExtended<_Element>::read (std::istream &is)
+ {
+ is >> _p;
+ return is;
+ }
+
+ template<typename _Element>
+ inline
+ std::ostream &ModularExtended<_Element>::write (std::ostream &os, const Element &x) const
+ {
+ return os << static_cast<int64_t>(x);
+ }
+
+ template<typename _Element>
+ inline
+ std::istream &ModularExtended<_Element>::read (std::istream &is, Element &x) const
+ {
+ int64_t tmp;
+ is >> tmp;
+ init(x,tmp);
+ return is;
+ }
+
+ // --------------------
+ // ----- Initialisation de Modular<double>
+
+ template<>
+ template<>
+ inline
+ typename ModularExtended<double>::Element&
+ ModularExtended<double>::init<const int64_t> (Element& x, const int64_t y) const
+ {
+ x = static_cast<Element>(std::abs(y) % _lp);
+ if (y < 0) negin(x);
+ return x;
+ }
+
+ template<>
+ template<>
+ inline
+ typename ModularExtended<double>::Element&
+ ModularExtended<double>::init<const uint64_t> (Element& x, const uint64_t y) const
+ {
+ return x = static_cast<Element>(y % (uint64_t)(_lp));
+ }
+
+ template<>
+ template<>
+ inline
+ typename ModularExtended<double>::Element&
+ ModularExtended<double>::init<const Integer &> (Element& x, const Integer& y) const
+ {
+ x = static_cast<Element>(y % _lp);
+ if (x < 0) x += _p;
+ return x;
+ }
+
+ template<>
+ template<>
+ inline
+ typename ModularExtended<double>::Element&
+ ModularExtended<double>::init<const typename ModularExtended<double>::Element &> (Element& x, const Element& y) const
+ {
+ return x = y;
+ }
+
+ // --------------------
+ // ----- Initialisation de ModularExtended<float>
+
+ template<>
+ template<>
+ inline
+ typename ModularExtended<float>::Element&
+ ModularExtended<float>::init(typename ModularExtended<float>::Element& r, const double a) const
+ {
+ r = static_cast<float>(std::fmod(a, _p));
+ if (r < 0.f) r += _p;
+ return r;
+ }
+
+ template<>
+ template<>
+ inline
+ typename ModularExtended<float>::Element&
+ ModularExtended<float>::init(typename ModularExtended<float>::Element& r, const int32_t a) const
+ {
+ r = static_cast<Element>(std::abs(a) % _lp);
+ if (a < 0) negin(r);
+ return r;
+ }
+
+ template<>
+ template<>
+ inline
+ typename ModularExtended<float>::Element&
+ ModularExtended<float>::init(typename ModularExtended<float>::Element& r, const uint32_t a) const
+ {
+ return r = static_cast<Element>(a % uint32_t(_lp));
+ }
+
+ template<>
+ template<>
+ inline
+ typename ModularExtended<float>::Element&
+ ModularExtended<float>::init(typename ModularExtended<float>::Element& r, const int64_t a) const
+ {
+ r = static_cast<Element>(std::abs(a) % int64_t(_lp));
+ if (a < 0) negin(r);
+ return r;
+ }
+
+ template<>
+ template<>
+ inline
+ typename ModularExtended<float>::Element&
+ ModularExtended<float>::init(typename ModularExtended<float>::Element& r, const uint64_t a) const
+ {
+ return r = static_cast<Element>(a % uint64_t(_lp));
+ }
+
+ template<>
+ template<>
+ inline
+ typename ModularExtended<float>::Element&
+ ModularExtended<float>::init(typename ModularExtended<float>::Element& r, const Integer& a) const
+ {
+ r = static_cast<Element>(a % _lp);
+ if (a < 0) negin(r);
+ return r;
+ }
+
+ // --------------
+ // Multiplication
+ // --------------
+ template<>
+ inline
+ typename ModularExtended<float>::Element&
+ ModularExtended<float>::mul(typename ModularExtended<float>::Element& r,
+ const typename ModularExtended<float>::Element& a,
+ const typename ModularExtended<float>::Element& b) const {
+#ifdef FP_FAST_FMAF
+ Element abh, abl, pql, q;
+ abh = a * b;
+ abl = fma(a, b, -abh);
+ q = std::floor(abh*_invp);
+ pql = fma (-q, _p, abh);
+ r = abl + pql;
+#else
+ Element abh, abl, pql, pqh, q;
+ mult_dekker(a, b, abh, abl);
+ q = std::floor(abh*_invp);
+ mult_dekker(-q, _p, pqh, pql);
+ r = (abh + pqh) + (abl + pql);
+#endif
+ if(r >= _p)
+ r-= _p;
+ else if(r < 0)
+ r += _p;
+#ifndef NDEBUG
+ assert((r < _p) && (r >= 0));
+#endif
+ return r;
+ }
+
+ template<>
+ inline
+ typename ModularExtended<double>::Element&
+ ModularExtended<double>::mul(typename ModularExtended<double>::Element& r,
+ const typename ModularExtended<double>::Element& a,
+ const typename ModularExtended<double>::Element& b) const {
+#ifdef FP_FAST_FMA
+ Element abh, abl, pql, q;
+ abh = a * b;
+ abl = fma(a, b, -abh);
+ q = std::floor(abh*_invp);
+ pql = fma (-q, _p, abh);
+ r = abl + pql;
+#else
+ Element abh, abl, pql, pqh, q;
+ mult_dekker(a, b, abh, abl);
+ q = std::floor(abh*_invp);
+ mult_dekker(-q, _p, pqh, pql);
+ r = (abh + pqh) + (abl + pql);
+#endif
+ if(r >= _p)
+ r-= _p;
+ else if(r < 0)
+ r += _p;
+#ifndef NDEBUG
+ assert((r < _p) && (r >= 0));
+#endif
+ return r;
+ }
+
+ // --------------
+ // Reduction
+ // --------------
+ template<>
+ inline
+ typename ModularExtended<float>::Element&
+ ModularExtended<float>::reduce (typename ModularExtended<float>::Element& a) const{
+#ifdef FP_FAST_FMAF
+ Element pql, q;
+ q = std::floor(a*_invp);
+ pql = fma (-q, _p, a);
+ a = pql;
+#else
+ Element pql, pqh, q;
+ q = std::floor(a*_invp);
+ mult_dekker(-q, _p, pqh, pql);
+ a = (a + pqh) + pql;
+#endif
+ if(a >= _p)
+ a-= _p;
+ else if(a < 0)
+ a += _p;
+#ifndef NDEBUG
+ assert((a < _p) && (a >= 0));
+#endif
+ return a;
+ }
+
+ template<>
+ inline
+ typename ModularExtended<double>::Element&
+ ModularExtended<double>::reduce (typename ModularExtended<double>::Element& a) const{
+#ifdef FP_FAST_FMA
+ Element pql, q;
+ q = std::floor(a*_invp);
+ pql = fma (-q, _p, a);
+ a = pql;
+#else
+ Element pql, pqh, q;
+ q = std::floor(a*_invp);
+ mult_dekker(-q, _p, pqh, pql);
+ a = (a + pqh) + pql;
+#endif
+ if(a >= _p)
+ a-= _p;
+ else if(a < 0)
+ a += _p;
+#ifndef NDEBUG
+ assert((a < _p) && (a >= 0));
+#endif
+ return a;
+ }
+
+} // Givaro
+
+#endif // __GIVARO_MODULAR_EXTENDED_INL
diff --git a/src/kernel/ring/modular-float.h b/src/kernel/ring/modular-float.h
index bd3d069..792de61 100644
--- a/src/kernel/ring/modular-float.h
+++ b/src/kernel/ring/modular-float.h
@@ -1,3 +1,5 @@
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
// ==========================================================================
// Copyright(c)'1994-2015 by The Givaro group
// This file is part of Givaro.
@@ -35,6 +37,7 @@ namespace Givaro
// ----- Exported Types and constantes
typedef Modular<float> Self_t;
typedef uint32_t Residu_t;
+ using Compute_t = float;
enum { size_rep = sizeof(Residu_t) };
// ----- Constantes
@@ -70,7 +73,10 @@ namespace Givaro
template<class T> inline T& characteristic(T& p) const { return p = _lp; }
inline Residu_t cardinality() const { return _lp; }
template<class T> inline T& cardinality(T& p) const { return p = _lp; }
- static inline Residu_t maxCardinality() { return 2896; }
+ static inline Residu_t maxCardinality() {
+ // 2896 = max { p / 2*p^2 < 2^24 }
+ return 2896;
+ }
static inline Residu_t minCardinality() { return 2; }
// ----- Checkers
diff --git a/src/kernel/ring/modular-int16.h b/src/kernel/ring/modular-int16.h
index 1b98615..099b2ca 100644
--- a/src/kernel/ring/modular-int16.h
+++ b/src/kernel/ring/modular-int16.h
@@ -1,3 +1,5 @@
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
// ==========================================================================
// $Source: /var/lib/cvs/Givaro/src/kernel/zpz/givzpz16std.h,v $
// Copyright(c)'1994-2009 by The Givaro group
@@ -11,16 +13,15 @@
//
// Modified by Pascal Giorgi on 2002/02/13 (pascal.giorgi at ens-lyon.fr)
// Modified by Alexis Breust on 2015/01/06 (alexis.breust at imag.fr)
+// Modified by Romain Lebreton on 2016/06/10 (romain.lebreton at lirmm.fr)
-/*! @file givzpz16std.h
- * @ingroup zpz
- * @brief Arithmetic on Z/pZ, with p a prime number less than 2^14.
- * Modulo typedef is a signed long number. In case it was modified
- * then Bézout algorithm must be changed (coefficient can be negative).
+/*! @file ring/modular-int16.h
+ * @ingroup ring
+ * @brief representation of <code>Z/mZ</code> over \c int16_t .
*/
-
-#ifndef __GIVARO_zpz16std_H
-#define __GIVARO_zpz16std_H
+
+#ifndef __GIVARO_modular_int16_H
+#define __GIVARO_modular_int16_H
#include "givaro/givinteger.h"
#include "givaro/givbasictype.h"
@@ -32,196 +33,215 @@
namespace Givaro {
- /*! @brief This class implement the standard arithmetic with Modulo Elements.
- * - The representation of an integer a in Zpz is the value a % p
- * - m max is 32768
- * - p max is 32749
- * .
- */
- template<typename COMP>
- class Modular<int16_t, COMP> : public virtual FiniteFieldInterface<int16_t>
- {
- public:
-
- // ----- Exported Types and constantes
- using Self_t = Modular<int16_t, COMP>;
- using Residu_t = uint16_t;
- using Compute_t = typename std::make_unsigned<COMP>::type;
- enum { size_rep = sizeof(Residu_t) };
-
- // ----- Representation of vector of the Element
- typedef Element* Array;
- typedef const Element* constArray;
-
- // ----- Constantes
- const Element zero;
- const Element one;
- const Element mOne;
-
- // ----- Constructors
- Modular()
- : zero(static_cast<Element>(0))
- , one(static_cast<Element>(1))
- , mOne(static_cast<Element>(-1))
- , _p(static_cast<Residu_t>(0)) {}
-
- Modular(const Residu_t p)
- : zero(static_cast<Element>(0))
- , one(static_cast<Element>(1))
- , mOne(static_cast<Element>(p-1))
- , _p(static_cast<Residu_t>(p))
- {
- assert(_p >= minCardinality());
- assert(_p <= maxCardinality());
- }
-
- Modular(const Self_t& F)
- : zero(F.zero), one(F.one), mOne(F.mOne), _p(F._p) {}
-
- // ----- Accessors
- inline Element minElement() const override { return zero; }
- inline Element maxElement() const override { return mOne; }
-
- // ----- Access to the modulus
- inline Residu_t residu() const { return _p; }
- inline Residu_t size() const { return _p; }
- inline Residu_t characteristic() const { return _p; }
- inline Residu_t cardinality() const { return _p; }
- template<class T> inline T& characteristic(T& p) const { return p = _p; }
- template<class T> inline T& cardinality(T& p) const { return p = _p; }
-
- static inline Residu_t maxCardinality();
- static inline Residu_t minCardinality() { return 2; }
-
- // ----- Checkers
- inline bool isZero(const Element& a) const override { return a == zero; }
- inline bool isOne (const Element& a) const override { return a == one; }
- inline bool isMOne(const Element& a) const override { return a == mOne; }
- inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
- inline size_t length(const Element a) const { return size_rep; }
-
- // ----- Ring-wise operators
- inline bool operator==(const Self_t& F) const { return _p == F._p; }
- inline bool operator!=(const Self_t& F) const { return _p != F._p; }
- inline Self_t& operator=(const Self_t& F)
- {
- F.assign(const_cast<Element&>(one), F.one);
- F.assign(const_cast<Element&>(zero), F.zero);
- F.assign(const_cast<Element&>(mOne), F.mOne);
- _p = F._p;
- return *this;
- }
-
- // ----- Initialisation
- Element& init (Element& x) const;
- Element& init (Element& x, const float y) const;
- Element& init (Element& x, const double y) const;
- Element& init (Element& x, const int32_t y) const;
- Element& init (Element& x, const uint32_t y) const;
- Element& init (Element& x, const int64_t y) const;
- Element& init (Element& x, const uint64_t y) const;
- Element& init (Element& x, const Integer& y) const;
- template<typename T> Element& init(Element& r, const T& a) const
- { r = Caster<Element>(a); return reduce(r); }
- void init(const size_t, Array a, constArray b) const;
-
- Element& assign (Element& x, const Element& y) const;
- void assign(const size_t sz, Array r, constArray a ) const;
-
- // ----- Convert and reduce
- template<typename T> T& convert(T& r, const Element& a) const
- { return r = static_cast<T>(a); }
-
- Element& reduce (Element& x, const Element& y) const;
- Element& reduce (Element& x) const;
-
- // ----- Classic arithmetic
- Element& mul(Element& r, const Element& a, const Element& b) const override;
- Element& div(Element& r, const Element& a, const Element& b) const override;
- Element& add(Element& r, const Element& a, const Element& b) const override;
- Element& sub(Element& r, const Element& a, const Element& b) const override;
- Element& neg(Element& r, const Element& a) const override;
- Element& inv(Element& r, const Element& a) const override;
-
- Element& mulin(Element& r, const Element& a) const override;
- Element& divin(Element& r, const Element& a) const override;
- Element& addin(Element& r, const Element& a) const override;
- Element& subin(Element& r, const Element& a) const override;
- Element& negin(Element& r) const override;
- Element& invin(Element& r) const override;
-
- // -- axpy: r <- a * x + y
- // -- axpyin: r <- a * x + r
- Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& axpyin(Element& r, const Element& a, const Element& x) const override;
-
- // -- axmy: r <- a * x - y
- // -- axmyin: r <- a * x - r
- Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& axmyin(Element& r, const Element& a, const Element& x) const override;
-
- // -- maxpy: r <- y - a * x
- // -- maxpyin: r <- r - a * x
- Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& maxpyin(Element& r, const Element& a, const Element& x) const override;
-
- // ----- Classic arithmetic on arrays
- void mul(const size_t sz, Array r, constArray a, constArray b) const;
- void mul(const size_t sz, Array r, constArray a, Element b) const;
- void div(const size_t sz, Array r, constArray a, constArray b) const;
- void div(const size_t sz, Array r, constArray a, Element b) const;
- void add(const size_t sz, Array r, constArray a, constArray b) const;
- void add(const size_t sz, Array r, constArray a, Element b) const;
- void sub(const size_t sz, Array r, constArray a, constArray b) const;
- void sub(const size_t sz, Array r, constArray a, Element b) const;
- void neg(const size_t sz, Array r, constArray a) const;
- void inv(const size_t sz, Array r, constArray a) const;
-
- void axpy (const size_t sz, Array r, constArray a, constArray x, constArray c) const;
- void axpyin (const size_t sz, Array r, constArray a, constArray x) const;
- void axmy (const size_t sz, Array r, constArray a, constArray x, constArray c) const;
- void maxpyin (const size_t sz, Array r, constArray a, constArray x) const;
-
- // <- \sum_i a[i], return 1 if a.size() ==0,
- Element& reduceadd ( Element& r, const size_t sz, constArray a ) const;
-
- // <- \prod_i a[i], return 1 if a.size() ==0,
- Element& reducemul ( Element& r, const size_t sz, constArray a ) const;
-
- // <- \sum_i a[i] * b[i]
- Element& dotprod ( Element& r, const size_t sz, constArray a, constArray b ) const;
- Element& dotprod ( Element& r, const int bound, const size_t sz, constArray a, constArray b ) const;
-
- // a -> r: uint32_t to double
- void i2d ( const size_t sz, double* r, constArray a ) const;
-
- // a -> r % p: double to uint32_t % p
- void d2i ( const size_t sz, Array r, const double* a ) const;
-
- // ----- Random generators
- typedef ModularRandIter<Self_t> RandIter;
- typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
- template< class Random > Element& random(Random& g, Element& r) const
- { return init(r, g()); }
- template< class Random > Element& nonzerorandom(Random& g, Element& a) const
- { while (isZero(init(a, g())))
- ;
- return a; }
-
- // --- IO methods
- std::istream& read (std::istream& s);
- std::ostream& write(std::ostream& s) const;
- std::istream& read (std::istream& s, Element& a) const;
- std::ostream& write(std::ostream& s, const Element a) const;
-
- protected:
- // -- data representation of the domain:
- Residu_t _p;
- };
+ /*! @brief This class implement the standard arithmetic with Modulo Elements.
+ * - The representation of an integer a in Zpz is the value a % p
+ * - m max is 32768
+ * - p max is 32749
+ * .
+ */
+ template<typename COMP>
+ class Modular<int16_t, COMP> : public virtual FiniteFieldInterface<int16_t>
+ {
+ public:
+
+ // ----- Exported Types and constantes
+ using Self_t = Modular<int16_t, COMP>;
+ using Residu_t = uint16_t;
+ using Compute_t = typename std::make_unsigned<COMP>::type;
+ enum { size_rep = sizeof(Residu_t) };
+
+ // ----- Representation of vector of the Element
+ typedef Element* Array;
+ typedef const Element* constArray;
+
+ // ----- Constantes
+ const Element zero;
+ const Element one;
+ const Element mOne;
+
+ // ----- Constructors
+ Modular()
+ : zero(static_cast<Element>(0))
+ , one(static_cast<Element>(1))
+ , mOne(static_cast<Element>(-1))
+ , _p(static_cast<Residu_t>(0))
+ , _bitsizep(0) {}
+
+ Modular(const Residu_t p)
+ : zero(static_cast<Element>(0))
+ , one(static_cast<Element>(1))
+ , mOne(static_cast<Element>(p-1))
+ , _p(static_cast<Residu_t>(p))
+ , _bitsizep(0)
+ {
+ assert(_p >= minCardinality());
+ assert(_p <= maxCardinality());
+ Residu_t __p = _p;
+ while (__p != 0) {
+ _bitsizep++;
+ __p >>= 1;
+ }
+ }
+
+ Modular(const Self_t& F)
+ : zero(F.zero), one(F.one), mOne(F.mOne), _p(F._p), _bitsizep(F._bitsizep) {}
+
+ // ----- Accessors
+ inline Element minElement() const override { return zero; }
+ inline Element maxElement() const override { return mOne; }
+
+ // ----- Access to the modulus
+ inline Residu_t residu() const { return _p; }
+ inline Residu_t size() const { return _p; }
+ inline Residu_t characteristic() const { return _p; }
+ inline Residu_t cardinality() const { return _p; }
+ template<class T> inline T& characteristic(T& p) const { return p = _p; }
+ template<class T> inline T& cardinality(T& p) const { return p = _p; }
+
+ static inline Residu_t maxCardinality();
+ static inline Residu_t minCardinality() { return 2; }
+
+ // ----- Checkers
+ inline bool isZero(const Element& a) const override { return a == zero; }
+ inline bool isOne (const Element& a) const override { return a == one; }
+ inline bool isMOne(const Element& a) const override { return a == mOne; }
+ inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
+ inline size_t length(const Element a) const { return size_rep; }
+
+ // ----- Ring-wise operators
+ inline bool operator==(const Self_t& F) const { return _p == F._p; }
+ inline bool operator!=(const Self_t& F) const { return _p != F._p; }
+ inline Self_t& operator=(const Self_t& F)
+ {
+ F.assign(const_cast<Element&>(one), F.one);
+ F.assign(const_cast<Element&>(zero), F.zero);
+ F.assign(const_cast<Element&>(mOne), F.mOne);
+ _p = F._p;
+ _bitsizep = F._bitsizep;
+ return *this;
+ }
+
+ // ----- Initialisation
+ Element& init (Element& x) const;
+ Element& init (Element& x, const float y) const;
+ Element& init (Element& x, const double y) const;
+ Element& init (Element& x, const int32_t y) const;
+ Element& init (Element& x, const uint32_t y) const;
+ Element& init (Element& x, const int64_t y) const;
+ Element& init (Element& x, const uint64_t y) const;
+ Element& init (Element& x, const Integer& y) const;
+ template<typename T> Element& init(Element& r, const T& a) const
+ { r = Caster<Element>(a); return reduce(r); }
+ void init(const size_t, Array a, constArray b) const;
+
+ Element& assign (Element& x, const Element& y) const;
+ void assign(const size_t sz, Array r, constArray a ) const;
+
+ // ----- Convert and reduce
+ template<typename T> T& convert(T& r, const Element& a) const
+ { return r = static_cast<T>(a); }
+
+ Element& reduce (Element& x, const Element& y) const;
+ Element& reduce (Element& x) const;
+
+ // ----- Classic arithmetic
+ Element& mul(Element& r, const Element& a, const Element& b) const override;
+ Element& div(Element& r, const Element& a, const Element& b) const override;
+ Element& add(Element& r, const Element& a, const Element& b) const override;
+ Element& sub(Element& r, const Element& a, const Element& b) const override;
+ Element& neg(Element& r, const Element& a) const override;
+ Element& inv(Element& r, const Element& a) const override;
+
+ Element& mulin(Element& r, const Element& a) const override;
+ Element& divin(Element& r, const Element& a) const override;
+ Element& addin(Element& r, const Element& a) const override;
+ Element& subin(Element& r, const Element& a) const override;
+ Element& negin(Element& r) const override;
+ Element& invin(Element& r) const override;
+
+ // Functions defined in modular-mulprecomp
+ //
+ // void precomp_p (Compute_t& invp) const
+ // Element& mul_precomp_p (Element& r, const Element& a, const Element& b, const Compute_t& invp) const
+ //
+ // void precomp_b (Compute_t& invb, const Element& b) const
+ // void precomp_b (Compute_t& invb, const Element& b, const Compute_t& invp) const
+ // Element& mul_precomp_b (Element& r, const Element& a, const Element& b, const Compute_t& invb) const
+
+#include "modular-mulprecomp.inl"
+
+ // -- axpy: r <- a * x + y
+ // -- axpyin: r <- a * x + r
+ Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& axpyin(Element& r, const Element& a, const Element& x) const override;
+
+ // -- axmy: r <- a * x - y
+ // -- axmyin: r <- a * x - r
+ Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& axmyin(Element& r, const Element& a, const Element& x) const override;
+
+ // -- maxpy: r <- y - a * x
+ // -- maxpyin: r <- r - a * x
+ Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& maxpyin(Element& r, const Element& a, const Element& x) const override;
+
+ // ----- Classic arithmetic on arrays
+ void mul(const size_t sz, Array r, constArray a, constArray b) const;
+ void mul(const size_t sz, Array r, constArray a, Element b) const;
+ void div(const size_t sz, Array r, constArray a, constArray b) const;
+ void div(const size_t sz, Array r, constArray a, Element b) const;
+ void add(const size_t sz, Array r, constArray a, constArray b) const;
+ void add(const size_t sz, Array r, constArray a, Element b) const;
+ void sub(const size_t sz, Array r, constArray a, constArray b) const;
+ void sub(const size_t sz, Array r, constArray a, Element b) const;
+ void neg(const size_t sz, Array r, constArray a) const;
+ void inv(const size_t sz, Array r, constArray a) const;
+
+ void axpy (const size_t sz, Array r, constArray a, constArray x, constArray c) const;
+ void axpyin (const size_t sz, Array r, constArray a, constArray x) const;
+ void axmy (const size_t sz, Array r, constArray a, constArray x, constArray c) const;
+ void maxpyin (const size_t sz, Array r, constArray a, constArray x) const;
+
+ // <- \sum_i a[i], return 1 if a.size() ==0,
+ Element& reduceadd ( Element& r, const size_t sz, constArray a ) const;
+
+ // <- \prod_i a[i], return 1 if a.size() ==0,
+ Element& reducemul ( Element& r, const size_t sz, constArray a ) const;
+
+ // <- \sum_i a[i] * b[i]
+ Element& dotprod ( Element& r, const size_t sz, constArray a, constArray b ) const;
+ Element& dotprod ( Element& r, const int bound, const size_t sz, constArray a, constArray b ) const;
+
+ // a -> r: uint32_t to double
+ void i2d ( const size_t sz, double* r, constArray a ) const;
+
+ // a -> r % p: double to uint32_t % p
+ void d2i ( const size_t sz, Array r, const double* a ) const;
+
+ // ----- Random generators
+ typedef ModularRandIter<Self_t> RandIter;
+ typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
+ template< class Random > Element& random(Random& g, Element& r) const
+ { return init(r, g()); }
+ template< class Random > Element& nonzerorandom(Random& g, Element& a) const
+ { while (isZero(init(a, g())))
+ ;
+ return a; }
+
+ // --- IO methods
+ std::istream& read (std::istream& s);
+ std::ostream& write(std::ostream& s) const;
+ std::istream& read (std::istream& s, Element& a) const;
+ std::ostream& write(std::ostream& s, const Element a) const;
+
+ protected:
+ // -- data representation of the domain:
+ Residu_t _p;
+ size_t _bitsizep;
+ };
} // namespace Givaro
#include "givaro/modular-int16.inl"
-#endif // __GIVARO_zpz16std_H
-// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
+#endif // __GIVARO_modular_int16_H
diff --git a/src/kernel/ring/modular-int32.h b/src/kernel/ring/modular-int32.h
index c724876..f356167 100644
--- a/src/kernel/ring/modular-int32.h
+++ b/src/kernel/ring/modular-int32.h
@@ -1,3 +1,5 @@
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
// ==========================================================================
// $Source: /var/lib/cvs/Givaro/src/kernel/zpz/givzpz32std.h,v $
// Copyright(c)'1994-2009 by The Givaro group
@@ -11,12 +13,11 @@
//
// Modified by Pascal Giorgi on 2002/02/13 (pascal.giorgi at ens-lyon.fr)
// Modified by Alexis Breust on 2014/12/16 (alexis.breust at imag.fr)
+// Modified by Romain Lebreton on 2016/06/10 (romain.lebreton at lirmm.fr)
-/*! @file givzpz32std.h
- * @ingroup zpz
- * @brief Arithmetic on \c Z/pZ, with p a prime number less than \f$2^32\f$.
- * Modulo typedef is a signed long number. In case it was modified
- * then Bézout algorithm must be changed (coefficient can be negative).
+/*! @file ring/modular-int32.h
+ * @ingroup ring
+ * @brief representation of <code>Z/mZ</code> over \c int32_t .
*/
#ifndef __GIVARO_modular_int32_H
@@ -30,154 +31,179 @@
namespace Givaro {
- /*! @brief This class implement the standard arithmetic with Modulo Elements.
- * - The representation of an integer a in Zpz is the value a % p
- * - m max is 46341
- * - p max is 46337
- * .
- */
- template<typename COMP>
- class Modular<int32_t, COMP> : public virtual FiniteFieldInterface<int32_t>
- {
- public:
-
- // ----- Exported Types and constantes
- using Self_t = Modular<int32_t, COMP>;
- using Residu_t = uint32_t;
- using Compute_t = typename std::make_unsigned<COMP>::type;
- enum { size_rep = sizeof(Residu_t) };
-
- // ----- Representation of vector of the Element
- typedef Element* Array;
- typedef const Element* constArray;
-
- // ----- Constantes
- const Element zero;
- const Element one;
- const Element mOne;
-
- // ----- Constructors
- Modular()
- : zero(0), one(1), mOne(-1), _p(0), _dp(0.0) {}
-
- Modular(Residu_t p)
- : zero(0), one(1), mOne((Element)p-1), _p(p), _dp((double)p)
- {
- assert(_p >= minCardinality());
- assert(_p <= maxCardinality());
- }
-
- Modular(const Self_t& F)
- : zero(F.zero), one(F.one), mOne(F.mOne), _p(F._p), _dp(F._dp) {}
-
- // ----- Accessors
- inline Element minElement() const override { return zero; }
- inline Element maxElement() const override { return mOne; }
-
- // ----- Access to the modulus
- inline Residu_t residu() const { return _p; }
- inline Residu_t size() const { return _p; }
- inline Residu_t characteristic() const { return _p; }
- inline Residu_t cardinality() const { return _p; }
- template<class T> inline T& characteristic(T& p) const { return p = _p; }
- template<class T> inline T& cardinality(T& p) const { return p = _p; }
-
- static inline Residu_t maxCardinality();
- static inline Residu_t minCardinality() { return 2; }
-
- // ----- Checkers
- inline bool isZero(const Element& a) const override { return a == zero; }
- inline bool isOne (const Element& a) const override { return a == one; }
- inline bool isMOne(const Element& a) const override { return a == mOne; }
- inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
- inline size_t length(const Element a) const { return size_rep; }
-
- // ----- Ring-wise operators
- inline bool operator==(const Self_t& F) const { return _p == F._p; }
- inline bool operator!=(const Self_t& F) const { return _p != F._p; }
- inline Self_t& operator=(const Self_t& F)
- {
- F.assign(const_cast<Element&>(one), F.one);
- F.assign(const_cast<Element&>(zero), F.zero);
- F.assign(const_cast<Element&>(mOne), F.mOne);
- _p = F._p;
- _dp = F._dp;
- return *this;
- }
-
- // ----- Initialisation
- Element& init (Element& x) const;
- Element& init (Element& x, const float y) const;
- Element& init (Element& x, const double y) const;
- Element& init (Element& x, const int64_t y) const;
- Element& init (Element& x, const uint64_t y) const;
- Element& init (Element& x, const Integer& y) const;
- template<typename T> Element& init(Element& r, const T& a) const
- { return reduce(Caster<Element,T>(r,a)); }
-
- Element& assign (Element& x, const Element& y) const;
-
- // ----- Convert and reduce
- template<typename T> T& convert(T& r, const Element& a) const
- { return Caster<T,Element>(r,a); }
-
- Element& reduce (Element& x, const Element& y) const;
- Element& reduce (Element& x) const;
-
- // ----- Classic arithmetic
- Element& mul(Element& r, const Element& a, const Element& b) const override;
- Element& div(Element& r, const Element& a, const Element& b) const override;
- Element& add(Element& r, const Element& a, const Element& b) const override;
- Element& sub(Element& r, const Element& a, const Element& b) const override;
- Element& neg(Element& r, const Element& a) const override;
- Element& inv(Element& r, const Element& a) const override;
-
- Element& mulin(Element& r, const Element& a) const override;
- Element& divin(Element& r, const Element& a) const override;
- Element& addin(Element& r, const Element& a) const override;
- Element& subin(Element& r, const Element& a) const override;
- Element& negin(Element& r) const override;
- Element& invin(Element& r) const override;
-
- // -- axpy: r <- a * x + y
- // -- axpyin: r <- a * x + r
- Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& axpyin(Element& r, const Element& a, const Element& x) const override;
-
- // -- axmy: r <- a * x - y
- // -- axmyin: r <- a * x - r
- Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& axmyin(Element& r, const Element& a, const Element& x) const override;
-
- // -- maxpy: r <- y - a * x
- // -- maxpyin: r <- r - a * x
- Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& maxpyin(Element& r, const Element& a, const Element& x) const override;
-
- // ----- Random generators
- typedef ModularRandIter<Self_t> RandIter;
- typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
- template< class Random > Element& random(Random& g, Element& r) const
- { return init(r, g()); }
- template< class Random > Element& nonzerorandom(Random& g, Element& a) const
- { while (isZero(init(a, g())))
- ;
- return a; }
-
- // --- IO methods
- std::istream& read (std::istream& s);
- std::ostream& write(std::ostream& s) const;
- std::istream& read (std::istream& s, Element& a) const;
- std::ostream& write(std::ostream& s, const Element a) const;
-
- protected:
- // -- data representation of the domain:
- Residu_t _p;
- double _dp;
- };
+ /*! @brief This class implement the standard arithmetic with Modulo Elements.
+ * - The representation of an integer a in Zpz is the value a % p
+ * - m max is 46341
+ * - p max is 46337
+ * .
+ */
+ template<typename COMP>
+ class Modular<int32_t, COMP> : public virtual FiniteFieldInterface<int32_t>
+ {
+ public:
+
+ // ----- Exported Types and constantes
+ using Self_t = Modular<int32_t, COMP>;
+ using Residu_t = uint32_t;
+ using Compute_t = typename std::make_unsigned<COMP>::type;
+ enum { size_rep = sizeof(Residu_t) };
+
+ // ----- Representation of vector of the Element
+ typedef Element* Array;
+ typedef const Element* constArray;
+
+ // ----- Constantes
+ const Element zero;
+ const Element one;
+ const Element mOne;
+
+ // ----- Constructors
+ Modular()
+ : zero(static_cast<Element>(0))
+ , one(static_cast<Element>(1))
+ , mOne(static_cast<Element>(-1))
+ , _p(static_cast<Residu_t>(0))
+ , _bitsizep(0) {}
+
+ Modular(const Residu_t p)
+ : zero(static_cast<Element>(0))
+ , one(static_cast<Element>(1))
+ , mOne(static_cast<Element>(p-1))
+ , _p(static_cast<Residu_t>(p))
+ , _bitsizep(0)
+ {
+ assert(_p >= minCardinality());
+ assert(_p <= maxCardinality());
+ Residu_t __p = _p;
+ while (__p != 0) {
+ _bitsizep++;
+ __p >>= 1;
+ }
+ }
+
+ Modular(const Self_t& F)
+ : zero(F.zero), one(F.one), mOne(F.mOne), _p(F._p), _dp(F._dp), _bitsizep(F._bitsizep) {}
+
+ // ----- Accessors
+ inline Element minElement() const override { return zero; }
+ inline Element maxElement() const override { return mOne; }
+
+ // ----- Access to the modulus
+ inline Residu_t residu() const { return _p; }
+ inline Residu_t size() const { return _p; }
+ inline Residu_t characteristic() const { return _p; }
+ inline Residu_t cardinality() const { return _p; }
+ template<class T> inline T& characteristic(T& p) const { return p = _p; }
+ template<class T> inline T& cardinality(T& p) const { return p = _p; }
+
+ static inline Residu_t maxCardinality();
+ static inline Residu_t minCardinality() { return 2; }
+
+ // ----- Checkers
+ inline bool isZero(const Element& a) const override { return a == zero; }
+ inline bool isOne (const Element& a) const override { return a == one; }
+ inline bool isMOne(const Element& a) const override { return a == mOne; }
+ inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
+ inline size_t length(const Element a) const { return size_rep; }
+
+ // ----- Ring-wise operators
+ inline bool operator==(const Self_t& F) const { return _p == F._p; }
+ inline bool operator!=(const Self_t& F) const { return _p != F._p; }
+ inline Self_t& operator=(const Self_t& F)
+ {
+ F.assign(const_cast<Element&>(one), F.one);
+ F.assign(const_cast<Element&>(zero), F.zero);
+ F.assign(const_cast<Element&>(mOne), F.mOne);
+ _p = F._p;
+ _dp = F._dp;
+ _bitsizep = F._bitsizep;
+ return *this;
+ }
+
+ // ----- Initialisation
+ Element& init (Element& x) const;
+ Element& init (Element& x, const float y) const;
+ Element& init (Element& x, const double y) const;
+ Element& init (Element& x, const int64_t y) const;
+ Element& init (Element& x, const uint64_t y) const;
+ Element& init (Element& x, const Integer& y) const;
+ template<typename T> Element& init(Element& r, const T& a) const
+ { return reduce(Caster<Element,T>(r,a)); }
+
+ Element& assign (Element& x, const Element& y) const;
+
+ // ----- Convert and reduce
+ template<typename T> T& convert(T& r, const Element& a) const
+ { return Caster<T,Element>(r,a); }
+
+ Element& reduce (Element& x, const Element& y) const;
+ Element& reduce (Element& x) const;
+
+ // ----- Classic arithmetic
+ Element& mul(Element& r, const Element& a, const Element& b) const override;
+ Element& div(Element& r, const Element& a, const Element& b) const override;
+ Element& add(Element& r, const Element& a, const Element& b) const override;
+ Element& sub(Element& r, const Element& a, const Element& b) const override;
+ Element& neg(Element& r, const Element& a) const override;
+ Element& inv(Element& r, const Element& a) const override;
+
+ Element& mulin(Element& r, const Element& a) const override;
+ Element& divin(Element& r, const Element& a) const override;
+ Element& addin(Element& r, const Element& a) const override;
+ Element& subin(Element& r, const Element& a) const override;
+ Element& negin(Element& r) const override;
+ Element& invin(Element& r) const override;
+
+ // Functions defined in modular-mulprecomp
+ //
+ // void precomp_p (Compute_t& invp) const
+ // Element& mul_precomp_p (Element& r, const Element& a, const Element& b, const Compute_t& invp) const
+ //
+ // void precomp_b (Compute_t& invb, const Element& b) const
+ // void precomp_b (Compute_t& invb, const Element& b, const Compute_t& invp) const
+ // Element& mul_precomp_b (Element& r, const Element& a, const Element& b, const Compute_t& invb) const
+
+#include"modular-mulprecomp.inl"
+
+ // -- axpy: r <- a * x + y
+ // -- axpyin: r <- a * x + r
+ Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& axpyin(Element& r, const Element& a, const Element& x) const override;
+
+ // -- axmy: r <- a * x - y
+ // -- axmyin: r <- a * x - r
+ Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& axmyin(Element& r, const Element& a, const Element& x) const override;
+
+ // -- maxpy: r <- y - a * x
+ // -- maxpyin: r <- r - a * x
+ Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& maxpyin(Element& r, const Element& a, const Element& x) const override;
+
+ // ----- Random generators
+ typedef ModularRandIter<Self_t> RandIter;
+ typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
+ template< class Random > Element& random(Random& g, Element& r) const
+ { return init(r, g()); }
+ template< class Random > Element& nonzerorandom(Random& g, Element& a) const
+ { while (isZero(init(a, g())))
+ ;
+ return a; }
+
+ // --- IO methods
+ std::istream& read (std::istream& s);
+ std::ostream& write(std::ostream& s) const;
+ std::istream& read (std::istream& s, Element& a) const;
+ std::ostream& write(std::ostream& s, const Element a) const;
+
+ protected:
+ // -- data representation of the domain:
+ Residu_t _p;
+ double _dp;
+ size_t _bitsizep;
+ };
}
#include "givaro/modular-int32.inl"
-#endif
-
+#endif // __GIVARO_modular_int32_H
diff --git a/src/kernel/ring/modular-int32.inl b/src/kernel/ring/modular-int32.inl
index fad6e82..2e8d532 100644
--- a/src/kernel/ring/modular-int32.inl
+++ b/src/kernel/ring/modular-int32.inl
@@ -44,7 +44,7 @@ namespace Givaro {
inline typename Modular<int32_t, COMP>::Element& Modular<int32_t, COMP>::mul
(Element& r, const Element& a, const Element& b) const
{
- return __GIVARO_MODULAR_INTEGER_MUL(r,_p,a,b);
+ return __GIVARO_MODULAR_INTEGER_MUL(r,_p,a,b);
}
template<typename COMP>
@@ -86,9 +86,8 @@ namespace Givaro {
template<typename COMP>
inline typename Modular<int32_t, COMP>::Element& Modular<int32_t, COMP>::mulin
(Element& r, const Element& a) const
- {
- r = Element(Compute_t(r)*Compute_t(a) % Compute_t(_p));
- return r;
+ {
+ return __GIVARO_MODULAR_INTEGER_MULIN(r, _p, a);
}
template<typename COMP>
@@ -129,8 +128,8 @@ namespace Givaro {
(Element& r) const
{
return inv(r, r);
- }
-
+ }
+
template<typename COMP>
inline typename Modular<int32_t, COMP>::Element& Modular<int32_t, COMP>::axpy
(Element& r, const Element& a, const Element& b, const Element& c) const
diff --git a/src/kernel/ring/modular-int64.h b/src/kernel/ring/modular-int64.h
index 376fdd9..256a1dc 100644
--- a/src/kernel/ring/modular-int64.h
+++ b/src/kernel/ring/modular-int64.h
@@ -1,3 +1,5 @@
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
// ==========================================================================
// $Source: /var/lib/cvs/Givaro/src/kernel/zpz/givzpz64std.h,v $
// Copyright(c)'1994-2009 by The Givaro group
@@ -8,17 +10,15 @@
// Authors: T. Gautier
// $Id: givzpz64std.h,v 1.21 2011-02-04 14:11:46 jgdumas Exp $
// ==========================================================================
+// Modified by Romain Lebreton on 2016/06/10 (romain.lebreton at lirmm.fr)
-/*! @file givzpz64std.h
- * @ingroup zpz
- * @brief Zpz on 64bit words
- * Arithmetic on Z/pZ, with p a prime number less than 2^64
- * Modulo typedef is a signed long number. In case it was modified
- * then Bézout algorithm must be changed (coefficient can be negative).
+/*! @file ring/modular-int64.h
+ * @ingroup ring
+ * @brief representation of <code>Z/mZ</code> over \c int64_t .
*/
-#ifndef __GIVARO_zpz64std_H
-#define __GIVARO_zpz64std_H
+#ifndef __GIVARO_modular_int64_H
+#define __GIVARO_modular_int64_H
#include "givaro/givinteger.h"
#include "givaro/givcaster.h"
@@ -29,183 +29,208 @@
namespace Givaro {
- /*! @brief This class implement the standard arithmetic with Modulo Elements.
- * - The representation of an integer a in Zpz is the value a % p
- * - m max is 3037000499
- * - p max is 3037000493
- * .
- */
- template<typename COMP>
- class Modular<int64_t, COMP> : public virtual FiniteFieldInterface<int64_t>
- {
- public:
-
- // ----- Exported Types and constantes
- using Self_t = Modular<int64_t, COMP>;
- using Residu_t = uint64_t;
- using Compute_t = typename std::make_unsigned<COMP>::type;
- enum { size_rep = sizeof(Residu_t) };
-
- // ----- Representation of vector of the Element
- typedef Element* Array;
- typedef const Element* constArray;
-
- // ----- Constantes
- const Element zero;
- const Element one;
- const Element mOne;
-
- // ----- Constructors
- Modular()
- : zero(0), one(1), mOne(-1), _p(0) {}
-
- Modular(Residu_t p)
- : zero(0), one(1), mOne((Element)p-1), _p(p)
- {
- assert(_p >= minCardinality());
- assert(_p <= maxCardinality());
- }
-
- Modular(const Self_t& F)
- : zero(F.zero), one(F.one), mOne(F.mOne), _p(F._p) {}
-
- // ----- Accessors
- inline Element minElement() const override { return zero; }
- inline Element maxElement() const override { return mOne; }
-
- // ----- Access to the modulus
- inline Residu_t residu() const { return _p; }
- inline Residu_t size() const { return _p; }
- inline Residu_t characteristic() const { return _p; }
- inline Residu_t cardinality() const { return _p; }
- template<class T> inline T& characteristic(T& p) const { return p = _p; }
- template<class T> inline T& cardinality(T& p) const { return p = _p; }
- static inline Residu_t maxCardinality() { return 3037000499_ui64; }
- static inline Residu_t minCardinality() { return 2_ui64; }
-
- // ----- Checkers
- inline bool isZero(const Element& a) const override { return a == zero; }
- inline bool isOne (const Element& a) const override { return a == one; }
- inline bool isMOne(const Element& a) const override { return a == mOne; }
- inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
- inline size_t length(const Element a) const { return size_rep; }
-
- // ----- Ring-wise operators
- inline bool operator==(const Self_t& F) const { return _p == F._p; }
- inline bool operator!=(const Self_t& F) const { return _p != F._p; }
- inline Self_t& operator=(const Self_t& F)
- {
- F.assign(const_cast<Element&>(one), F.one);
- F.assign(const_cast<Element&>(zero), F.zero);
- F.assign(const_cast<Element&>(mOne), F.mOne);
- _p = F._p;
- return *this;
- }
-
- // ----- Initialisation
- Element& init (Element& x) const;
- Element& init (Element& x, const Integer& y) const;
- template<typename T> Element& init(Element& r, const T& a) const
- { r = Caster<Element>(a); return reduce(r); }
- void init(const size_t, Array a, constArray b) const;
-
- Element& assign (Element& x, const Element& y) const;
- void assign(const size_t sz, Array r, constArray a ) const;
-
- // ----- Convert and reduce
- template<typename T> T& convert(T& r, const Element& a) const
- { return r = static_cast<T>(a); }
-
- Element& reduce (Element& x, const Element& y) const;
- Element& reduce (Element& x) const;
-
- // ----- Classic arithmetic
- Element& mul(Element& r, const Element& a, const Element& b) const override;
- Element& div(Element& r, const Element& a, const Element& b) const override;
- Element& add(Element& r, const Element& a, const Element& b) const override;
- Element& sub(Element& r, const Element& a, const Element& b) const override;
- Element& neg(Element& r, const Element& a) const override;
- Element& inv(Element& r, const Element& a) const override;
-
- Element& mulin(Element& r, const Element& a) const override;
- Element& divin(Element& r, const Element& a) const override;
- Element& addin(Element& r, const Element& a) const override;
- Element& subin(Element& r, const Element& a) const override;
- Element& negin(Element& r) const override;
- Element& invin(Element& r) const override;
-
- // -- axpy: r <- a * x + y
- // -- axpyin: r <- a * x + r
- Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& axpyin(Element& r, const Element& a, const Element& x) const override;
-
- // -- axmy: r <- a * x - y
- // -- axmyin: r <- a * x - r
- Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& axmyin(Element& r, const Element& a, const Element& x) const override;
-
- // -- maxpy: r <- y - a * x
- // -- maxpyin: r <- r - a * x
- Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& maxpyin(Element& r, const Element& a, const Element& x) const override;
-
- // ----- Classic arithmetic on arrays
- void mul(const size_t sz, Array r, constArray a, constArray b) const;
- void mul(const size_t sz, Array r, constArray a, Element b) const;
- void div(const size_t sz, Array r, constArray a, constArray b) const;
- void div(const size_t sz, Array r, constArray a, Element b) const;
- void add(const size_t sz, Array r, constArray a, constArray b) const;
- void add(const size_t sz, Array r, constArray a, Element b) const;
- void sub(const size_t sz, Array r, constArray a, constArray b) const;
- void sub(const size_t sz, Array r, constArray a, Element b) const;
- void neg(const size_t sz, Array r, constArray a) const;
- void inv(const size_t sz, Array r, constArray a) const;
-
- void axpy (const size_t sz, Array r, constArray a, constArray x, constArray c) const;
- void axpyin (const size_t sz, Array r, constArray a, constArray x) const;
- void axmy (const size_t sz, Array r, constArray a, constArray x, constArray c) const;
- void maxpyin (const size_t sz, Array r, constArray a, constArray x) const;
-
- // <- \sum_i a[i], return 1 if a.size() ==0,
- Element& reduceadd ( Element& r, const size_t sz, constArray a ) const;
-
- // <- \prod_i a[i], return 1 if a.size() ==0,
- Element& reducemul ( Element& r, const size_t sz, constArray a ) const;
-
- // <- \sum_i a[i] * b[i]
- Element& dotprod ( Element& r, const size_t sz, constArray a, constArray b ) const;
- Element& dotprod ( Element& r, const int bound, const size_t sz, constArray a, constArray b ) const;
-
- // a -> r: uint32_t to double
- void i2d ( const size_t sz, double* r, constArray a ) const;
-
- // a -> r % p: double to uint32_t % p
- void d2i ( const size_t sz, Array r, const double* a ) const;
-
- // ----- Random generators
- typedef ModularRandIter<Self_t> RandIter;
- typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
- template< class Random > Element& random(Random& g, Element& r) const
- { return init(r, g()); }
- template< class Random > Element& nonzerorandom(Random& g, Element& a) const
- { while (isZero(init(a, g())))
- ;
- return a; }
-
- // --- IO methods
- std::istream& read (std::istream& s);
- std::ostream& write(std::ostream& s) const;
- std::istream& read (std::istream& s, Element& a) const;
- std::ostream& write(std::ostream& s, const Element a) const;
-
- protected:
- // -- data representation of the domain:
- Residu_t _p;
- };
+ /*! @brief This class implement the standard arithmetic with Modulo Elements.
+ * - The representation of an integer a in Zpz is the value a % p
+ * - m max is 3037000499
+ * - p max is 3037000493
+ * .
+ */
+ template<typename COMP>
+ class Modular<int64_t, COMP> : public virtual FiniteFieldInterface<int64_t>
+ {
+ public:
+
+ // ----- Exported Types and constantes
+ using Self_t = Modular<int64_t, COMP>;
+ using Residu_t = uint64_t;
+ using Compute_t = typename std::make_unsigned<COMP>::type;
+ enum { size_rep = sizeof(Residu_t) };
+
+ // ----- Representation of vector of the Element
+ typedef Element* Array;
+ typedef const Element* constArray;
+
+ // ----- Constantes
+ const Element zero;
+ const Element one;
+ const Element mOne;
+
+ // ----- Constructors
+ Modular()
+ : zero(static_cast<Element>(0))
+ , one(static_cast<Element>(1))
+ , mOne(static_cast<Element>(-1))
+ , _p(static_cast<Residu_t>(0))
+ , _bitsizep(0) {}
+
+ Modular(const Residu_t p)
+ : zero(static_cast<Element>(0))
+ , one(static_cast<Element>(1))
+ , mOne(static_cast<Element>(p-1))
+ , _p(static_cast<Residu_t>(p))
+ , _bitsizep(0)
+ {
+ assert(_p >= minCardinality());
+ assert(_p <= maxCardinality());
+ Residu_t __p = _p;
+ while (__p != 0) {
+ _bitsizep++;
+ __p >>= 1;
+ }
+ }
+
+ Modular(const Self_t& F)
+ : zero(F.zero), one(F.one), mOne(F.mOne), _p(F._p), _bitsizep(F._bitsizep) {}
+
+ // ----- Accessors
+ inline Element minElement() const override { return zero; }
+ inline Element maxElement() const override { return mOne; }
+
+ // ----- Access to the modulus
+ inline Residu_t residu() const { return _p; }
+ inline Residu_t size() const { return _p; }
+ inline Residu_t characteristic() const { return _p; }
+ inline Residu_t cardinality() const { return _p; }
+ template<class T> inline T& characteristic(T& p) const { return p = _p; }
+ template<class T> inline T& cardinality(T& p) const { return p = _p; }
+ static inline Residu_t maxCardinality();
+ static inline Residu_t minCardinality() { return 2_ui64; }
+
+ // ----- Checkers
+ inline bool isZero(const Element& a) const override { return a == zero; }
+ inline bool isOne (const Element& a) const override { return a == one; }
+ inline bool isMOne(const Element& a) const override { return a == mOne; }
+ inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
+ inline size_t length(const Element a) const { return size_rep; }
+
+ // ----- Ring-wise operators
+ inline bool operator==(const Self_t& F) const { return _p == F._p; }
+ inline bool operator!=(const Self_t& F) const { return _p != F._p; }
+ inline Self_t& operator=(const Self_t& F)
+ {
+ F.assign(const_cast<Element&>(one), F.one);
+ F.assign(const_cast<Element&>(zero), F.zero);
+ F.assign(const_cast<Element&>(mOne), F.mOne);
+ _p = F._p;
+ _bitsizep = F._bitsizep;
+ return *this;
+ }
+
+ // ----- Initialisation
+ Element& init (Element& x) const;
+ Element& init (Element& x, const Integer& y) const;
+ template<typename T> Element& init(Element& r, const T& a) const
+ { r = Caster<Element>(a); return reduce(r); }
+ void init(const size_t, Array a, constArray b) const;
+
+ Element& assign (Element& x, const Element& y) const;
+ void assign(const size_t sz, Array r, constArray a ) const;
+
+ // ----- Convert and reduce
+ template<typename T> T& convert(T& r, const Element& a) const
+ { return r = static_cast<T>(a); }
+
+ Element& reduce (Element& x, const Element& y) const;
+ Element& reduce (Element& x) const;
+
+ // ----- Classic arithmetic
+ Element& mul(Element& r, const Element& a, const Element& b) const override;
+ Element& div(Element& r, const Element& a, const Element& b) const override;
+ Element& add(Element& r, const Element& a, const Element& b) const override;
+ Element& sub(Element& r, const Element& a, const Element& b) const override;
+ Element& neg(Element& r, const Element& a) const override;
+ Element& inv(Element& r, const Element& a) const override;
+
+ Element& mulin(Element& r, const Element& a) const override;
+ Element& divin(Element& r, const Element& a) const override;
+ Element& addin(Element& r, const Element& a) const override;
+ Element& subin(Element& r, const Element& a) const override;
+ Element& negin(Element& r) const override;
+ Element& invin(Element& r) const override;
+
+ // Functions defined in modular-mulprecomp
+ //
+ // void precomp_p (Compute_t& invp) const
+ // Element& mul_precomp_p (Element& r, const Element& a, const Element& b, const Compute_t& invp) const
+ //
+ // void precomp_b (Compute_t& invb, const Element& b) const
+ // void precomp_b (Compute_t& invb, const Element& b, const Compute_t& invp) const
+ // Element& mul_precomp_b (Element& r, const Element& a, const Element& b, const Compute_t& invb) const
+
+#include "modular-mulprecomp.inl"
+
+ // -- axpy: r <- a * x + y
+ // -- axpyin: r <- a * x + r
+ Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& axpyin(Element& r, const Element& a, const Element& x) const override;
+
+ // -- axmy: r <- a * x - y
+ // -- axmyin: r <- a * x - r
+ Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& axmyin(Element& r, const Element& a, const Element& x) const override;
+
+ // -- maxpy: r <- y - a * x
+ // -- maxpyin: r <- r - a * x
+ Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& maxpyin(Element& r, const Element& a, const Element& x) const override;
+
+ // ----- Classic arithmetic on arrays
+ void mul(const size_t sz, Array r, constArray a, constArray b) const;
+ void mul(const size_t sz, Array r, constArray a, Element b) const;
+ void div(const size_t sz, Array r, constArray a, constArray b) const;
+ void div(const size_t sz, Array r, constArray a, Element b) const;
+ void add(const size_t sz, Array r, constArray a, constArray b) const;
+ void add(const size_t sz, Array r, constArray a, Element b) const;
+ void sub(const size_t sz, Array r, constArray a, constArray b) const;
+ void sub(const size_t sz, Array r, constArray a, Element b) const;
+ void neg(const size_t sz, Array r, constArray a) const;
+ void inv(const size_t sz, Array r, constArray a) const;
+
+ void axpy (const size_t sz, Array r, constArray a, constArray x, constArray c) const;
+ void axpyin (const size_t sz, Array r, constArray a, constArray x) const;
+ void axmy (const size_t sz, Array r, constArray a, constArray x, constArray c) const;
+ void maxpyin (const size_t sz, Array r, constArray a, constArray x) const;
+
+ // <- \sum_i a[i], return 1 if a.size() ==0,
+ Element& reduceadd ( Element& r, const size_t sz, constArray a ) const;
+
+ // <- \prod_i a[i], return 1 if a.size() ==0,
+ Element& reducemul ( Element& r, const size_t sz, constArray a ) const;
+
+ // <- \sum_i a[i] * b[i]
+ Element& dotprod ( Element& r, const size_t sz, constArray a, constArray b ) const;
+ Element& dotprod ( Element& r, const int bound, const size_t sz, constArray a, constArray b ) const;
+
+ // a -> r: uint32_t to double
+ void i2d ( const size_t sz, double* r, constArray a ) const;
+
+ // a -> r % p: double to uint32_t % p
+ void d2i ( const size_t sz, Array r, const double* a ) const;
+
+ // ----- Random generators
+ typedef ModularRandIter<Self_t> RandIter;
+ typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
+ template< class Random > Element& random(Random& g, Element& r) const
+ { return init(r, g()); }
+ template< class Random > Element& nonzerorandom(Random& g, Element& a) const
+ { while (isZero(init(a, g())))
+ ;
+ return a; }
+
+ // --- IO methods
+ std::istream& read (std::istream& s);
+ std::ostream& write(std::ostream& s) const;
+ std::istream& read (std::istream& s, Element& a) const;
+ std::ostream& write(std::ostream& s, const Element a) const;
+
+ protected:
+ // -- data representation of the domain:
+ Residu_t _p;
+ size_t _bitsizep;
+ };
} // namespace Givaro
#include "givaro/modular-int64.inl"
-#endif // __GIVARO_zpz64std_H
-// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
+#endif // __GIVARO_modular_uint64_H
diff --git a/src/kernel/ring/modular-int64.inl b/src/kernel/ring/modular-int64.inl
index 8375b6b..ae31362 100644
--- a/src/kernel/ring/modular-int64.inl
+++ b/src/kernel/ring/modular-int64.inl
@@ -21,11 +21,21 @@ namespace Givaro {
template<>
inline Modular<int64_t, uint64_t>::Residu_t
- Modular<int64_t, uint64_t>::maxCardinality() { return 2147483647u; } // 2^31 - 1
+ Modular<int64_t, uint64_t>::maxCardinality() { return 4294967295_ui64; } // 2^32 - 1
template<>
inline Modular<int64_t, int64_t>::Residu_t
- Modular<int64_t, int64_t>::maxCardinality() { return 2147483647u; }
+ Modular<int64_t, int64_t>::maxCardinality() { return 4294967295_ui64; }
+
+#ifdef __GIVARO_HAVE_INT128
+ template<>
+ inline Modular<int64_t, uint128_t>::Residu_t
+ Modular<int64_t, uint128_t>::maxCardinality() { return 9223372036854775807_ui64; } // 2^63 - 1
+
+ template<>
+ inline Modular<int64_t, int128_t>::Residu_t
+ Modular<int64_t, int128_t>::maxCardinality() { return 9223372036854775807_ui64; }
+#endif
// ------------------------
// ----- Classic arithmetic
@@ -513,7 +523,7 @@ namespace Givaro {
template<>
inline std::ostream& Modular<int64_t, int64_t>::write (std::ostream& s ) const
{
- return s << "Modular<int64_t, uint64_t> modulo " << residu();
+ return s << "Modular<int64_t, int64_t> modulo " << residu();
}
template<>
@@ -522,6 +532,20 @@ namespace Givaro {
return s << "Modular<int64_t, uint64_t> modulo " << residu();
}
+#ifdef __GIVARO_HAVE_INT128
+ template<>
+ inline std::ostream& Modular<int64_t, int128_t>::write (std::ostream& s ) const
+ {
+ return s << "Modular<int64_t, int128_t> modulo " << residu();
+ }
+
+ template<>
+ inline std::ostream& Modular<int64_t, uint128_t>::write (std::ostream& s ) const
+ {
+ return s << "Modular<int64_t, uint128_t> modulo " << residu();
+ }
+#endif
+
template<typename COMP>
inline std::istream& Modular<int64_t, COMP>::read (std::istream& s, Element& a) const
{
diff --git a/src/kernel/ring/modular-int8.h b/src/kernel/ring/modular-int8.h
index 1a4531b..0f2db8c 100644
--- a/src/kernel/ring/modular-int8.h
+++ b/src/kernel/ring/modular-int8.h
@@ -1,3 +1,5 @@
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
// ==========================================================================
// $Source: /var/lib/cvs/Givaro/src/kernel/zpz/givzpz16std.h,v $
// Copyright(c)'1994-2009 by The Givaro group
@@ -11,14 +13,13 @@
//
// Modified by Pascal Giorgi on 2002/02/13 (pascal.giorgi at ens-lyon.fr)
// Modified by Alexis Breust on 2015/01/06 (alexis.breust at imag.fr)
+// Modified by Romain Lebreton on 2016/06/10 (romain.lebreton at lirmm.fr)
-/*! @file givzpz16std.h
- * @ingroup zpz
- * @brief Arithmetic on Z/pZ, with p a prime number less than 2^14.
- * Modulo typedef is a signed long number. In case it was modified
- * then Bézout algorithm must be changed (coefficient can be negative).
+/*! @file ring/modular-int8.h
+ * @ingroup ring
+ * @brief representation of <code>Z/mZ</code> over \c int8_t .
*/
-
+
#ifndef __GIVARO_modular_int8_H
#define __GIVARO_modular_int8_H
@@ -32,195 +33,217 @@
namespace Givaro {
- /*! @brief This class implement the standard arithmetic with Modulo Elements.
- * - The representation of an integer a in Zpz is the value a % p
- * - m max is 32768
- * - p max is 32749
- * .
- */
- template<typename COMP>
- class Modular<int8_t, COMP> : public virtual FiniteFieldInterface<int8_t>
- {
- public:
-
- // ----- Exported Types and constantes
- using Self_t = Modular<int8_t, COMP>;
- using Residu_t = uint8_t;
- using Compute_t = typename std::make_unsigned<COMP>::type;
- enum { size_rep = sizeof(Residu_t) };
-
-
- // ----- Representation of vector of the Element
- typedef Element* Array;
- typedef const Element* constArray;
-
- // ----- Constantes
- const Element zero;
- const Element one;
- const Element mOne;
-
- // ----- Constructors
- Modular()
- : _p(static_cast<Residu_t>(0)) {}
-
- Modular(const Residu_t p)
- : zero(static_cast<Element>(0))
- , one(static_cast<Element>(1))
- , mOne(static_cast<Element>(p-1))
- , _p(static_cast<Residu_t>(p))
- {
- assert(_p >= minCardinality());
- assert(_p <= maxCardinality());
- }
-
- Modular(const Self_t& F)
- : zero(F.zero), one(F.one), mOne(F.mOne), _p(F._p) {}
-
- // ----- Accessors
- inline Element minElement() const override { return zero; }
- inline Element maxElement() const override { return mOne; }
-
- // ----- Access to the modulus
- inline Residu_t residu() const { return _p; }
- inline Residu_t size() const { return _p; }
- inline Residu_t characteristic() const { return _p; }
- inline Residu_t cardinality() const { return _p; }
- template<class T> inline T& characteristic(T& p) const { return p = _p; }
- template<class T> inline T& cardinality(T& p) const { return p = _p; }
- static inline Residu_t maxCardinality();
- static inline Residu_t minCardinality() { return 2; }
-
- // ----- Checkers
- inline bool isZero(const Element& a) const override { return a == zero; }
- inline bool isOne (const Element& a) const override { return a == one; }
- inline bool isMOne(const Element& a) const override { return a == mOne; }
- inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
- inline size_t length(const Element a) const { return size_rep; }
-
- // ----- Ring-wise operators
- inline bool operator==(const Self_t& F) const { return _p == F._p; }
- inline bool operator!=(const Self_t& F) const { return _p != F._p; }
- inline Self_t& operator=(const Self_t& F)
- {
- F.assign(const_cast<Element&>(one), F.one);
- F.assign(const_cast<Element&>(zero), F.zero);
- F.assign(const_cast<Element&>(mOne), F.mOne);
- _p = F._p;
- return *this;
- }
-
- // ----- Initialisation
- Element& init (Element& x) const;
- Element& init (Element& x, const float y) const;
- Element& init (Element& x, const double y) const;
- Element& init (Element& x, const int16_t y) const;
- Element& init (Element& x, const uint16_t y) const;
- Element& init (Element& x, const int32_t y) const;
- Element& init (Element& x, const uint32_t y) const;
- Element& init (Element& x, const int64_t y) const;
- Element& init (Element& x, const uint64_t y) const;
- Element& init (Element& x, const Integer& y) const;
- template<typename T> Element& init(Element& r, const T& a) const
- { r = Caster<Element>(a); return reduce(r); }
- void init(const size_t, Array a, constArray b) const;
-
- Element& assign (Element& x, const Element& y) const;
- void assign(const size_t sz, Array r, constArray a ) const;
-
- // ----- Convert and reduce
- template<typename T> T& convert(T& r, const Element& a) const
- { return r = static_cast<T>(a); }
-
- Element& reduce (Element& x, const Element& y) const;
- Element& reduce (Element& x) const;
-
- // ----- Classic arithmetic
- Element& mul(Element& r, const Element& a, const Element& b) const override;
- Element& div(Element& r, const Element& a, const Element& b) const override;
- Element& add(Element& r, const Element& a, const Element& b) const override;
- Element& sub(Element& r, const Element& a, const Element& b) const override;
- Element& neg(Element& r, const Element& a) const override;
- Element& inv(Element& r, const Element& a) const override;
-
- Element& mulin(Element& r, const Element& a) const override;
- Element& divin(Element& r, const Element& a) const override;
- Element& addin(Element& r, const Element& a) const override;
- Element& subin(Element& r, const Element& a) const override;
- Element& negin(Element& r) const override;
- Element& invin(Element& r) const override;
-
- // -- axpy: r <- a * x + y
- // -- axpyin: r <- a * x + r
- Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& axpyin(Element& r, const Element& a, const Element& x) const override;
-
- // -- axmy: r <- a * x - y
- // -- axmyin: r <- a * x - r
- Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& axmyin(Element& r, const Element& a, const Element& x) const override;
-
- // -- maxpy: r <- y - a * x
- // -- maxpyin: r <- r - a * x
- Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& maxpyin(Element& r, const Element& a, const Element& x) const override;
-
- // ----- Classic arithmetic on arrays
- void mul(const size_t sz, Array r, constArray a, constArray b) const;
- void mul(const size_t sz, Array r, constArray a, Element b) const;
- void div(const size_t sz, Array r, constArray a, constArray b) const;
- void div(const size_t sz, Array r, constArray a, Element b) const;
- void add(const size_t sz, Array r, constArray a, constArray b) const;
- void add(const size_t sz, Array r, constArray a, Element b) const;
- void sub(const size_t sz, Array r, constArray a, constArray b) const;
- void sub(const size_t sz, Array r, constArray a, Element b) const;
- void neg(const size_t sz, Array r, constArray a) const;
- void inv(const size_t sz, Array r, constArray a) const;
-
- void axpy (const size_t sz, Array r, constArray a, constArray x, constArray c) const;
- void axpyin (const size_t sz, Array r, constArray a, constArray x) const;
- void axmy (const size_t sz, Array r, constArray a, constArray x, constArray c) const;
- void maxpyin (const size_t sz, Array r, constArray a, constArray x) const;
-
- // <- \sum_i a[i], return 1 if a.size() ==0,
- Element& reduceadd ( Element& r, const size_t sz, constArray a ) const;
-
- // <- \prod_i a[i], return 1 if a.size() ==0,
- Element& reducemul ( Element& r, const size_t sz, constArray a ) const;
-
- // <- \sum_i a[i] * b[i]
- Element& dotprod ( Element& r, const size_t sz, constArray a, constArray b ) const;
- Element& dotprod ( Element& r, const int bound, const size_t sz, constArray a, constArray b ) const;
-
- // a -> r: uint32_t to double
- void i2d ( const size_t sz, double* r, constArray a ) const;
-
- // a -> r % p: double to uint32_t % p
- void d2i ( const size_t sz, Array r, const double* a ) const;
-
- // ----- Random generators
- typedef ModularRandIter<Self_t> RandIter;
- typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
- template< class Random > Element& random(Random& g, Element& r) const
- { return init(r, g()); }
- template< class Random > Element& nonzerorandom(Random& g, Element& a) const
- { while (isZero(init(a, g())))
- ;
- return a; }
-
- // --- IO methods
- std::istream& read (std::istream& s);
- std::ostream& write(std::ostream& s) const;
- std::istream& read (std::istream& s, Element& a) const;
- std::ostream& write(std::ostream& s, const Element a) const;
-
- protected:
- // -- data representation of the domain:
- Residu_t _p;
- };
+ /*! @brief This class implement the standard arithmetic with Modulo Elements.
+ * - The representation of an integer a in Zpz is the value a % p
+ * - m max is 32768
+ * - p max is 32749
+ * .
+ */
+ template<typename COMP>
+ class Modular<int8_t, COMP> : public virtual FiniteFieldInterface<int8_t>
+ {
+ public:
+
+ // ----- Exported Types and constantes
+ using Self_t = Modular<int8_t, COMP>;
+ using Residu_t = uint8_t;
+ using Compute_t = typename std::make_unsigned<COMP>::type;
+ enum { size_rep = sizeof(Residu_t) };
+
+
+ // ----- Representation of vector of the Element
+ typedef Element* Array;
+ typedef const Element* constArray;
+
+ // ----- Constantes
+ const Element zero;
+ const Element one;
+ const Element mOne;
+
+ // ----- Constructors
+ Modular()
+ : zero(static_cast<Element>(0))
+ , one(static_cast<Element>(1))
+ , mOne(static_cast<Element>(-1))
+ , _p(static_cast<Residu_t>(0))
+ , _bitsizep(0) {}
+
+ Modular(const Residu_t p)
+ : zero(static_cast<Element>(0))
+ , one(static_cast<Element>(1))
+ , mOne(static_cast<Element>(p-1))
+ , _p(static_cast<Residu_t>(p))
+ , _bitsizep(0)
+ {
+ assert(_p >= minCardinality());
+ assert(_p <= maxCardinality());
+ Residu_t __p = _p;
+ while (__p != 0) {
+ _bitsizep++;
+ __p >>= 1;
+ }
+ }
+
+ Modular(const Self_t& F)
+ : zero(F.zero), one(F.one), mOne(F.mOne), _p(F._p), _bitsizep(F._bitsizep) {}
+
+ // ----- Accessors
+ inline Element minElement() const override { return zero; }
+ inline Element maxElement() const override { return mOne; }
+
+ // ----- Access to the modulus
+ inline Residu_t residu() const { return _p; }
+ inline Residu_t size() const { return _p; }
+ inline Residu_t characteristic() const { return _p; }
+ inline Residu_t cardinality() const { return _p; }
+ template<class T> inline T& characteristic(T& p) const { return p = _p; }
+ template<class T> inline T& cardinality(T& p) const { return p = _p; }
+ static inline Residu_t maxCardinality();
+ static inline Residu_t minCardinality() { return 2; }
+
+ // ----- Checkers
+ inline bool isZero(const Element& a) const override { return a == zero; }
+ inline bool isOne (const Element& a) const override { return a == one; }
+ inline bool isMOne(const Element& a) const override { return a == mOne; }
+ inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
+ inline size_t length(const Element a) const { return size_rep; }
+
+ // ----- Ring-wise operators
+ inline bool operator==(const Self_t& F) const { return _p == F._p; }
+ inline bool operator!=(const Self_t& F) const { return _p != F._p; }
+ inline Self_t& operator=(const Self_t& F)
+ {
+ F.assign(const_cast<Element&>(one), F.one);
+ F.assign(const_cast<Element&>(zero), F.zero);
+ F.assign(const_cast<Element&>(mOne), F.mOne);
+ _p = F._p;
+ _bitsizep = F._bitsizep;
+ return *this;
+ }
+
+ // ----- Initialisation
+ Element& init (Element& x) const;
+ Element& init (Element& x, const float y) const;
+ Element& init (Element& x, const double y) const;
+ Element& init (Element& x, const int16_t y) const;
+ Element& init (Element& x, const uint16_t y) const;
+ Element& init (Element& x, const int32_t y) const;
+ Element& init (Element& x, const uint32_t y) const;
+ Element& init (Element& x, const int64_t y) const;
+ Element& init (Element& x, const uint64_t y) const;
+ Element& init (Element& x, const Integer& y) const;
+ template<typename T> Element& init(Element& r, const T& a) const
+ { r = Caster<Element>(a); return reduce(r); }
+ void init(const size_t, Array a, constArray b) const;
+
+ Element& assign (Element& x, const Element& y) const;
+ void assign(const size_t sz, Array r, constArray a ) const;
+
+ // ----- Convert and reduce
+ template<typename T> T& convert(T& r, const Element& a) const
+ { return r = static_cast<T>(a); }
+
+ Element& reduce (Element& x, const Element& y) const;
+ Element& reduce (Element& x) const;
+
+ // ----- Classic arithmetic
+ Element& mul(Element& r, const Element& a, const Element& b) const override;
+ Element& div(Element& r, const Element& a, const Element& b) const override;
+ Element& add(Element& r, const Element& a, const Element& b) const override;
+ Element& sub(Element& r, const Element& a, const Element& b) const override;
+ Element& neg(Element& r, const Element& a) const override;
+ Element& inv(Element& r, const Element& a) const override;
+
+ Element& mulin(Element& r, const Element& a) const override;
+ Element& divin(Element& r, const Element& a) const override;
+ Element& addin(Element& r, const Element& a) const override;
+ Element& subin(Element& r, const Element& a) const override;
+ Element& negin(Element& r) const override;
+ Element& invin(Element& r) const override;
+
+ // Functions defined in modular-mulprecomp
+ //
+ // void precomp_p (Compute_t& invp) const
+ // Element& mul_precomp_p (Element& r, const Element& a, const Element& b, const Compute_t& invp) const
+ //
+ // void precomp_b (Compute_t& invb, const Element& b) const
+ // void precomp_b (Compute_t& invb, const Element& b, const Compute_t& invp) const
+ // Element& mul_precomp_b (Element& r, const Element& a, const Element& b, const Compute_t& invb) const
+
+#include "modular-mulprecomp.inl"
+
+ // -- axpy: r <- a * x + y
+ // -- axpyin: r <- a * x + r
+ Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& axpyin(Element& r, const Element& a, const Element& x) const override;
+
+ // -- axmy: r <- a * x - y
+ // -- axmyin: r <- a * x - r
+ Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& axmyin(Element& r, const Element& a, const Element& x) const override;
+
+ // -- maxpy: r <- y - a * x
+ // -- maxpyin: r <- r - a * x
+ Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& maxpyin(Element& r, const Element& a, const Element& x) const override;
+
+ // ----- Classic arithmetic on arrays
+ void mul(const size_t sz, Array r, constArray a, constArray b) const;
+ void mul(const size_t sz, Array r, constArray a, Element b) const;
+ void div(const size_t sz, Array r, constArray a, constArray b) const;
+ void div(const size_t sz, Array r, constArray a, Element b) const;
+ void add(const size_t sz, Array r, constArray a, constArray b) const;
+ void add(const size_t sz, Array r, constArray a, Element b) const;
+ void sub(const size_t sz, Array r, constArray a, constArray b) const;
+ void sub(const size_t sz, Array r, constArray a, Element b) const;
+ void neg(const size_t sz, Array r, constArray a) const;
+ void inv(const size_t sz, Array r, constArray a) const;
+
+ void axpy (const size_t sz, Array r, constArray a, constArray x, constArray c) const;
+ void axpyin (const size_t sz, Array r, constArray a, constArray x) const;
+ void axmy (const size_t sz, Array r, constArray a, constArray x, constArray c) const;
+ void maxpyin (const size_t sz, Array r, constArray a, constArray x) const;
+
+ // <- \sum_i a[i], return 1 if a.size() ==0,
+ Element& reduceadd ( Element& r, const size_t sz, constArray a ) const;
+
+ // <- \prod_i a[i], return 1 if a.size() ==0,
+ Element& reducemul ( Element& r, const size_t sz, constArray a ) const;
+
+ // <- \sum_i a[i] * b[i]
+ Element& dotprod ( Element& r, const size_t sz, constArray a, constArray b ) const;
+ Element& dotprod ( Element& r, const int bound, const size_t sz, constArray a, constArray b ) const;
+
+ // a -> r: uint32_t to double
+ void i2d ( const size_t sz, double* r, constArray a ) const;
+
+ // a -> r % p: double to uint32_t % p
+ void d2i ( const size_t sz, Array r, const double* a ) const;
+
+ // ----- Random generators
+ typedef ModularRandIter<Self_t> RandIter;
+ typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
+ template< class Random > Element& random(Random& g, Element& r) const
+ { return init(r, g()); }
+ template< class Random > Element& nonzerorandom(Random& g, Element& a) const
+ { while (isZero(init(a, g())))
+ ;
+ return a; }
+
+ // --- IO methods
+ std::istream& read (std::istream& s);
+ std::ostream& write(std::ostream& s) const;
+ std::istream& read (std::istream& s, Element& a) const;
+ std::ostream& write(std::ostream& s, const Element a) const;
+
+ protected:
+ // -- data representation of the domain:
+ Residu_t _p;
+ size_t _bitsizep;
+ };
} // namespace Givaro
#include "givaro/modular-int8.inl"
#endif // __GIVARO_modular_int8_H
-// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
diff --git a/src/kernel/ring/modular-integer.h b/src/kernel/ring/modular-integer.h
index 9a43d5b..adfd769 100644
--- a/src/kernel/ring/modular-integer.h
+++ b/src/kernel/ring/modular-integer.h
@@ -36,7 +36,7 @@ namespace Givaro
{
public:
// ----- Exported Types and constantes
- typedef Modular<Integer> Self_t;
+ typedef Modular<Integer> Self_t;
typedef Integer Residu_t; // - type to store residue
enum { size_rep = sizeof(Residu_t) }; // - size of the storage type
@@ -167,10 +167,10 @@ namespace Givaro
/* Specialisation for Modular<integer> field*/
template <>
- class ModularRandIter<Modular<Integer, Integer> >
+ class ModularRandIter<Modular<Integer> >
{
public:
- typedef Modular<Integer, Integer> Ring;
+ typedef Modular<Integer> Ring;
typedef Ring::Element Element;
ModularRandIter(const Ring& R, const size_t& size = 0, const size_t& seed = 0)
diff --git a/src/kernel/ring/modular-inttype.h b/src/kernel/ring/modular-inttype.h
index 394ac8d..ac70282 100644
--- a/src/kernel/ring/modular-inttype.h
+++ b/src/kernel/ring/modular-inttype.h
@@ -36,7 +36,7 @@ namespace Givaro
{
public:
// ----- Exported Types and constantes
- typedef Modular<IntType> Self_t;
+ typedef Modular<IntType,COMP> Self_t;
typedef IntType Residu_t;
enum { size_rep = sizeof(Residu_t) };
using Element = typename FiniteFieldInterface<IntType>::Element;
@@ -58,7 +58,7 @@ namespace Givaro
, one(1)
, mOne(p-static_cast<Residu_t>(1)), _p(p)
{
- assert(_p >= minCardinality());
+ assert(_p >= minCardinality());
}
template<class IntConvType>
@@ -68,7 +68,7 @@ namespace Givaro
, mOne( Caster<Residu_t>(p-1) )
, _p( Caster<Residu_t>(p) )
{
- assert(_p >= minCardinality());
+ assert(_p >= minCardinality());
}
Modular(const Integer& p, const Integer& e=Integer::one)
@@ -77,7 +77,7 @@ namespace Givaro
, mOne( Caster<Residu_t>(p-1) )
, _p( Caster<Residu_t>(p) )
{
- assert(_p >= minCardinality());
+ assert(_p >= minCardinality());
}
Modular(const Self_t& F)
diff --git a/src/kernel/ring/modular-mulprecomp.inl b/src/kernel/ring/modular-mulprecomp.inl
new file mode 100644
index 0000000..330a048
--- /dev/null
+++ b/src/kernel/ring/modular-mulprecomp.inl
@@ -0,0 +1,95 @@
+#define SUITABLE_INTEGRAL_TYPES(ELEMENT, COMPUTE_T, s) \
+ (2 * sizeof(ELEMENT) >= sizeof(COMPUTE_T)) && (sizeof(COMPUTE_T) == s) \
+ && (std::is_integral<ELEMENT>::value) && (std::is_integral<COMPUTE_T>::value)
+
+#define SUITABLE_INTEGRAL_TYPES_AND_HALF_SIZE(ELEMENT, COMPUTE_T, s) \
+ (2 * sizeof(ELEMENT) == sizeof(COMPUTE_T)) && (sizeof(COMPUTE_T) == s) \
+ && (std::is_integral<ELEMENT>::value) && (std::is_integral<COMPUTE_T>::value)
+
+#define SUITABLE_INTEGRAL_TYPES_AND_FULL_SIZE(ELEMENT, COMPUTE_T, s) \
+ (sizeof(ELEMENT) == sizeof(COMPUTE_T)) && (sizeof(COMPUTE_T) == s) \
+ && (std::is_integral<ELEMENT>::value) && (std::is_integral<COMPUTE_T>::value)
+
+template<typename ELEMENT = Element, typename COMPUTE_T = Compute_t, int s = sizeof(COMPUTE_T)>
+inline
+typename std::enable_if<SUITABLE_INTEGRAL_TYPES (ELEMENT, COMPUTE_T, s)>::type
+precomp_p
+(Compute_t& invp) const
+{
+ assert( _bitsizep <= (4*s-2));
+ invp = (static_cast<Compute_t>(1) << (4*s + _bitsizep - 1)) / static_cast<Compute_t>(_p);
+}
+
+template<typename ELEMENT = Element, typename COMPUTE_T = Compute_t, int s = sizeof(COMPUTE_T)>
+inline
+typename std::enable_if<SUITABLE_INTEGRAL_TYPES (ELEMENT, COMPUTE_T, s), Element&>::type
+mul_precomp_p
+(Element& r, const Element& a, const Element& b, const Compute_t& invp) const
+{
+ Compute_t prod = static_cast<Compute_t>(static_cast<Residu_t>(a))*static_cast<Compute_t>(static_cast<Residu_t>(b));
+ Compute_t prodhi = (prod >> (_bitsizep - 2)); // Could fit into an Element but no use
+ Residu_t c = static_cast<Residu_t>((prodhi * invp) >> (4*s+1));
+ Residu_t rr = static_cast<Residu_t>(prod) - c * _p;
+ rr -= (rr >= _p)?_p:0;
+ return r = static_cast<Element>(rr);
+}
+
+template<typename ELEMENT = Element, typename COMPUTE_T = Compute_t, int s = sizeof(COMPUTE_T)>
+inline
+typename std::enable_if<SUITABLE_INTEGRAL_TYPES (ELEMENT, COMPUTE_T, s)>::type
+precomp_b
+(Compute_t& invb, const Element& b) const{
+ assert( _bitsizep <= (4*s-1));
+ invb = (static_cast<Compute_t>(1) << (4*s)) * static_cast<Compute_t>(static_cast<Residu_t>(b)) / static_cast<Compute_t>(_p);
+}
+
+template<typename ELEMENT = Element, typename COMPUTE_T = Compute_t, int s = sizeof(COMPUTE_T)>
+inline
+typename std::enable_if<SUITABLE_INTEGRAL_TYPES_AND_HALF_SIZE (ELEMENT, COMPUTE_T, s)>::type
+precomp_b
+(Compute_t& invb, const Element& b, const Compute_t& invp) const{
+ assert( _bitsizep <= (4*s-2));
+
+ invb = (static_cast<Compute_t>(static_cast<Residu_t>(b)) * invp) >> (_bitsizep - 1);
+ Residu_t r = - static_cast<Residu_t>(invb) * _p;
+
+ // Quotient can only be off by two since (2^n * b) / 2^(r-1) is an exact division
+ bool flag = (r >= (_p));
+ invb += (flag)?1:0;
+ r -= (flag)?_p:0;
+ invb += (r >= (_p))?1:0;
+}
+
+template<typename ELEMENT = Element, typename COMPUTE_T = Compute_t, int s = sizeof(COMPUTE_T)>
+inline
+typename std::enable_if<SUITABLE_INTEGRAL_TYPES_AND_FULL_SIZE (ELEMENT, COMPUTE_T, s)>::type
+precomp_b
+(Compute_t& invb, const Element& b, const Element& invp) const{
+ assert( _bitsizep <= (4*s-2));
+ // Prod = 2^(4*s) * b
+ invb = (static_cast<Compute_t>(b) * static_cast<Compute_t>(static_cast<Residu_t>(invp))) >> (_bitsizep - 1);
+
+ // The problem is that Residu_t should be unsigned half size
+ // Warning : left shift count >= width of type when Element is half size of Compute_t
+ // and then << does nothing instead of mapping to zero
+ Residu_t r = (static_cast<Residu_t>(b) << (4*s)) - static_cast<Residu_t>(invb) * _p;
+
+ // Quotient can only be off by two since (2^n * b) / 2^(r-1) is an exact division
+ bool flag = (r >= (_p));
+ invb += (flag)?1:0;
+ r -= (flag)?_p:0;
+ invb += (r >= (_p))?1:0;
+}
+
+template<typename ELEMENT = Element, typename COMPUTE_T = Compute_t, int s = sizeof(COMPUTE_T)>
+inline
+typename std::enable_if<SUITABLE_INTEGRAL_TYPES (ELEMENT, COMPUTE_T, s), Element&>::type
+mul_precomp_b
+(Element& r, const Element& a, const Element& b, const Compute_t& invb) const {
+ Residu_t c = static_cast<Residu_t> ((static_cast<Compute_t>(a) * invb) >> (4*s));
+ Residu_t rr = static_cast<Residu_t>(a) * static_cast<Residu_t>(b) - c * _p;
+ rr -= (rr >= _p)?_p:0;
+ return r = static_cast<Residu_t>(rr);
+}
+
+#undef SUITABLE_INTEGRAL_TYPES
diff --git a/src/kernel/ring/modular-ruint.h b/src/kernel/ring/modular-ruint.h
index a13b0ae..5092b1a 100644
--- a/src/kernel/ring/modular-ruint.h
+++ b/src/kernel/ring/modular-ruint.h
@@ -175,9 +175,7 @@ namespace Givaro
template< class Random > Element& random(Random& g, Element& r) const
{ RecInt::rand(r); mod_n(r, _p); return r; }
template< class Random > Element& nonzerorandom(Random& g, Element& a) const
- { while (isZero(random(g, a)))
- ;
- return a; }
+ { while (isZero(random(g, a))) { } return a; }
// --- IO methods
std::istream& read (std::istream& s);
diff --git a/src/kernel/ring/modular-uint16.h b/src/kernel/ring/modular-uint16.h
index cb05319..e4dc103 100644
--- a/src/kernel/ring/modular-uint16.h
+++ b/src/kernel/ring/modular-uint16.h
@@ -1,3 +1,5 @@
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
// ==========================================================================
// Copyright(c)'1994-2015 by The Givaro group
// This file is part of Givaro.
@@ -8,9 +10,9 @@
// A. Breust (adapted)
// ==========================================================================
-/*! @file field/modular-uint64.h
- * @ingroup field
- * @brief representation of <code>Z/mZ</code> over \c uint64_t .
+/*! @file ring/modular-uint16.h
+ * @ingroup ring
+ * @brief representation of <code>Z/mZ</code> over \c uint16_t .
*/
#ifndef __GIVARO_modular_uint16_H
@@ -24,163 +26,187 @@
namespace Givaro {
- /** \brief Specialization of Modular to uint64_t element type with efficient dot product.
- *
- * Efficient element operations for dot product, mul, axpy, by using floating point
- * inverse of modulus (borrowed from NTL) and some use of non-normalized intermediate values.
- *
- * For some uses this is the most efficient field for primes in the range from half word
- * to 2^30.
- *
- * Requires: Modulus < 2^30.
- * Intended use: 2^15 < prime modulus < 2^30.
- * \ingroup field
- */
-
- template <typename COMP>
- class Modular<uint16_t, COMP> : public virtual FiniteFieldInterface<uint16_t>
- {
- public:
- // ----- Exported Types and constantes
- using Self_t = Modular<uint16_t, COMP>;
- using Compute_t = typename std::make_unsigned<COMP>::type;
- using Residu_t = uint16_t;
- enum { size_rep = sizeof(Residu_t) };
-
- // ----- Constantes
- const Element zero;
- const Element one;
- const Element mOne;
-
- // ----- Constructors
- Modular()
- : zero(Element(0)), one(Element(1)), mOne(Element(0)), _p(Element(0)) {}
-
- Modular(Residu_t p)
- : zero(Element(0)), one(Element(1)), mOne(Element(p-1)), _p(p)
- {
- assert(_p >= minCardinality());
- assert(_p <= maxCardinality());
- }
-
- Modular(const Self_t& F)
- : zero(F.zero), one(F.one), mOne(F.mOne), _p(F._p) {}
-
- // ----- Accessors
- inline Element minElement() const override { return zero; }
- inline Element maxElement() const override { return mOne; }
-
- // ----- Access to the modulus
- inline Residu_t residu() const { return _p; }
- inline Residu_t size() const { return _p; }
- inline Residu_t characteristic() const { return _p; }
- inline Residu_t cardinality() const { return _p; }
- template<class T> inline T& characteristic(T& p) const { return p = _p; }
- template<class T> inline T& cardinality(T& p) const { return p = _p; }
- static inline Residu_t maxCardinality();
- static inline Residu_t minCardinality() { return 2; }
-
- // ----- Checkers
- inline bool isZero(const Element& a) const override { return a == zero; }
- inline bool isOne (const Element& a) const override { return a == one; }
- inline bool isMOne(const Element& a) const override { return a == mOne; }
- inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
- inline size_t length(const Element a) const { return size_rep; }
-
- // ----- Ring-wise operators
- bool operator==(const Self_t& F) const { return _p == F._p; }
- bool operator!=(const Self_t& F) const { return _p != F._p; }
- Self_t& operator=(const Self_t& F)
- {
- F.assign(const_cast<Element&>(one), F.one);
- F.assign(const_cast<Element&>(zero), F.zero);
- F.assign(const_cast<Element&>(mOne), F.mOne);
- _p = F._p;
- return *this;
- }
-
- // ----- Initialisation
- Element& init (Element& x) const
- { return x = 0; }
- Element& init (Element& x, const float a) const;
- Element& init (Element& x, const double a) const;
- Element& init (Element& x, const int32_t a) const;
- Element& init (Element& x, const uint32_t a) const;
- Element& init (Element& x, const int64_t a) const;
- Element& init (Element& x, const uint64_t a) const;
- Element& init (Element& x, const Integer& a) const;
- template<typename T> Element& init(Element& r, const T& a) const
- {
- reduce(r, Caster<Element>((a < 0)? -a : a));
- if (a < 0) negin(r);
- return r;
- }
-
- Element& assign (Element& x, const Element& y) const
- { return x = y; }
-
- // ----- Convert and reduce
- template<typename T> T& convert(T& r, const Element& a) const
- { return r = Caster<T>(a); }
-
- Element& reduce (Element& x, const Element& y) const
- { x = y % _p; return x; }
- Element& reduce (Element& x) const
- { x %= _p; return x; }
-
- // ----- Classic arithmetic
- Element& mul(Element& r, const Element& a, const Element& b) const override;
- Element& div(Element& r, const Element& a, const Element& b) const override;
- Element& add(Element& r, const Element& a, const Element& b) const override;
- Element& sub(Element& r, const Element& a, const Element& b) const override;
- Element& neg(Element& r, const Element& a) const override;
- Element& inv(Element& r, const Element& a) const override;
-
- Element& mulin(Element& r, const Element& a) const override;
- Element& divin(Element& r, const Element& a) const override;
- Element& addin(Element& r, const Element& a) const override;
- Element& subin(Element& r, const Element& a) const override;
- Element& negin(Element& r) const override;
- Element& invin(Element& r) const override;
-
- // -- axpy: r <- a * x + y
- // -- axpyin: r <- a * x + r
- Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& axpyin(Element& r, const Element& a, const Element& x) const override;
-
- // -- axmy: r <- a * x - y
- // -- axmyin: r <- a * x - r
- Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& axmyin(Element& r, const Element& a, const Element& x) const override;
-
- // -- maxpy: r <- y - a * x
- // -- maxpyin: r <- r - a * x
- Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& maxpyin(Element& r, const Element& a, const Element& x) const override;
-
- // ----- Random generators
- typedef ModularRandIter<Self_t> RandIter;
- typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
- template< class Random > Element& random(Random& g, Element& r) const
- { return init(r, g()); }
- template< class Random > Element& nonzerorandom(Random& g, Element& a) const
- { while (isZero(init(a, g())))
- ;
- return a; }
-
- // --- IO methods
- std::ostream& write(std::ostream& s) const;
- std::istream& read (std::istream& s, Element& a) const;
- std::ostream& write(std::ostream& s, const Element& a) const;
-
- protected:
- // -- data representation of the domain:
- Residu_t _p;
- };
+ /** \brief Specialization of Modular to uint64_t element type with efficient dot product.
+ *
+ * Efficient element operations for dot product, mul, axpy, by using floating point
+ * inverse of modulus (borrowed from NTL) and some use of non-normalized intermediate values.
+ *
+ * For some uses this is the most efficient field for primes in the range from half word
+ * to 2^30.
+ *
+ * Requires: Modulus < 2^30.
+ * Intended use: 2^15 < prime modulus < 2^30.
+ * \ingroup field
+ */
+
+ template <typename COMP>
+ class Modular<uint16_t, COMP> : public virtual FiniteFieldInterface<uint16_t>
+ {
+ public:
+ // ----- Exported Types and constantes
+ using Self_t = Modular<uint16_t, COMP>;
+ using Compute_t = typename std::make_unsigned<COMP>::type;
+ using Residu_t = uint16_t;
+ enum { size_rep = sizeof(Residu_t) };
+
+ // ----- Constantes
+ const Element zero;
+ const Element one;
+ const Element mOne;
+ // ----- Constructors
+ Modular()
+ : zero(static_cast<Element>(0))
+ , one(static_cast<Element>(1))
+ , mOne(static_cast<Element>(-1))
+ , _p(static_cast<Residu_t>(0))
+ , _bitsizep(0) {}
+
+ Modular(const Residu_t p)
+ : zero(static_cast<Element>(0))
+ , one(static_cast<Element>(1))
+ , mOne(static_cast<Element>(p-1))
+ , _p(static_cast<Residu_t>(p))
+ , _bitsizep(0)
+ {
+ assert(_p >= minCardinality());
+ assert(_p <= maxCardinality());
+ Residu_t __p = _p;
+ while (__p != 0) {
+ _bitsizep++;
+ __p >>= 1;
+ }
+ }
+
+ Modular(const Self_t& F)
+ : zero(F.zero), one(F.one), mOne(F.mOne), _p(F._p), _bitsizep(F._bitsizep) {}
+
+ // ----- Accessors
+ inline Element minElement() const override { return zero; }
+ inline Element maxElement() const override { return mOne; }
+
+ // ----- Access to the modulus
+ inline Residu_t residu() const { return _p; }
+ inline Residu_t size() const { return _p; }
+ inline Residu_t characteristic() const { return _p; }
+ inline Residu_t cardinality() const { return _p; }
+ template<class T> inline T& characteristic(T& p) const { return p = _p; }
+ template<class T> inline T& cardinality(T& p) const { return p = _p; }
+ static inline Residu_t maxCardinality();
+ static inline Residu_t minCardinality() { return 2; }
+
+ // ----- Checkers
+ inline bool isZero(const Element& a) const override { return a == zero; }
+ inline bool isOne (const Element& a) const override { return a == one; }
+ inline bool isMOne(const Element& a) const override { return a == mOne; }
+ inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
+ inline size_t length(const Element a) const { return size_rep; }
+
+ // ----- Ring-wise operators
+ bool operator==(const Self_t& F) const { return _p == F._p; }
+ bool operator!=(const Self_t& F) const { return _p != F._p; }
+ Self_t& operator=(const Self_t& F)
+ {
+ F.assign(const_cast<Element&>(one), F.one);
+ F.assign(const_cast<Element&>(zero), F.zero);
+ F.assign(const_cast<Element&>(mOne), F.mOne);
+ _p = F._p;
+ _bitsizep = F._bitsizep;
+ return *this;
+ }
+
+ // ----- Initialisation
+ Element& init (Element& x) const
+ { return x = 0; }
+ Element& init (Element& x, const float a) const;
+ Element& init (Element& x, const double a) const;
+ Element& init (Element& x, const int32_t a) const;
+ Element& init (Element& x, const uint32_t a) const;
+ Element& init (Element& x, const int64_t a) const;
+ Element& init (Element& x, const uint64_t a) const;
+ Element& init (Element& x, const Integer& a) const;
+ template<typename T> Element& init(Element& r, const T& a) const
+ {
+ reduce(r, Caster<Element>((a < 0)? -a : a));
+ if (a < 0) negin(r);
+ return r;
+ }
+
+ Element& assign (Element& x, const Element& y) const
+ { return x = y; }
+
+ // ----- Convert and reduce
+ template<typename T> T& convert(T& r, const Element& a) const
+ { return r = Caster<T>(a); }
+
+ Element& reduce (Element& x, const Element& y) const
+ { x = y % _p; return x; }
+ Element& reduce (Element& x) const
+ { x %= _p; return x; }
+
+ // ----- Classic arithmetic
+ Element& mul(Element& r, const Element& a, const Element& b) const override;
+ Element& div(Element& r, const Element& a, const Element& b) const override;
+ Element& add(Element& r, const Element& a, const Element& b) const override;
+ Element& sub(Element& r, const Element& a, const Element& b) const override;
+ Element& neg(Element& r, const Element& a) const override;
+ Element& inv(Element& r, const Element& a) const override;
+
+ Element& mulin(Element& r, const Element& a) const override;
+ Element& divin(Element& r, const Element& a) const override;
+ Element& addin(Element& r, const Element& a) const override;
+ Element& subin(Element& r, const Element& a) const override;
+ Element& negin(Element& r) const override;
+ Element& invin(Element& r) const override;
+
+ // Functions defined in modular-mulprecomp
+ //
+ // void precomp_p (Compute_t& invp) const
+ // Element& mul_precomp_p (Element& r, const Element& a, const Element& b, const Compute_t& invp) const
+ //
+ // void precomp_b (Compute_t& invb, const Element& b) const
+ // void precomp_b (Compute_t& invb, const Element& b, const Compute_t& invp) const
+ // Element& mul_precomp_b (Element& r, const Element& a, const Element& b, const Compute_t& invb) const
+
+#include "modular-mulprecomp.inl"
+
+ // -- axpy: r <- a * x + y
+ // -- axpyin: r <- a * x + r
+ Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& axpyin(Element& r, const Element& a, const Element& x) const override;
+
+ // -- axmy: r <- a * x - y
+ // -- axmyin: r <- a * x - r
+ Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& axmyin(Element& r, const Element& a, const Element& x) const override;
+
+ // -- maxpy: r <- y - a * x
+ // -- maxpyin: r <- r - a * x
+ Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& maxpyin(Element& r, const Element& a, const Element& x) const override;
+
+ // ----- Random generators
+ typedef ModularRandIter<Self_t> RandIter;
+ typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
+ template< class Random > Element& random(Random& g, Element& r) const
+ { return init(r, g()); }
+ template< class Random > Element& nonzerorandom(Random& g, Element& a) const
+ { while (isZero(init(a, g())))
+ ;
+ return a; }
+
+ // --- IO methods
+ std::ostream& write(std::ostream& s) const;
+ std::istream& read (std::istream& s, Element& a) const;
+ std::ostream& write(std::ostream& s, const Element& a) const;
+
+ protected:
+ // -- data representation of the domain:
+ Residu_t _p;
+ size_t _bitsizep;
+ };
}
#include "givaro/modular-uint16.inl"
#endif //__GIVARO_modular_uint16_H
-
diff --git a/src/kernel/ring/modular-uint32.h b/src/kernel/ring/modular-uint32.h
index 946b8c0..f56c26a 100644
--- a/src/kernel/ring/modular-uint32.h
+++ b/src/kernel/ring/modular-uint32.h
@@ -1,3 +1,5 @@
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
// ==========================================================================
// $Source: /var/lib/cvs/Givaro/src/kernel/zpz/givzpz32uns.h,v $
// Copyright(c)'1994-2009 by The Givaro group
@@ -10,15 +12,15 @@
// ==========================================================================
//
// Modified by Pascal Giorgi on 2002/02/13 (pascal.giorgi at ens-lyon.fr)
+// Modified by Romain Lebreton on 2016/06/10 (romain.lebreton at lirmm.fr)
-/*! @file givzpz32uns.h
- * @ingroup zpz
- * @brief Arithmetic on Z/pZ, with p a prime number less than 2^32.
- * Modulo typedef is a signed long number. In case it was modified
- * then Bézout algorithm must be changed (coefficient can be negative).
+/*! @file ring/modular-uint32.h
+ * @ingroup ring
+ * @brief representation of <code>Z/mZ</code> over \c uint32_t .
*/
-#ifndef __GIVARO_zpz32unsigned_H
-#define __GIVARO_zpz32unsigned_H
+
+#ifndef __GIVARO_modular_uint32_H
+#define __GIVARO_modular_uint32_H
#include "givaro/givinteger.h"
#include "givaro/ring-interface.h"
@@ -28,153 +30,180 @@
namespace Givaro
{
- /*! @brief This class implement the standard arithmetic with Modulo Elements.
- * - The representation of an integer a in Zpz is the value a % p
- * - m max is 65536
- * - p max is 65521
- * .
- */
- template<typename COMP>
- class Modular<uint32_t, COMP> : public virtual FiniteFieldInterface<uint32_t>
- {
- public:
- // ----- Exported Types and constantes
- using Self_t = Modular<uint32_t, COMP>;
- using Compute_t = typename std::make_unsigned<COMP>::type;
- using Residu_t = uint32_t;
- enum { size_rep = sizeof(Residu_t) };
-
- // ----- Constantes
- const Element zero;
- const Element one;
- const Element mOne;
-
- // ----- Constructors
- Modular() : zero(0), one(1), mOne(0), _p(0), _dp(0.0) {}
-
- Modular(Residu_t p)
- : zero(0), one(1), mOne((Element)p-1), _p(p), _dp((double)p)
- {
- assert(_p >= minCardinality());
- assert(_p <= maxCardinality());
- }
-
- Modular(const Self_t& F)
- : zero(F.zero), one(F.one), mOne(F.mOne), _p(F._p), _dp(F._dp) {}
-
- // ----- Accessors
- inline Element minElement() const override { return zero; }
- inline Element maxElement() const override { return mOne; }
-
- // ----- Access to the modulus
- inline Residu_t residu() const { return _p; }
- inline Residu_t size() const { return _p; }
- inline Residu_t characteristic() const { return _p; }
- inline Residu_t cardinality() const { return _p; }
- template<class T> inline T& characteristic(T& p) const { return p = _p; }
- template<class T> inline T& cardinality(T& p) const { return p = _p; }
- static inline Residu_t maxCardinality();
- static inline Residu_t minCardinality() { return 2; }
-
- // ----- Checkers
- inline bool isZero(const Element& a) const override { return a == zero; }
- inline bool isOne (const Element& a) const override { return a == one; }
- inline bool isMOne(const Element& a) const override { return a == mOne; }
- inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
- inline size_t length(const Element a) const { return size_rep; }
-
- // ----- Ring-wise operators
- bool operator==(const Self_t& F) const { return _p == F._p; }
- bool operator!=(const Self_t& F) const { return _p != F._p; }
- Self_t& operator=(const Self_t& F)
- {
- F.assign(const_cast<Element&>(one), F.one);
- F.assign(const_cast<Element&>(zero), F.zero);
- F.assign(const_cast<Element&>(mOne), F.mOne);
- _p = F._p;
- _dp = F._dp;
- return *this;
- }
-
- // ----- Initialisation
- Element& init (Element& x) const
- { return x = 0; }
- Element& init (Element& x, const double a) const;
- Element& init (Element& x, const int64_t a) const;
- Element& init (Element& x, const uint64_t a) const;
- Element& init (Element& x, const Integer& a) const;
- template<typename T> Element& init(Element& r, const T& a) const
- {
- reduce( Caster<Element,T>(r, (a < 0)? -a : a) );
- return (a<0 ? negin(r) : r) ;
- }
-
- Element& assign (Element& x, const Element& y) const
- { return x = y; }
-
- // ----- Convert and reduce
- template<typename T> T& convert(T& r, const Element& a) const
- { return Caster<T,Element>(r,a); }
-
- Element& reduce (Element& x, const Element& y) const
- { return x = y % _p;}
- Element& reduce (Element& x) const
- { return x %= _p; }
-
- // ----- Classic arithmetic
- Element& mul(Element& r, const Element& a, const Element& b) const override;
- Element& div(Element& r, const Element& a, const Element& b) const override;
- Element& add(Element& r, const Element& a, const Element& b) const override;
- Element& sub(Element& r, const Element& a, const Element& b) const override;
- Element& neg(Element& r, const Element& a) const override;
- Element& inv(Element& r, const Element& a) const override;
-
- Element& mulin(Element& r, const Element& a) const override;
- Element& divin(Element& r, const Element& a) const override;
- Element& addin(Element& r, const Element& a) const override;
- Element& subin(Element& r, const Element& a) const override;
- Element& negin(Element& r) const override;
- Element& invin(Element& r) const override;
-
- // -- axpy: r <- a * x + y
- // -- axpyin: r <- a * x + r
- Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& axpyin(Element& r, const Element& a, const Element& x) const override;
-
- // -- axmy: r <- a * x - y
- // -- axmyin: r <- a * x - r
- Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& axmyin(Element& r, const Element& a, const Element& x) const override;
-
- // -- maxpy: r <- y - a * x
- // -- maxpyin: r <- r - a * x
- Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& maxpyin(Element& r, const Element& a, const Element& x) const override;
-
- // ----- Random generators
- typedef ModularRandIter<Self_t> RandIter;
- typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
- template< class Random > Element& random(Random& g, Element& r) const
- { return init(r, g()); }
- template< class Random > Element& nonzerorandom(Random& g, Element& a) const
- { while (isZero(init(a, g())))
- ;
- return a; }
-
- // --- IO methods
- std::ostream& write(std::ostream& s) const;
- std::istream& read (std::istream& s, Element& a) const;
- std::ostream& write(std::ostream& s, const Element a) const;
-
- protected:
-
- Residu_t _p;
- double _dp;
-
- };
+ /*! @brief This class implement the standard arithmetic with Modulo Elements.
+ * - The representation of an integer a in Zpz is the value a % p
+ * - m max is 65536
+ * - p max is 65521
+ * .
+ */
+ template<typename COMP>
+ class Modular<uint32_t, COMP> : public virtual FiniteFieldInterface<uint32_t>
+ {
+ public:
+ // ----- Exported Types and constantes
+ using Self_t = Modular<uint32_t, COMP>;
+ using Compute_t = typename std::make_unsigned<COMP>::type;
+ using Residu_t = uint32_t;
+ enum { size_rep = sizeof(Residu_t) };
+
+ // ----- Constantes
+ const Element zero;
+ const Element one;
+ const Element mOne;
+
+ // ----- Constructors
+ // ----- Constructors
+ Modular()
+ : zero(static_cast<Element>(0))
+ , one(static_cast<Element>(1))
+ , mOne(static_cast<Element>(-1))
+ , _p(static_cast<Residu_t>(0))
+ , _bitsizep(0) {}
+
+ Modular(const Residu_t p)
+ : zero(static_cast<Element>(0))
+ , one(static_cast<Element>(1))
+ , mOne(static_cast<Element>(p-1))
+ , _p(static_cast<Residu_t>(p))
+ , _bitsizep(0)
+ {
+ assert(_p >= minCardinality());
+ assert(_p <= maxCardinality());
+ Residu_t __p = _p;
+ while (__p != 0) {
+ _bitsizep++;
+ __p >>= 1;
+ }
+ }
+
+ Modular(const Self_t& F)
+ : zero(F.zero), one(F.one), mOne(F.mOne), _p(F._p), _dp(F._dp), _bitsizep(F._bitsizep) {}
+
+ // ----- Accessors
+ inline Element minElement() const override { return zero; }
+ inline Element maxElement() const override { return mOne; }
+
+ // ----- Access to the modulus
+ inline Residu_t residu() const { return _p; }
+ inline Residu_t size() const { return _p; }
+ inline Residu_t characteristic() const { return _p; }
+ inline Residu_t cardinality() const { return _p; }
+ template<class T> inline T& characteristic(T& p) const { return p = _p; }
+ template<class T> inline T& cardinality(T& p) const { return p = _p; }
+ static inline Residu_t maxCardinality();
+ static inline Residu_t minCardinality() { return 2; }
+
+ // ----- Checkers
+ inline bool isZero(const Element& a) const override { return a == zero; }
+ inline bool isOne (const Element& a) const override { return a == one; }
+ inline bool isMOne(const Element& a) const override { return a == mOne; }
+ inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
+ inline size_t length(const Element a) const { return size_rep; }
+
+ // ----- Ring-wise operators
+ bool operator==(const Self_t& F) const { return _p == F._p; }
+ bool operator!=(const Self_t& F) const { return _p != F._p; }
+ Self_t& operator=(const Self_t& F)
+ {
+ F.assign(const_cast<Element&>(one), F.one);
+ F.assign(const_cast<Element&>(zero), F.zero);
+ F.assign(const_cast<Element&>(mOne), F.mOne);
+ _p = F._p;
+ _dp = F._dp;
+ _bitsizep = F._bitsizep;
+ return *this;
+ }
+
+ // ----- Initialisation
+ Element& init (Element& x) const
+ { return x = 0; }
+ Element& init (Element& x, const double a) const;
+ Element& init (Element& x, const int64_t a) const;
+ Element& init (Element& x, const uint64_t a) const;
+ Element& init (Element& x, const Integer& a) const;
+ template<typename T> Element& init(Element& r, const T& a) const
+ {
+ reduce( Caster<Element,T>(r, (a < 0)? -a : a) );
+ return (a<0 ? negin(r) : r) ;
+ }
+
+ Element& assign (Element& x, const Element& y) const
+ { return x = y; }
+
+ // ----- Convert and reduce
+ template<typename T> T& convert(T& r, const Element& a) const
+ { return Caster<T,Element>(r,a); }
+
+ Element& reduce (Element& x, const Element& y) const
+ { return x = y % _p;}
+ Element& reduce (Element& x) const
+ { return x %= _p; }
+
+ // ----- Classic arithmetic
+ Element& mul(Element& r, const Element& a, const Element& b) const override;
+ Element& div(Element& r, const Element& a, const Element& b) const override;
+ Element& add(Element& r, const Element& a, const Element& b) const override;
+ Element& sub(Element& r, const Element& a, const Element& b) const override;
+ Element& neg(Element& r, const Element& a) const override;
+ Element& inv(Element& r, const Element& a) const override;
+
+ Element& mulin(Element& r, const Element& a) const override;
+ Element& divin(Element& r, const Element& a) const override;
+ Element& addin(Element& r, const Element& a) const override;
+ Element& subin(Element& r, const Element& a) const override;
+ Element& negin(Element& r) const override;
+ Element& invin(Element& r) const override;
+
+ // Functions defined in modular-mulprecomp
+ //
+ // void precomp_p (Compute_t& invp) const
+ // Element& mul_precomp_p (Element& r, const Element& a, const Element& b, const Compute_t& invp) const
+ //
+ // void precomp_b (Compute_t& invb, const Element& b) const
+ // void precomp_b (Compute_t& invb, const Element& b, const Compute_t& invp) const
+ // Element& mul_precomp_b (Element& r, const Element& a, const Element& b, const Compute_t& invb) const
+
+#include "modular-mulprecomp.inl"
+
+ // -- axpy: r <- a * x + y
+ // -- axpyin: r <- a * x + r
+ Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& axpyin(Element& r, const Element& a, const Element& x) const override;
+
+ // -- axmy: r <- a * x - y
+ // -- axmyin: r <- a * x - r
+ Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& axmyin(Element& r, const Element& a, const Element& x) const override;
+
+ // -- maxpy: r <- y - a * x
+ // -- maxpyin: r <- r - a * x
+ Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& maxpyin(Element& r, const Element& a, const Element& x) const override;
+
+ // ----- Random generators
+ typedef ModularRandIter<Self_t> RandIter;
+ typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
+ template< class Random > Element& random(Random& g, Element& r) const
+ { return init(r, g()); }
+ template< class Random > Element& nonzerorandom(Random& g, Element& a) const
+ { while (isZero(init(a, g())))
+ ;
+ return a; }
+
+ // --- IO methods
+ std::ostream& write(std::ostream& s) const;
+ std::istream& read (std::istream& s, Element& a) const;
+ std::ostream& write(std::ostream& s, const Element a) const;
+
+ protected:
+
+ Residu_t _p;
+ double _dp;
+ size_t _bitsizep;
+
+ };
}
#include "givaro/modular-uint32.inl"
-#endif
-
+#endif // __GIVARO_modular_uint32_H
diff --git a/src/kernel/ring/modular-uint64.h b/src/kernel/ring/modular-uint64.h
index f9d4dbf..dc2d7b7 100644
--- a/src/kernel/ring/modular-uint64.h
+++ b/src/kernel/ring/modular-uint64.h
@@ -1,3 +1,5 @@
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
// ==========================================================================
// Copyright(c)'1994-2015 by The Givaro group
// This file is part of Givaro.
@@ -24,157 +26,182 @@
namespace Givaro
{
- /** \brief Specialization of Modular to uint64_t element type with efficient dot product.
- *
- * Efficient element operations for dot product, mul, axpy, by using floating point
- * inverse of modulus (borrowed from NTL) and some use of non-normalized intermediate values.
- *
- * For some uses this is the most efficient field for primes in the range from half word
- * to 2^30.
- *
- * Requires: Modulus < 2^30.
- * Intended use: 2^15 < prime modulus < 2^30.
- * \ingroup field
- */
- template<typename COMP>
- class Modular<uint64_t, COMP> : public virtual FiniteFieldInterface<uint64_t>
- {
- public:
- // ----- Exported Types and constantes
- using Self_t = Modular<uint64_t, COMP>;
- using Compute_t = typename std::make_unsigned<COMP>::type;
- using Residu_t = uint64_t;
-
- enum { size_rep = sizeof(Residu_t) };
-
- // ----- Constantes
- const Element zero;
- const Element one;
- const Element mOne;
-
- // ----- Constructors
- Modular()
- : zero(0), one(1), mOne(0), _p(0) {}
-
- Modular(Residu_t p)
- : zero(0), one(1), mOne((Element)p-1), _p(p)
- {
- assert(_p >= minCardinality());
- assert(_p <= maxCardinality());
- }
-
- Modular(const Self_t& F)
- : zero(F.zero), one(F.one), mOne(F.mOne), _p(F._p) {}
-
- // ----- Accessors
- inline Element minElement() const override { return zero; }
- inline Element maxElement() const override { return mOne; }
-
- // ----- Access to the modulus
- inline Residu_t residu() const { return _p; }
- inline Residu_t size() const { return _p; }
- inline Residu_t characteristic() const { return _p; }
- inline Residu_t cardinality() const { return _p; }
- template<class T> inline T& characteristic(T& p) const { return p = _p; }
- template<class T> inline T& cardinality(T& p) const { return p = _p; }
- static inline Residu_t maxCardinality();
- static inline Residu_t minCardinality() { return 2; }
-
- // ----- Checkers
- inline bool isZero(const Element& a) const override { return a == zero; }
- inline bool isOne (const Element& a) const override { return a == one; }
- inline bool isMOne(const Element& a) const override { return a == mOne; }
- inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
- inline size_t length(const Element a) const { return size_rep; }
-
- // ----- Ring-wise operators
- bool operator==(const Self_t& F) const { return _p == F._p; }
- bool operator!=(const Self_t& F) const { return _p != F._p; }
- Self_t& operator=(const Self_t& F)
- {
- F.assign(const_cast<Element&>(one), F.one);
- F.assign(const_cast<Element&>(zero), F.zero);
- F.assign(const_cast<Element&>(mOne), F.mOne);
- _p = F._p;
- return *this;
- }
-
- // ----- Initialisation
- Element& init (Element& x) const
- { return x = 0; }
- Element& init (Element& x, const Integer& a) const;
- template<typename T> Element& init(Element& r, const T& a) const
- {
- reduce(r, Caster<Element>((a < 0)? -a : a));
- if (a < 0) negin(r);
- return r;
- }
-
- Element& assign (Element& x, const Element& y) const
- { return x = y; }
-
- // ----- Convert and reduce
- template<typename T> T& convert(T& r, const Element& a) const
- { return r = Caster<T>(a); }
-
- Element& reduce (Element& x, const Element& y) const
- { x = y % _p; return x; }
- Element& reduce (Element& x) const
- { x %= _p; return x; }
-
- // ----- Classic arithmetic
- Element& mul(Element& r, const Element& a, const Element& b) const override;
- Element& div(Element& r, const Element& a, const Element& b) const override;
- Element& add(Element& r, const Element& a, const Element& b) const override;
- Element& sub(Element& r, const Element& a, const Element& b) const override;
- Element& neg(Element& r, const Element& a) const override;
- Element& inv(Element& r, const Element& a) const override;
-
- Element& mulin(Element& r, const Element& a) const override;
- Element& divin(Element& r, const Element& a) const override;
- Element& addin(Element& r, const Element& a) const override;
- Element& subin(Element& r, const Element& a) const override;
- Element& negin(Element& r) const override;
- Element& invin(Element& r) const override;
-
- // -- axpy: r <- a * x + y
- // -- axpyin: r <- a * x + r
- Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& axpyin(Element& r, const Element& a, const Element& x) const override;
-
- // -- axmy: r <- a * x - y
- // -- axmyin: r <- a * x - r
- Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& axmyin(Element& r, const Element& a, const Element& x) const override;
-
- // -- maxpy: r <- y - a * x
- // -- maxpyin: r <- r - a * x
- Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& maxpyin(Element& r, const Element& a, const Element& x) const override;
-
- // ----- Random generators
- typedef ModularRandIter<Self_t> RandIter;
- typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
- template< class Random > Element& random(Random& g, Element& r) const
- { return init(r, g()); }
- template< class Random > Element& nonzerorandom(Random& g, Element& a) const
- { while (isZero(init(a, g())))
- ;
- return a; }
-
- // --- IO methods
- std::istream& read (std::istream& s);
- std::ostream& write(std::ostream& s) const;
- std::istream& read (std::istream& s, Element& a) const;
- std::ostream& write(std::ostream& s, const Element& a) const;
-
- protected:
-
- Residu_t _p;
- };
+ /** \brief Specialization of Modular to uint64_t element type with efficient dot product.
+ *
+ * Efficient element operations for dot product, mul, axpy, by using floating point
+ * inverse of modulus (borrowed from NTL) and some use of non-normalized intermediate values.
+ *
+ * For some uses this is the most efficient field for primes in the range from half word
+ * to 2^30.
+ *
+ * Requires: Modulus < 2^30.
+ * Intended use: 2^15 < prime modulus < 2^30.
+ * \ingroup field
+ */
+ template<typename COMP>
+ class Modular<uint64_t, COMP> : public virtual FiniteFieldInterface<uint64_t>
+ {
+ public:
+ // ----- Exported Types and constantes
+ using Self_t = Modular<uint64_t, COMP>;
+ using Compute_t = typename std::make_unsigned<COMP>::type;
+ using Residu_t = uint64_t;
+
+ enum { size_rep = sizeof(Residu_t) };
+
+ // ----- Constantes
+ const Element zero;
+ const Element one;
+ const Element mOne;
+
+ // ----- Constructors
+ Modular()
+ : zero(static_cast<Element>(0))
+ , one(static_cast<Element>(1))
+ , mOne(static_cast<Element>(-1))
+ , _p(static_cast<Residu_t>(0))
+ , _bitsizep(0) {}
+
+ Modular(const Residu_t p)
+ : zero(static_cast<Element>(0))
+ , one(static_cast<Element>(1))
+ , mOne(static_cast<Element>(p-1))
+ , _p(static_cast<Residu_t>(p))
+ , _bitsizep(0)
+ {
+ assert(_p >= minCardinality());
+ assert(_p <= maxCardinality());
+ Residu_t __p = _p;
+ while (__p != 0) {
+ _bitsizep++;
+ __p >>= 1;
+ }
+ }
+
+ Modular(const Self_t& F)
+ : zero(F.zero), one(F.one), mOne(F.mOne), _p(F._p), _bitsizep(F._bitsizep) {}
+
+ // ----- Accessors
+ inline Element minElement() const override { return zero; }
+ inline Element maxElement() const override { return mOne; }
+
+ // ----- Access to the modulus
+ inline Residu_t residu() const { return _p; }
+ inline Residu_t size() const { return _p; }
+ inline Residu_t characteristic() const { return _p; }
+ inline Residu_t cardinality() const { return _p; }
+ template<class T> inline T& characteristic(T& p) const { return p = _p; }
+ template<class T> inline T& cardinality(T& p) const { return p = _p; }
+ static inline Residu_t maxCardinality();
+ static inline Residu_t minCardinality() { return 2_ui64; }
+
+ // ----- Checkers
+ inline bool isZero(const Element& a) const override { return a == zero; }
+ inline bool isOne (const Element& a) const override { return a == one; }
+ inline bool isMOne(const Element& a) const override { return a == mOne; }
+ inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
+ inline size_t length(const Element a) const { return size_rep; }
+
+ // ----- Ring-wise operators
+ bool operator==(const Self_t& F) const { return _p == F._p; }
+ bool operator!=(const Self_t& F) const { return _p != F._p; }
+ Self_t& operator=(const Self_t& F)
+ {
+ F.assign(const_cast<Element&>(one), F.one);
+ F.assign(const_cast<Element&>(zero), F.zero);
+ F.assign(const_cast<Element&>(mOne), F.mOne);
+ _p = F._p;
+ _bitsizep = F._bitsizep;
+ return *this;
+ }
+
+ // ----- Initialisation
+ Element& init (Element& x) const
+ { return x = 0; }
+ Element& init (Element& x, const Integer& a) const;
+ template<typename T> Element& init(Element& r, const T& a) const
+ {
+ reduce(r, Caster<Element>((a < 0)? -a : a));
+ if (a < 0) negin(r);
+ return r;
+ }
+
+ Element& assign (Element& x, const Element& y) const
+ { return x = y; }
+
+ // ----- Convert and reduce
+ template<typename T> T& convert(T& r, const Element& a) const
+ { return r = Caster<T>(a); }
+
+ Element& reduce (Element& x, const Element& y) const
+ { x = y % _p; return x; }
+ Element& reduce (Element& x) const
+ { x %= _p; return x; }
+
+ // ----- Classic arithmetic
+ Element& mul(Element& r, const Element& a, const Element& b) const override;
+ Element& div(Element& r, const Element& a, const Element& b) const override;
+ Element& add(Element& r, const Element& a, const Element& b) const override;
+ Element& sub(Element& r, const Element& a, const Element& b) const override;
+ Element& neg(Element& r, const Element& a) const override;
+ Element& inv(Element& r, const Element& a) const override;
+
+ Element& mulin(Element& r, const Element& a) const override;
+ Element& divin(Element& r, const Element& a) const override;
+ Element& addin(Element& r, const Element& a) const override;
+ Element& subin(Element& r, const Element& a) const override;
+ Element& negin(Element& r) const override;
+ Element& invin(Element& r) const override;
+
+ // Functions defined in modular-mulprecomp
+ //
+ // void precomp_p (Compute_t& invp) const
+ // Element& mul_precomp_p (Element& r, const Element& a, const Element& b, const Compute_t& invp) const
+ //
+ // void precomp_b (Compute_t& invb, const Element& b) const
+ // void precomp_b (Compute_t& invb, const Element& b, const Compute_t& invp) const
+ // Element& mul_precomp_b (Element& r, const Element& a, const Element& b, const Compute_t& invb) const
+
+#include "modular-mulprecomp.inl"
+
+ // -- axpy: r <- a * x + y
+ // -- axpyin: r <- a * x + r
+ Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& axpyin(Element& r, const Element& a, const Element& x) const override;
+
+ // -- axmy: r <- a * x - y
+ // -- axmyin: r <- a * x - r
+ Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& axmyin(Element& r, const Element& a, const Element& x) const override;
+
+ // -- maxpy: r <- y - a * x
+ // -- maxpyin: r <- r - a * x
+ Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& maxpyin(Element& r, const Element& a, const Element& x) const override;
+
+ // ----- Random generators
+ typedef ModularRandIter<Self_t> RandIter;
+ typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
+ template< class Random > Element& random(Random& g, Element& r) const
+ { return init(r, g()); }
+ template< class Random > Element& nonzerorandom(Random& g, Element& a) const
+ { while (isZero(init(a, g())))
+ ;
+ return a; }
+
+ // --- IO methods
+ std::istream& read (std::istream& s);
+ std::ostream& write(std::ostream& s) const;
+ std::istream& read (std::istream& s, Element& a) const;
+ std::ostream& write(std::ostream& s, const Element& a) const;
+
+ protected:
+
+ Residu_t _p;
+ size_t _bitsizep;
+ };
}
#include "givaro/modular-uint64.inl"
-#endif
-
+#endif // __GIVARO_modular_uint64_H
diff --git a/src/kernel/ring/modular-uint64.inl b/src/kernel/ring/modular-uint64.inl
index c3a65e1..32be263 100644
--- a/src/kernel/ring/modular-uint64.inl
+++ b/src/kernel/ring/modular-uint64.inl
@@ -21,11 +21,21 @@ namespace Givaro
template<>
inline Modular<uint64_t, uint64_t>::Residu_t
- Modular<uint64_t, uint64_t>::maxCardinality() { return 4294967295u; } // 2^32 - 1
+ Modular<uint64_t, uint64_t>::maxCardinality() { return 4294967295_ui64; } // 2^32 - 1
template<>
inline Modular<uint64_t, int64_t>::Residu_t
- Modular<uint64_t, int64_t>::maxCardinality() { return 4294967295u; }
+ Modular<uint64_t, int64_t>::maxCardinality() { return 4294967295_ui64; }
+
+#ifdef __GIVARO_HAVE_INT128
+ template<>
+ inline Modular<uint64_t, uint128_t>::Residu_t
+ Modular<uint64_t, uint128_t>::maxCardinality() { return 9223372036854775807_ui64; } // 2^63 - 1
+
+ template<>
+ inline Modular<uint64_t, int128_t>::Residu_t
+ Modular<uint64_t, int128_t>::maxCardinality() { return 9223372036854775807_ui64; }
+#endif
// ------------------------
// ----- Classic arithmetic
@@ -215,7 +225,22 @@ namespace Givaro
inline std::ostream &Modular<uint64_t, uint64_t>::write (std::ostream &os) const
{
return os << "Modular<uint64_t, uint64_t> modulo " << _p;
- }
+ }
+
+#ifdef __GIVARO_HAVE_INT128
+ template<>
+ inline std::ostream& Modular<uint64_t, int128_t>::write (std::ostream& s ) const
+ {
+ return s << "Modular<uint64_t, int128_t> modulo " << residu();
+ }
+
+ template<>
+ inline std::ostream& Modular<uint64_t, uint128_t>::write (std::ostream& s ) const
+ {
+ return s << "Modular<uint64_t, uint128_t> modulo " << residu();
+ }
+#endif
+
template<typename COMP>
inline std::istream &Modular<uint64_t, COMP>::read (std::istream &is)
diff --git a/src/kernel/ring/modular-uint8.h b/src/kernel/ring/modular-uint8.h
index 20cf058..a727c4b 100644
--- a/src/kernel/ring/modular-uint8.h
+++ b/src/kernel/ring/modular-uint8.h
@@ -1,3 +1,5 @@
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
// ==========================================================================
// Copyright(c)'1994-2015 by The Givaro group
// This file is part of Givaro.
@@ -6,11 +8,12 @@
// see the COPYRIGHT file for more details.
// Authors: Brice Boyer (briceboyer) <boyer.brice at gmail.com>
// A. Breust (taken from FFLAS-FFPACK)
+// R. Lebreton
// ==========================================================================
-/*! @file field/modular-uint64.h
- * @ingroup field
- * @brief representation of <code>Z/mZ</code> over \c uint64_t .
+/*! @file ring/modular-uint8.h
+ * @ingroup ring
+ * @brief representation of <code>Z/mZ</code> over \c uint8_t .
*/
#ifndef __GIVARO_modular_uint8_H
@@ -23,157 +26,179 @@
namespace Givaro {
- template <typename COMP>
- class Modular<uint8_t, COMP> : public virtual FiniteFieldInterface<uint8_t>
- {
- public:
- // ----- Exported Types and constantes
- using Self_t = Modular<uint8_t, COMP>;
- using Compute_t = typename std::make_unsigned<COMP>::type;
- using Residu_t = uint8_t;
-
- enum { size_rep = sizeof(Residu_t) };
-
- // ----- Constantes
- const Element zero;
- const Element one;
- const Element mOne;
-
- // ----- Constructors
- Modular()
- : _p(static_cast<Residu_t>(0)) {}
-
- Modular(const Residu_t p)
- : zero(static_cast<Element>(0))
- , one(static_cast<Element>(1))
- , mOne(static_cast<Element>(p-1))
- , _p(static_cast<Residu_t>(p))
+ template <typename COMP>
+ class Modular<uint8_t, COMP> : public virtual FiniteFieldInterface<uint8_t>
{
- assert(_p >= minCardinality());
- assert(_p <= maxCardinality());
- }
-
- Modular(const Self_t& F)
- : zero(F.zero), one(F.one), mOne(F.mOne), _p(F._p) {}
-
- // ----- Accessors
- inline Element minElement() const override { return zero; }
- inline Element maxElement() const override { return mOne; }
-
- // ----- Access to the modulus
- inline Residu_t residu() const { return _p; }
- inline Residu_t size() const { return _p; }
- inline Residu_t characteristic() const { return _p; }
- inline Residu_t cardinality() const { return _p; }
- template<class T> inline T& characteristic(T& p) const { return p = _p; }
- template<class T> inline T& cardinality(T& p) const { return p = _p; }
- static inline Residu_t maxCardinality();
- static inline Residu_t minCardinality() { return 2; }
-
- // ----- Checkers
- inline bool isZero(const Element& a) const override { return a == zero; }
- inline bool isOne (const Element& a) const override { return a == one; }
- inline bool isMOne(const Element& a) const override { return a == mOne; }
- inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
- inline size_t length(const Element a) const { return size_rep; }
-
- // ----- Ring-wise operators
- bool operator==(const Self_t& F) const { return _p == F._p; }
- bool operator!=(const Self_t& F) const { return _p != F._p; }
- Self_t& operator=(const Self_t& F)
- {
- F.assign(const_cast<Element&>(one), F.one);
- F.assign(const_cast<Element&>(zero), F.zero);
- F.assign(const_cast<Element&>(mOne), F.mOne);
- _p = F._p;
- return *this;
- }
-
- // ----- Initialisation
- Element& init (Element& x) const
- { return x = 0; }
- Element& init (Element& x, const float a) const;
- Element& init (Element& x, const double a) const;
- Element& init (Element& x, const int16_t a) const;
- Element& init (Element& x, const uint16_t a) const;
- Element& init (Element& x, const int32_t a) const;
- Element& init (Element& x, const uint32_t a) const;
- Element& init (Element& x, const int64_t a) const;
- Element& init (Element& x, const uint64_t a) const;
- Element& init (Element& x, const Integer& a) const;
- template<typename T> Element& init(Element& r, const T& a) const
- {
- reduce(r, Caster<Element>((a < 0)? -a : a));
- if (a < 0) negin(r);
- return r;
- }
-
- Element& assign (Element& x, const Element& y) const
- { return x = y; }
-
- // ----- Convert and reduce
- template<typename T> T& convert(T& r, const Element& a) const
- { return r = Caster<T>(a); }
-
- Element& reduce (Element& x, const Element& y) const
- { x = y % _p; return x; }
- Element& reduce (Element& x) const
- { x %= _p; return x; }
-
- // ----- Classic arithmetic
- Element& mul(Element& r, const Element& a, const Element& b) const override;
- Element& div(Element& r, const Element& a, const Element& b) const override;
- Element& add(Element& r, const Element& a, const Element& b) const override;
- Element& sub(Element& r, const Element& a, const Element& b) const override;
- Element& neg(Element& r, const Element& a) const override;
- Element& inv(Element& r, const Element& a) const override;
-
- Element& mulin(Element& r, const Element& a) const override;
- Element& divin(Element& r, const Element& a) const override;
- Element& addin(Element& r, const Element& a) const override;
- Element& subin(Element& r, const Element& a) const override;
- Element& negin(Element& r) const override;
- Element& invin(Element& r) const override;
-
- // -- axpy: r <- a * x + y
- // -- axpyin: r <- a * x + r
- Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& axpyin(Element& r, const Element& a, const Element& x) const override;
-
- // -- axmy: r <- a * x - y
- // -- axmyin: r <- a * x - r
- Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& axmyin(Element& r, const Element& a, const Element& x) const override;
-
- // -- maxpy: r <- y - a * x
- // -- maxpyin: r <- r - a * x
- Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
- Element& maxpyin(Element& r, const Element& a, const Element& x) const override;
-
- // ----- Random generators
- typedef ModularRandIter<Self_t> RandIter;
- typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
- template< class Random > Element& random(Random& g, Element& r) const
- { return init(r, g()); }
- template< class Random > Element& nonzerorandom(Random& g, Element& a) const
- { while (isZero(init(a, g())))
- ;
- return a; }
-
- // --- IO methods
- std::istream& read (std::istream& s);
- std::ostream& write(std::ostream& s) const;
- std::istream& read (std::istream& s, Element& a) const;
- std::ostream& write(std::ostream& s, const Element& a) const;
-
- protected:
-
- Residu_t _p;
- };
+ public:
+ // ----- Exported Types and constantes
+ using Self_t = Modular<uint8_t, COMP>;
+ using Compute_t = typename std::make_unsigned<COMP>::type;
+ using Residu_t = uint8_t;
+
+ enum { size_rep = sizeof(Residu_t) };
+
+ // ----- Constantes
+ const Element zero;
+ const Element one;
+ const Element mOne;
+
+ // ----- Constructors
+ Modular()
+ : zero(static_cast<Element>(0))
+ , one(static_cast<Element>(1))
+ , mOne(static_cast<Element>(-1))
+ , _p(static_cast<Residu_t>(0))
+ , _bitsizep(0) {}
+
+ Modular(const Residu_t p)
+ : zero(static_cast<Element>(0))
+ , one(static_cast<Element>(1))
+ , mOne(static_cast<Element>(p-1))
+ , _p(static_cast<Residu_t>(p))
+ , _bitsizep(0)
+ {
+ assert(_p >= minCardinality());
+ assert(_p <= maxCardinality());
+ Residu_t __p = _p;
+ while (__p != 0) {
+ _bitsizep++;
+ __p >>= 1;
+ }
+ }
+
+ Modular(const Self_t& F)
+ : zero(F.zero), one(F.one), mOne(F.mOne), _p(F._p), _bitsizep(F._bitsizep) {}
+
+ // ----- Accessors
+ inline Element minElement() const override { return zero; }
+ inline Element maxElement() const override { return mOne; }
+
+ // ----- Access to the modulus
+ inline Residu_t residu() const { return _p; }
+ inline Residu_t size() const { return _p; }
+ inline Residu_t characteristic() const { return _p; }
+ inline Residu_t cardinality() const { return _p; }
+ template<class T> inline T& characteristic(T& p) const { return p = _p; }
+ template<class T> inline T& cardinality(T& p) const { return p = _p; }
+ static inline Residu_t maxCardinality();
+ static inline Residu_t minCardinality() { return 2; }
+
+ // ----- Checkers
+ inline bool isZero(const Element& a) const override { return a == zero; }
+ inline bool isOne (const Element& a) const override { return a == one; }
+ inline bool isMOne(const Element& a) const override { return a == mOne; }
+ inline bool areEqual(const Element& a, const Element& b) const override { return a == b; }
+ inline size_t length(const Element a) const { return size_rep; }
+
+ // ----- Ring-wise operators
+ bool operator==(const Self_t& F) const { return _p == F._p; }
+ bool operator!=(const Self_t& F) const { return _p != F._p; }
+ Self_t& operator=(const Self_t& F)
+ {
+ F.assign(const_cast<Element&>(one), F.one);
+ F.assign(const_cast<Element&>(zero), F.zero);
+ F.assign(const_cast<Element&>(mOne), F.mOne);
+ _p = F._p;
+ _bitsizep = F._bitsizep;
+ return *this;
+ }
+
+ // ----- Initialisation
+ Element& init (Element& x) const
+ { return x = 0; }
+ Element& init (Element& x, const float a) const;
+ Element& init (Element& x, const double a) const;
+ Element& init (Element& x, const int16_t a) const;
+ Element& init (Element& x, const uint16_t a) const;
+ Element& init (Element& x, const int32_t a) const;
+ Element& init (Element& x, const uint32_t a) const;
+ Element& init (Element& x, const int64_t a) const;
+ Element& init (Element& x, const uint64_t a) const;
+ Element& init (Element& x, const Integer& a) const;
+ template<typename T> Element& init(Element& r, const T& a) const
+ {
+ reduce(r, Caster<Element>((a < 0)? -a : a));
+ if (a < 0) negin(r);
+ return r;
+ }
+
+ Element& assign (Element& x, const Element& y) const
+ { return x = y; }
+
+ // ----- Convert and reduce
+ template<typename T> T& convert(T& r, const Element& a) const
+ { return r = Caster<T>(a); }
+
+ Element& reduce (Element& x, const Element& y) const
+ { x = y % _p; return x; }
+ Element& reduce (Element& x) const
+ { x %= _p; return x; }
+
+ // ----- Classic arithmetic
+ Element& mul(Element& r, const Element& a, const Element& b) const override;
+ Element& div(Element& r, const Element& a, const Element& b) const override;
+ Element& add(Element& r, const Element& a, const Element& b) const override;
+ Element& sub(Element& r, const Element& a, const Element& b) const override;
+ Element& neg(Element& r, const Element& a) const override;
+ Element& inv(Element& r, const Element& a) const override;
+
+ Element& mulin(Element& r, const Element& a) const override;
+ Element& divin(Element& r, const Element& a) const override;
+ Element& addin(Element& r, const Element& a) const override;
+ Element& subin(Element& r, const Element& a) const override;
+ Element& negin(Element& r) const override;
+ Element& invin(Element& r) const override;
+
+ // Functions defined in modular-mulprecomp
+ //
+ // void precomp_p (Compute_t& invp) const
+ // Element& mul_precomp_p (Element& r, const Element& a, const Element& b, const Compute_t& invp) const
+ //
+ // void precomp_b (Compute_t& invb, const Element& b) const
+ // void precomp_b (Compute_t& invb, const Element& b, const Compute_t& invp) const
+ // Element& mul_precomp_b (Element& r, const Element& a, const Element& b, const Compute_t& invb) const
+
+#include "modular-mulprecomp.inl"
+
+ // -- axpy: r <- a * x + y
+ // -- axpyin: r <- a * x + r
+ Element& axpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& axpyin(Element& r, const Element& a, const Element& x) const override;
+
+ // -- axmy: r <- a * x - y
+ // -- axmyin: r <- a * x - r
+ Element& axmy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& axmyin(Element& r, const Element& a, const Element& x) const override;
+
+ // -- maxpy: r <- y - a * x
+ // -- maxpyin: r <- r - a * x
+ Element& maxpy (Element& r, const Element& a, const Element& x, const Element& y) const override;
+ Element& maxpyin(Element& r, const Element& a, const Element& x) const override;
+
+ // ----- Random generators
+ typedef ModularRandIter<Self_t> RandIter;
+ typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
+ template< class Random > Element& random(Random& g, Element& r) const
+ { return init(r, g()); }
+ template< class Random > Element& nonzerorandom(Random& g, Element& a) const
+ { while (isZero(init(a, g())))
+ ;
+ return a; }
+
+ // --- IO methods
+ std::istream& read (std::istream& s);
+ std::ostream& write(std::ostream& s) const;
+ std::istream& read (std::istream& s, Element& a) const;
+ std::ostream& write(std::ostream& s, const Element& a) const;
+
+ protected:
+
+ Residu_t _p;
+ size_t _bitsizep;
+ };
}
#include "givaro/modular-uint8.inl"
-#endif
-
+#endif // __GIVARO_modular_uint8_H
diff --git a/src/kernel/ring/modular-uint8.inl b/src/kernel/ring/modular-uint8.inl
index b3299b7..d736bbb 100644
--- a/src/kernel/ring/modular-uint8.inl
+++ b/src/kernel/ring/modular-uint8.inl
@@ -1,3 +1,5 @@
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
// ==========================================================================
// Copyright(c)'1994-2015 by The Givaro group
// This file is part of Givaro.
@@ -300,5 +302,4 @@ namespace Givaro {
} // namespace Givaro
-#endif
-
+#endif // __GIVARO_modular_uint8_INL
diff --git a/src/kernel/ring/modular.h b/src/kernel/ring/modular.h
index 221af44..b199751 100644
--- a/src/kernel/ring/modular.h
+++ b/src/kernel/ring/modular.h
@@ -24,13 +24,13 @@
#include "givaro/modular-uint16.h"
#include "givaro/modular-int32.h"
#include "givaro/modular-uint32.h"
+#include "givaro/modular-int64.h"
+#include "givaro/modular-uint64.h"
#include "givaro/modular-float.h"
#include "givaro/modular-double.h"
#include "givaro/modular-integer.h"
#include "givaro/modular-inttype.h"
#include "givaro/modular-log16.h"
-#include "givaro/modular-int64.h"
-#include "givaro/modular-uint64.h"
#include "givaro/modular-ruint.h"
#endif
diff --git a/src/kernel/ring/montgomery-int32.h b/src/kernel/ring/montgomery-int32.h
index 019be20..b841a74 100644
--- a/src/kernel/ring/montgomery-int32.h
+++ b/src/kernel/ring/montgomery-int32.h
@@ -183,8 +183,7 @@ namespace Givaro
template< class Random > Element& random(Random& g, Element& r) const
{ return init(r, g()); }
template< class Random > Element& nonzerorandom(Random& g, Element& a) const
- { while (isZero(init(a, g())));
- return a; }
+ { while (isZero(init(a, g()))) {} return a; }
// --- IO methods
std::ostream& write(std::ostream& s) const;
diff --git a/src/kernel/ring/montgomery-ruint.h b/src/kernel/ring/montgomery-ruint.h
index 4296bfb..feee3d1 100755
--- a/src/kernel/ring/montgomery-ruint.h
+++ b/src/kernel/ring/montgomery-ruint.h
@@ -173,10 +173,10 @@ namespace Givaro
// ----- Random generators
typedef ModularRandIter<Self_t> RandIter;
typedef GeneralRingNonZeroRandIter<Self_t> NonZeroRandIter;
- template< class Random > Element& random(Random& g, Element& r) const { RecInt::rand(r); mod_n(r, _p); return r; }
- template< class Random > Element& nonzerorandom(Random& g, Element& a) const
- { while (isZero(random(g, a)));
- return a; }
+ template< class Random > Element& random(Random& g, Element& r) const
+ { RecInt::rand(r); mod_n(r, _p); return r; }
+ template< class Random > Element& nonzerorandom(Random& g, Element& a) const
+ { while (isZero(random(g, a))) {} return a; }
// --- IO methods
std::istream& read (std::istream& s);
diff --git a/src/kernel/system/givconfig.h b/src/kernel/system/givconfig.h
index e3278fc..cc1bbdd 100644
--- a/src/kernel/system/givconfig.h
+++ b/src/kernel/system/givconfig.h
@@ -66,8 +66,8 @@
// - zz: revision number
#define GIVARO_MAJOR_VERSION 4
#define GIVARO_MINOR_VERSION 0
-#define GIVARO_REVISION_VERSION 1
-#define GIVARO_VERSION 40001
+#define GIVARO_REVISION_VERSION 2
+#define GIVARO_VERSION 40002
// -- Defines this value both to compile the library of user program
// value: integer that defines debug level trace information (not well defined)
@@ -427,7 +427,6 @@ template<> struct Signed_Trait<unsigned long> : public GIVARO_numeric_limits<un
#endif
-
#if defined(_OPENMP) || defined(OMP_H) || defined(__OMP_H) || defined(__pmp_omp_h)
# ifndef __GIVARO_USE_OPENMP
# define __GIVARO_USE_OPENMP 1
@@ -436,6 +435,13 @@ template<> struct Signed_Trait<unsigned long> : public GIVARO_numeric_limits<un
//# undef __GIVARO_USE_OPENMP
#endif
+#ifdef __GIVARO_HAVE_INT128
+/* Define int128 type */
+#define int128_t __int128_t
+
+/* Define uint128 type */
+#define uint128_t __uint128_t
+#endif
#endif
// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 81afb19..36c961e 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -33,7 +33,8 @@ BASIC_TESTS = \
test-geom \
test-integer \
test-conversion \
- test-ratrecon
+ test-ratrecon \
+ test-modularmulprecomp
# *_mg tests accept -DMG_DEFAULT=MG_ACTIVE to enable Montgomery representation
RECINT_TESTS= \
@@ -97,6 +98,7 @@ test_random_SOURCES = test-random.C
test_ifactor_SOURCES = test-ifactor.C
test_ffarith_SOURCES = test-ffarith.C
test_ringarith_SOURCES= test-ringarith.C
+test_modularmulprecomp_SOURCES= test-modularmulprecomp.C
test_mod_SOURCES = test-mod.C
test_modsqroot_SOURCES= test-modsqroot.C
test_trunc_SOURCES = test-trunc.C
diff --git a/tests/jenkins-maker.sh b/tests/jenkins-maker.sh
new file mode 100755
index 0000000..c6c2d29
--- /dev/null
+++ b/tests/jenkins-maker.sh
@@ -0,0 +1,85 @@
+#!/bin/bash
+# This file is part of the Givaro library.
+# It is distributed under the terms of the CeCILL-B licence
+# (see COPYING)
+# Created by AB - 2014/12/03
+# Modified by AC - 2016/04/08
+# Modified by CP - 2016/06/22
+
+# Some influential environment variables:
+# CXX C++ compiler command
+# CXXFLAGS C++ compiler flags
+
+# Note: This script is intended to be launched
+# by the Jenkins web interface whenever it needs
+# to compile the project.
+# But should be stored in /<slave_jenkins_path>/makers/
+# It is launched from the svn:trunk root directory.
+
+SOURCE_DIRECTORY=$( cd "$( dirname "$0" )" && pwd )
+
+
+#=============================#
+# Change only these variables #
+#=============================#
+CXX=`pwd | awk -F/ '{print $(NF)}'`
+
+JENKINS_DIR=${SOURCE_DIRECTORY%%/workspace/*}
+LOCAL_DIR="$JENKINS_DIR"/local
+
+# Where to install givaro binaries
+# Keep default for local installation.
+PREFIX_INSTALL="$LOCAL_DIR/$CXX"
+
+# Add path to compilers (if needed)
+export PATH="$PATH":"/usr/local/bin":"$LOCAL_DIR/$CXX/bin"
+# Add specific locations (if needed)
+export LD_LIBRARY_PATH="$LD_LIBRARY_PATH":"/usr/local/lib":"$LOCAL_DIR/$CXX/lib"
+
+# Where is GMP installed (compiled with cxx interface)
+# Keep empty if in usual folders (i.e. /usr or /usr/local)
+GMP_PATH="$LOCAL_DIR/$CXX"
+
+# /!\ Warning /!\ This could be an issue if you changed
+# the local installation directory
+rm -rf "$PREFIX_INSTALL"/bin/givaro* "$PREFIX_INSTALL"/include/givaro* "$PREFIX_INSTALL"/include/gmp++ "$PREFIX_INSTALL"/include/recint "$PREFIX_INSTALL"/lib/libgivaro*
+
+
+#================#
+# Setup Variables#
+#================#
+
+if [ "$CXX" == "icpc" ]; then
+ distribution=`uname -m`
+ if [ "$distribution" == "i686" ]; then
+ source /usr/local/bin/compilervars.sh ia32
+ else
+ source /usr/local/bin/compilervars.sh intel64
+ fi
+fi
+
+# Particular case for Fedora23: g++=g++-5.3
+vm_name=`uname -n | cut -d"-" -f1`
+
+if [[ "$vm_name" == "fedora" && "$CXX" == "g++-5.3" ]]; then
+ CXX="g++"
+fi
+
+#==================================#
+# Automated installation and tests #
+#==================================#
+
+echo "|=== JENKINS AUTOMATED SCRIPT ===| ./autogen.sh CXX=$CXX CXXFLAGS=$CXXFLAGS --prefix=$PREFIX_INSTALL --with-gmp=$GMP_PATH"
+./autogen.sh CXX=$CXX CXXFLAGS=$CXXFLAGS --prefix="$PREFIX_INSTALL" --with-gmp="$GMP_PATH"
+V="$?"; if test "x$V" != "x0"; then exit "$V"; fi
+
+echo "|=== JENKINS AUTOMATED SCRIPT ===| make install"
+make install
+V="$?"; if test "x$V" != "x0"; then exit "$V"; fi
+
+echo "|=== JENKINS AUTOMATED SCRIPT ===| make perfpublisher"
+make perfpublisher
+
+echo "|=== JENKINS AUTOMATED SCRIPT ===| make examples"
+make examples
+V="$?"; if test "x$V" != "x0"; then exit "$V"; fi
diff --git a/tests/perfpublisher.sh b/tests/perfpublisher.sh
index 2c3c452..15dba49 100755
--- a/tests/perfpublisher.sh
+++ b/tests/perfpublisher.sh
@@ -8,12 +8,23 @@ XMLFILE=$1
tests=$2
COMPILER=$3
+# choose gdate on OS X
+if command -v "gdate" >/dev/null; then
+ DATE=gdate
+else
+ DATE=date
+fi
#=================#
# Plateform infos #
#=================#
COMPILERVERSION=$($COMPILER --version 2>&1 | head -1)
-CPUFREQ=$(lscpu | grep "MHz" | rev | cut -f1 -d' ' | rev)
+
+if command -v "lscpu" >/dev/null; then
+ CPUFREQ=$(lscpu | grep "MHz" | rev | cut -f1 -d' ' | rev)
+else
+ CPUFREQ=$((`sysctl -n hw.cpufrequency`/1000000))
+fi
ARCH=$(uname -m)
OSNAME=$(uname -s)
OSVERSION=$(uname -r)
@@ -45,8 +56,8 @@ echo '<report name="tests-report" categ="tests">' >> $XMLFILE
#=======#
echo '<start>' >> $XMLFILE
-echo '<date format="YYYYMMDD" val="'$(date +%Y%m%d)'" />' >> $XMLFILE
-echo '<time format="HHMMSS" val="'$(date +%H%M%S)'" />' >> $XMLFILE
+echo '<date format="YYYYMMDD" val="'$($DATE +%Y%m%d)'" />' >> $XMLFILE
+echo '<time format="HHMMSS" val="'$($DATE +%H%M%S)'" />' >> $XMLFILE
echo '</start>' >> $XMLFILE
#=======#
@@ -59,9 +70,9 @@ do
then
#File does not exist: compile it
echo '[Compiling]' $test
- COMPILESTART=$(date +%s%3N)
+ COMPILESTART=$($DATE +%s%3N)
COMPILELOG=$(make $test 2>&1; echo 'Returned state: '$?)
- COMPILEEND=$(date +%s%3N)
+ COMPILEEND=$($DATE +%s%3N)
COMPILETIME=$(($COMPILEEND - $COMPILESTART))
COMPILECHECK=$(echo $COMPILELOG | grep -o '[^ ]*$')
COMPILETIMERELEVANT='true'
@@ -92,9 +103,9 @@ do
#Compilation success
echo '[Executing]' $test
EXECUTED='yes'
- EXECUTIONSTART=$(date +%s%3N)
+ EXECUTIONSTART=$($DATE +%s%3N)
EXECUTIONLOG=$(./$test 2>&1; echo 'Returned state: '$?)
- EXECUTIONEND=$(date +%s%3N)
+ EXECUTIONEND=$($DATE +%s%3N)
EXECUTIONTIME=$(($EXECUTIONEND - $EXECUTIONSTART))
EXECUTIONCHECK=$(echo $EXECUTIONLOG | grep -o '[^ ]*$')
diff --git a/tests/test-ffarith.C b/tests/test-ffarith.C
index ef59427..554df2c 100755
--- a/tests/test-ffarith.C
+++ b/tests/test-ffarith.C
@@ -7,6 +7,7 @@
#include <iostream>
+#include <givaro/givinteger.h>
#include <givaro/modular.h>
#include <givaro/modular-balanced.h>
#include <givaro/montgomery.h>
@@ -43,8 +44,11 @@ using namespace Givaro;
template<class Field>
bool invertible(const Field& F, const typename Field::Element& a)
{
- auto ai(a);
- return F.mulin(F.inv(ai, a), a) == F.one;
+// auto ai(a);
+// return F.mulin(F.inv(ai, a), a) == F.one;
+ Integer ai;
+ F.convert(ai,a);
+ return (gcd(ai,Integer(F.characteristic()))==1);
}
template<class Field>
@@ -93,7 +97,7 @@ int TestOneField(const Field& F, const typename Field::Element& first)
F.assign(a, first);
typename Field::RandIter g(F);
- while (!invertible(F, g.random(b)));
+ while (!invertible(F, g.random(b))) {}
F.init(c); // empty constructor
F.init(d); // empty constructor
@@ -217,7 +221,7 @@ int TestField(const Field& F, const uint64_t seed)
JEONETESTE(F,x);
for (size_t i = 0; i< NBITER; ++i) {
- while (F.isZero(g.random(x)));
+ while (F.isZero(g.random(x))) {}
JEONETESTE(F,x);
}
@@ -368,12 +372,19 @@ int main(int argc, char ** argv)
// Zech log finite field with 256 elements
// and prescribed 1 + x +x^3 +x^4 +x^8 irreducible polynomial
- std::vector< GFqDom<int64_t>::Residu_t > Irred(9);
- Irred[0] = 1; Irred[1] = 1; Irred[2] = 0; Irred[3] = 1;
- Irred[4] = 1; Irred[5] = 0; Irred[6] = 0; Irred[7] = 0;
- Irred[8] = 1;
+// std::vector< GFqDom<int64_t>::Residu_t > Irred(9);
+// Irred[0] = 1; Irred[1] = 1; Irred[2] = 0; Irred[3] = 1;
+// Irred[4] = 1; Irred[5] = 0; Irred[6] = 0; Irred[7] = 0;
+// Irred[8] = 1;
+ std::vector< int64_t > Irred {1,1,0,1,1,0,0,0,1};
TEST_SPECIFIC(GFqDom<int64_t>, GF256, 2, 8, Irred);
+ // Zech log finite field with 343 elements
+ // and prescribed 3 +x^3 irreducible polynomial
+ // and prescribed 5 +3x +4x^2 generator polynomial
+ std::vector< int64_t > Irred3 {3,0,0,1}, Gene3 {5,3,4};
+ TEST_SPECIFIC(GFqDom<int64_t>, GF343, 7, 3, Irred3, Gene3);
+
TEST_SPECIFIC(GFqDom<int32_t>, GF625, 5, 4);
TEST_SPECIFIC(GFqExt<int32_t>, GF81, 3, 4);
diff --git a/tests/test-modularmulprecomp.C b/tests/test-modularmulprecomp.C
new file mode 100644
index 0000000..e75bb18
--- /dev/null
+++ b/tests/test-modularmulprecomp.C
@@ -0,0 +1,211 @@
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: t; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
+// Copyright(c)'1994-2009,2010 by The Givaro group
+// This file is part of Givaro.
+// Copyright (C) 2016 Romain Lebreton, adapted from test-ringarith.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/givinteger.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; \
+ throw std::string( "Message d'erreur" ); \
+ return -1; \
+ }
+
+#define NBITER 100
+
+template<class Ring>
+int TestOneMulPrecomp(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;
+ typename Ring::Compute_t invp, invb, invb2;
+
+ F.assign(a, x);
+ F.assign(b, y);
+ if (a < 0){
+ std::cout << "x : " << x << "\n a : " << a << std::endl;
+ }
+
+ F.mul(c,a,b);
+
+ F.precomp_p(invp);
+ F.mul_precomp_p(d,a,b,invp);
+ TESTE_EG(c,d);
+
+ F.precomp_b(invb, b);
+ F.mul_precomp_b(d,a,b,invb);
+ TESTE_EG(c,d);
+
+ F.precomp_b(invb2, b, invp);
+ TESTE_EG(invb,invb2);
+ F.mul_precomp_b(d,a,b,invb);
+ TESTE_EG(c,d);
+
+ return 0;
+}
+
+#define JETESTE( F, x, y ) \
+ if (TestOneMulPrecomp(F,x,y)) { \
+ std::cout << #x << " " << #y << " failed !" << std::endl; \
+ return -1; \
+ }
+
+
+template<class Ring>
+int TestMulPrecomp(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);
+ JETESTE(F,x,y);
+
+ F.init(x, Ring::maxCardinality()-1U);
+ F.init(y, Ring::maxCardinality()-1U);
+ JETESTE(F,x,y);
+
+ F.assign(x, F.maxElement());
+ F.assign(y, F.maxElement());
+ JETESTE(F,x,y);
+
+ F.assign(x, F.minElement());
+ F.assign(y, F.maxElement());
+ JETESTE(F,x,y);
+
+ F.assign(x, F.minElement());
+ F.assign(y, F.minElement());
+ JETESTE(F,x,y);
+
+ for (size_t i = 0; i< NBITER; ++i) {
+ g.random(x); g.random(y);
+ JETESTE(F,x,y);
+ }
+
+ return 0;
+}
+
+#undef JETESTE
+
+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);
+
+ //--------------------------------//
+ //----- Test mult with precomp ---//
+
+ /**********************/
+ /* Modular<8, ...> */
+ /**********************/
+ Modular<int8_t, int8_t> M88p (3); //2-bit prime
+ TestMulPrecomp(M88p,seed);
+ Modular<int8_t, uint8_t> M8u8p (3);
+ TestMulPrecomp(M8u8p,seed);
+ Modular<uint8_t, int8_t> Mu88p (3);
+ TestMulPrecomp(Mu88p,seed);
+ Modular<uint8_t, uint8_t> Mu8u8p (3);
+ TestMulPrecomp(Mu8u8p,seed);
+
+ Modular<int8_t, int16_t> M816p (61); //6-bit prime
+ TestMulPrecomp(M816p,seed);
+ Modular<int8_t, uint16_t> M8u16p (61);
+ TestMulPrecomp(M8u16p,seed);
+ Modular<uint8_t, int16_t> Mu816p (61);
+ TestMulPrecomp(Mu816p,seed);
+ Modular<uint8_t, uint16_t> Mu8u16p (61);
+ TestMulPrecomp(Mu8u16p,seed);
+
+ /**********************/
+ /* Modular<16, ...> */
+ /**********************/
+ Modular<int16_t, int16_t> M1616p (61); //6-bit prime
+ TestMulPrecomp(M1616p,seed);
+ Modular<int16_t, uint16_t> M16u16p (61);
+ TestMulPrecomp(M16u16p,seed);
+ Modular<uint16_t, int16_t> Mu1616p (61);
+ TestMulPrecomp(Mu1616p,seed);
+ Modular<uint16_t, uint16_t> Mu16u16p (61);
+ TestMulPrecomp(Mu16u16p,seed);
+
+ Modular<int16_t, int32_t> M1632p (16381); //14-bit prime
+ TestMulPrecomp(M1632p,seed);
+ Modular<int16_t, uint32_t> M16u32p (16381);
+ TestMulPrecomp(M16u32p,seed);
+ Modular<uint16_t, int32_t> Mu1632p (16381);
+ TestMulPrecomp(Mu1632p,seed);
+ Modular<uint16_t, uint32_t> Mu16u32p (16381);
+ TestMulPrecomp(Mu16u32p,seed);
+
+ /**********************/
+ /* Modular<32, ...> */
+ /**********************/
+ Modular<int32_t, int32_t> M3232p (16381); //14-bit prime
+ TestMulPrecomp(M3232p,seed);
+ Modular<int32_t, uint32_t> M32u32p (16381);
+ TestMulPrecomp(M32u32p,seed);
+ Modular<uint32_t, int32_t> Mu3232p (16381);
+ TestMulPrecomp(Mu3232p,seed);
+ Modular<uint32_t, uint32_t> Mu32u32p (16381);
+ TestMulPrecomp(Mu32u32p,seed);
+
+
+ Modular<int32_t, int64_t> M3264p (1073741789); //30-bit prime
+ TestMulPrecomp(M3264p,seed);
+ Modular<int32_t, uint64_t> M32u64p (1073741789);
+ TestMulPrecomp(M32u64p,seed);
+ Modular<uint32_t, int64_t> Mu3264p (1073741789);
+ TestMulPrecomp(Mu3264p,seed);
+ Modular<uint32_t, uint64_t> Mu32u64p (1073741789);
+ TestMulPrecomp(Mu32u64p,seed);
+
+ Modular<uint32_t, uint64_t> Mu32u64p2 (143513);
+ TestMulPrecomp(Mu32u64p2,seed);
+
+
+ /**********************/
+ /* Modular<64, ...> */
+ /**********************/
+ Modular<int64_t, int64_t> M6464p (1073741789); //30-bit prime
+ TestMulPrecomp(M6464p,seed);
+ Modular<int64_t, uint64_t> M64u64p (1073741789);
+ TestMulPrecomp(M64u64p,seed);
+ Modular<uint64_t, int64_t> Mu6464p (1073741789);
+ TestMulPrecomp(Mu6464p,seed);
+ Modular<uint64_t, uint64_t> Mu64u64p (1073741789);
+ TestMulPrecomp(Mu64u64p,seed);
+
+#ifdef __GIVARO_HAVE_INT128
+ Modular<int64_t, int128_t> M64128p (4611686018427387847ul); //62-bit prime
+ TestMulPrecomp(M64128p,seed);
+ Modular<int64_t, uint128_t> M64u128p (4611686018427387847ul);
+ TestMulPrecomp(M64u128p,seed);
+ Modular<uint64_t, int128_t> Mu64128p (4611686018427387847ul);
+ TestMulPrecomp(Mu64128p,seed);
+ Modular<uint64_t, uint128_t> Mu64u128p (4611686018427387847ul);
+ TestMulPrecomp(Mu64u128p,seed);
+#endif
+
+ return 0;
+}
diff --git a/tests/test-ringarith.C b/tests/test-ringarith.C
index d67e005..f6d8cc6 100644
--- a/tests/test-ringarith.C
+++ b/tests/test-ringarith.C
@@ -13,17 +13,19 @@
#include <givaro/givinteger.h>
#include <givaro/zring.h>
#include <givaro/gfq.h>
+#include <givaro/modular-extended.h>
#include <recint/recint.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; \
+ if (!F.areEqual((a),(b))) { \
+ F.write( F.write(std::cout,a) << "!=",b) \
+ << " failed (at line " << __LINE__ << ")" << std::endl; \
return -1; \
- }
+ }
#define JETESTE( a, s ) \
if (TestRing( (a), (s)) ) { \
@@ -62,6 +64,9 @@ int TestOneRing(const Ring& F, const typename Ring::Element& x, const typename R
// F.write(std::cerr << "1: ", F.one) << std::endl;
TESTE_EG(a, F.one);
+ F.init(a, -1);
+ TESTE_EG(a, F.mOne);
+
F.assign(a, x);
F.assign(b, y);
F.init(c); // empty constructor
@@ -106,15 +111,15 @@ int TestOneRing(const Ring& F, const typename Ring::Element& x, const typename R
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;
+// 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
@@ -141,6 +146,23 @@ int TestOneRing(const Ring& F, const typename Ring::Element& x, const typename R
// 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.init(a, -3);
+ F.assign(b , y);
+ F.mul(c,a,b);
+ F.addin(c,b);
+ F.addin(c,b);
+ F.addin(c,b);
+ F.init(d, 0U);
+ TESTE_EG(c, d);
#ifdef GIVARO_DEBUG
F.write(std::cerr );
@@ -157,8 +179,8 @@ int TestOneRing(const Ring& F, const typename Ring::Element& x, const typename R
template<class Ring>
int TestRing(const Ring& F, const uint64_t seed)
{
- typename Ring::Element x, y, a;
- typename Ring::RandIter g(F, 0_ui64, seed);
+ typename Ring::Element x, y;
+ typename Ring::RandIter g(F, 0_ui64, seed);
F.init(x, 7U);
F.init(y, -29.0);
@@ -180,10 +202,10 @@ int TestRing(const Ring& F, const uint64_t seed)
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);
- }
+ for (size_t i = 0; i< NBITER; ++i) {
+ g.random(x); g.random(y);
+ JEONETESTE(F,x,y);
+ }
return 0;
}
@@ -222,6 +244,7 @@ int TestPolRing(const Ring& F, const uint64_t seed)
return 0;
}
+
int main(int argc, char ** argv)
{
auto seed = static_cast<uint64_t>(argc>1?atoi(argv[1]):BaseTimer::seed());
@@ -239,6 +262,11 @@ int main(int argc, char ** argv)
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...) \
Ring Name(Modulus); \
JETESTE(Name, seed);
@@ -260,7 +288,12 @@ int main(int argc, char ** argv)
TEST_SPECIFIC(ModularUCUS, UCUS4, 4);
TEST_SPECIFIC(ModularUSUZ, USUZ4, 4);
TEST_SPECIFIC(ModularUZULL, UZULL4, 4);
- //TEST_SPECIFIC(ModularFD, FD4, 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);
@@ -271,11 +304,11 @@ int main(int argc, char ** argv)
//--------------//
//----- 75 -----//
- TEST_SPECIFIC(Modular<int8_t>, C75, 75);
+ 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, 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);
@@ -286,6 +319,12 @@ int main(int argc, char ** argv)
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);
@@ -298,6 +337,8 @@ int main(int argc, char ** argv)
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"); \
@@ -353,7 +394,12 @@ int main(int argc, char ** argv)
TEST_LAST(ModularZULL, ZULLmax);
TEST_LAST(ModularUCUS, UCUSmax);
TEST_LAST(ModularUSUZ, USUZmax);
- TEST_LAST(ModularUZULL, UZULLmax);
+ 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);
@@ -366,6 +412,8 @@ int main(int argc, char ** argv)
TEST_LAST(Montgomery<int32_t>, MZmax);
TEST_LAST(Montgomery<RecInt::ruint128>, MRUmax);
+ TEST_LAST(ModularExtended<float>, MEFmax);
+ TEST_LAST(ModularExtended<double>, MEDmax);
return 0;
}
diff --git a/tests/test-ruint_operators.C b/tests/test-ruint_operators.C
index 74dc1db..56b53ff 100644
--- a/tests/test-ruint_operators.C
+++ b/tests/test-ruint_operators.C
@@ -92,7 +92,8 @@ int main(void)
while ((r = USItype(rand())) < 2);
z %= r; gz %= r; ruint_to_mpz(gcmp, z);
if (gcmp != gz) return 5;
- while (s == 0) rand(s); ruint_to_mpz(gs, s); s %= s; gs %= gs; ruint_to_mpz(gcmp, s);
+ while (s == 0) rand(s);
+ ruint_to_mpz(gs, s); s %= s; gs %= gs; ruint_to_mpz(gcmp, s);
if (gs != gcmp) return 5;
x *= y; gx *= gy; gx %= size; ruint_to_mpz(gcmp, x);
@@ -109,7 +110,8 @@ int main(void)
while ((r = USItype(rand())) < 2);
y /= r; gy /= r; ruint_to_mpz(gcmp, y);
if (gcmp != gy) return 7;
- while (s == 0) rand(s); ruint_to_mpz(gs, s); s /= s; gs /= gs; ruint_to_mpz(gcmp, s);
+ while (s == 0) rand(s);
+ ruint_to_mpz(gs, s); s /= s; gs /= gs; ruint_to_mpz(gcmp, s);
if (gs != gcmp) return 7;
// Refresh
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/givaro.git
More information about the debian-science-commits
mailing list