diff --git a/autotune/winograd.C b/autotune/winograd.C
index 303f5d7e5..2ab303ebb 100644
--- a/autotune/winograd.C
+++ b/autotune/winograd.C
@@ -65,7 +65,7 @@ int main () {
     typedef FIELD Field;
     Field F(17);
     typedef Field::Element Element ;
-    size_t n=512, nmax=4000, prec=512, nbest=0, count=0;
+    size_t n=512, nmax=10000, prec=512, nbest=0, count=0;
     TTimer chrono;
     bool bound=false;
 
diff --git a/benchmarks/Makefile.am b/benchmarks/Makefile.am
index 43324f644..41c7e734c 100755
--- a/benchmarks/Makefile.am
+++ b/benchmarks/Makefile.am
@@ -33,7 +33,7 @@ endif
 
 PERFPUBLISHERFILE=benchmarks-report.xml
 
-FFLA_BENCH =    benchmark-fgemm benchmark-fgemm-rns benchmark-wino benchmark-ftrsm  benchmark-fgesv benchmark-ftrsv benchmark-ftrtri  benchmark-inverse benchmark-fsytrf benchmark-fsyrk  benchmark-lqup benchmark-fsyr2k benchmark-pluq benchmark-charpoly benchmark-charpoly-mp benchmark-fgemm-mp benchmark-fgemv-mp benchmark-ftrsm-mp benchmark-lqup-mp benchmark-checkers benchmark-fadd-lvl2 benchmark-fdot benchmark-fgemv benchmark-quasisep benchmark-sss benchmark-storage-transpose benchmark-qscomp
+FFLA_BENCH =    benchmark-fgemm benchmark-fgemm-rns benchmark-wino benchmark-ftrsm  benchmark-fgesv benchmark-ftrsv benchmark-ftrtri  benchmark-inverse benchmark-fsytrf benchmark-fsyrk  benchmark-lqup benchmark-fsyr2k benchmark-pluq benchmark-charpoly benchmark-charpoly-mp benchmark-fgemm-mp benchmark-fgemv-mp benchmark-ftrsm-mp benchmark-lqup-mp benchmark-checkers benchmark-fadd-lvl2 benchmark-fdot benchmark-fgemv benchmark-quasisep benchmark-sss benchmark-storage-transpose benchmark-qscomp benchmark-rns-double
 
 BLAS_BENCH =    benchmark-sgemm$(EXEEXT) benchmark-dgemm benchmark-dtrsm
 LAPA_BENCH =    benchmark-dtrtri benchmark-dgetri benchmark-dgetrf benchmark-dsytrf
@@ -87,6 +87,7 @@ benchmark_fadd_lvl2_SOURCES = benchmark-fadd-lvl2.C
 benchmark_fdot_SOURCES = benchmark-fdot.C
 benchmark_fgemv_SOURCES = benchmark-fgemv.C
 benchmark_storage_transpose_SOURCES = benchmark-storage-transpose.C
+benchmark_rns_double_SOURCES = benchmark-rns-double.C
 
 benchmark_sgemm_CXXFLAGS = $(AM_CXXFLAGS) -D__SGEMM__
 
diff --git a/benchmarks/benchmark-dtrsm.C b/benchmarks/benchmark-dtrsm.C
index 576835957..5ffd80ce6 100644
--- a/benchmarks/benchmark-dtrsm.C
+++ b/benchmarks/benchmark-dtrsm.C
@@ -104,7 +104,7 @@ int main(int argc, char** argv) {
     // -----------
     // Standard output for benchmark - Alexis Breust 2014/11/14
     std::cout << "Time: " << time / double(iter)
-    << " Gfops: " << (2.*double(n)/1000.*double(n)/1000.*double(n)/1000.0) / time * double(iter) / 3.;
+    << " Gfops: " << (1.*double(n)/1000.*double(n)/1000.*double(n)/1000.0) / time * double(iter) ;
     FFLAS::writeCommandString(std::cout, as) << std::endl;
 
     return 0;
diff --git a/benchmarks/benchmark-fgemm-mp.C b/benchmarks/benchmark-fgemm-mp.C
index 4d349cc06..90058dcab 100644
--- a/benchmarks/benchmark-fgemm-mp.C
+++ b/benchmarks/benchmark-fgemm-mp.C
@@ -28,7 +28,8 @@
 // everywhere in the call stack
 #define __FFLASFFPACK_OPENBLAS_NT_ALREADY_SET 1
 
-#define INTEGER_NO_RNS 1
+//#define INTEGER_NO_RNS 1
+#define PROFILE_FGEMM_MP
 
 #if not defined(MG_DEFAULT)
 #define MG_DEFAULT MG_ACTIVE
@@ -88,7 +89,7 @@ static Argument as[] = {
 };
 
 template<typename Ints>
-int tmain(){
+int tmain(){ 
     srand( (int)seed);
     srand48(seed);
     Givaro::Integer::seeding(seed);
@@ -183,17 +184,16 @@ int tmain(){
         fmpz_mat_clear(BB);
 #endif
         //END FLINT CODE //
-        using  FFLAS::CuttingStrategy::Recursive;
-        using  FFLAS::StrategyParameter::TwoDAdaptive;
+        //using  FFLAS::CuttingStrategy::Recursive;
+        //using  FFLAS::StrategyParameter::TwoDAdaptive;
         // RNS MUL_LA
         chrono.clear();chrono.start();
         // 		PAR_BLOCK{
         //             FFLAS::fgemm(F,FFLAS::FflasNoTrans,FFLAS::FflasNoTrans,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc, SPLITTER(NUM_THREADS,Recursive,TwoDAdaptive) );
         // 		}
         {
-            FFLAS::fgemm(F,FFLAS::FflasNoTrans,FFLAS::FflasNoTrans,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc,FFLAS::ParSeqHelper::Sequential());
+          FFLAS::fgemm(F,FFLAS::FflasNoTrans,FFLAS::FflasNoTrans,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc,FFLAS::ParSeqHelper::Sequential());
         }
-
         chrono.stop();
         time+=chrono.realtime();
 
@@ -234,7 +234,7 @@ int tmain(){
 int main(int argc, char** argv){
   
 #ifdef __FFLASFFPACK_OPENBLAS_NUM_THREADS
-    openblas_set_num_threads(__FFLASFFPACK_OPENBLAS_NUM_THREADS);
+   openblas_set_num_threads(__FFLASFFPACK_OPENBLAS_NUM_THREADS);
 #endif
 
     FFLAS::parseArguments(argc,argv,as);
diff --git a/benchmarks/benchmark-fgemm.C b/benchmarks/benchmark-fgemm.C
index 91ec24f94..3378b1537 100644
--- a/benchmarks/benchmark-fgemm.C
+++ b/benchmarks/benchmark-fgemm.C
@@ -28,8 +28,8 @@
 //#define WINO_PARALLEL_TMPS
 //#define __FFLASFFPACK_FORCE_SEQ
 //#define PFGEMM_WINO_SEQ 32
-//#define CLASSIC_SEQ
-#define CLASSIC_HYBRID
+#define CLASSIC_SEQ
+//#define CLASSIC_HYBRID
 //#define WINO_SEQ
 //#define FFT_PROFILER
 //#define PROFILE_FGEMM_MP
@@ -82,15 +82,18 @@ int main(int argc, char** argv) {
     //  typedef Givaro::Modular<Givaro::Integer> Field;
     // typedef Givaro::Modular<int64_t> Field;
     typedef Givaro::Modular<double> Field;
-    //  typedef Givaro::Modular<float> Field;
+    //typedef Givaro::Modular<float> Field;
     //  typedef Givaro::ModularBalanced<float> Field;
-    //  typedef Givaro::ModularBalanced<double> Field;
+    //typedef Givaro::ModularBalanced<double> Field;
     // typedef Givaro::ModularBalanced<int64_t> Field;
     //  typedef Givaro::Modular<Givaro::Integer> Field;
+    //typedef Givaro::DoubleDomain Field;
+
+
     typedef Field::Element Element;
 
     Field F(q);
-    if (q > F.maxCardinality()) return 1;
+    //if (q > F.maxCardinality()) return 1;
 
     Timer chrono, TimFreivalds;
     double timev=0.0;
@@ -99,6 +102,8 @@ int main(int argc, char** argv) {
     Element * A, * B, * C;
 
     Field::RandIter G(F);
+    //G.setBitsize(20);
+
     A = fflas_new(F,m,k,Alignment::CACHE_PAGESIZE);
     //#pragma omp parallel for collapse(2) schedule(runtime)
     PAR_BLOCK { pfrand(F,G, m,k,A,m/size_t(NBK)); }
@@ -112,6 +117,7 @@ int main(int argc, char** argv) {
     PAR_BLOCK { pfzero(F, m,n,C,m/NBK); }
 
 
+
     for (size_t i=0;i<=iter;++i){
 
         chrono.clear();
@@ -205,7 +211,7 @@ int main(int argc, char** argv) {
             }
             else{
 
-                MMHelper<Field,MMHelperAlgo::Winograd>//,
+              MMHelper<Field,MMHelperAlgo::Winograd>//,
                 //typename FieldTraits<Field>::value,
                 //ParSeqHelper::Sequential>
                 WH (F, nbw, ParSeqHelper::Sequential());
diff --git a/benchmarks/benchmark-ftrsm.C b/benchmarks/benchmark-ftrsm.C
index 44d85255e..e239b955e 100644
--- a/benchmarks/benchmark-ftrsm.C
+++ b/benchmarks/benchmark-ftrsm.C
@@ -81,7 +81,7 @@ int main(int argc, char** argv) {
         FFLAS::ReadMatrix (file1.c_str(),F,m,m,A);
     }
     else{
-        A = FFLAS::fflas_new (F,m,m,Alignment::CACHE_PAGESIZE);
+      A = FFLAS::fflas_new (F,m,m,Alignment::CACHE_PAGESIZE);
         PAR_BLOCK{ FFLAS::pfrand(F,G,m,m,A,m/NBK); }
 
         for (size_t k=0;k<(size_t)m;++k)
@@ -102,10 +102,10 @@ for (size_t i=0;i<=iter;++i){
     if (i) chrono.start();
 
     if (!p){
-        FFLAS::ParSeqHelper::Sequential H;
-        FFLAS::ftrsm (F, FFLAS::FflasLeft, FFLAS::FflasLower,
-                      FFLAS::FflasNoTrans, FFLAS::FflasNonUnit,
-                      m,n, F.one, A, m, B, n, H);
+      FFLAS::ParSeqHelper::Sequential H;
+      FFLAS::ftrsm (F, FFLAS::FflasLeft, FFLAS::FflasLower,
+                    FFLAS::FflasNoTrans, FFLAS::FflasNonUnit,
+                    m,n, F.one, A, m, B, n, H);
     }
     else{
         FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Block,FFLAS::StrategyParameter::Threads> PSH(t);
diff --git a/benchmarks/benchmark-pluq.C b/benchmarks/benchmark-pluq.C
index 07f19979a..19516a4b1 100644
--- a/benchmarks/benchmark-pluq.C
+++ b/benchmarks/benchmark-pluq.C
@@ -269,7 +269,7 @@ int main(int argc, char** argv) {
             if (slab)
                 R = FFPACK::LUdivine (F, diag, FFLAS::FflasNoTrans, m, n, A, n, P, Q);
             else
-                R = FFPACK::PLUQ(F, diag, m, n, A, n, P, Q);
+              R = FFPACK::PLUQ(F, diag, m, n, A, n, P, Q);
         }
         if (i) {chrono.stop(); time[i-1]=chrono.realtime();}
 
diff --git a/benchmarks/benchmark-rns-double.C b/benchmarks/benchmark-rns-double.C
new file mode 100644
index 000000000..20be8f688
--- /dev/null
+++ b/benchmarks/benchmark-rns-double.C
@@ -0,0 +1,335 @@
+/*
+ * Copyright (C) 2023 the FFLAS-FFPACK group
+ *
+ * Written by Pascal Giorgi <pascal.giorgi@lirmm.fr>
+ *
+ * This file is Free Software and part of FFLAS-FFPACK.
+ *
+ * ========LICENCE========
+ * This file is part of the library FFLAS-FFPACK.
+ *
+ * FFLAS-FFPACK is free software: you can redistribute it and/or modify
+ * it under the terms of the  GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * This library 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
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+ * ========LICENCE========
+ *.
+ */
+
+//#define CHECK_RNS
+//#define RNS_DEBUG
+//#define BENCH_RNS
+ 
+#define __FFLASFFPACK_SEQUENTIAL
+
+#include "fflas-ffpack/fflas-ffpack-config.h"
+
+
+#include <iostream>
+#include <vector>
+#include <string>
+#include <chrono>
+#include <algorithm>
+using namespace std; 
+
+#include "fflas-ffpack/utils/timer.h"
+#include "fflas-ffpack/fflas/fflas.h"
+#include "fflas-ffpack/field/rns-double.h"
+#include "fflas-ffpack/field/rns-double-extended.h"
+#include "fflas-ffpack/field/rns-integer.h"
+#include "fflas-ffpack/utils/args-parser.h"
+#include "givaro/givinteger.h"
+#include "givaro/modular-integer.h"
+#include "givaro/zring.h"
+#include <givaro/givrns.h>
+
+//#define BENCH_FLINT
+#ifdef BENCH_FLINT
+#define __GMP_BITS_PER_MP_LIMB 64
+extern "C" {
+#include "flint/longlong.h"
+#include "flint/long_extras.h" 
+#include "flint/fmpz_mat.h"
+#include "flint/fmpz.h"
+#include "flint/flint.h"
+}
+#endif
+
+
+#define LOOPS_TIME 1
+
+template <typename Field>
+void run_bench(size_t n, size_t primes_bits, size_t b, size_t iters, size_t seed, bool matmul=false){
+	Givaro::Integer p;
+	FFLAS::Timer chrono;
+	double timeFFLASToRNS=0.,timeFFLASFromRNS=0.,timeFFLASPrecomp=0.;
+	double timeFFLASToRNS_ext=0.,timeFFLASFromRNS_ext=0.,timeFFLASPrecomp_ext=0.;
+#ifdef BENCH_FLINT
+	double timeFlintToRNS=0., timeFlintToRNSnaive=0.,timeFlintFromRNS=0.,timeFlintFromRNSnaive=0.,timeFlintPrecomp=0.;
+#endif
+	size_t loop=0;
+	size_t bits=b;
+	int logb= log((double)n)/log(2.0);
+
+	size_t fflas_primes_bits       = (primes_bits ? primes_bits : 41-logb);
+	size_t fflas_primes_bits_ext   = 2* fflas_primes_bits;
+	// std::cout<<"FFLAS     primes_bitsize= "<<fflas_primes_bits<<std::endl;
+	// std::cout<<"FFLAS EXT primes_bitsize= "<<fflas_primes_bits_ext<<std::endl;
+
+#ifdef BENCH_FLINT
+	const size_t flint_primes_bits = (primes_bits ? primes_bits :NMOD_MAT_OPTIMAL_MODULUS_BITS);
+	//std::cout<<"FLINT prime bit is "<<flint_primes_bits<<std::endl;
+#endif
+
+	auto start = std::chrono::system_clock::now();
+	auto end = std::chrono::system_clock::now();
+
+	size_t rns_bitsize;
+	Givaro::Integer rns_bound;
+	for (;( loop<iters || std::chrono::duration<double>(end-start) < std::chrono::seconds(LOOPS_TIME)) ;loop++){
+
+		Givaro::Integer::random_exact_2exp(p, bits);
+		//nextprime(p,p); avoid p to b a prime has it is not needed here
+		Field F(p);
+		typename Field::RandIter Rand(F,0,(uint64_t)seed);
+		typename Field::Element_ptr A;
+		A= FFLAS::fflas_new(F,n,n);
+
+		FFLAS::frand(F,Rand,n,n,A,n);
+
+		if (matmul){
+			rns_bound=p*p*(uint64_t)n;
+			rns_bitsize=rns_bound.bitsize();
+		}else {
+			rns_bound=p;
+			rns_bitsize=rns_bound.bitsize();
+		}
+
+		int logb= log(double(rns_bitsize))/log(2.); // log of the integer bitsize of entry in A
+
+		fflas_primes_bits       = (primes_bits ? primes_bits : std::min(25,41-logb));
+		fflas_primes_bits_ext   =  2* fflas_primes_bits;
+
+
+
+#ifdef BENCH_FLINT
+		// FLINT RNS CODE //
+		{
+			fmpz_mat_t AA;
+			fmpz_mat_init (AA, n, n);
+			flint_rand_t randstate;
+			flint_randinit(randstate);
+			fmpz_mat_randbits(AA, randstate, bits);
+
+			fmpz_comb_t comb;
+			fmpz_comb_temp_t comb_temp;
+			nmod_mat_t * mod_A;
+			mp_limb_t * primes;
+			mp_limb_t * residues;    
+			size_t num_primes = (rns_bitsize + flint_primes_bits - 1) / flint_primes_bits;
+
+			/* FLINT RNS initialization */
+			residues = (mp_limb_t *)flint_malloc(sizeof(mp_limb_t) * num_primes);			
+			mod_A    = (nmod_mat_t *)flint_malloc(sizeof(nmod_mat_t) * num_primes);
+			primes   = (mp_limb_t *)flint_malloc(sizeof(mp_limb_t) * num_primes);
+			primes[0] = n_nextprime(UWORD(1) << flint_primes_bits, 0);
+			nmod_mat_init(mod_A[0], AA->r, AA->c, primes[0]);
+			for (size_t i = 1; i < num_primes; i++)
+				{
+					primes[i] = n_nextprime(primes[i-1], 0);
+					nmod_mat_init(mod_A[i], AA->r, AA->c, primes[i]);
+				}
+
+			chrono.clear();chrono.start();
+			fmpz_comb_init(comb, primes, num_primes);
+			fmpz_comb_temp_init(comb_temp, comb);
+			chrono.stop();
+			timeFlintPrecomp+=chrono.usertime();
+			chrono.clear();chrono.start();
+
+			/* Calculate residues of AA */
+			fmpz_mat_multi_mod_ui_precomp(mod_A, num_primes, AA, comb, comb_temp); 
+
+			chrono.stop();
+			timeFlintToRNS+=chrono.usertime();
+			chrono.clear();chrono.start();
+
+
+			/* Calculate residue without COMB approach (NO FAST RNS) */
+			//fmpz_mat_multi_mod_ui(mod_A, num_primes, AA);
+			// for (long i = 0; i < AA->r; i++)
+			// 	{   for (long j = 0; j < AA->c; j++)
+			// 			for (long l = 0; l < num_primes; l++)
+			// 				nmod_mat_entry(mod_A[l],i,j) =  fmpz_fdiv_ui(fmpz_mat_entry(AA,i,j), primes[l]);
+			// 		}
+
+			chrono.stop();
+			timeFlintToRNSnaive+=chrono.usertime();
+			chrono.clear();chrono.start();
+			/* Chinese remaindering */
+			for (long i = 0; i < AA->r; i++)
+				{   for (long j = 0; j < AA->c; j++)
+						{   for (size_t k = 0; k < num_primes; k++)
+								residues[k] = mod_A[k]->rows[i][j];
+							fmpz_multi_CRT_ui(&AA->rows[i][j], residues, comb, comb_temp, 1);
+						}
+				}
+			chrono.stop();
+			timeFlintFromRNS+=chrono.usertime();
+
+			fmpz_mat_clear(AA);
+			flint_free(mod_A);
+			fmpz_comb_temp_clear(comb_temp);
+			fmpz_comb_clear(comb);
+			flint_free(residues);
+			flint_free(primes);
+
+		}
+		//END FLINT CODE //
+#endif
+
+		// FFLAS RNS DOUBLE CODE
+		if (fflas_primes_bits<=41-logb){
+			chrono.clear();	chrono.start();
+			// construct an RNS structure and its associated Domain
+			FFPACK::rns_double RNS(rns_bound, fflas_primes_bits,false,seed);
+			typedef FFPACK::RNSInteger<FFPACK::rns_double> RnsDomain;
+			RnsDomain Zrns(RNS);
+			typename RnsDomain::Element_ptr mod_A = FFLAS::fflas_new(Zrns,n,n);
+			chrono.stop();
+			timeFFLASPrecomp+=chrono.usertime();
+			chrono.clear();	chrono.start();
+			RNS.init(n,n,mod_A._ptr,mod_A._stride,A,n,p);
+			chrono.stop();
+			timeFFLASToRNS+=chrono.usertime();
+			chrono.clear();	chrono.start();
+			RNS.convert(n,n,Givaro::Integer(0), A,n, mod_A._ptr,mod_A._stride);
+			//FFLAS::fconvert_rns(Zrns,n,n,Givaro::Integer(0),A,n,mod_A);
+			chrono.stop();
+			timeFFLASFromRNS+=chrono.usertime();
+			FFLAS::fflas_delete(mod_A);
+		}
+
+		// FFLAS RNS DOUBLE EXTENDED CODE
+		if (fflas_primes_bits_ext<=2*(41-logb)){
+			chrono.clear();	chrono.start();	
+			// construct an RNS structure and its associated Domain
+			typedef FFPACK::rns_double_extended RNSExt;
+			RNSExt RNS(rns_bound, fflas_primes_bits_ext,false,seed);
+			typedef FFPACK::RNSInteger<RNSExt> RnsDomain;
+			RnsDomain Zrns(RNS);
+			typename RnsDomain::Element_ptr mod_A = FFLAS::fflas_new(Zrns,n,n);
+			chrono.stop();
+			timeFFLASPrecomp_ext+=chrono.usertime();
+			chrono.clear();	chrono.start();
+
+			RNS.init(n,n,mod_A._ptr,mod_A._stride,A,n,p);
+			chrono.stop();
+			timeFFLASToRNS_ext+=chrono.usertime();
+			chrono.clear();	chrono.start();	
+			//FFLAS::fconvert_rns(Zrns,n,n,Givaro::Integer(0),A,n,mod_A);
+			RNS.convert(n,n,Givaro::Integer(0), A,n, mod_A._ptr,mod_A._stride);
+			chrono.stop();
+			timeFFLASFromRNS_ext+=chrono.usertime();
+			FFLAS::fflas_delete(mod_A);
+		}
+		FFLAS::fflas_delete(A);
+		end = std::chrono::system_clock::now();
+	}
+
+#define SPC1 10
+#define SPC2 12
+#define SPC3 33
+#define PREC 2
+
+	cout<<setw(SPC1+5)
+	    <<loop<<" |"
+	    <<setw(SPC1)
+	    <<n<<" |"
+	    <<setw(SPC1)
+	    <<bits<<" |";
+	   
+
+	//#define TIMING(x,y,z) cout<<setw(SPC2)<<std::fixed << std::setprecision(PREC)<<(x+y+z)/loop<<setw(SPC3)<<"("+to_string(x/loop)+","+to_string(y/loop)+","+to_string(z/loop)+")"<<" |";
+
+#define TIMING(b,x,y,z) cout<<setw(SPC1)<<b<<std::fixed << std::setprecision(PREC) \
+			  <<setw(SPC3/3)<<((x/loop)*1000000.)\
+			  <<setw(SPC3/3)<<((y/loop)*(1000000./(n*n)))\
+			  <<setw(SPC3/3)<<((z/loop)*(1000000./(n*n)))<<" |";
+
+
+
+	TIMING(fflas_primes_bits, timeFFLASPrecomp,timeFFLASFromRNS,timeFFLASToRNS);
+	TIMING(fflas_primes_bits_ext, timeFFLASPrecomp_ext,timeFFLASFromRNS_ext,timeFFLASToRNS_ext);
+#ifdef BECNH_FLINT
+	TIMING(flint_primes_bits,timeFlintPrecomp,timeFlintFromRNS,timeFlintToRNS);
+#endif
+	cout<<endl;
+}
+
+
+
+int main(int argc, char** argv){
+	static size_t iters = 10 ;
+	static unsigned long b = 512 ;
+	static unsigned long p = 0 ;
+	static size_t m = 128 ;
+	static bool longbench=false;
+	static bool matmul=false;
+	static size_t seed= time(NULL);	
+	static Argument as[] = {
+		{ 'b', "-b B", "Set the bitsize of the matrix entries.",         TYPE_INT , &b },
+		{ 'p', "-p B", "Set the bitsize of the RNS Prime.",              TYPE_INT , &p },
+		{ 'm', "-m M", "Set the dimension m of the matrix.",             TYPE_INT , &m },
+		{ 'i', "-i R", "Set minimal number of repetitions.",             TYPE_INT , &iters },
+		{ 'l', "-l Y", "Running a long benchmark .",                     TYPE_BOOL, &longbench },
+		{ 't', "-t Y", "Running a matmul size benchmark .",              TYPE_BOOL, &matmul },
+		{ 's', "-s S", "Set the seeding value",                          TYPE_INT , &seed },
+		END_OF_ARGUMENTS
+	};
+	FFLAS::parseArguments(argc,argv,as);
+
+
+
+
+	typedef Givaro::Modular<Givaro::Integer> Field;
+	
+	cout<<"### running RNS conversions benchmark ###"<<endl;
+	if (matmul)
+		cout<<"### matrix multiplication parameters are applied: rns bitsize is  (log n + 2 bitsize)"<<endl; 
+	cout<<"### timing in microseconds (precomp,from, to) [from/to are given per element]"<<endl;
+	cout<<"### |"
+	    <<setw(SPC1)<<"loop"<<" |"
+	    <<setw(SPC1)<<"matrix dim"<<" |"
+	    <<setw(SPC1)<<"bitsize"<<" |"
+	    <<setw(SPC3+SPC1)<<"FFLAS (log2(mi), precomp, from, to )"<<" |"
+	    <<setw(SPC3+SPC1)<<"FFLAS_EXT (log2(mi), precomp, from, to)"<<" |"
+#ifdef BECNH_FLINT
+	    <<setw(SPC3)<<"FLINT (log2(mi), precomp, from, to)"<<" |"
+#endif
+	    <<endl;
+
+	if (!longbench)
+		run_bench<Field>(m,p,b,iters,seed,matmul);
+	else {
+		size_t matdim=m;	
+		for (size_t bits=128; bits<(1<<20);bits<<=1)
+			run_bench<Field>(matdim, p, bits, iters,seed++,matmul);
+
+	}
+	return 0;
+}
+
+
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
+
diff --git a/examples/matmul.C b/examples/matmul.C
index 598cb957b..6dc076ef6 100644
--- a/examples/matmul.C
+++ b/examples/matmul.C
@@ -31,7 +31,7 @@
 int main(int argc, char** argv)
 {
     if (argc != 4) {
-        std::cerr << "Usage: matmul <p> <matrixA> <matrixB>" << std::endl;
+      std::cerr << "Usage: matmul <p> <matrixA> <matrixB>" << std::endl;
         return -1;
     }
 
diff --git a/fflas-ffpack/config-blas.h b/fflas-ffpack/config-blas.h
index e503ab382..ce4932b19 100644
--- a/fflas-ffpack/config-blas.h
+++ b/fflas-ffpack/config-blas.h
@@ -425,7 +425,8 @@ static const char* EXT_BLAS_DIAG         (CBLAS_DIAG t)      { if (t == CblasUni
 // define external link to LAPACK routines
 extern "C" {
     //!@bug we should also allow lapacke from MLK
-    void dgetrf_ (const CBLAS_INT *, const CBLAS_INT *, double *, const CBLAS_INT *, CBLAS_INT *, CBLAS_INT *);
+  void dgetrf_ (const CBLAS_INT *, const CBLAS_INT *, double *, const CBLAS_INT *, CBLAS_INT *, CBLAS_INT *);
+  void sgetrf_ (const CBLAS_INT *, const CBLAS_INT *, float *, const CBLAS_INT *, CBLAS_INT *, CBLAS_INT *);
     void dgetri_ (const CBLAS_INT *, double *, const CBLAS_INT *, const CBLAS_INT *, double *, const CBLAS_INT *, CBLAS_INT *);
     void dtrtri_ (const char *, const char *, const CBLAS_INT *, double *, const CBLAS_INT *, CBLAS_INT *);
     void dswap_ (const CBLAS_INT *, double *, const CBLAS_INT *, double *, const CBLAS_INT *);
@@ -447,13 +448,23 @@ extern "C" {
     // return A=P.L.U (L unitary) with ColMajor
     // return A=L.U.P (U unitary) with RowMajor
     //! @bug Order is not used. we should use ATLAS/interfaces/lapack/C/src/clapack_dgetrf.c or similar
-    inline CBLAS_INT clapack_dgetrf(const blas_enum CBLAS_ORDER, const CBLAS_INT M, const CBLAS_INT N,
-                                    double *A, const CBLAS_INT lda, CBLAS_INT *ipiv)
-    {
-        CBLAS_INT info;
-        dgetrf_ ( &M, &N, A, &lda, ipiv, &info);
-        return info;
-    }
+  inline CBLAS_INT clapack_dgetrf(const blas_enum CBLAS_ORDER, const CBLAS_INT M, const CBLAS_INT N,
+                                  double *A, const CBLAS_INT lda, CBLAS_INT *ipiv)
+  {
+    CBLAS_INT info;
+    dgetrf_ ( &M, &N, A, &lda, ipiv, &info);
+    return info;
+  }
+
+  inline CBLAS_INT clapack_sgetrf(const blas_enum CBLAS_ORDER, const CBLAS_INT M, const CBLAS_INT N,
+                                  float *A, const CBLAS_INT lda, CBLAS_INT *ipiv)
+  {
+    CBLAS_INT info;
+    sgetrf_ ( &M, &N, A, &lda, ipiv, &info);
+    return info;
+  }
+
+
 
     inline CBLAS_INT clapack_dgetri(const blas_enum CBLAS_ORDER, const CBLAS_INT N, double *A,
                                     const CBLAS_INT lda, const CBLAS_INT *ipiv)
@@ -517,13 +528,12 @@ extern "C" {
 extern "C" {
     // LAPACK routines
 
-    CBLAS_INT clapack_dgetrf(const blas_enum CBLAS_ORDER Order, const CBLAS_INT M, const CBLAS_INT N,
-                             double *A, const CBLAS_INT lda, CBLAS_INT *ipiv);
+  CBLAS_INT clapack_dgetrf(const blas_enum CBLAS_ORDER Order, const CBLAS_INT M, const CBLAS_INT N,
+                           double *A, const CBLAS_INT lda, CBLAS_INT *ipiv);
     CBLAS_INT clapack_dgetri(const blas_enum CBLAS_ORDER Order, const CBLAS_INT N, double *A,
                              const CBLAS_INT lda, const CBLAS_INT *ipiv);
     CBLAS_INT clapack_dtrtri(const blas_enum CBLAS_ORDER Order,const blas_enum CBLAS_UPLO Uplo,
                              const blas_enum CBLAS_DIAG Diag,const CBLAS_INT N, double *A, const CBLAS_INT lda);
-
 }
 #endif // CLAPACK ?
 
diff --git a/fflas-ffpack/fflas/fflas.h b/fflas-ffpack/fflas/fflas.h
index efd71b7e2..f545d5a84 100644
--- a/fflas-ffpack/fflas/fflas.h
+++ b/fflas-ffpack/fflas/fflas.h
@@ -103,6 +103,7 @@
 #include "fflas_fadd.h"
 #include "fflas_fscal.h"
 #include "fflas_fassign.h"
+#include "fflas_transpose.h"
 
 #include "fflas_fgemm.inl"
 #include "fflas_pfgemm.inl"
@@ -121,6 +122,7 @@
 #include "fflas_faxpy.inl"
 #include "fflas_fdot.inl"
 
+
 //---------------------------------------------------------------------
 // MultiPrecision routines
 //---------------------------------------------------------------------
diff --git a/fflas-ffpack/fflas/fflas_fadd.h b/fflas-ffpack/fflas/fflas_fadd.h
index b693ff13a..cb0416706 100644
--- a/fflas-ffpack/fflas/fflas_fadd.h
+++ b/fflas-ffpack/fflas/fflas_fadd.h
@@ -64,6 +64,7 @@ namespace FFLAS {
           typename Field::ConstElement_ptr B, const size_t incb,
           typename Field::Element_ptr C, const size_t incc)
     {
+
         details::fadd<Field, true>(F,N,A,inca,B,incb,C,incc
                                    , typename FieldTraits<Field>::category() );
     }
diff --git a/fflas-ffpack/fflas/fflas_fadd.inl b/fflas-ffpack/fflas/fflas_fadd.inl
index 917cb0207..b7ebb9e67 100644
--- a/fflas-ffpack/fflas/fflas_fadd.inl
+++ b/fflas-ffpack/fflas/fflas_fadd.inl
@@ -353,7 +353,7 @@ namespace FFLAS { namespace details {
          )
     {
         if (inca == 1 && incb == 1 && incc == 1) {
-            for (size_t i=0; i<N; i++) {
+          for (size_t i=0; i<N; i++) {
                 if (ADD)
                     F.add (C[i], A[i], B[i]);
                 else
@@ -367,7 +367,7 @@ namespace FFLAS { namespace details {
                 if (ADD)
                     F.add (*Ci, *Ai, *Bi);
                 else
-                    F.add (*Ci, *Ai, *Bi);
+                    F.add (*Ci, *Ai, *Bi);            
         }
     }
 
@@ -397,7 +397,7 @@ namespace FFLAS { namespace details {
           typename Field::Element_ptr C, const size_t incc
           , FieldCategories::UnparametricTag
          )
-    {
+    {      
         if (inca == 1 && incb == 1 && incc == 1) {
             if (ADD)
                 FFLAS::vectorised::add(C,A,B,N);
diff --git a/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical_mp.inl b/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical_mp.inl
index 14a2ac3a1..2cd44decb 100644
--- a/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical_mp.inl
+++ b/fflas-ffpack/fflas/fflas_fgemm/fgemm_classical_mp.inl
@@ -322,7 +322,8 @@ namespace FFLAS {
            Givaro::Integer* C, const size_t ldc,
            MMHelper<Givaro::ZRing<Givaro::Integer>, MMHelperAlgo::Classic, ModeCategories::ConvertTo<ElementCategories::RNSElementTag>, ParSeq >  & H)
     {
-        //std::cerr<<"Entering fgemm<ZRing<Integer>> ParSeq"<<std::endl;
+      //std::cerr<<"Entering fgemm<ZRing<Integer>> ParSeq"<<std::endl;
+      //std::cerr<<A<<" "<<B<<std::endl;
 #ifdef PROFILE_FGEMM_MP
         Timer chrono;
         chrono.start();
@@ -376,7 +377,7 @@ namespace FFLAS {
 #ifdef PROFILE_FGEMM_MP
         chrono.stop();
         std::cout<<"-------------------------------"<<std::endl;
-        std::cout<<"FGEMM_MP: nb prime: "<<RNS._size<<std::endl;
+        std::cout<<"FGEMM_MP: nb prime: "<<RNS._size<< " (bitsize<"<<prime_bitsize<<")"<<std::endl;
         std::cout<<"FGEMM_MP:     init: "<<uint64_t(chrono.realtime()*1000)<<"ms"<<std::endl;
         chrono.start();
 #endif
@@ -510,7 +511,7 @@ namespace FFLAS {
                                    MMHelper<Givaro::Modular<Givaro::Integer>, MMHelperAlgo::Classic, ModeCategories::ConvertTo<ElementCategories::RNSElementTag> > & H)
     {
         // compute the product over Z
-        //std::cerr<<"Entering fgemm<Modular<Integer>>"<<std::endl;
+      //std::cerr<<"Entering fgemm<Modular<Integer>>"<<std::endl;
         typedef Givaro::ZRing<Givaro::Integer> IntegerDomain;
         Givaro::Integer p;
         F.cardinality(p);
diff --git a/fflas-ffpack/fflas/fflas_fscal.inl b/fflas-ffpack/fflas/fflas_fscal.inl
index 257b3734c..4beac1378 100644
--- a/fflas-ffpack/fflas/fflas_fscal.inl
+++ b/fflas-ffpack/fflas/fflas_fscal.inl
@@ -273,7 +273,7 @@ namespace FFLAS { namespace details {
     inline void
     fscalin (const Field& F, const size_t n, const typename Field::Element a,
              typename Field::Element_ptr X, const size_t incX, FC)
-    {
+    {  
         if (F.isOne(a))
             return ;
 
@@ -344,7 +344,7 @@ namespace FFLAS {
     fscalin (const Field& F, const size_t n, const typename Field::Element a,
              typename Field::Element_ptr X, const size_t incX)
     {
-        details::fscalin(F,n,a,X,incX,typename FieldTraits<Field>::category());
+      details::fscalin(F,n,a,X,incX,typename FieldTraits<Field>::category());
     }
 
     template<class Field>
diff --git a/fflas-ffpack/fflas/fflas_helpers.inl b/fflas-ffpack/fflas/fflas_helpers.inl
index 85478d4d2..15c56a31d 100644
--- a/fflas-ffpack/fflas/fflas_helpers.inl
+++ b/fflas-ffpack/fflas/fflas_helpers.inl
@@ -195,7 +195,7 @@ namespace FFLAS {
         }
         bool Aunfit(){ return Protected::unfit(std::max(static_cast<const DFElt&>(-Amin),Amax));}
         bool Bunfit(){ return Protected::unfit(std::max(static_cast<const DFElt&>(-Bmin),Bmax));}
-        void setOutBounds(const size_t k, const DFElt alpha, const DFElt beta)
+        void setOutBounds(const uint64_t k, const DFElt alpha, const DFElt beta)
         {
             if (beta<0){
                 Outmin = beta*Cmax;
diff --git a/fflas-ffpack/fflas/fflas_level1.inl b/fflas-ffpack/fflas/fflas_level1.inl
index 0e3c94f4c..f1af3bf5f 100644
--- a/fflas-ffpack/fflas/fflas_level1.inl
+++ b/fflas-ffpack/fflas/fflas_level1.inl
@@ -117,7 +117,7 @@ namespace FFLAS {
         typename Field::ConstElement_ptr Yi = Y ;
         for (; Xi < X+n*incX; Xi+=incX, Yi += incY )
             F.convert( *Xi , *Yi);
-    }
+    } 
 
     /** fnegin
      * \f$x \gets - x\f$.
diff --git a/fflas-ffpack/fflas/fflas_simd/simd128.inl b/fflas-ffpack/fflas/fflas_simd/simd128.inl
index df87a9ae3..09c14904a 100644
--- a/fflas-ffpack/fflas/fflas_simd/simd128.inl
+++ b/fflas-ffpack/fflas/fflas_simd/simd128.inl
@@ -105,7 +105,7 @@ Simd128_impl<std::is_arithmetic<T>::value, std::is_integral<T>::value, std::is_s
 
 #include "simd128_int16.inl"
 #include "simd128_int32.inl"
-#ifdef __x86_64__
+#if defined( __x86_64__) or defined(__aarch64__)
 #include "simd128_int64.inl"
 #endif
 #endif //#ifdef SIMD_INT
diff --git a/fflas-ffpack/fflas/fflas_transpose.h b/fflas-ffpack/fflas/fflas_transpose.h
index 3ff5b264a..5c4891655 100644
--- a/fflas-ffpack/fflas/fflas_transpose.h
+++ b/fflas-ffpack/fflas/fflas_transpose.h
@@ -291,7 +291,7 @@ namespace FFLAS {
         }
     }
 
-    /**************************************************************************/
+  /**************************************************************************/
     /*
      * Perfom transposition on the matrix A and store the result in matrix B.
      *      B[j,i] = A[i,j]     for 0 <= i < m, 0 <= j < n
diff --git a/fflas-ffpack/field/Makefile.am b/fflas-ffpack/field/Makefile.am
index a1ac51f67..798740a41 100644
--- a/fflas-ffpack/field/Makefile.am
+++ b/fflas-ffpack/field/Makefile.am
@@ -27,9 +27,12 @@ pkgincludesubdir=$(pkgincludedir)/field
 RNS=rns.h			        \
 	rns.inl				\
 	rns-double.h			\
+	rns-common.h\
+rns-double-extended.h			\
 	rns-double-elt.h		\
 	rns-double.inl			\
 	rns-double-recint.inl		\
+	rns-double-extended.inl			\
 	rns-integer.h			\
 	rns-integer-mod.h
 
diff --git a/fflas-ffpack/field/rns-double-extended.h b/fflas-ffpack/field/rns-double-extended.h
new file mode 100644
index 000000000..87e881dca
--- /dev/null
+++ b/fflas-ffpack/field/rns-double-extended.h
@@ -0,0 +1,431 @@
+/*
+ * Copyright (C) 2014 the FFLAS-FFPACK group
+ *
+ * Written by Pascal Giorgi <pascal.giorgi@lirmm.fr>
+ *
+ *
+ * ========LICENCE========
+ * This file is part of the library FFLAS-FFPACK.
+ *
+ * FFLAS-FFPACK is free software: you can redistribute it and/or modify
+ * it under the terms of the  GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * This library 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
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+ * ========LICENCE========
+ *.
+ */
+
+/*! @file field/rns-double-extended.h
+ * @ingroup field
+ * @brief  rns structure with double support (full precision)
+ */
+
+#ifndef __FFPACK_rns_double_extended_H
+#define __FFPACK_rns_double_extended_H
+
+// Bigger multiple of s lesser or equal than x, s must be a power of two
+#ifndef ROUND_DOWN
+#define ROUND_DOWN(x, s) ((x) & ~((s)-1))
+#endif
+
+#include <vector>
+#include <givaro/modular-floating.h>
+#include <givaro/givinteger.h>
+#include <givaro/givintprime.h>
+#include "givaro/modular-extended.h"
+#include <recint/ruint.h>
+#include "fflas-ffpack/config-blas.h"
+#include "fflas-ffpack/utils/fflas_memory.h"
+#include "fflas-ffpack/utils/align-allocator.h"
+#include "fflas-ffpack/field/rns-double-elt.h"
+
+namespace FFPACK {
+
+  // forward declaration of different variants
+  struct  rns_double_extended_V1 ;
+  struct  rns_double_extended_V2 ;
+
+  // set default value to one variant
+  using rns_double_extended = rns_double_extended_V2;
+
+
+#define RALIGN setw(40)<<right
+	/* Structure that handles rns representation given a bound and bitsize for prime moduli, allow large moduli
+	 * support sign representation (i.e. the bound must be twice larger then ||A||)
+	 */
+	struct rns_double_extended_V1 {
+		typedef Givaro::Integer integer;
+		typedef Givaro::ModularExtended<double> ModField;
+		
+		std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _basis; // the rns moduli (mi)
+		std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _basisMax; // (mi-1)
+		std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _negbasis; // (-mi)
+		std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _invbasis; // the inverse of rns moduli (1/mi)
+		std::vector<ModField> _field_rns; // the associated prime field for each mi
+		integer                  _M; // the product of the mi's
+		std::vector<integer>         _Mi; // _M/mi
+		std::vector<double>         _MMi; // (_Mi)^(-1) mod mi
+		std::vector<double>      _crt_in[6]; //  2^(16*j) mod mi 
+		std::vector<double>     _crt_out[6]; //  (_Mi._MMi) written in base 2^16
+		size_t                _size; // the size of the rns basis (number of mi's)
+		size_t               _pbits; // the size in bit of the mi's
+		size_t                 _ldm; // log[2^16](_M)
+
+		typedef double                        BasisElement;
+		typedef rns_double_elt                     Element;
+		typedef rns_double_elt_ptr             Element_ptr;
+		typedef rns_double_elt_cstptr     ConstElement_ptr;
+
+		rns_double_extended_V1(const integer& bound, size_t pbits, bool rnsmod=false, long seed=time(NULL))
+			:  _M(1), _size(0), _pbits(pbits)
+		{
+			integer::seeding(seed);
+			integer prime; Givaro::IntPrimeDom IPD;
+			integer sum=1;
+			while (_M < bound*sum) {
+				_basis.resize(_size+1);
+				do {
+					integer::random_exact_2exp(prime, _pbits-1);
+					IPD.nextprimein(prime);
+				} while (_M%prime == 0);
+				_basis[_size]=prime;
+				_size++;
+				_M*=prime;
+				if (rnsmod) sum+=prime;
+			}
+			precompute_cst();
+		}
+
+		rns_double_extended_V1(size_t pbits, size_t size, long seed=time(NULL))
+			:  _M(1), _size(size), _pbits(pbits)
+		{
+			integer::seeding(seed);
+			integer prime; Givaro::IntPrimeDom IPD;
+			integer sum=1;
+			_basis.resize(size);
+			_negbasis.resize(size);
+			_basisMax.resize(size);
+			for(size_t i = 0 ; i < _size ; ++i){
+				integer::random_exact_2exp(prime, _pbits-1);
+				IPD.nextprimein(prime);
+				_basis[i]=prime;
+				_basisMax[i] = prime-1;
+				_negbasis[i] = 0-prime;
+				_M*=prime;
+			}
+			precompute_cst();
+		}
+
+		template<typename Vect>
+		rns_double_extended_V1(const Vect& basis, bool rnsmod=false, long seed=time(NULL))
+			:  _basis(basis.begin(),basis.end()), _basisMax(basis.size()), _negbasis(basis.size()), _M(1), _size(basis.size()), _pbits(0)
+		{
+			for(size_t i=0;i<_size;i++){
+				_M*=_basis[i];
+				_pbits=std::max(_pbits, integer(_basis[i]).bitsize());
+			}
+			precompute_cst();
+		}
+
+
+		void precompute_cst(){
+
+      // Check that _pbits <= 48 
+      if ( _pbits > 48 ){
+        std::cout<<"FFLAS Error in rns_double_extended: primes bitsize "<<_pbits<<" is too large ... aborting\n";
+        std::terminate();
+      }
+
+			_ldm = (_M.bitsize()/48) + ((_M.bitsize()%48)?1:0) ;
+			_invbasis.resize(_size);
+			_basisMax.resize(_size);
+			_negbasis.resize(_size);
+			_field_rns.resize(_size);
+			_Mi.resize(_size);
+			_MMi.resize(_size);
+			for (size_t i=0;i<6;i++){
+				_crt_in[i] .resize(_size*_ldm);
+				_crt_out[i].resize(_size*_ldm);
+			}
+			const unsigned int MASK=0xFFFF;
+			for (size_t i=0;i<_size;i++){
+				_invbasis[i]  = 1./_basis[i];
+				_basisMax[i] = _basis[i]-1;
+				_negbasis[i] = 0-_basis[i];
+				_field_rns[i] = ModField(_basis[i]);
+				_Mi[i]        = _M/(uint64_t)_basis[i];
+				_field_rns[i].init(_MMi[i], _Mi[i] % (double)_basis[i]);
+				_field_rns[i].invin(_MMi[i]);
+				integer tmp= _Mi[i]*(uint64_t)_MMi[i];
+				double a0,a1,a2;
+				for(size_t j=0;j<_ldm;j++){
+					uint64_t limb= tmp[0];
+					a0= double(limb&MASK);
+					a1= double((limb>>16)&MASK);
+					a2= double((limb>>32)&MASK);					
+					_crt_out[0][j+i*_ldm]=a0;
+					_crt_out[1][j+i*_ldm]=a1;
+					_crt_out[2][j+i*_ldm]=a2;
+					_crt_out[3][j+i*_ldm]=a0+a1;
+					_crt_out[4][j+i*_ldm]=a1+a2;
+					_crt_out[5][j+i*_ldm]=a0+a1+a2;					
+					tmp>>=48;
+				}
+				double beta=double(1UL<<48);
+				double  acc=1;
+				for(size_t j=0;j<_ldm;j++){
+					uint64_t limb= acc;
+					a0= double(limb&MASK);
+					a1= double((limb>>16)&MASK);
+					a2= double((limb>>32)&MASK);					
+					_crt_in[0][j+i*_ldm]=a0;
+					_crt_in[1][j+i*_ldm]=a1;
+					_crt_in[2][j+i*_ldm]=a2;
+					_crt_in[3][j+i*_ldm]=a0+a1;
+					_crt_in[4][j+i*_ldm]=a1+a2;
+					_crt_in[5][j+i*_ldm]=a0+a1+a2;					
+					_field_rns[i].mulin(acc,beta);
+				}
+			}
+#ifdef RNS_DEBUG
+			std::cout<<"basis= [";
+			for(size_t i=0;i<_size;i++)
+				std::cout<<(int64_t)_basis[i]<<(i!=_size-1?",":"]\n");
+      Givaro::ModularExtended<double> ZZ(2UL<<48);;
+			for(size_t l=0;l<6;l++){
+        std::cout<<"crtin"<<l<<"=\n";
+        FFLAS::WriteMatrix(std::cout, ZZ, _size, _ldm, _crt_in[l].data(), _ldm);
+        std::cout<<std::endl;				
+			}
+			for(size_t l=0;l<6;l++){
+				std::cout<<"crtout"<<l<<"=\n";
+        FFLAS::WriteMatrix(std::cout, ZZ, _size, _ldm, _crt_out[l].data(), _ldm);
+        std::cout<<std::endl;				
+			}
+#endif
+		}
+
+		// Arns must be an array of m*n*_size
+		// abs(||A||) <= maxA
+		void init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, const integer& maxA, bool RNS_MAJOR=false) const
+		{
+			init(m,n,Arns,rda,A,lda, uint64_t(maxA.bitsize()/48 + (maxA.bitsize()%48?1:0)),RNS_MAJOR);
+		}
+	  
+    void init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false) const;
+		void convert(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false) const;
+
+		// reduce entries of Arns to be less than the rns basis elements
+		void reduce(size_t n, double* Arns, size_t rda, bool RNS_MAJOR=false) const;
+
+	}; // end of struct rns_double_extended
+
+
+
+
+
+	struct rns_double_extended_V2 {
+		typedef Givaro::Integer integer;
+		typedef Givaro::ModularExtended<double> ModField;
+		
+		std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _basis; // the rns moduli (mi)
+		std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _basisMax; // (mi-1)
+		std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _negbasis; // (-mi)
+		std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _invbasis; // the inverse of rns moduli (1/mi)
+		std::vector<ModField> _field_rns; // the associated prime field for each mi
+		integer                  _M; // the product of the mi's
+		std::vector<integer>         _Mi; // _M/mi
+		std::vector<double>         _MMi; // (_Mi)^(-1) mod mi
+		std::vector<double>      _crt_in; //  2^(16*j) mod mi -> each entry are splitted into two line (a0 + 2^27 a1)
+		std::vector<double>     _crt_out; //  (_Mi._MMi) written in base 2^16
+		size_t                _size; // the size of the rns basis (number of mi's)
+		size_t               _pbits; // the size in bit of the mi's
+		size_t                 _ldm; // log[2^16](_M)
+		uint64_t    _shift;
+
+		typedef double                        BasisElement;
+		typedef rns_double_elt                     Element;
+		typedef rns_double_elt_ptr             Element_ptr;
+		typedef rns_double_elt_cstptr     ConstElement_ptr;
+
+		rns_double_extended_V2 (const integer& bound, size_t pbits, bool rnsmod=false, long seed=time(NULL))
+			:  _M(1), _size(0), _pbits(pbits)
+		{
+			if (pbits > 52){
+				std::cerr<<"FFLAS- RNS EXTENDED ERROR: prime bitsize above 52 bits ... log(p)="<<pbits<<std::endl;
+				std::terminate();
+			}
+			_shift=pbits/2;
+
+			if( (bound.bitsize()/uint64_t(pbits)) * (integer(1)<<_shift) >= (integer(1)<<37) ){
+				std::cerr<<"FFLAS- RNS EXTENDED ERROR: prime bitsize ("<<pbits<<") too large to handle RNS basis with "<<bound.bitsize()<<" bits"<<std::endl;
+				std::terminate();				
+      }
+
+
+			integer::seeding(seed);
+			Givaro::IntPrimeDom IPD;
+			integer prime;
+			integer sum=1;
+			while (_M < bound*sum) {
+				_basis.resize(_size+1);
+				do {
+					integer::random_exact_2exp(prime, _pbits-1);
+					IPD.nextprimein(prime);
+				} while (_M%prime == 0);
+				_basis[_size]=prime;
+				_size++;
+				_M*=prime;
+				if (rnsmod) sum+=prime;
+			}
+			precompute_cst();
+		}
+		// can force to reduce integer entries larger than M
+		void precompute_cst(size_t K=0){
+
+      // Check that _pbits satisfies log(_M)/16* 2^(_pbits/2) * 2^16 < 2^53 (Required for correctness)
+      // => pbits <= 2*(41 - loglog(_M))
+      if ( _pbits > 2*(41 - log(double(_M.bitsize()))/log(2.))){
+        std::cout<<"FFLAS Error in rns_double_extended: primes bitsize "<<_pbits<<" is too large ... aborting\n";
+        std::terminate();
+      }
+
+			if (K!=0)
+				_ldm=K;
+			else
+				_ldm = (_M.bitsize()/16) + ((_M.bitsize()%16)?1:0) ;
+			_invbasis.resize(_size);
+			_field_rns.resize(_size);
+			_Mi.resize(_size);
+			_MMi.resize(_size);
+			_basisMax.resize(_size);
+			_negbasis.resize(_size);			
+			_crt_in.resize(2*_size*_ldm);
+			_crt_out.resize(_size*_ldm);
+			//const unsigned int MASK=0xFFFF;
+#ifdef BENCH_RNS_PRECOMP
+			Givaro::Timer chrono;
+			double t1=0.,t2=0.,t3=0.;
+#endif
+			for (size_t i=0;i<_size;i++){
+#ifdef BENCH_RNS_PRECOMP
+				chrono.start();
+#endif
+				_invbasis[i]  = 1./_basis[i];
+				_basisMax[i] = _basis[i]-1;
+				_negbasis[i] = 0-_basis[i];
+				_field_rns[i] = ModField(_basis[i]);
+				_Mi[i]        = _M/(uint64_t)_basis[i];
+				_field_rns[i].init(_MMi[i], _Mi[i] % (double)_basis[i]);				
+				_field_rns[i].invin(_MMi[i]);
+				integer tmp= _Mi[i]*(uint64_t)_MMi[i];
+				const mpz_t*    m0     = reinterpret_cast<const mpz_t*>(&tmp);
+				const uint16_t* m0_ptr = reinterpret_cast<const uint16_t*>(m0[0]->_mp_d);
+				size_t maxs=std::min(_ldm,(tmp.size())*sizeof(mp_limb_t)/2);// to ensure 32 bits portability
+#ifdef BENCH_RNS_PRECOMP
+				chrono.stop();
+				t1+=chrono.usertime();
+				chrono.start();
+#endif 
+				size_t l=0;
+#ifdef __FFLASFFPACK_HAVE_LITTLE_ENDIAN
+        for(;l<maxs;l++)
+          _crt_out[l+i*_ldm]=m0_ptr[l];
+#else
+        for(;l<maxs;l++)
+          _crt_out[l+i*_ldm]=m0_ptr[l ^ ((sizeof(mp_limb_t)/2U) - 1U)];
+#endif
+        for(;l<_ldm;l++)
+					_crt_out[l+i*_ldm]=0.;;
+#ifdef BENCH_RNS_PRECOMP
+				chrono.stop();
+        t2+=chrono.usertime();
+        chrono.start();
+#endif
+        double beta=double(1<<16);
+        double  acc=1;
+
+        for(size_t j=0;j<_ldm;j++){
+          uint64_t acci= (uint64_t)acc;
+          _crt_in[j+i*_ldm]=acci  & ((1<<_shift)-1);
+          _crt_in[j+(i+_size)*_ldm]=(acci >> _shift);
+          _field_rns[i].mulin(acc,beta);					
+
+        }
+#ifdef BENCH_RNS_PRECOMP
+				chrono.stop();
+				t3+=chrono.usertime();
+#endif
+			}
+#ifdef BENCH_RNS_PRECOMP
+			std::cout<<"RNS precomp t1="<<t1<<std::endl;
+			std::cout<<"RNS precomp t2="<<t2<<std::endl;
+			std::cout<<"RNS precomp t3="<<t3<<std::endl;
+#endif
+
+
+#ifdef RNS_DEBUG
+			std::cout<<"RNS double ext - basis:= [";
+			for(size_t i=0;i<_size;i++)
+				std::cout<<(int64_t)_basis[i]<<(i!=_size-1?",":"];\n");
+			Givaro::ModularExtended<double> ZZ(2UL<<48);;
+			std::cout<<"CRTmat:=";
+      //write_field(ZZ,std::cout,_crt_out.data(), _size, _ldm,_ldm,true);
+      FFLAS::WriteMatrix(std::cout, ZZ, _size, _ldm, _crt_out.data(), _ldm);
+#endif			
+		}
+
+    // Arns must be an array of m*n*_size
+    // abs(||A||) <= maxA    
+
+
+
+    void init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false, const FFLAS::FFLAS_TRANSPOSE ta=FFLAS::FflasNoTrans) const;
+    void convert(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false, const FFLAS::FFLAS_TRANSPOSE ta=FFLAS::FflasNoTrans) const;
+    void convert_bis(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false, const FFLAS::FFLAS_TRANSPOSE ta=FFLAS::FflasNoTrans) const;
+
+    // Arns must be an array of m*n*_size with  abs(||A||) <= maxA    
+    void init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, const integer& maxA, bool RNS_MAJOR=false) const{
+      init(m,n,Arns,rda,A,lda, maxA.bitsize()/16 + (maxA.bitsize()%16?1:0),RNS_MAJOR, FFLAS::FflasNoTrans);
+    }
+    void init_transpose(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, const integer& maxA, bool RNS_MAJOR=false) const{
+      init(m,n,Arns,rda,A,lda, maxA.bitsize()/16 + (maxA.bitsize()%16?1:0),RNS_MAJOR, FFLAS::FflasTrans);
+    }
+    void init_transpose(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false) const{
+      init(m,n,Arns,rda,A,lda,k,RNS_MAJOR, FFLAS::FflasTrans);
+    }
+    void convert_transpose(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false) const {
+      convert(m,n,gamma,A,lda,Arns,rda,RNS_MAJOR,FFLAS::FflasTrans);
+    }
+
+
+
+    // reduce entries of Arns to be less than the rns basis elements
+    void reduce(size_t n, double* Arns, size_t rda, bool RNS_MAJOR=false) const;
+
+    template<size_t K>
+    void init(size_t m, size_t n, double* Arns, size_t rda, const RecInt::ruint<K>* A, size_t lda, size_t k, bool RNS_MAJOR=false) const;
+    template<size_t K>
+    void convert(size_t m, size_t n, integer gamma, RecInt::ruint<K>* A, size_t lda, const double* Arns, size_t rda, integer p=0,bool RNS_MAJOR=false) const;
+    
+	};
+
+
+   
+	
+} // end of namespace FFPACK
+
+#include "rns-double-extended.inl"
+
+#endif //__FFPACK_rns_double_extended_H
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
diff --git a/fflas-ffpack/field/rns-double-extended.inl b/fflas-ffpack/field/rns-double-extended.inl
new file mode 100644
index 000000000..67392ffd1
--- /dev/null
+++ b/fflas-ffpack/field/rns-double-extended.inl
@@ -0,0 +1,1103 @@
+/*
+ * Copyright (C) 2016 the FFLAS-FFPACK group
+ *
+ * Written by Pascal Giorgi <pascal.giorgi@lirmm.fr>
+ *
+ *
+ * ========LICENCE========
+ * This file is part of the library FFLAS-FFPACK.
+ *
+ * FFLAS-FFPACK is free software: you can redistribute it and/or modify
+ * it under the terms of the  GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * This library 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
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+ * ========LICENCE========
+ *.
+ */
+
+
+#ifndef __FFLASFFPACK_field_rns_double_extended_INL
+#define __FFLASFFPACK_field_rns_double_extended_INL
+#include "fflas-ffpack/utils/Matio.h"
+#include "fflas-ffpack/fflas/fflas_freduce.h"
+
+namespace FFLAS {
+  template<>
+  void fscalin(const Givaro::ModularExtended<double> & F, const size_t n, double a,		
+               double* X, const size_t incX)
+  {
+    for(size_t i=0;i<n;i+=incX){
+      F.mulin(X[i],a);
+    }
+  }
+}
+
+namespace FFPACK {
+
+  /**************************************
+   * V1 Variant of RNS DOUBLE EXTENDED
+   **************************************/
+
+
+  // reduce entries of Arns to be less than the rns basis elements
+  inline void rns_double_extended_V1::reduce(size_t n, double* Arns, size_t rda, bool RNS_MAJOR) const{
+
+#ifdef __FFLASFFPACK_HAVE_SSE4_1_INSTRUCTIONS
+    using simd = Simd<double>;
+    using vect_t = typename simd::vect_t;
+
+    if(_size % simd::vect_size == 0){
+      //#pragma omp parallel for schedule(static, 256)			  
+      for(size_t i = 0 ; i < n ; i++){
+        vect_t tmp1, tmp2, v, min, max, basis, inv, neg;
+        min = simd::set1(0.);
+        for(size_t j = 0 ; j < _size ; j+=simd::vect_size){
+          basis = simd::load(_basis.data()+j);
+          inv   = simd::load(_invbasis.data()+j);
+          max   = simd::load(_basisMax.data()+j);
+          neg   = simd::load(_negbasis.data()+j);
+          v     = simd::load(Arns+i*_size+j);
+          simd::mod(v, basis, inv, neg, min, max, tmp1,tmp2);
+          simd::store(Arns+i*_size+j, v);
+        }
+      }
+    } else{
+      //#pragma omp parallel for schedule(static, 256)			  
+      for(size_t i = 0 ; i < n ; i++){
+        vect_t tmp1, tmp2, tmp3, v, min, max, basis, inv, neg;
+        size_t j = 0;
+        for( ; j < ROUND_DOWN(_size, simd::vect_size) ; j+=simd::vect_size){
+          basis = simd::load(_basis.data()+j);
+          inv   = simd::load(_invbasis.data()+j);
+          max   = simd::load(_basisMax.data()+j);
+          neg   = simd::load(_negbasis.data()+j);
+          v     = simd::load(Arns+i*_size+j);
+          simd::mod(v, basis, inv, neg, min, max, tmp1,tmp2);
+          simd::store(Arns+i*_size+j, v);
+        }
+        for( ; j < _size ; ++j){
+          _field_rns[j].reduce(Arns[i*_size+j]);
+        }
+      }
+    }
+#else
+    for(size_t i=0;i<_size;i++)
+      FFLAS::freduce (_field_rns[i],n,Arns+i*rda,1);
+
+#endif
+
+  }
+
+
+
+  void rns_double_extended_V1::init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR) const
+  {
+    if (_size>= (1<<16)){
+      std::cerr<<"RNS EXTENDED DOUBLE: init Error -> the nbr of moduli in RNS basis is > 2^16, not implemented. aborting\n";std::terminate();
+    }
+#ifdef BENCH_RNS
+    if (m!=1 && n!=1){
+      std::cerr<<RALIGN<<"RNS double ext (To) --> rns size ("<<_size<<") kronecker size ("<<k<<") data dim ("<<m*n<<")"<<std::endl;
+      std::cerr<<"RNS double ext -> Numbit(M)="<<_M.bitsize()<<std::endl;
+    }
+#endif
+    //init(m*n,Arns,A,lda);
+    if (k>_ldm){
+      FFPACK::failure()(__func__,__FILE__,__LINE__,"rns_struct: init (too large entry)");
+      std::cerr<<"k="<<k<<" _ldm="<<_ldm<<std::endl;
+    }
+    size_t mn=m*n;
+    size_t mnk=mn*k;
+    double *A_beta = FFLAS::fflas_new<double >(mnk*6);
+    double *A_beta0=A_beta;
+    double *A_beta1=A_beta+1*mnk;
+    double *A_beta2=A_beta+2*mnk;
+    double *A_beta3=A_beta+3*mnk;
+    double *A_beta4=A_beta+4*mnk;
+    double *A_beta5=A_beta+5*mnk;
+    size_t mnsize=mn*_size;
+    double *A_rns_tmp = FFLAS::fflas_new<double >(mnsize*6);
+    double *A_rns0=A_rns_tmp;
+    double *A_rns1=A_rns_tmp+1*mnsize;
+    double *A_rns2=A_rns_tmp+2*mnsize;
+    double *A_rns3=A_rns_tmp+3*mnsize;
+    double *A_rns4=A_rns_tmp+4*mnsize;
+    double *A_rns5=A_rns_tmp+5*mnsize;
+
+		
+    const integer* Aiter=A;
+    // split A into A_beta according to a Kronecker transform in base 2^48
+    Givaro::Timer tkr; tkr.start();
+    PARFOR1D(i,m,SPLITTER(NUM_THREADS),
+             for(size_t j=0;j<n;j++){
+               size_t idx=j+i*n;
+               const mpz_t*    m0     = reinterpret_cast<const mpz_t*>(Aiter+j+i*lda);
+               const uint16_t* m0_ptr = reinterpret_cast<const uint16_t*>(m0[0]->_mp_d);
+               size_t l=0,h=0;
+               // size in base 2^16
+               size_t k16=((Aiter[j+i*lda].size())*sizeof(mp_limb_t)/2);
+               size_t maxs=std::min(k,((Aiter[j+i*lda].size())*sizeof(mp_limb_t)/2)/3);// to ensure 32 bits portability
+               double a0,a1,a2;
+               if (m0[0]->_mp_size >= 0)
+                 for (;l<maxs;l++){
+                   a0= m0_ptr[h++];//std::cout<<"a0="<<(uint64_t)a0<<std::endl;
+                   a1= m0_ptr[h++];//std::cout<<"a1="<<(uint64_t)a1<<std::endl;
+                   a2= m0_ptr[h++];//std::cout<<"a2="<<(uint64_t)a2<<std::endl;
+                   A_beta0[l+idx*k]= a0;
+                   A_beta1[l+idx*k]= a1;
+                   A_beta2[l+idx*k]= a2;
+                   A_beta3[l+idx*k]= a0+a1;
+                   A_beta4[l+idx*k]= a1+a2;
+                   A_beta5[l+idx*k]= a0+a1+a2;							  							  
+                 }
+               else
+                 for (;l<maxs;l++){
+                   a0= -double(m0_ptr[h++]);
+                   a1= -double(m0_ptr[h++]);
+                   a2= -double(m0_ptr[h++]);
+                   A_beta0[l+idx*k]= a0;
+                   A_beta1[l+idx*k]= a1;
+                   A_beta2[l+idx*k]= a2;
+                   A_beta3[l+idx*k]= a0+a1;
+                   A_beta4[l+idx*k]= a1+a2;
+                   A_beta5[l+idx*k]= a0+a1+a2;							  							  
+                 }
+               for (;l<k;l++){
+                 a0= (h<k16)?m0_ptr[h++]:0.;
+                 a1= (h<k16)?m0_ptr[h++]:0.;
+                 a2= (h<k16)?m0_ptr[h++]:0.;
+                 A_beta0[l+idx*k]= a0;
+                 A_beta1[l+idx*k]= a1;
+                 A_beta2[l+idx*k]= a2;
+                 A_beta3[l+idx*k]= a0+a1;
+                 A_beta4[l+idx*k]= a1+a2;
+                 A_beta5[l+idx*k]= a0+a1+a2;							  							  							  
+               }					
+             }
+             );
+
+#ifdef CHECK_RNS
+    for (size_t i=0;i<m;i++)
+      for (size_t j=0;j<n;j++){
+        //std::cout<<"A="<<A[i*lda+j]<<std::endl;
+				
+				
+        integer tmp=0;
+        int idx=j+i*n;
+        for(int64_t l=(int64_t)k-1;l>=0;l--){
+          int64_t limb,c0,c1,c2,c3,c4;
+          c0= A_beta0[l+idx*k]; // 1
+          c1= A_beta3[l+idx*k] - A_beta1[l+idx*k] - A_beta0[l+idx*k] ; // X
+          c2= A_beta5[l+idx*k] - A_beta3[l+idx*k] - A_beta4[l+idx*k] + 2*A_beta1[l+idx*k]; // X^2
+          c3= A_beta4[l+idx*k] - A_beta1[l+idx*k] - A_beta2[l+idx*k]; // X^3
+          c4= A_beta2[l+idx*k]; // X^4
+          if (c1!=0. || c3!=0.) std::cout<<"RNS EXTENDED: error in splitting entries (linear form)\n";
+          limb= c0+(c2<<16)+(c4<<32);
+          tmp=(tmp<<48)+limb;
+        }
+        if (tmp!=A[i*lda+j]){
+          std::cout<<"RNS EXTENDED: error in splitting entries (16-adic) --> "<<tmp<<" != "<<A[i*lda+j]<<std::endl;;
+          std::cout<<"MAX="<<std::min(k,((Aiter[j+i*lda].size())*sizeof(mp_limb_t)/2)/3)<<std::endl;
+          std::cout<<"MAX1="<<k<<std::endl;
+          std::cout<<"MAX2="<<((Aiter[j+i*lda].size())*sizeof(mp_limb_t)/2)/3<<std::endl;
+          std::cout<<"bisize: "<<Aiter[j+i*lda].bitsize()<<std::endl;
+
+          tmp=A[i*lda+j];
+          for(size_t l=0;l<k;l++){
+            std::cout<<"l="<<l<<"  -->  "<<((tmp[0] & 0xFFFF))             << " != "<<int64_t(A_beta0[l+idx*k])<<std::endl;
+            std::cout<<"l="<<l<<"  -->  "<<((tmp[0] & 0xFFFFFFFF)>>16)     << " != "<<int64_t(A_beta1[l+idx*k])<<std::endl;
+            std::cout<<"l="<<l<<"  -->  "<<((tmp[0] & 0xFFFFFFFFFFFF)>>32) << " != "<<int64_t(A_beta2[l+idx*k])<<std::endl;
+            tmp>>=48;
+          }
+        }
+      }													
+#endif
+    tkr.stop();
+#ifdef BENCH_RNS
+    if (m!=1 && n!=1)
+      std::cerr<<RALIGN<<"RNS double ext (To) - Kronecker : "<<tkr.usertime()<<std::endl;
+#endif
+
+
+    Givaro::ZRing<double> ZD;
+
+    // Using Helper for potential parallelism -> need to be activated by hand
+    FFLAS::MMHelper<Givaro::ZRing<double>, FFLAS::MMHelperAlgo::Winograd>  MMH (ZD, -1, FFLAS::ParSeqHelper::Sequential());
+
+    if (RNS_MAJOR==false) {
+      // A_rns = _crt_in x A_beta^T
+      Givaro::Timer tfgemm; tfgemm.start();
+      FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in[0].data(),_ldm,A_beta0,k,0.,A_rns0, mn, MMH);
+      FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in[1].data(),_ldm,A_beta1,k,0.,A_rns1, mn, MMH);
+      FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in[2].data(),_ldm,A_beta2,k,0.,A_rns2, mn, MMH);
+      FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in[3].data(),_ldm,A_beta3,k,0.,A_rns3, mn, MMH);
+      FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in[4].data(),_ldm,A_beta4,k,0.,A_rns4, mn, MMH);
+      FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in[5].data(),_ldm,A_beta5,k,0.,A_rns5, mn, MMH);
+
+      tfgemm.stop();
+#ifdef BENCH_RNS
+      if (m!=1 && n!=1)
+        std::cerr<<RALIGN<<"RNS double ext (To)  - fgemm : "<<tfgemm.usertime()<<std::endl;
+#endif
+    }
+    else {
+      // Arns =  A_beta x _crt_in^T
+      std::cerr<<"NOT YET IMPLEMENTED .... aborting\n"; std::terminate();
+    }
+    Givaro::Timer tred; tred.start();
+
+#ifdef RNS_DEBUG			
+    std::cout<<"A:=Matrix("<<m<<","<<n<<",[";
+    for(size_t i=0;i<m;i++){
+      std::cout<<"[";
+      for(size_t j=0;j<n;j++)
+        std::cout<<(int64_t)A[j+i*lda]<<(j!=n-1?",":(i!=m-1?"],":"]"));
+    }
+    std::cout<<"]);\n";			
+
+    for(size_t l=0;l<6;l++){
+      std::cout<<"Abeta"<<l<<":=Matrix("<<_size<<","<<mn<<",[";
+      for(size_t i=0;i<k;i++){
+        std::cout<<"[";
+        for(size_t j=0;j<mn;j++)
+          std::cout<<(int64_t)A_beta[l*mnk+j+i*mn]<<(j!=mn-1?",":(i!=k-1?"],":"]"));
+      }
+      std::cout<<"]);\n";
+    }
+
+    for(size_t l=0;l<6;l++){
+      std::cout<<"Arns"<<l<<":=Matrix("<<_size<<","<<mn<<",[";
+      for(size_t i=0;i<_size;i++){
+        std::cout<<"[";
+        for(size_t j=0;j<mn;j++)
+          std::cout<<(int64_t)A_rns_tmp[l*mnsize+j+i*mn]<<(j!=mn-1?",":(i!=_size-1?"],":"]"));
+      }
+      std::cout<<"]);\n";
+    }
+#endif
+			
+			
+    double c0,c1,c2,c3,c4;
+    //double C,C1,C3;
+    //int64_t c0,c1,c2,c3,c4;
+
+    const double two16=1UL<<16,
+      two32=1UL<<32,
+      two48=1UL<<48,
+      two64=double(1UL<<32)*double(1UL<<32); 
+			
+    PARFOR1D(i,_size,SPLITTER(NUM_THREADS),
+             double two16_mod_mi,two48_mod_mi;
+             _field_rns[i].init(two16_mod_mi,(1<<16));
+             _field_rns[i].mul(two48_mod_mi,two16_mod_mi,two16_mod_mi);
+             _field_rns[i].mulin(two48_mod_mi,two16_mod_mi);
+					 
+             double two_16_invmi = two16 * _invbasis[i];
+             double two_32_invmi = two32 * _invbasis[i];
+             double two_48_invmi = two48 * _invbasis[i]; 
+             double two_64_invmi = two64 * _invbasis[i];
+             double q1,q2,q3,q4;
+
+             // std::cout<<"M1:="<<two_16_invmi<<";\n";
+             // std::cout<<"M2:="<<two_32_invmi<<";\n";
+             // std::cout<<"M3:="<<two_48_invmi<<";\n";
+             // std::cout<<"M4:="<<two_64_invmi<<";\n";
+					 
+					 
+             for(size_t j=0;j<mn;j++)
+               {
+                 size_t h=j+i*mn;
+                 c0= A_rns0[h]; // 1
+                 c1= A_rns3[h] - A_rns1[h] - A_rns0[h] ; // X
+                 c2= A_rns5[h] - A_rns3[h] - A_rns4[h] + 2*A_rns1[h]; // X^2
+                 c3= A_rns4[h] - A_rns1[h] - A_rns2[h]; // X^3
+                 c4= A_rns2[h]; // X^4
+#ifdef RNS_DEBUG
+                 std::cout<<(int64_t)c0
+                          <<"+2^16*"<<(int64_t)c1
+                          <<"+2^32*"<<(int64_t)c2
+                          <<"+2^48*"<<(int64_t)c3
+                          <<"+2^64*"<<(int64_t)c4<<";\n";
+#endif				 
+                 // compute c=c0+c1.2^16+c2.2^32+c3.2^48+c4.2^64 mod mi							 
+                 // _field_rns[i].axpy(c,c4,two16_mod_mi,c3);
+                 // _field_rns[i].axpy(c,c ,two16_mod_mi,c2);
+                 // _field_rns[i].axpy(c,c ,two16_mod_mi,c1);
+                 // _field_rns[i].axpy(c,c ,two16_mod_mi,c0);
+                 // Arns[j+i*rda]= c+(c>=0?0:_basis[i]);
+
+                 q1= std::floor(c1*two_16_invmi);//std::cout<<"Q1:="<<(int64_t)q1<<";\n";
+                 q2= std::floor(c2*two_32_invmi);//std::cout<<"Q2:="<<(int64_t)q2<<";\n";
+                 q3= std::floor(c3*two_48_invmi);//std::cout<<"Q3:="<<(int64_t)q3<<";\n";
+                 q4= std::floor(c4*two_64_invmi);//std::cout<<"Q4:="<<(int64_t)q4<<";\n";
+
+                 c1*=two16;
+                 c2*=two32;
+                 c3*=two48;
+                 c4*=two64;
+							 
+                 c1=fma(q1,_negbasis[i],c1); //std::cout<<"D1:="<<(int64_t)c1<<";\n";
+                 c2=fma(q2,_negbasis[i],c2); //std::cout<<"D2:="<<(int64_t)c2<<";\n";
+                 c3=fma(q3,_negbasis[i],c3); //std::cout<<"D3:="<<(int64_t)c3<<";\n";
+                 c4=fma(q4,_negbasis[i],c4); //std::cout<<"D4:="<<(int64_t)c4<<";\n";
+
+                 c0+=c1;
+                 c2+=c3;
+                 c0+=c4+c2;
+                 while(c0<0.)        c0+=_basis[i];
+                 while(c0>_basis[i]) c0-=_basis[i];
+                 Arns[j+i*rda]= c0;
+							 
+                 // c1+=c2<<16;
+                 // c3+=c4<<16;
+                 // _field_rns[i].reduce(C1,c1);
+                 // _field_rns[i].reduce(C3,c3);
+                 // _field_rns[i].axpy(C1,C1,two16_mod_mi,c0);
+                 // _field_rns[i].mul(C3,C3,two48_mod_mi);
+                 // _field_rns[i].add(C,C1,C3);
+                 // Arns[j+i*rda]= C+(C>=0?0:_basis[i]);
+               }
+             );
+    //reduce(mn,Arns,rda,RNS_MAJOR);
+			
+    tred.stop();
+#ifdef BENCH_RNS
+    if (m!=1 && n!=1)
+      std::cerr<<RALIGN<<"RNS double ext (To)  - Reduce : "<<tred.usertime()<<std::endl;
+#endif
+	
+    FFLAS::fflas_delete(A_beta);
+    FFLAS::fflas_delete(A_rns_tmp);
+
+#ifdef CHECK_RNS
+    bool ok=true;
+    for (size_t i=0;i<m;i++)
+      for(size_t j=0;j<n;j++)
+        for(size_t k=0;k<_size;k++){
+          ok&= (((A[i*lda+j] % (int64_t) _basis[k])+(A[i*lda+j]<0?(int64_t)_basis[k]:0)) == (int64_t) Arns[i*n+j+k*rda]);
+          if (((A[i*lda+j] % (int64_t) _basis[k])+(A[i*lda+j]<0?(int64_t)_basis[k]:0))
+              != (int64_t) Arns[i*n+j+k*rda])
+            {
+              std::cout<<((A[i*lda+j] % (int64_t) _basis[k])+(A[i*lda+j]<0?(int64_t)_basis[k]:0))
+                       <<" != "
+                       <<(int64_t) Arns[i*n+j+k*rda]
+                       <<" --> "<<A[i*lda+j]<<"  mod ("<<(uint64_t)_basis[k]<<")"<<std::endl;
+            }
+          else{
+            // std::cout<<"OK -> "<<((A[i*lda+j] % (int64_t) _basis[k])+(A[i*lda+j]<0?(int64_t)_basis[k]:0))
+            // 		 <<" == "
+            // 		 <<(int64_t) Arns[i*n+j+k*rda]
+            // 		 <<" --> "<<A[i*lda+j]<<"  mod ("<<(uint64_t)_basis[k]<<")"<<std::endl;
+						
+          }
+        }
+    std::cout<<"RNS EXTENDED Double (To)  ... "<<(ok?"OK":"ERROR")<<std::endl;
+    if (!ok) std::terminate();
+#endif
+
+
+
+			
+  }
+
+
+
+
+  struct uint48_t {
+    uint64_t data :48;
+    uint48_t () =default;
+    uint48_t (uint64_t x) :data(x){}
+  } __attribute__((packed));
+		
+
+		
+  void rns_double_extended_V1::convert(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR) const {
+    // assume the number of moduli is less than 2^16
+    if (_size>= (1<<16)){
+      std::cerr<<"RNS EXTENDED DOUBLE: convert Error -> the nbr of moduli in RNS basis is > 2^16, not implemented. aborting\n";std::terminate();
+    }
+#ifdef BENCH_RNS
+    if (m!=1 && n!=1)
+      std::cerr<<RALIGN<<"RNS double ext (From) --> rns size ("<<_size<<") kronecker size ("<<_ldm<<") data dim ("<<m*n<<")"<<std::endl;
+#endif
+			
+#ifdef CHECK_RNS
+    integer* Acopy=new integer[m*n];
+    for(size_t i=0;i<m;i++)
+      for(size_t j=0;j<n;j++)
+        Acopy[i*n+j]=A[i*lda+j];			
+#endif
+
+    integer hM= (_M-1)>>1;
+    size_t  mn= m*n;
+    size_t mnldm=mn*_ldm;
+    double *A_beta= FFLAS::fflas_new<double>(6*mnldm);			
+    double *A_beta0=A_beta;
+    double *A_beta1=A_beta+1*mnldm;
+    double *A_beta2=A_beta+2*mnldm;
+    double *A_beta3=A_beta+3*mnldm;
+    double *A_beta4=A_beta+4*mnldm;
+    double *A_beta5=A_beta+5*mnldm;
+    size_t mnsize=mn*_size;
+    double *A_rns_tmp = FFLAS::fflas_new<double >(mnsize*6);
+    double *A_rns0=A_rns_tmp;
+    double *A_rns1=A_rns_tmp+1*mnsize;
+    double *A_rns2=A_rns_tmp+2*mnsize;
+    double *A_rns3=A_rns_tmp+3*mnsize;
+    double *A_rns4=A_rns_tmp+4*mnsize;
+    double *A_rns5=A_rns_tmp+5*mnsize;
+			
+    Givaro::Timer tsplit;
+    tsplit.start();			
+    uint64_t tmp,idx;
+    double aa0,aa1,aa2;
+    for (size_t i=0;i<_size;i++)
+      for(size_t j=0;j<mn;j++){
+        idx=i*mn+j;
+        tmp=Arns[idx];
+        aa0=tmp&0xFFFF;
+        aa1=(tmp>>16)&0xFFFF;
+        aa2=tmp>>32;
+        A_rns0[idx]=aa0;
+        A_rns1[idx]=aa1;
+        A_rns2[idx]=aa2;
+        A_rns3[idx]=aa0+aa1;
+        A_rns4[idx]=aa1+aa2;
+        A_rns5[idx]=aa0+aa1+aa2;
+      }
+    tsplit.stop();
+#ifdef BENCH_RNS
+    if (m!=1 && n!=1)				
+      std::cerr<<RALIGN<<"RNS double ext (From) -  split : "<<tsplit.usertime()<<std::endl;
+#endif
+    Givaro::Timer tfgemmc;tfgemmc.start();
+    if (RNS_MAJOR==false){
+      // compute A_beta = Ap^T x M_beta
+      FFLAS::fgemm(Givaro::ZRing<double>(),FFLAS::FflasTrans, FFLAS::FflasNoTrans,(int) mn,(int) _ldm,(int) _size, 1.0 , A_rns0,(int) rda, _crt_out[0].data(),(int) _ldm, 0., A_beta0,(int)_ldm, FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive >());
+      FFLAS::fgemm(Givaro::ZRing<double>(),FFLAS::FflasTrans, FFLAS::FflasNoTrans,(int) mn,(int) _ldm,(int) _size, 1.0 , A_rns1,(int) rda, _crt_out[1].data(),(int) _ldm, 0., A_beta1,(int)_ldm, FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive >());
+      FFLAS::fgemm(Givaro::ZRing<double>(),FFLAS::FflasTrans, FFLAS::FflasNoTrans,(int) mn,(int) _ldm,(int) _size, 1.0 , A_rns2,(int) rda, _crt_out[2].data(),(int) _ldm, 0., A_beta2,(int)_ldm, FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive >());
+      FFLAS::fgemm(Givaro::ZRing<double>(),FFLAS::FflasTrans, FFLAS::FflasNoTrans,(int) mn,(int) _ldm,(int) _size, 1.0 , A_rns3,(int) rda, _crt_out[3].data(),(int) _ldm, 0., A_beta3,(int)_ldm, FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive >());
+      FFLAS::fgemm(Givaro::ZRing<double>(),FFLAS::FflasTrans, FFLAS::FflasNoTrans,(int) mn,(int) _ldm,(int) _size, 1.0 , A_rns4,(int) rda, _crt_out[4].data(),(int) _ldm, 0., A_beta4,(int)_ldm, FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive >());
+      FFLAS::fgemm(Givaro::ZRing<double>(),FFLAS::FflasTrans, FFLAS::FflasNoTrans,(int) mn,(int) _ldm,(int) _size, 1.0 , A_rns5,(int) rda, _crt_out[5].data(),(int) _ldm, 0., A_beta5,(int)_ldm, FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive >());
+    }
+    else {
+      // compute A_beta = Ap x M_Beta 
+      std::cerr<<"NOT YET IMPLEMENTED .... aborting\n"; std::terminate();
+    }
+    tfgemmc.stop();
+#ifdef BENCH_RNS
+    if (m!=1 && n!=1)
+      std::cerr<<RALIGN<<"RNS double ext (From) -  fgemm : "<<tfgemmc.usertime()<<std::endl;
+#endif
+    // compute A using inverse Kronecker transform of A_beta expressed in base 2^48
+
+    FFLAS::fflas_delete( A_rns_tmp);
+#ifdef RNS_DEBUG			
+    std::cout<<"basis:= [";
+    for(size_t i=0;i<_size;i++)
+      std::cout<<(int64_t)_basis[i]<<(i!=_size-1?",":"];\n");
+
+    std::cout<<"Arns:=Matrix("<<_size<<","<<mn<<",[";
+    for(size_t i=0;i<_size;i++){
+      std::cout<<"[";
+      for(size_t j=0;j<mn;j++)
+        std::cout<<(int64_t)Arns[j+i*mn]<<(j!=mn-1?",":(i!=_size-1?"],":"]"));
+    }
+    std::cout<<"]);\n";
+			
+    for(size_t l=0;l<6;l++){
+      std::cout<<"Arns"<<l<<":=Matrix("<<_size<<","<<mn<<",[";
+      for(size_t i=0;i<_size;i++){
+        std::cout<<"[";
+        for(size_t j=0;j<mn;j++)
+          std::cout<<(int64_t)A_rns_tmp[l*mnsize+j+i*mn]<<(j!=mn-1?",":(i!=_size-1?"],":"]"));
+      }
+      std::cout<<"]);\n";
+    }
+			
+    for(size_t l=0;l<6;l++){
+      std::cout<<"Abeta"<<l<<":=Matrix("<<mn<<","<<_size<<",[";
+      for(size_t i=0;i<mn;i++){
+        std::cout<<"[";
+        for(size_t j=0;j<_ldm;j++)
+          std::cout<<(int64_t)A_beta[l*mnldm+j+i*_ldm]<<(j!=_ldm-1?",":(i!=mn-1?"],":"]"));
+      }
+      std::cout<<"]);\n";
+    }
+
+#endif
+
+    Givaro::Timer tkroc;
+    tkroc.start();
+			
+    integer* Aiter= A;
+    size_t k=_ldm;
+    size_t k64= (k*48+64)/64+1;
+#ifdef RNS_DEBUG
+    std::cout<<"kron base 48: "<<k<<std::endl;
+    std::cout<<"kron base 64: "<<k64<<std::endl;
+#endif
+    std::vector<uint16_t> A0_tmp(k64<<2,0),A1_tmp(k64<<2,0),A2_tmp(k64<<2,0),A3_tmp(k64<<2,0),A4_tmp(k64<<2,0);
+    uint48_t *A0,*A1,*A2,*A3,*A4;
+    A0= reinterpret_cast<uint48_t*>(A0_tmp.data());
+    A1= reinterpret_cast<uint48_t*>(A1_tmp.data()+1);
+    A2= reinterpret_cast<uint48_t*>(A2_tmp.data()+2);
+    A3= reinterpret_cast<uint48_t*>(A3_tmp.data()+3);
+    A4= reinterpret_cast<uint48_t*>(A4_tmp.data()+4);
+    integer a0,a1,a2,a3,a4,res;
+    mpz_t *m0,*m1,*m2,*m3,*m4;
+    m0= reinterpret_cast<mpz_t*>(&a0);
+    m1= reinterpret_cast<mpz_t*>(&a1);
+    m2= reinterpret_cast<mpz_t*>(&a2);
+    m3= reinterpret_cast<mpz_t*>(&a3);
+    m4= reinterpret_cast<mpz_t*>(&a4);
+    mp_limb_t *m0_d,*m1_d,*m2_d,*m3_d,*m4_d;
+    m0_d = m0[0]->_mp_d;
+    m1_d = m1[0]->_mp_d;
+    m2_d = m2[0]->_mp_d;
+    m3_d = m3[0]->_mp_d;
+    m4_d = m4[0]->_mp_d;
+    m0[0]->_mp_alloc = m1[0]->_mp_alloc = m2[0]->_mp_alloc = m3[0]->_mp_alloc = m4[0]->_mp_alloc = (int) (k64*8/sizeof(mp_limb_t)); // to ensure 32 bits portability
+    m0[0]->_mp_size  = m1[0]->_mp_size  = m2[0]->_mp_size  = m3[0]->_mp_size  = m4[0]->_mp_size  = (int) (k64*8/sizeof(mp_limb_t)); // to ensure 32 bits portability
+    //		auto sp=SPLITTER();
+    //		PARFOR1D(i,m,sp,
+#ifdef RNS_DEBUG
+    std::cout<<"M:="<<_M<<std::endl;
+#endif
+    uint64_t c0,c1,c2,c3,c4;
+    for(size_t i=0;i<m;i++)
+      for (size_t j=0;j<n;j++){
+        size_t idx=i*n+j;
+#ifdef RNS_DEBUG
+        std::cout<<"c0:=0;";
+        std::cout<<"c1:=0;";
+        std::cout<<"c2:=0;";
+        std::cout<<"c3:=0;";
+        std::cout<<"c4:=0;";
+#endif	
+        for (size_t l=0;l<k;l++){
+          size_t idxl=l+idx*k;
+          c0= A_beta0[idxl]; // 1
+          c1= A_beta3[idxl] - A_beta1[idxl] - A_beta0[idxl] ; // X
+          c2= A_beta5[idxl] - A_beta3[idxl] - A_beta4[idxl] + 2*A_beta1[idxl]; // X^2
+          c3= A_beta4[idxl] - A_beta1[idxl] - A_beta2[idxl]; // X^3
+          c4= A_beta2[idxl]; // X^4
+#ifdef RNS_DEBUG
+          std::cout<<"CCC2"<<l<<":="<<c2<<"\n;";
+						
+          std::cout<<"c0:=c0+2^"<<l*48<<"*"<<c0<<";";
+          std::cout<<"c1:=c1+2^"<<l*48+16<<"*"<<c1<<";";
+          std::cout<<"c2:=c2+2^"<<l*48+32<<"*"<<c2<<";";
+          std::cout<<"c3:=c3+2^"<<l*48+48<<"*"<<c3<<";";
+          std::cout<<"c4:=c4+2^"<<l*48+64<<"*"<<c4<<";";
+#endif
+          A0[l]= c0;
+          A1[l]= c1;
+          A2[l]= c2;
+          A3[l]= c3;
+          A4[l]= c4;
+        }
+#ifdef RNS_DEBUG
+        std::cout<<"0 mod "<<_M<<";\n";
+#endif
+        // see A0,A1,A2,A3 as a the gmp integers a0,a1,a2,a3
+        m0[0]->_mp_d= reinterpret_cast<mp_limb_t*>(A0_tmp.data());
+        m1[0]->_mp_d= reinterpret_cast<mp_limb_t*>(A1_tmp.data());
+        m2[0]->_mp_d= reinterpret_cast<mp_limb_t*>(A2_tmp.data());
+        m3[0]->_mp_d= reinterpret_cast<mp_limb_t*>(A3_tmp.data());
+        m4[0]->_mp_d= reinterpret_cast<mp_limb_t*>(A4_tmp.data());
+        res = a0;res+= a1;res+= a2;res+= a3;res+=a4;
+        res%=_M;
+
+#ifdef RNS_DEBUG
+        std::cout<<"a0:="<<a0<<";\n";
+        std::cout<<"a1:="<<a1<<";\n";
+        std::cout<<"a2:="<<a2<<";\n";
+        std::cout<<"a3:="<<a3<<";\n";
+        std::cout<<"a4:="<<a4<<";\n";
+        std::cout<<"res:="<<res<<";\n";
+#endif	
+        // get the correct result according to the expected sign of A
+        if (res>hM)
+          res-=_M;
+        if (gamma==0)
+          Aiter[j+i*lda]=res;
+        else
+          if (gamma==integer(1))
+            Aiter[j+i*lda]+=res;
+          else
+            if (gamma==integer(-1))
+              Aiter[j+i*lda]=res-Aiter[j+i*lda];
+            else{
+              Aiter[j+i*lda]*=gamma;
+              Aiter[j+i*lda]+=res;
+            }
+					
+      }
+    tkroc.stop();
+#ifdef BENCH_RNS
+    if (m!=1 && n!=1)
+      std::cerr<<RALIGN<<"RNS double ext (From) - Kronecker : "<<tkroc.usertime()<<std::endl;
+#endif
+			
+    m0[0]->_mp_d = m0_d;
+    m1[0]->_mp_d = m1_d;
+    m2[0]->_mp_d = m2_d;
+    m3[0]->_mp_d = m3_d;
+    m4[0]->_mp_d = m4_d;
+    m0[0]->_mp_alloc = m1[0]->_mp_alloc = m2[0]->_mp_alloc= m3[0]->_mp_alloc = m4[0]->_mp_alloc = 1;
+    m0[0]->_mp_size  = m1[0]->_mp_size  = m2[0]->_mp_size = m3[0]->_mp_size  = m4[0]->_mp_size  = 0;
+    FFLAS::fflas_delete( A_beta);
+				 
+#ifdef CHECK_RNS
+    bool ok=true;
+    for (size_t i=0;i<m;i++)
+      for(size_t j=0;j<n;j++)
+        for(size_t k=0;k<_size;k++){
+          int64_t _p =(int64_t) _basis[k];
+          integer curr=A[i*lda+j] - gamma*Acopy[i*n+j];
+          ok&= ( curr% _p +(curr%_p<0?_p:0) == (int64_t) Arns[i*n+j+k*rda]);
+          if (( curr% _p +(curr%_p<0?_p:0)) != (int64_t) Arns[i*n+j+k*rda])
+            std::cout<<"A("<<i<<","<<j<<") -> "<<A[i*lda+j]<<" mod "<<(int64_t) _basis[k]<<"!="<<(int64_t) Arns[i*n+j+k*rda]<<";"<<std::endl;
+        }
+    std::cout<<"RNS EXTENDED Double (From)  ... "<<(ok?"OK":"ERROR")<<std::endl;
+    if (!ok) std::terminate();
+#endif
+  }
+
+
+
+  /**************************************
+   * V2 Variant of RNS DOUBLE EXTENDED
+   **************************************/
+
+  // reduce entries of Arns to be less than the rns basis elements
+  inline void rns_double_extended_V2::reduce(size_t n, double* Arns, size_t rda, bool RNS_MAJOR) const{
+
+#ifdef __FFLASFFPACK_HAVE_SSE4_1_INSTRUCTIONS
+    using simd = Simd<double>;
+    using vect_t = typename simd::vect_t;
+
+    if(_size % simd::vect_size == 0){
+      //#pragma omp parallel for schedule(static, 256)			  
+      for(size_t i = 0 ; i < n ; i++){
+        vect_t tmp1, tmp2, v, min, max, basis, inv, neg;
+        min = simd::set1(0.);
+        for(size_t j = 0 ; j < _size ; j+=simd::vect_size){
+          basis = simd::load(_basis.data()+j);
+          inv   = simd::load(_invbasis.data()+j);
+          max   = simd::load(_basisMax.data()+j);
+          neg   = simd::load(_negbasis.data()+j);
+          v     = simd::load(Arns+i*_size+j);
+          simd::mod(v, basis, inv, neg, min, max, tmp1,tmp2);
+          simd::store(Arns+i*_size+j, v);
+        }
+      }
+    } else{
+      //#pragma omp parallel for schedule(static, 256)			  
+      for(size_t i = 0 ; i < n ; i++){
+        vect_t tmp1, tmp2, tmp3, v, min, max, basis, inv, neg;
+        size_t j = 0;
+        for( ; j < ROUND_DOWN(_size, simd::vect_size) ; j+=simd::vect_size){
+          basis = simd::load(_basis.data()+j);
+          inv   = simd::load(_invbasis.data()+j);
+          max   = simd::load(_basisMax.data()+j);
+          neg   = simd::load(_negbasis.data()+j);
+          v     = simd::load(Arns+i*_size+j);
+          simd::mod(v, basis, inv, neg, min, max, tmp1,tmp2);
+          simd::store(Arns+i*_size+j, v);
+        }
+        for( ; j < _size ; ++j){
+          _field_rns[j].reduce(Arns[i*_size+j]);
+        }
+      }
+    }
+#else
+    for(size_t i=0;i<_size;i++)
+      FFLAS::freduce (_field_rns[i],n,Arns+i*rda,1);
+
+#endif
+
+  }
+
+
+  inline void rns_double_extended_V2::init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR, const FFLAS::FFLAS_TRANSPOSE trans) const
+  {
+    if (k>_ldm){
+      FFPACK::failure()(__func__,__FILE__,__LINE__,"rns_double_extended [init] -> rns basis is too small to handle integers with 2^(16*k) values ");
+      std::cerr<<"with k="<<k<<" _ldm="<<_ldm<<std::endl;
+    }
+    const size_t mn=m*n;
+    Givaro::ZRing<double> ZD;
+
+    if (mn) {
+      // split A into A_beta according to a Kronecker transform in base 2^16
+      double *A_beta = FFLAS::fflas_new<double >(mn*k);
+      RNS_COMMON::Kronecker_base16(trans,m,n,A_beta,A,lda,k);
+
+      // Using Helper for potential parallelism -> need to be activated by hand
+      //using ParallelStrategy= FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive>;
+      using ParallelStrategy= FFLAS::ParSeqHelper::Sequential;
+      FFLAS::MMHelper<Givaro::ZRing<double>, FFLAS::MMHelperAlgo::Winograd>  MatMulHelper (ZD, -1, ParallelStrategy());
+
+      // Pseudo reduction of the matrix entries
+      double *Arns_tmp = FFLAS::fflas_new<double >(mn*_size*2); 
+      if (RNS_MAJOR==false) // Arns = _crt_in x A_beta^T        
+        FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,2*_size,mn,k,1.0,_crt_in.data(),_ldm,A_beta,k,0.,Arns_tmp,mn,MatMulHelper);
+      else // Arns =  A_beta x _crt_in^T
+        FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,mn,2*_size,k,1.0,A_beta, k, _crt_in.data(),_ldm,0.,Arns_tmp,2*_size,MatMulHelper);
+
+
+      // Reduce the pseudo reduce matrix modulo RNS basis
+      // reduce(mn,Arns,rda,RNS_MAJOR);
+      /*
+        NEED TO MODIFY THE CODE FOR TAKING INTO ACCOUNT RNS MAJOR MODE       
+       */
+      size_t inc            = (RNS_MAJOR) ? 2*_size : 1;
+      size_t inc_input      = (RNS_MAJOR) ? _size : 1;
+      size_t real_rda       = (RNS_MAJOR) ? 1 : mn;
+      double *Arns_tmp_high = (RNS_MAJOR) ? Arns_tmp+_size: Arns_tmp+ mn*_size;
+
+
+      // Givaro::ZRing<int64_t> ZZ;
+      // int64_t * ArnsConverted = new int64_t[2*mn*_size];
+
+     
+
+      // std::cout<<"basis=[";
+      // for(size_t i=0;i<_size-1;i++)
+      //   std::cout<<(int64_t)_basis[i]<<",";
+      // std::cout<<(int64_t)_basis[_size-1]<<"]\n";
+
+      // std::cout<<"shift="<<uint64_t(1<<_shift)<<std::endl;
+
+      // FFLAS::fconvert(Givaro::ZRing<double>(),2*_size*m*n,ArnsConverted,1,Arns_tmp,1);
+      // FFLAS::WriteMatrix(std::cout<<"AA=",ZZ,m*n,2*_size, ArnsConverted, 2*_size,FFLAS::FflasSageMath)<<std::endl;
+
+
+
+      //FFLAS::fadd(_field_rns[i],mn,Arns_tmp+i*real_rda,inc,Arns_tmp+(_size+i)*real_rda,inc, Arns+i*real_rda,inc); //  -> scaling by 2^shift allow a proper modular reduction
+
+
+
+      for(size_t i=0;i<_size;i++){
+        // PG: it seems that the fscalin is buggy when called below -> use fscal + faddin instead of fscalin + fadd
+        //FFLAS::fscalin (_field_rns[i],mn, double(1<<_shift), Arns_tmp_high+i*real_rda, inc);
+        //FFLAS::fadd(_field_rns[i],mn,Arns_tmp+i*real_rda,inc, Arns_tmp_high+i*real_rda,inc, Arns+i*real_rda,inc_input);
+
+        FFLAS::fscal (_field_rns[i],mn, double(1<<_shift), Arns_tmp_high+i*real_rda, inc, Arns+i*real_rda, inc_input); // This two lines cannot be replaced with an axpy
+        FFLAS::faddin(_field_rns[i],mn, Arns_tmp+i*real_rda,inc, Arns+i*real_rda,inc_input);
+
+        // FFLAS::fconvert(Givaro::ZRing<double>(),_size*m*n,ArnsConverted,1,Arns,1);
+        // FFLAS::WriteMatrix(std::cout<<"AA"<<i<<"=",ZZ,m*n,_size, ArnsConverted, _size,FFLAS::FflasSageMath)<<std::endl;
+
+
+
+
+        FFLAS::finit (_field_rns[i],mn, Arns+i*real_rda,inc_input, Arns+i*real_rda,inc_input); //  needed to remove negative entry (i.e. kronecker keeps negative value)
+
+        //FFLAS::freduce (_field_rns[i],mn, Arns+i*real_rda,inc); // not needed this is already handled within the fscalin
+      }
+
+
+      // FFLAS::fconvert(Givaro::ZRing<double>(),_size*m*n,ArnsConverted,1,Arns,1);
+      // FFLAS::WriteMatrix(std::cout<<"A=",ZZ,m*n,_size, ArnsConverted, _size,FFLAS::FflasSageMath)<<std::endl;
+      // delete[] ArnsConverted;
+
+
+
+#ifdef CHECK_RNS // check the result
+      RNS_COMMON::check_integers_to_rns(trans, m, n, Arns, rda, A, lda, _basis, RNS_MAJOR);
+#endif
+
+      FFLAS::fflas_delete(A_beta);
+      FFLAS::fflas_delete(Arns_tmp);
+    }
+  }
+
+  inline void rns_double_extended_V2::convert(size_t m, size_t n, integer gamma, integer* A, size_t lda,
+                                              const double* Arns, size_t rda, bool RNS_MAJOR, const FFLAS::FFLAS_TRANSPOSE trans) const
+  {
+    const size_t  mn= m*n, mn2=2*mn;
+    uint64_t _mask= ((1<<_shift)-1);
+    Givaro::ZRing<double> ZD;
+    Givaro::ZRing<Givaro::Integer> ZZ;
+    if (mn) {
+#ifdef CHECK_RNS
+      integer* Acopy=new integer[m*n];
+      FFLAS::fassign(ZZ,m,n,Acopy,n,A,lda);
+#endif
+      double *A_beta   = FFLAS::fflas_new<double>(mn2*_ldm);
+
+      // Split Arns into two matrix (lower half bit and upper half bit)
+      double *Arns_low = FFLAS::fflas_new<double >(mn2*_size);
+   
+      // Using Helper for potential parallelism -> need to be activated by hand
+      //using ParallelStrategy= FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive>;
+      using ParallelStrategy= FFLAS::ParSeqHelper::Sequential;
+      FFLAS::MMHelper<Givaro::ZRing<double>, FFLAS::MMHelperAlgo::Winograd>  MatMulHelper (ZD, -1, ParallelStrategy());
+
+      if (RNS_MAJOR==false) {// compute A_beta = Arns^T x M_beta
+        double* Arns_high = Arns_low+mn;
+        for (size_t i=0;i<_size;i++)
+          for (size_t j=0;j<mn;j++){
+            int64_t acci= (int64_t)Arns[i*mn+j];
+            Arns_low [i*mn2+j] = acci  &  _mask;
+            Arns_high[i*mn2+j] = (acci >> _shift);
+          }
+        FFLAS::fgemm(Givaro::ZRing<double>(),FFLAS::FflasTrans, FFLAS::FflasNoTrans, mn2, _ldm, _size, 1.0 , Arns_low, mn2, _crt_out.data(), _ldm, 0., A_beta,_ldm, MatMulHelper);
+      }
+      else { // compute A_beta = Arns x M_Beta
+        double* Arns_high = Arns_low+mn*_size;
+        for (size_t i=0;i<mn;i++)
+          for (size_t j=0;j<_size;j++){
+            uint64_t acci= (uint64_t)Arns[i*_size+j];;
+            Arns_low [i*_size+j] = acci  &  _mask;
+            Arns_high[i*_size+j] = (acci >> _shift);
+          }
+        FFLAS::fgemm(Givaro::ZRing<double>(),FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, mn2, _ldm, _size, 1.0 , Arns_low, _size, _crt_out.data(), _ldm, 0., A_beta, _ldm, MatMulHelper);
+      }
+
+      // std::cout<<"****\n RNS partial reconstruction:\n";
+      // std::cout<<"shift="<<_shift<<"\n";
+      // Givaro::ZRing<int64_t> ZZ;
+      // int64_t * ArnsInt = new int64_t[mn2*_size];
+      // FFLAS::fconvert(Givaro::ZRing<double>(),_size,mn2,ArnsInt,mn2,Arns_low,mn2 );
+      // std::cout<<"Arns=";
+      // FFLAS::WriteMatrix(std::cout,ZZ,_size, mn2, ArnsInt, mn2,FFLAS::FflasSageMath)<<std::endl;
+      // std::cout<<"M=";
+      // FFLAS::WriteMatrix(std::cout,Givaro::ZRing<double>(), _size,_ldm,_crt_out.data(), _ldm,FFLAS::FflasSageMath)<<std::endl;
+
+
+      // int64_t * AbetaInt = new int64_t[mn2*_ldm];
+      // FFLAS::fconvert(Givaro::ZRing<double>(),mn2,_ldm,AbetaInt,_ldm,A_beta,_ldm );
+
+      // FFLAS::WriteMatrix(std::cout,Givaro::ZRing<int64_t>(), mn2,_ldm, AbetaInt, _ldm)<<std::endl;//FFLAS::FflasSageMath)<<std::endl;
+
+      // std::cout<<"******\n";
+
+
+      // inverse kronecker of the higher bit part
+      RNS_COMMON::inverse_Kronecker_base16(trans,m,n,A_beta+mn*_ldm ,A,lda,_ldm, _M, gamma,  (1LL<<_shift));
+      //FFLAS::WriteMatrix(std::cout,Givaro::ZRing<Givaro::Integer>(),m, n, A, n)<<std::endl;//FFLAS::FflasSageMath)<<std::endl;
+
+      // inverse kronecker of the lower bit part
+      RNS_COMMON::inverse_Kronecker_base16(trans,m,n,A_beta         ,A,lda,_ldm, _M, 1);
+
+      //FFLAS::WriteMatrix(std::cout,Givaro::ZRing<Givaro::Integer>(),m, n, A, n)<<std::endl;//FFLAS::FflasSageMath)<<std::endl;
+
+
+      FFLAS::fflas_delete(A_beta);
+      FFLAS::fflas_delete(Arns_low);
+#ifdef CHECK_RNS
+      RNS_COMMON::check_rns_to_integers(trans, m, n, Arns, rda, A, lda, gamma, Acopy, n, _basis, RNS_MAJOR);
+      FFLAS::fflas_delete(Acopy);
+#endif
+
+
+
+// #ifdef CHECK_RNS
+//       bool ok=true;
+//       for (size_t i=0;i<m;i++)
+//         for(size_t j=0;j<n;j++)
+//           for(size_t k=0;k<_size;k++){
+//             int64_t _p =(int64_t) _basis[k];
+//             integer curr=A[i*lda+j] - gamma*Acopy[i*n+j];
+//             ok&= ( curr% _p +(curr%_p<0?_p:0) == (int64_t) Arns[i*n+j+k*rda]);
+//           }
+//       std::cout<<"RNS convert ... "<<(ok?"OK":"ERROR")<<std::endl;
+// #endif
+    }
+  
+  }
+
+
+  inline void rns_double_extended_V2::convert_bis(size_t m, size_t n, integer gamma, integer* A, size_t lda,
+                                              const double* Arns, size_t rda, bool RNS_MAJOR, const FFLAS::FFLAS_TRANSPOSE trans) const
+	{
+#ifdef BENCH_RNS
+	  if (m!=1 && n!=1)
+	    std::cerr<<"RNS double ext (From) --> rns size ("<<_size<<") kronecker size ("<<_ldm<<") data dim ("<<m*n<<")"<<std::endl;
+#endif
+
+#ifdef CHECK_RNS
+		integer* Acopy=new integer[m*n];
+		for(size_t i=0;i<m;i++)
+			for(size_t j=0;j<n;j++)
+				Acopy[i*n+j]=A[i*lda+j];
+
+#endif
+		Givaro::Timer tsplit;tsplit.start();
+		size_t mn=m*n;
+		double *Arns_tmp = FFLAS::fflas_new<double >(mn*_size*2);
+		double* Arns_tmp2= Arns_tmp+mn;
+		for (size_t i=0;i<_size;i++)
+		  for (size_t j=0;j<mn;j++){
+		    uint64_t acci= (uint64_t)Arns[i*rda+j];;
+        Arns_tmp [i*2*mn+j] = acci  & ((1<<_shift)-1);
+        Arns_tmp2[i*2*mn+j] = (acci >> _shift);
+		  }
+		
+    
+#ifdef RNS_DEBUG
+		Givaro::ModularExtended<double> ZZ(2UL<<48);;
+		std::cout<<"Arns:=";
+    FFLAS::WriteMatrix(std::cout, ZZ, _size, mn, Arns, mn);
+		std::cout<<"Arns1:=";
+    FFLAS::WriteMatrix(std::cout, ZZ, _size, mn, Arns_tmp, 2*mn);
+		std::cout<<"Arns2:=";
+    FFLAS::WriteMatrix(std::cout, ZZ, _size, mn, Arns_tmp2, 2*mn);
+#endif
+		
+		integer hM= (_M-1)>>1;
+		double *A_beta= FFLAS::fflas_new<double>(2*mn*_ldm);
+		double *A_beta2 = A_beta+mn*_ldm;
+		tsplit.stop();
+#ifdef BENCH_RNS			
+		if(m>1 && n>1) std::cerr<<"RNS EXTENDED double (From) - split : "<<tsplit.usertime()<<std::endl;
+#endif
+
+		
+		Givaro::Timer tfgemmc;tfgemmc.start();
+		if (RNS_MAJOR==false)
+      // compute A_beta = Ap^T x M_beta
+			FFLAS::fgemm(Givaro::ZRing<double>(),FFLAS::FflasTrans, FFLAS::FflasNoTrans,(int) mn*2,(int) _ldm,(int) _size, 1.0 , Arns_tmp,(int) mn*2, _crt_out.data(),(int) _ldm, 0., A_beta,(int)_ldm,
+                   FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive >());
+		else // compute A_beta = Ap x M_Beta
+			cblas_dgemm(CblasRowMajor,CblasNoTrans, CblasNoTrans, (int)mn*2, (int)_ldm, (int)_size, 1.0 , Arns_tmp, (int)_size, _crt_out.data(), (int)_ldm, 0., A_beta,(int)_ldm);
+
+		FFLAS::fflas_delete( Arns_tmp);
+		tfgemmc.stop();
+#ifdef RNS_DEBUG
+		std::cout<<"Abeta1:=";write_field(ZZ,std::cout,A_beta, _ldm, mn,mn,true);
+		std::cout<<"Abeta2:=";write_field(ZZ,std::cout,A_beta2, _ldm, mn,mn,true);		
+#endif
+
+
+#ifdef BENCH_RNS
+		if(m>1 && n>1) std::cerr<<"RNS EXTENDED double (From) - fgemm : "<<tfgemmc.usertime()<<std::endl;
+#endif
+    // compute A using inverse Kronecker transform of A_beta expressed in base 2^log_beta
+		integer* Aiter= A;
+		size_t k=_ldm;
+		size_t k4=((k+3)>>2)+ (((k+3)%4==0)?0:1);
+		std::vector<uint16_t> A0(k4<<2,0),A1(k4<<2,0),A2(k4<<2,0),A3(k4<<2,0);
+		integer a0,a1,a2,a3,res,res2;
+		mpz_t *m0,*m1,*m2,*m3;
+		m0= reinterpret_cast<mpz_t*>(&a0);
+		m1= reinterpret_cast<mpz_t*>(&a1);
+		m2= reinterpret_cast<mpz_t*>(&a2);
+		m3= reinterpret_cast<mpz_t*>(&a3);
+		mp_limb_t *m0_d,*m1_d,*m2_d,*m3_d;
+		m0_d = m0[0]->_mp_d;
+		m1_d = m1[0]->_mp_d;
+		m2_d = m2[0]->_mp_d;
+		m3_d = m3[0]->_mp_d;
+		m0[0]->_mp_alloc = m1[0]->_mp_alloc = m2[0]->_mp_alloc = m3[0]->_mp_alloc = (int) (k4*8/sizeof(mp_limb_t)); // to ensure 32 bits portability
+		m0[0]->_mp_size  = m1[0]->_mp_size  = m2[0]->_mp_size  = m3[0]->_mp_size  = (int) (k4*8/sizeof(mp_limb_t)); // to ensure 32 bits portability
+		Givaro::Timer tkroc;
+		tkroc.start();
+    //		auto sp=SPLITTER();
+    //		PARFOR1D(i,m,sp,
+		for(size_t i=0;i<m;i++)
+			for (size_t j=0;j<n;j++){
+				size_t idx=i*n+j;
+				for (size_t l=0;l<k;l++){
+					uint64_t tmp=(uint64_t)A_beta[l+idx*k];
+					uint16_t* tptr= reinterpret_cast<uint16_t*>(&tmp);
+					A0[l  ]= tptr[0];
+					A1[l+1]= tptr[1];
+					A2[l+2]= tptr[2];
+					A3[l+3]= tptr[3];
+				}
+        // see A0,A1,A2,A3 as a the gmp integers a0,a1,a2,a3
+				m0[0]->_mp_d= reinterpret_cast<mp_limb_t*>(&A0[0]);
+				m1[0]->_mp_d= reinterpret_cast<mp_limb_t*>(&A1[0]);
+				m2[0]->_mp_d= reinterpret_cast<mp_limb_t*>(&A2[0]);
+				m3[0]->_mp_d= reinterpret_cast<mp_limb_t*>(&A3[0]);
+				res = a0;res+= a1;res+= a2;res+= a3;
+
+		
+				for (size_t l=0;l<k;l++){
+				  uint64_t tmp=(uint64_t)A_beta2[l+idx*k];
+					uint16_t* tptr= reinterpret_cast<uint16_t*>(&tmp);
+					A0[l  ]= tptr[0];
+					A1[l+1]= tptr[1];
+					A2[l+2]= tptr[2];
+					A3[l+3]= tptr[3];
+				}
+        // see A0,A1,A2,A3 as a the gmp integers a0,a1,a2,a3
+				m0[0]->_mp_d= reinterpret_cast<mp_limb_t*>(&A0[0]);
+				m1[0]->_mp_d= reinterpret_cast<mp_limb_t*>(&A1[0]);
+				m2[0]->_mp_d= reinterpret_cast<mp_limb_t*>(&A2[0]);
+				m3[0]->_mp_d= reinterpret_cast<mp_limb_t*>(&A3[0]);
+				res2 = a0;res2+= a1;res2+= a2;res2+= a3;
+#ifdef RNS_DEBUG
+				std::cout<<"res1:="<<res<<";\n";
+				std::cout<<"res2:="<<res2<<";\n";
+#endif
+        res+=(res2<<_shift);
+				res%=_M;
+				
+				// get the correct result according to the expected sign of A
+				if (res>hM)
+					res-=_M;
+				if (gamma==0)
+					Aiter[j+i*lda]=res;
+				else
+					if (gamma==integer(1))
+						Aiter[j+i*lda]+=res;
+					else
+						if (gamma==integer(-1))
+							Aiter[j+i*lda]=res-Aiter[j+i*lda];
+						else{
+							Aiter[j+i*lda]*=gamma;
+							Aiter[j+i*lda]+=res;
+						}
+
+			}
+    tkroc.stop();
+#ifdef BENCH_RNS			
+    if(m>1 && n>1) std::cerr<<"RNS EXTENDED double (From) -  Convert : "<<tkroc.usertime()<<std::endl;
+#endif
+		m0[0]->_mp_d = m0_d;
+		m1[0]->_mp_d = m1_d;
+		m2[0]->_mp_d = m2_d;
+		m3[0]->_mp_d = m3_d;
+		m0[0]->_mp_alloc = m1[0]->_mp_alloc = m2[0]->_mp_alloc= m3[0]->_mp_alloc = 1;
+		m0[0]->_mp_size  = m1[0]->_mp_size  = m2[0]->_mp_size = m3[0]->_mp_size  = 0;
+		FFLAS::fflas_delete( A_beta);
+
+
+#ifdef CHECK_RNS
+		bool ok=true;
+		for (size_t i=0;i<m;i++)
+		  for(size_t j=0;j<n;j++)
+		    for(size_t k=0;k<_size;k++){
+		      int64_t _p =(int64_t) _basis[k];
+		      integer curr=A[i*lda+j] - gamma*Acopy[i*n+j];
+		      ok&= ( curr% _p +(curr%_p<0?_p:0) == (int64_t) Arns[i*n+j+k*rda]);
+		      if (!ok) std::cout<<A[i*lda+j]<<" mod "<<(int64_t) _basis[k]<<"="<<(int64_t) Arns[i*n+j+k*rda]<<";"<<std::endl;
+				}
+		std::cout<<"RNS EXTENDED double (From) ... "<<(ok?"OK":"ERROR")<<std::endl;
+#endif
+
+	}
+
+  
+
+}// end of namespace FFPACK
+
+
+#endif // __FFLASFFPACK_field_rns_double_extended_INL
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
diff --git a/fflas-ffpack/field/rns-double.h b/fflas-ffpack/field/rns-double.h
index 8286c8768..447f45796 100644
--- a/fflas-ffpack/field/rns-double.h
+++ b/fflas-ffpack/field/rns-double.h
@@ -52,413 +52,219 @@
 
 namespace FFPACK {
 
-    /* Structure that handles rns representation given a bound and bitsize for prime moduli
-     * support sign representation (i.e. the bound must be twice larger then ||A||)
-     */
-    struct rns_double {
-        typedef Givaro::Integer integer;
-        typedef Givaro::Modular<double> ModField;
-
-        std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _basis; // the rns moduli (mi)
-        std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _basisMax; // (mi-1)
-        std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _negbasis; // (-mi)
-        std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _invbasis; // the inverse of rns moduli (1/mi)
-        std::vector<ModField> _field_rns; // the associated prime field for each mi
-        integer                  _M; // the product of the mi's
-        std::vector<integer>         _Mi; // _M/mi
-        std::vector<double>         _MMi; // (_Mi)^(-1) mod mi
-        std::vector<double>      _crt_in; //  2^(16*j) mod mi
-        std::vector<double>     _crt_out; //  (_Mi._MMi) written in base 2^16
-        size_t                _size; // the size of the rns basis (number of mi's)
-        size_t               _pbits; // the size in bit of the mi's
-        size_t                 _ldm; // log[2^16](_M)
-        integer                  _mi_sum; // the product of the mi's
-
-        typedef double                        BasisElement;
-        typedef rns_double_elt                     Element;
-        typedef rns_double_elt_ptr             Element_ptr;
-        typedef rns_double_elt_cstptr     ConstElement_ptr;
-
-        rns_double(const integer& bound, size_t pbits, bool rnsmod=false, long seed=time(NULL))
-        :  _M(1), _size(0), _pbits(pbits), _mi_sum(1)
-        {
-            integer::seeding(seed);
-            Givaro::IntPrimeDom IPD;
-            integer prime;
-            while (_M < bound*_mi_sum) {
-                _basis.resize(_size+1);
-                do {
-                    integer::random_exact_2exp(prime, _pbits-1);
-                    IPD.nextprimein(prime);
-                } while (_M%prime == 0);
-                _basis[_size]=prime;
-                _size++;
-                _M*=prime;
-                if (rnsmod) _mi_sum+=prime;
-            }
-            // std::ostream_iterator<uint64_t> out_it (std::cout,", ");
-            // std::cout<<"RNS basis =";
-            // std::copy ( _basis.begin(), _basis.end(), out_it );
-            // std::cout<<std::endl;
-            // std::cout<<"RNS sum Mi ="<<_mi_sum<<"\n";
-            precompute_cst();
-        }
-
-        rns_double(size_t pbits, size_t size, long seed=time(NULL))
-        :  _M(1), _size(size), _pbits(pbits), _mi_sum(1)
-        {
-            integer::seeding(seed);
-            Givaro::IntPrimeDom IPD;
-            integer prime;
-            _basis.resize(size);
-            _negbasis.resize(size);
-            _basisMax.resize(size);
-            for(size_t i = 0 ; i < _size ; ++i){
-                integer::random_exact_2exp(prime, _pbits-1);
-                IPD.nextprimein(prime);
-                _basis[i]=prime;
-                _basisMax[i] = prime-1;
-                _negbasis[i] = 0-prime;
-                _M*=prime;
-            }
-            precompute_cst();
-        }
-
-        template<typename Vect>
-        rns_double(const Vect& basis, bool rnsmod=false, long seed=time(NULL))
-        :  _basis(basis.begin(),basis.end()), _basisMax(basis.size()), _negbasis(basis.size()), _M(1), _size(basis.size()), _pbits(0), _mi_sum(1)
-        {
-            for(size_t i=0;i<_size;i++){
-                //std::cout<<"basis["<<i<<"]="<<_basis[i]<<std::endl;
-                _M*=_basis[i];
-                _pbits=std::max(_pbits, integer(_basis[i]).bitsize());
-            }
-            //std::cout<<"M="<<_M<<std::endl;
-            precompute_cst();
-        }
-
-        rns_double(const RNSIntegerMod<rns_double>& basis, bool rnsmod=false, long seed=time(NULL)) {
-
-        }
-
-        // can force to reduce integer entries larger than M
-        void precompute_cst(size_t K=0){
-            if (K!=0)
-                _ldm=K;
-            else
-                _ldm = (_M.bitsize()/16) + ((_M.bitsize()%16)?1:0) ;
-            _invbasis.resize(_size);
-            _field_rns.resize(_size);
-            _Mi.resize(_size);
-            _MMi.resize(_size);
-            _basisMax.resize(_size);
-            _negbasis.resize(_size);
-            _crt_in.resize(_size*_ldm);
-            _crt_out.resize(_size*_ldm);
-            //const unsigned int MASK=0xFFFF;
-            //Givaro::Timer chrono;
-            //double t1=0.,t2=0.,t3=0.;
-
-            for (size_t i=0;i<_size;i++){
-                //chrono.start();
-                _invbasis[i]  = 1./_basis[i];
-                _basisMax[i] = _basis[i]-1;
-                _negbasis[i] = 0-_basis[i];
-                _field_rns[i] = ModField(_basis[i]);
-                _Mi[i]        = _M/(uint64_t)_basis[i];
-                _field_rns[i].init(_MMi[i], _Mi[i] % (double)_basis[i]);
-                _field_rns[i].invin(_MMi[i]);
-                integer tmp= _Mi[i]*(uint64_t)_MMi[i];
-                const mpz_t*    m0     = reinterpret_cast<const mpz_t*>(&tmp);
-                const uint16_t* m0_ptr = reinterpret_cast<const uint16_t*>(m0[0]->_mp_d);
-                size_t maxs=std::min(_ldm,(tmp.size())*sizeof(mp_limb_t)/2);// to ensure 32 bits portability
-                //chrono.stop();
-                //t1+=chrono.usertime();
-                //chrono.start();
-                /*
-                   for(size_t j=0;j<_ldm;j++){
-                   _crt_out[j+i*_ldm]=double(tmp[0]&MASK);
-                   tmp>>=16; // Bad idea -> too slow (must get the lowest limb of the integer)
-
-                   }
-                   */
-                size_t l=0;
+  /* Structure that handles rns representation with moduli stored in double FP numbers
+   * support sign representation (i.e. the bound must be twice larger then ||A||)
+   */
+  struct rns_double {
+    typedef Givaro::Integer integer;
+    typedef Givaro::Modular<double> ModField;
+
+    std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _basis; // the rns moduli (mi)
+    std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _basisMax; // (mi-1)
+    std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _negbasis; // (-mi)
+    std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _invbasis; // the inverse of rns moduli (1/mi)
+    std::vector<ModField> _field_rns; // the associated prime field for each mi
+    integer                  _M; // the product of the mi's
+    std::vector<integer>         _Mi; // _M/mi
+    std::vector<double>         _MMi; // (_Mi)^(-1) mod mi
+    std::vector<double>      _crt_in; //  2^(16*j) mod mi
+    std::vector<double>     _crt_out; //  (_Mi._MMi) written in base 2^16
+    size_t                _size; // the size of the rns basis (number of mi's)
+    size_t               _pbits; // the size in bit of the mi's
+    size_t                 _ldm; // log[2^16](_M)
+    integer                  _mi_sum; // the product of the mi's
+
+    typedef double                        BasisElement;
+    typedef rns_double_elt                     Element;
+    typedef rns_double_elt_ptr             Element_ptr;
+    typedef rns_double_elt_cstptr     ConstElement_ptr;
+
+
+    // construct an RNS basis with primes of bit-length (pbits) ensuring to represent integers lying in [0, bound[
+    // Rmk: when rnsmod is set to true the RNS basis (m1,m2, ..., mk) satisfies that : m1*m2*...*mk >= bound * (m1+m2+...+mk)
+    rns_double(const integer& bound, size_t pbits, bool rnsmod=false, long seed=time(NULL))
+      :  _M(1), _size(0), _pbits(pbits), _mi_sum(1)
+    {
+      integer::seeding(seed);
+      Givaro::IntPrimeDom IPD;
+      integer prime;
+      while (_M < bound*_mi_sum) {
+        _basis.resize(_size+1);
+        do {
+          integer::random_exact_2exp(prime, _pbits-1);
+          IPD.nextprimein(prime);
+        } while (_M%prime == 0);
+        _basis[_size]=prime;
+        _size++;
+        _M*=prime;
+        if (rnsmod) _mi_sum+=prime;
+      }
+      precompute_cst();
+    }
+
+    // construct an RNS basis with (size) primes of bit-length (pbits) 
+    rns_double(size_t pbits, size_t size, long seed=time(NULL))
+      :  _M(1), _size(size), _pbits(pbits), _mi_sum(1)
+    {
+      integer::seeding(seed);
+      Givaro::IntPrimeDom IPD;
+      integer prime;
+      _basis.resize(size);
+      _negbasis.resize(size);
+      _basisMax.resize(size);
+      for(size_t i = 0 ; i < _size ; ++i){
+        integer::random_exact_2exp(prime, _pbits-1);
+        IPD.nextprimein(prime);
+        _basis[i]=prime;
+        _basisMax[i] = prime-1;
+        _negbasis[i] = 0-prime;
+        _M*=prime;
+      }
+      precompute_cst();
+    }
+
+    // construct an RNS basis from a vector of relatively prime numbers (basis)
+    template<typename Vect>
+    rns_double(const Vect& basis, long seed=time(NULL))
+      :  _basis(basis.begin(),basis.end()), _basisMax(basis.size()), _negbasis(basis.size()), _M(1), _size(basis.size()), _pbits(0), _mi_sum(1)
+    {
+      for(size_t i=0;i<_size;i++){
+        _M*=_basis[i];
+        _pbits=std::max(_pbits, integer(_basis[i]).bitsize());
+      }
+      precompute_cst();
+    }
+
+    rns_double(const RNSIntegerMod<rns_double>& basis, bool rnsmod=false, long seed=time(NULL)) {
+
+    }
+
+    // can force to reduce integer entries larger than M
+    void precompute_cst(size_t K=0){
+
+      // Check that _pbits satisfies log(_M)/16* 2^(_pbits) * 2^16 < 2^53 (Required for correctness)
+      // => pbits <= 41 - loglog(_M)
+      if ( _pbits > 41 - log(double(_M.bitsize()))/log(2.)){
+        std::cout<<"FFLAS Error in rns_double: primes bitsize is too large ... aborting\n";
+        std::terminate();
+      }
+
+      if (K!=0)
+        _ldm=K;
+      else
+        _ldm = (_M.bitsize()/16) + ((_M.bitsize()%16)?1:0) ;
+      _invbasis.resize(_size);
+      _field_rns.resize(_size);
+      _Mi.resize(_size);
+      _MMi.resize(_size);
+      _basisMax.resize(_size);
+      _negbasis.resize(_size);
+      _crt_in.resize(_size*_ldm);
+      _crt_out.resize(_size*_ldm);
+      for (size_t i=0;i<_size;i++){
+        _invbasis[i]  = 1./_basis[i];
+        _basisMax[i] = _basis[i]-1;
+        _negbasis[i] = 0-_basis[i];
+        _field_rns[i] = ModField(_basis[i]);
+        _Mi[i]        = _M/(uint64_t)_basis[i];
+        _field_rns[i].init(_MMi[i], _Mi[i] % (double)_basis[i]);
+        _field_rns[i].invin(_MMi[i]);
+        integer tmp= _Mi[i]*(uint64_t)_MMi[i];
+        const mpz_t*    m0     = reinterpret_cast<const mpz_t*>(&tmp);
+        const uint16_t* m0_ptr = reinterpret_cast<const uint16_t*>(m0[0]->_mp_d);
+        size_t maxs=std::min(_ldm,(tmp.size())*sizeof(mp_limb_t)/2);// to ensure 32 bits portability
+
+        size_t l=0;
 #ifdef __FFLASFFPACK_HAVE_LITTLE_ENDIAN
-                for(;l<maxs;l++)
-                    _crt_out[l+i*_ldm]=m0_ptr[l];
+        for(;l<maxs;l++)
+          _crt_out[l+i*_ldm]=m0_ptr[l];
 #else
-                for(;l<maxs;l++)
-                    _crt_out[l+i*_ldm]=m0_ptr[l ^ ((sizeof(mp_limb_t)/2U) - 1U)];
+        for(;l<maxs;l++)
+          _crt_out[l+i*_ldm]=m0_ptr[l ^ ((sizeof(mp_limb_t)/2U) - 1U)];
 #endif
-                for(;l<_ldm;l++)
-                    _crt_out[l+i*_ldm]=0.;;
-                // chrono.stop();
-                // t2+=chrono.usertime();
-                // chrono.start();
-                double beta=double(1<<16);
-                double  acc=1;
-                for(size_t j=0;j<_ldm;j++){
-                    _crt_in[j+i*_ldm]=acc;
-                    _field_rns[i].mulin(acc,beta);
-                }
-                // chrono.stop();
-                // t3+=chrono.usertime();
-
-            }
-            // std::cout<<"t1="<<t1<<std::endl;
-            // std::cout<<"t2="<<t2<<std::endl;
-            // std::cout<<"t3="<<t3<<std::endl;
+        for(;l<_ldm;l++)
+          _crt_out[l+i*_ldm]=0.;;
+        double beta=double(1<<16);
+        double  acc=1;
+        for(size_t j=0;j<_ldm;j++){
+          _crt_in[j+i*_ldm]=acc;
+          _field_rns[i].mulin(acc,beta);
         }
+      }
+    }
 
-        // Arns must be an array of m*n*_size
-        // abs(||A||) <= maxA
-        template<typename T>
-        void init(size_t m, size_t n, double* Arns, size_t rda, const T* A, size_t lda,
-                  const integer& maxA, bool RNS_MAJOR=false) const
-        {
-            init(m,n,Arns,rda,A,lda, maxA.bitsize()/16 + (maxA.bitsize()%16?1:0),RNS_MAJOR);
-        }
+    // rda is the distance between two consecutive residues : i.e. A_ij mod mk and A_ij mod m_(k+1) -> it is either 1 (RNS_MAJOR) or mxn (NOT RNS_MAJOR)
+    void init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false, const FFLAS::FFLAS_TRANSPOSE ta=FFLAS::FflasNoTrans) const;
+    void convert(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false, const FFLAS::FFLAS_TRANSPOSE ta=FFLAS::FflasNoTrans) const;
 
-        void init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false) const;
-        void init_transpose(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false) const;
-        void convert(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false) const;
-        void convert_transpose(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false) const;
+    // Arns must be an array of m*n*_size with  abs(||A||) <= maxA    
+    void init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, const integer& maxA, bool RNS_MAJOR=false) const{
+      init(m,n,Arns,rda,A,lda, maxA.bitsize()/16 + (maxA.bitsize()%16?1:0),RNS_MAJOR, FFLAS::FflasNoTrans);
+    }
+    void init_transpose(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, const integer& maxA, bool RNS_MAJOR=false) const{
+      init(m,n,Arns,rda,A,lda, maxA.bitsize()/16 + (maxA.bitsize()%16?1:0),RNS_MAJOR, FFLAS::FflasTrans);
+    }
+    void init_transpose(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false) const{
+      init(m,n,Arns,rda,A,lda,k,RNS_MAJOR, FFLAS::FflasTrans);
+    }
+    void convert_transpose(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false) const {
+      convert(m,n,gamma,A,lda,Arns,rda,RNS_MAJOR,FFLAS::FflasTrans);
+    }
 
-        // reduce entries of Arns to be less than the rns basis elements
-        void reduce(size_t n, double* Arns, size_t rda, bool RNS_MAJOR=false) const;
 
-        template<size_t K>
-        void init(size_t m, size_t n, double* Arns, size_t rda, const RecInt::ruint<K>* A, size_t lda, size_t k, bool RNS_MAJOR=false) const;
-        template<size_t K>
-        void convert(size_t m, size_t n, integer gamma, RecInt::ruint<K>* A, size_t lda, const double* Arns, size_t rda, integer p=0,bool RNS_MAJOR=false) const;
+    // reduce entries of Arns to be less than the rns basis elements
+    void reduce(size_t n, double* Arns, size_t rda, bool RNS_MAJOR=false) const;
 
+    template<size_t K>
+    void init(size_t m, size_t n, double* Arns, size_t rda, const RecInt::ruint<K>* A, size_t lda, size_t k, bool RNS_MAJOR=false) const;
+    template<size_t K>
+    void convert(size_t m, size_t n, integer gamma, RecInt::ruint<K>* A, size_t lda, const double* Arns, size_t rda, integer p=0,bool RNS_MAJOR=false) const;
 
-    }; // end of struct rns_double
 
-    /* Structure that handles rns representation given a bound and bitsize for prime moduli, allow large moduli
-     * support sign representation (i.e. the bound must be twice larger then ||A||)
-     */
-    struct rns_double_extended {
-        typedef Givaro::Integer integer;
-        typedef Givaro::ModularExtended<double> ModField;
-
-        std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _basis; // the rns moduli (mi)
-        std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _basisMax; // (mi-1)
-        std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _negbasis; // (-mi)
-        std::vector<double, AlignedAllocator<double, Alignment::CACHE_LINE>>       _invbasis; // the inverse of rns moduli (1/mi)
-        std::vector<ModField> _field_rns; // the associated prime field for each mi
-        integer                  _M; // the product of the mi's
-        std::vector<integer>         _Mi; // _M/mi
-        std::vector<double>         _MMi; // (_Mi)^(-1) mod mi
-        std::vector<double>      _crt_in; //  2^(16*j) mod mi
-        std::vector<double>     _crt_out; //  (_Mi._MMi) written in base 2^16
-        size_t                _size; // the size of the rns basis (number of mi's)
-        size_t               _pbits; // the size in bit of the mi's
-        size_t                 _ldm; // log[2^16](_M)
-
-        typedef double                        BasisElement;
-        typedef rns_double_elt                     Element;
-        typedef rns_double_elt_ptr             Element_ptr;
-        typedef rns_double_elt_cstptr     ConstElement_ptr;
-
-        rns_double_extended(const integer& bound, size_t pbits, bool rnsmod=false, long seed=time(NULL))
-        :  _M(1), _size(0), _pbits(pbits)
-        {
-            integer::seeding(seed);
-            integer prime; Givaro::IntPrimeDom IPD;
-            integer sum=1;
-            while (_M < bound*sum) {
-                _basis.resize(_size+1);
-                do {
-                    integer::random_exact_2exp(prime, _pbits-1);
-                    IPD.nextprimein(prime);
-                } while (_M%prime == 0);
-                _basis[_size]=prime;
-                _size++;
-                _M*=prime;
-                if (rnsmod) sum+=prime;
-            }
-            precompute_cst();
-        }
-
-        rns_double_extended(size_t pbits, size_t size, long seed=time(NULL))
-        :  _M(1), _size(size), _pbits(pbits)
-        {
-            integer::seeding(seed);
-            integer prime; Givaro::IntPrimeDom IPD;
-            integer sum=1;
-            _basis.resize(size);
-            _negbasis.resize(size);
-            _basisMax.resize(size);
-            for(size_t i = 0 ; i < _size ; ++i){
-                integer::random_exact_2exp(prime, _pbits-1);
-                IPD.nextprimein(prime);
-                _basis[i]=prime;
-                _basisMax[i] = prime-1;
-                _negbasis[i] = 0-prime;
-                _M*=prime;
-            }
-            precompute_cst();
-        }
+  }; // end of struct rns_double
 
-        template<typename Vect>
-        rns_double_extended(const Vect& basis, bool rnsmod=false, long seed=time(NULL))
-        :  _basis(basis.begin(),basis.end()), _basisMax(basis.size()), _negbasis(basis.size()), _M(1), _size(basis.size()), _pbits(0)
-        {
-            for(size_t i=0;i<_size;i++){
-                //std::cout<<"basis["<<i<<"]="<<_basis[i]<<std::endl;
-                _M*=_basis[i];
-                _pbits=std::max(_pbits, integer(_basis[i]).bitsize());
-            }
-            //std::cout<<"M="<<_M<<std::endl;
-            precompute_cst();
-        }
 
 
-        void precompute_cst(){
-            _ldm = (_M.bitsize()/16) + ((_M.bitsize()%16)?1:0) ;
-            _invbasis.resize(_size);
-            _basisMax.resize(_size);
-            _negbasis.resize(_size);
-            _field_rns.resize(_size);
-            _Mi.resize(_size);
-            _MMi.resize(_size);
-            _crt_in.resize(_size*_ldm);
-            _crt_out.resize(_size*_ldm);
-            const unsigned int MASK=0xFFFF;
-            for (size_t i=0;i<_size;i++){
-                _invbasis[i]  = 1./_basis[i];
-                _basisMax[i] = _basis[i]-1;
-                _negbasis[i] = 0-_basis[i];
-                _field_rns[i] = ModField(_basis[i]);
-                _Mi[i]        = _M/(uint64_t)_basis[i];
-                _field_rns[i].init(_MMi[i], _Mi[i] % (double)_basis[i]);
-                _field_rns[i].invin(_MMi[i]);
-                integer tmp= _Mi[i]*(uint64_t)_MMi[i];
-                for(size_t j=0;j<_ldm;j++){
-                    _crt_out[j+i*_ldm]=double(tmp&MASK);
-                    tmp>>=16;
-                }
-                double beta=double(1<<16);
-                double  acc=1;
-                for(size_t j=0;j<_ldm;j++){
-                    _crt_in[j+i*_ldm]=acc;
-                    _field_rns[i].mulin(acc,beta);
-                }
-            }
-        }
+  template<typename RNS>
+  class rnsRandIter {
+    std::vector<typename RNS::ModField::RandIter> _RNS_rand;
+    const RNS& _domain;
 
-        // Arns must be an array of m*n*_size
-        // abs(||A||) <= maxA
-        void init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda,
-                  const integer& maxA, bool RNS_MAJOR=false) const
-        {
-            init(m*n,Arns,A,lda);
-        }
-
-        void init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR=false){
-            init(m*n,Arns,A,lda);
-        }
-        void convert(size_t m, size_t n, integer gamma, integer* A, size_t lda, const double* Arns, size_t rda, bool RNS_MAJOR=false){
-            convert(m*n, A, Arns);
-        }
-        void init(size_t m, double* Arns, const integer* A, size_t lda) const;
-        void convert(size_t m, integer *A, const double *Arns) const;
-
-#if defined(__FFLASFFPACK_HAVE_SSE4_1_INSTRUCTIONS)
-
-        template<class SimdT>
-        inline void splitSimd(const SimdT x, SimdT & x_h, SimdT & x_l) const {
-            using simd = Simd<double>;
-            using vect_t = typename simd::vect_t;
-            vect_t vc = simd::set1((double)((1 << 27)+1));
-            vect_t tmp = simd::mul(vc, x);
-            x_h = simd::add(tmp, simd::sub(x, tmp));
-            x_l = simd::sub(x, x_h);
-        }
-
-        template<class SimdT>
-        inline void multSimd(const SimdT va, const SimdT vb, SimdT & vs, SimdT & vt) const{
-            using simd = Simd<double>;
-            using vect_t = typename simd::vect_t;
-            vect_t vah, val, vbh, vbl;
-            vs = simd::mul(va, vb);
-            //#ifdef __FMA__
-            vt = simd::fnmadd(va, vb, vs);
-            //#else
-            splitSimd(va, vah, val);
-            splitSimd(vb, vbh, vbl);
-            vt = simd::add(simd::add(simd::sub(simd::mul(vah, vbh), vs), simd::mul(vah, vbl)), simd::add(simd::mul(val, vbh), simd::mul(val, vbl)));
-            //#endif
-        }
-
-        template<class SimdT>
-        inline SimdT modSimd(const SimdT a, const SimdT p, const SimdT ip, const SimdT np) const{
-            using simd = Simd<double>;
-            using vect_t = typename simd::vect_t;
-            vect_t pqh, pql, abl, abh;
-            vect_t q = simd::floor(simd::mul(a, ip));
-            multSimd(p, q, pqh, pql);
-            vect_t r = simd::add(simd::sub(a, pqh), pql);
-            abh = simd::greater_eq(r, p);
-            abl = simd::lesser(r, simd::zero());
-            abh = simd::vand(abh, np);
-            abl = simd::vand(abl, p);
-            abh = simd::vor(abh, abl);
-            return r = simd::add(r, abh);
-        }
-
-#endif // __FFLASFFPACK_HAVE_SSE4_1_INSTRUCTIONS
-
-        // reduce entries of Arns to be less than the rns basis elements
-        void reduce(size_t n, double* Arns, size_t rda, bool RNS_MAJOR=false) const;
-
-
-
-    }; // end of struct rns_double_extended
-
-
-    template<typename RNS>
-    class rnsRandIter {
-        std::vector<typename RNS::ModField::RandIter> _RNS_rand;
-        const RNS& _domain;
-
-    public:
-        rnsRandIter(const RNS& R, uint64_t seed=0)
-        : _domain(R) {
-            for(const auto& F : R._field_rns)
-                _RNS_rand.emplace_back(F,seed);
-        }
-
-        /** RNS ring Element random assignement.
-         * Element is supposed to be initialized
-         * @return random ring Element
-         */
-        typename RNS::Element& random(typename RNS::Element& elt) const {
-            auto coefficient(elt._ptr);
-            for(auto & iter : _RNS_rand) {
-                iter.random( *coefficient );
-                coefficient += elt._stride;
-            }
-            return elt;
-        }
-
-        typename RNS::Element& operator()(typename RNS::Element& elt) const {
-            return this->random(elt);
-        }
-
-        typename RNS::Element operator()() const {
-            typename RNS::Element tmp; _domain.init(tmp);
-            return this->operator()(tmp);
-        }
-        typename RNS::Element random() const {
-            return this->operator()();
-        }
-
-        const RNS& ring() const { return _domain; }
-
-    };
+  public:
+    rnsRandIter(const RNS& R, uint64_t seed=0)
+      : _domain(R) {
+      for(const auto& F : R._field_rns)
+        _RNS_rand.emplace_back(F,seed);
+    }
 
+    /** RNS ring Element random assignement.
+     * Element is supposed to be initialized
+     * @return random ring Element
+     */
+    typename RNS::Element& random(typename RNS::Element& elt) const {
+      auto coefficient(elt._ptr);
+      for(auto & iter : _RNS_rand) {
+        iter.random( *coefficient );
+        coefficient += elt._stride;
+      }
+      return elt;
+    }
+
+    typename RNS::Element& operator()(typename RNS::Element& elt) const {
+      return this->random(elt);
+    }
+
+    typename RNS::Element operator()() const {
+      typename RNS::Element tmp; _domain.init(tmp);
+      return this->operator()(tmp);
+    }
+    typename RNS::Element random() const {
+      return this->operator()();
+    }
+
+    const RNS& ring() const { return _domain; }
+
+  };
 
 } // end of namespace FFPACK
 
@@ -466,10 +272,10 @@ namespace FFPACK {
 #include "rns-double-recint.inl"
 namespace FFLAS {
 
-    template<>
-    inline void fflas_delete (FFPACK::rns_double_elt_ptr A) {FFLAS::fflas_delete( A._ptr);}
-    template<>
-    inline void fflas_delete (FFPACK::rns_double_elt_cstptr A) {delete[] A._ptr;}
+  template<>
+  inline void fflas_delete (FFPACK::rns_double_elt_ptr A) {FFLAS::fflas_delete( A._ptr);}
+  template<>
+  inline void fflas_delete (FFPACK::rns_double_elt_cstptr A) {delete[] A._ptr;}
 
 }
 
diff --git a/fflas-ffpack/field/rns-double.inl b/fflas-ffpack/field/rns-double.inl
index 7713129d3..deec21b4f 100644
--- a/fflas-ffpack/field/rns-double.inl
+++ b/fflas-ffpack/field/rns-double.inl
@@ -29,605 +29,182 @@
 #define __FFLASFFPACK_field_rns_double_INL
 
 #include "fflas-ffpack/fflas/fflas_freduce.h"
+#include "fflas-ffpack/field/rns-common.h"
 
-namespace FFPACK {
 
-    // Arns must be an array of m*n*_size
-    // abs(||A||) < 2^(16k)
-    inline void rns_double::init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR) const
-    {
-        if (k>_ldm){
-            FFPACK::failure()(__func__,__FILE__,__LINE__,"rns_double [init] -> rns basis is too small to handle integers with 2^(16*k) values ");
-            std::cerr<<"with k="<<k<<" _ldm="<<_ldm<<std::endl;
-        }
-        const size_t mn=m*n;
-        if (mn) {
-        double *A_beta = FFLAS::fflas_new<double >(mn*k);
-        const integer* Aiter=A;
-        // split A into A_beta according to a Kronecker transform in base 2^16
-        //		auto sp=SPLITTER(MAX_THREADS,FFLAS::CuttingStrategy::Column,FFLAS::StrategyParameter::Threads);
-
-        Givaro::Timer tkr; tkr.start();
-        // #ifndef __FFLASFFPACK_SEQUENTIAL
-        // 			auto sp=SPLITTER(MAX_THREADS);
-        // #else
-        // 			auto sp=SPLITTER(1);
-        // #endif
-        // FOR2D(i,j,m,n,sp,
-        //       TASK(MODE(READ(Aiter[0]) READWRITE(A_beta[0])),
-        //for(size_t i=0;i<m;i++)
-        //PAR_BLOCK{
-        //			FOR1D(i,m,sp,
-        PARFOR1D(i,m,SPLITTER(NUM_THREADS),
-
-                 for(size_t j=0;j<n;j++){
-                 size_t idx=j+i*n;
-                 const mpz_t*    m0     = reinterpret_cast<const mpz_t*>(Aiter+j+i*lda);
-                 const uint16_t* m0_ptr = reinterpret_cast<const uint16_t*>(m0[0]->_mp_d);
-                 size_t l=0;
-                 //size_t maxs=std::min(k,(Aiter[j+i*lda].size())<<2);
-                 size_t maxs=std::min(k,(Aiter[j+i*lda].size())*sizeof(mp_limb_t)/2);// to ensure 32 bits portability
-
-#ifdef __FFLASFFPACK_HAVE_LITTLE_ENDIAN
-                 if (m0[0]->_mp_size >= 0)
-                 for (;l<maxs;l++)
-                 A_beta[l+idx*k]=  m0_ptr[l];
-                 else
-                 for (;l<maxs;l++)
-                 A_beta[l+idx*k]= - double(m0_ptr[l]);
-#else
-                 if (m0[0]->_mp_size >= 0)
-                 for (;l<maxs;l++) {
-                 size_t l2 = l ^ ((sizeof(mp_limb_t)/2U) - 1U);
-                 A_beta[l+idx*k]=  m0_ptr[l2];
-                 }
-                 else
-                     for (;l<maxs;l++) {
-                         size_t l2 = l ^ ((sizeof(mp_limb_t)/2U) - 1U);
-                         A_beta[l+idx*k]= - double(m0_ptr[l2]);
-                     }
-#endif
-                 for (;l<k;l++)
-                     A_beta[l+idx*k]=  0.;
-
-                 // 	   );
-                 }
-        );
-
-        tkr.stop();
-        //if(m>1 && n>1) std::cerr<<"Kronecker : "<<tkr.realtime()<<std::endl;
-        if (RNS_MAJOR==false) {
-            // Arns = _crt_in x A_beta^T
-            Givaro::Timer tfgemm; tfgemm.start();
-#ifndef ENABLE_CHECKER_fgemm
-            FFLAS::fgemm (Givaro::ZRing<double>(), FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in.data(),_ldm,A_beta,k,0.,Arns,rda,
-                          FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive>());
-#else
-            cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasTrans,(int)_size,(int)mn,(int)k,1.0,_crt_in.data(),(int)_ldm,A_beta,(int)k,0.,Arns,(int)rda);
-#endif
-            tfgemm.stop();
-            //if(m>1 && n>1) 	std::cerr<<"fgemm : "<<tfgemm.realtime()<<std::endl;
-        }
-        else {
-            // Arns =  A_beta x _crt_in^T
-#ifndef ENABLE_CHECKER_fgemm
-            FFLAS::fgemm (Givaro::ZRing<double>(), FFLAS::FflasNoTrans,FFLAS::FflasTrans,mn,_size,k,1.0,A_beta, k, _crt_in.data(),_ldm,0.,Arns,_size,
-                          FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive>());
-#else
-            cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasTrans,(int)mn,(int)_size,(int)k,1.0,A_beta,(int)k,_crt_in.data(),(int)_ldm,0.,Arns,(int)_size);
-#endif
-        }
-        Givaro::Timer tred; tred.start();
-
-        reduce(mn,Arns,rda,RNS_MAJOR);
-        tred.stop();
-        //if(m>1 && n>1) 			std::cerr<<"Reduce : "<<tred.realtime()<<std::endl;
-
-        FFLAS::fflas_delete( A_beta);
-
-#ifdef CHECK_RNS
-        bool ok=true;
-        for (size_t i=0;i<m;i++)
-            for(size_t j=0;j<n;j++)
-                for(size_t k=0;k<_size;k++){
-                    ok&= (((A[i*lda+j] % (int64_t) _basis[k])+(A[i*lda+j]<0?(int64_t)_basis[k]:0)) == (int64_t) Arns[i*n+j+k*rda]);
-                    if (((A[i*lda+j] % (int64_t) _basis[k])+(A[i*lda+j]<0?(int64_t)_basis[k]:0))
-                        != (int64_t) Arns[i*n+j+k*rda])
-                    {
-                        std::cout<<((A[i*lda+j] % (int64_t) _basis[k])+(A[i*lda+j]<0?(int64_t)_basis[k]:0))
-                        <<" != "
-                        <<(int64_t) Arns[i*n+j+k*rda]
-                        <<std::endl;
-                    }
-                }
-        std::cout<<"RNS freduce ... "<<(ok?"OK":"ERROR")<<std::endl;
-#endif
-        }
+  // Arns must be an array of m*n*_size
+  // abs(||A||) < 2^(16k)
+  inline void rns_double::init(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR, const FFLAS::FFLAS_TRANSPOSE trans) const
+  {
+    if (k>_ldm){
+      FFPACK::failure()(__func__,__FILE__,__LINE__,"rns_double [init] -> rns basis is too small to handle integers with 2^(16*k) values ");
+      std::cerr<<"with k="<<k<<" _ldm="<<_ldm<<std::endl;
     }
+    const size_t mn=m*n;
+    Givaro::ZRing<double> ZD;
+
+    if (mn) {
+      // split A into A_beta according to a Kronecker transform in base 2^16
+      double *A_beta = FFLAS::fflas_new<double >(mn*k);
+      RNS_COMMON::Kronecker_base16(trans,m,n,A_beta,A,lda,k);
+
+      // Using Helper for potential parallelism -> need to be activated by hand
+      //using ParallelStrategy= FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive>;
+      using ParallelStrategy= FFLAS::ParSeqHelper::Sequential;
+      FFLAS::MMHelper<Givaro::ZRing<double>, FFLAS::MMHelperAlgo::Winograd>  MatMulHelper (ZD, -1, ParallelStrategy());
+
+      // Pseudo reduction of the matrix entries
+      if (RNS_MAJOR==false) 
+        // Arns = _crt_in x A_beta^T
+        FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in.data(),_ldm,A_beta,k,0.,Arns,mn,MatMulHelper);
+      else 
+        // Arns =  A_beta x _crt_in^T
+        FFLAS::fgemm (ZD, FFLAS::FflasNoTrans,FFLAS::FflasTrans,mn,_size,k,1.0,A_beta, k, _crt_in.data(),_ldm,0.,Arns,_size,MatMulHelper);
+
+      // Reduce the pseudo reduce matrix modulo RNS basis
+      reduce(mn,Arns,rda,RNS_MAJOR);
+
+#ifdef CHECK_RNS // check the result
+      RNS_COMMON::check_integers_to_rns(trans, m, n, Arns, rda, A, lda, _basis, RNS_MAJOR);
+#endif
+      FFLAS::fflas_delete(A_beta);
+    }
+  }
 
-    // Arns must be an array of m*n*_size
-    // abs(||A||) < 2^(16k)
-    inline void rns_double::init_transpose(size_t m, size_t n, double* Arns, size_t rda, const integer* A, size_t lda, size_t k, bool RNS_MAJOR) const
-    {
-        if (k>_ldm)
-            FFPACK::failure()(__func__,__FILE__,__LINE__,"rns_struct: init (too large entry)");
-
-        const size_t mn=m*n;
-        if (mn) {
-        double *A_beta = FFLAS::fflas_new<double >(mn*k);
-        const integer* Aiter=A;
-        // split A into A_beta according to a Kronecker transform in base 2^16
-        for(size_t j=0;j<n;j++){
-            for(size_t i=0;i<m;i++){
-                size_t idx=i+j*m;
-                const mpz_t*    m0     = reinterpret_cast<const mpz_t*>(Aiter+j+i*lda);
-                const uint16_t* m0_ptr = reinterpret_cast<const uint16_t*>(m0[0]->_mp_d);
-                size_t l=0;
-                //size_t maxs=std::min(k,(Aiter[j+i*lda].size())<<2);
-                size_t maxs=std::min(k,(Aiter[j+i*lda].size())*sizeof(mp_limb_t)/2); // to ensure 32 bits portability
-#ifdef __FFLASFFPACK_HAVE_LITTLE_ENDIAN
-                if (m0[0]->_mp_size >= 0)
-                    for (;l<maxs;l++)
-                        A_beta[l+idx*k]=  m0_ptr[l];
-                else
-                    for (;l<maxs;l++)
-                        A_beta[l+idx*k]= - double(m0_ptr[l]);
-#else
-                if (m0[0]->_mp_size >= 0)
-                    for (;l<maxs;l++) {
-                        size_t l2 = l ^ ((sizeof(mp_limb_t)/2U) - 1U);
-                        A_beta[l+idx*k]=  m0_ptr[l2];
-                    }
-                else
-                    for (;l<maxs;l++) {
-                        size_t l2 = l ^ ((sizeof(mp_limb_t)/2U) - 1U);
-                        A_beta[l+idx*k]= - double(m0_ptr[l2]);
-                    }
-#endif
-                for (;l<k;l++)
-                    A_beta[l+idx*k]=  0.;
-            }
-        }
-        if (RNS_MAJOR==false) {
-            // Arns = _crt_in x A_beta^T
-#ifndef ENABLE_CHECKER_fgemm
-            FFLAS::fgemm (Givaro::ZRing<double>(), FFLAS::FflasNoTrans,FFLAS::FflasTrans,_size,mn,k,1.0,_crt_in.data(),_ldm, A_beta, k, 0.,Arns,rda,
-                          FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive>());
-#else
-            cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasTrans,(int)_size,(int)mn,(int)k,1.0,_crt_in.data(),(int)_ldm,A_beta,(int)k,0.,Arns,(int)rda);
-#endif
-        }
-        else {
-            // Arns =  A_beta x _crt_in^T
-#ifndef ENABLE_CHECKER_fgemm
-            FFLAS::fgemm (Givaro::ZRing<double>(), FFLAS::FflasNoTrans,FFLAS::FflasTrans,mn,_size,k,1.0,A_beta, k, _crt_in.data(),_ldm,0.,Arns,_size,
-                          FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive>());
-#else
-            cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasTrans,(int)mn,(int)_size,(int)k,1.0,A_beta,(int)k,_crt_in.data(),(int)_ldm,0.,Arns,(int)_size);
-#endif
-        }
-        reduce(mn,Arns,rda,RNS_MAJOR);
-
-        FFLAS::fflas_delete( A_beta);
 
+  inline void rns_double::convert(size_t m, size_t n, integer gamma, integer* A, size_t lda,
+                                  const double* Arns, size_t rda, bool RNS_MAJOR, const FFLAS::FFLAS_TRANSPOSE trans) const
+  {
+    const size_t  mn= m*n;
+    Givaro::ZRing<double> ZD;
+    Givaro::ZRing<Givaro::Integer> ZZ;
+    if (mn) {
 #ifdef CHECK_RNS
-        bool ok=true;
-        for (size_t i=0;i<m;i++)
-            for(size_t j=0;j<n;j++)
-                for(size_t k=0;k<_size;k++)
-                    ok&= (((A[i*lda+j] % (int64_t) _basis[k])+(A[i*lda+j]<0?(int64_t)_basis[k]:0))
-                          == (int64_t) Arns[j*m+i+k*rda]);
-        std::cout<<"RNS freduce ... "<<(ok?"OK":"ERROR")<<std::endl;
-#endif
-        }
-    }
+      integer* Acopy=new integer[m*n];
+      FFLAS::fassign(ZZ,m,n,Acopy,n,A,lda);
+#endif 
 
-    inline void rns_double::convert(size_t m, size_t n, integer gamma, integer* A, size_t lda,
-                                    const double* Arns, size_t rda, bool RNS_MAJOR) const
-    {
-        const size_t  mn= m*n;
-        if (mn) {
-#ifdef CHECK_RNS
-        integer* Acopy=new integer[m*n];
-        for(size_t i=0;i<m;i++)
-            for(size_t j=0;j<n;j++)
-                Acopy[i*n+j]=A[i*lda+j];
+      double *A_beta= FFLAS::fflas_new<double>(mn*_ldm);
 
-#endif
+      // Using Helper for potential parallelism -> need to be activated by hand
+      //using ParallelStrategy= FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive>;
+      using ParallelStrategy= FFLAS::ParSeqHelper::Sequential;
+      FFLAS::MMHelper<Givaro::ZRing<double>, FFLAS::MMHelperAlgo::Winograd>  MatMulHelper (ZD, -1, ParallelStrategy());
 
-        integer hM= (_M-1)>>1;
-        double *A_beta= FFLAS::fflas_new<double>(mn*_ldm);
-        Givaro::Timer tfgemmc;tfgemmc.start();
-        if (RNS_MAJOR==false) {// compute A_beta = Ap^T x M_beta
-#ifndef ENABLE_CHECKER_fgemm
-            FFLAS::fgemm(Givaro::ZRing<double>(),FFLAS::FflasTrans, FFLAS::FflasNoTrans, mn, _ldm, _size, 1.0 , Arns, rda, _crt_out.data(), _ldm, 0., A_beta,_ldm, FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive >());
-#else
-            cblas_dgemm(CblasRowMajor,CblasTrans, CblasNoTrans, (int)mn, (int)_ldm, (int)_size, 1.0 , Arns, (int)rda, _crt_out.data(), (int)_ldm, 0., A_beta,(int)_ldm);
-#endif
-        }
-        else  {// compute A_beta = Ap x M_Beta
-#ifndef ENABLE_CHECKER_fgemm
-            FFLAS::fgemm(Givaro::ZRing<double>(),FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, mn, _ldm, _size, 1.0 , Arns, _size, _crt_out.data(), _ldm, 0., A_beta, _ldm, FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive >());
-#else
-            cblas_dgemm(CblasRowMajor,CblasNoTrans, CblasNoTrans, (int)mn, (int)_ldm, (int)_size, 1.0 , Arns, (int)_size, _crt_out.data(), (int)_ldm, 0., A_beta,(int)_ldm);
-#endif
-        }
-        tfgemmc.stop();
-        //if(m>1 && n>1) std::cerr<<"fgemm Convert : "<<tfgemmc.realtime()<<std::endl;
-        // compute A using inverse Kronecker transform of A_beta expressed in base 2^log_beta
-        integer* Aiter= A;
-        size_t k=_ldm;
-        size_t k4=((k+3)>>2)+ (((k+3)%4==0)?0:1);
-        std::vector<uint16_t> A0(k4<<2,0),A1(k4<<2,0),A2(k4<<2,0),A3(k4<<2,0);
-        integer a0,a1,a2,a3,res;
-        mpz_t *m0,*m1,*m2,*m3;
-        m0= reinterpret_cast<mpz_t*>(&a0);
-        m1= reinterpret_cast<mpz_t*>(&a1);
-        m2= reinterpret_cast<mpz_t*>(&a2);
-        m3= reinterpret_cast<mpz_t*>(&a3);
-        mp_limb_t *m0_d,*m1_d,*m2_d,*m3_d;
-        m0_d = m0[0]->_mp_d;
-        m1_d = m1[0]->_mp_d;
-        m2_d = m2[0]->_mp_d;
-        m3_d = m3[0]->_mp_d;
-        m0[0]->_mp_alloc = m1[0]->_mp_alloc = m2[0]->_mp_alloc = m3[0]->_mp_alloc = (int) (k4*8/sizeof(mp_limb_t)); // to ensure 32 bits portability
-        m0[0]->_mp_size  = m1[0]->_mp_size  = m2[0]->_mp_size  = m3[0]->_mp_size  = (int) (k4*8/sizeof(mp_limb_t)); // to ensure 32 bits portability
-        Givaro::Timer tkroc;
-        tkroc.start();
-        //		auto sp=SPLITTER();
-        //		PARFOR1D(i,m,sp,
-        for(size_t i=0;i<m;i++)
-            for (size_t j=0;j<n;j++){
-                size_t idx=i*n+j;
-                for (size_t l=0;l<k;l++){
-                    uint64_t tmp=(uint64_t)A_beta[l+idx*k];
-                    uint16_t* tptr= reinterpret_cast<uint16_t*>(&tmp);
-#ifdef __FFLASFFPACK_HAVE_LITTLE_ENDIAN
-                    A0[l  ]= tptr[0];
-                    A1[l+1]= tptr[1];
-                    A2[l+2]= tptr[2];
-                    A3[l+3]= tptr[3];
-#else
-                    A0[l     ^ ((sizeof(mp_limb_t)/2U) - 1U)] = tptr[3];
-                    A1[(l+1) ^ ((sizeof(mp_limb_t)/2U) - 1U)] = tptr[2];
-                    A2[(l+2) ^ ((sizeof(mp_limb_t)/2U) - 1U)] = tptr[1];
-                    A3[(l+3) ^ ((sizeof(mp_limb_t)/2U) - 1U)] = tptr[0];
-#endif
-                }
-                // see A0,A1,A2,A3 as a the gmp integers a0,a1,a2,a3
-                m0[0]->_mp_d= reinterpret_cast<mp_limb_t*>(&A0[0]);
-                m1[0]->_mp_d= reinterpret_cast<mp_limb_t*>(&A1[0]);
-                m2[0]->_mp_d= reinterpret_cast<mp_limb_t*>(&A2[0]);
-                m3[0]->_mp_d= reinterpret_cast<mp_limb_t*>(&A3[0]);
-                res = a0;res+= a1;res+= a2;res+= a3;
-                res%=_M;
+      if (RNS_MAJOR==false) // compute A_beta = Ap^T x M_beta
+        FFLAS::fgemm(Givaro::ZRing<double>(),FFLAS::FflasTrans, FFLAS::FflasNoTrans, mn, _ldm, _size, 1.0 , Arns, rda, _crt_out.data(), _ldm, 0., A_beta,_ldm, MatMulHelper);
+      else  // compute A_beta = Ap x M_Beta
+        FFLAS::fgemm(Givaro::ZRing<double>(),FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, mn, _ldm, _size, 1.0 , Arns, _size, _crt_out.data(), _ldm, 0., A_beta, _ldm, MatMulHelper);
 
-                // get the correct result according to the expected sign of A
-                if (res>hM)
-                    res-=_M;
-                if (gamma==0)
-                    Aiter[j+i*lda]=res;
-                else
-                    if (gamma==integer(1))
-                        Aiter[j+i*lda]+=res;
-                    else
-                        if (gamma==integer(-1))
-                            Aiter[j+i*lda]=res-Aiter[j+i*lda];
-                        else{
-                            Aiter[j+i*lda]*=gamma;
-                            Aiter[j+i*lda]+=res;
-                        }
 
-            }
-        tkroc.stop();
-        //if(m>1 && n>1) std::cerr<<"Kronecker Convert : "<<tkroc.realtime()<<std::endl;
+      RNS_COMMON::inverse_Kronecker_base16(trans,m,n,A_beta,A,lda,_ldm, _M, gamma);
 
-        m0[0]->_mp_d = m0_d;
-        m1[0]->_mp_d = m1_d;
-        m2[0]->_mp_d = m2_d;
-        m3[0]->_mp_d = m3_d;
-        m0[0]->_mp_alloc = m1[0]->_mp_alloc = m2[0]->_mp_alloc= m3[0]->_mp_alloc = 1;
-        m0[0]->_mp_size  = m1[0]->_mp_size  = m2[0]->_mp_size = m3[0]->_mp_size  = 0;
-        FFLAS::fflas_delete( A_beta);
+      //FFLAS::WriteMatrix(std::cout,Givaro::ZRing<Givaro::Integer>(),m,n, A, lda, FFLAS::FflasSageMath)<<std::endl;
+      //FFLAS::WriteMatrix(std::cout,Givaro::ZRing<double>(),1,_size, _basis.data(), _size, FFLAS::FflasSageMath)<<std::endl;      
+      // if (RNS_MAJOR)
+      //   FFLAS::WriteMatrix(std::cout,Givaro::ZRing<double>(),mn,_size, Arns, rda,FFLAS::FflasSageMath)<<std::endl;
+      // else
+      //   FFLAS::WriteMatrix(std::cout,Givaro::ZRing<double>(),_size, mn, Arns, rda,FFLAS::FflasSageMath)<<std::endl;
 
+      FFLAS::fflas_delete(A_beta);
 #ifdef CHECK_RNS
-        bool ok=true;
-        for (size_t i=0;i<m;i++)
-            for(size_t j=0;j<n;j++)
-                for(size_t k=0;k<_size;k++){
-                    int64_t _p =(int64_t) _basis[k];
-                    integer curr=A[i*lda+j] - gamma*Acopy[i*n+j];
-                    ok&= ( curr% _p +(curr%_p<0?_p:0) == (int64_t) Arns[i*n+j+k*rda]);
-                    //std::cout<<A[i*lda+j]<<" mod "<<(int64_t) _basis[k]<<"="<<(int64_t) Arns[i*n+j+k*rda]<<";"<<std::endl;
-                }
-        std::cout<<"RNS convert ... "<<(ok?"OK":"ERROR")<<std::endl;
+      RNS_COMMON::check_rns_to_integers(trans, m, n, Arns, rda, A, lda, gamma, Acopy, n, _basis, RNS_MAJOR);
+      FFLAS::fflas_delete(Acopy);
 #endif
-        }
-    }
 
-    inline void rns_double::convert_transpose(size_t m, size_t n, integer gamma, integer* A, size_t lda,
-                                              const double* Arns, size_t rda, bool RNS_MAJOR) const
-    {
-        const size_t  mn= m*n;
-        if (mn) {
-        integer hM= (_M-1)>>1;
-        double *A_beta= FFLAS::fflas_new<double>(mn*_ldm);
-
-        if (RNS_MAJOR==false){
-            // compute A_beta = Ap^T x M_beta
-#ifndef ENABLE_CHECKER_fgemm
-            FFLAS::fgemm(Givaro::ZRing<double>(),FFLAS::FflasTrans, FFLAS::FflasNoTrans, mn, _ldm, _size, 1.0 , Arns, rda, _crt_out.data(), _ldm, 0., A_beta, _ldm, FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive >());
-#else
-            cblas_dgemm(CblasRowMajor,CblasTrans, CblasNoTrans,(int) mn,(int) _ldm,(int) _size, 1.0 , Arns,(int) rda, _crt_out.data(),(int) _ldm, 0., A_beta,(int)_ldm);
-#endif
-        }
-        else { // compute A_beta = Ap x M_Beta
-#ifndef ENABLE_CHECKER_fgemm
-            FFLAS::fgemm(Givaro::ZRing<double>(),FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, mn, _ldm, _size, 1.0 , Arns, rda, _crt_out.data(), _ldm, 0., A_beta, _ldm, FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::TwoDAdaptive >());
-#else
-            cblas_dgemm(CblasRowMajor,CblasNoTrans, CblasNoTrans, (int)mn, (int)_ldm, (int)_size, 1.0 , Arns, (int)rda, _crt_out.data(), (int)_ldm, 0., A_beta,(int)_ldm);
-#endif
-        }
-        // compute A using inverse Kronecker transform of A_beta expressed in base 2^log_beta
-        integer* Aiter= A;
-        size_t k=_ldm;
-        size_t k4=((k+3)>>2)+ (((k+3)%4==0)?0:1);
-        std::vector<uint16_t> A0(k4<<2,0),A1(k4<<2,0),A2(k4<<2,0),A3(k4<<2,0);
-        integer a0,a1,a2,a3,res;
-        mpz_t *m0,*m1,*m2,*m3;
-        m0= reinterpret_cast<mpz_t*>(&a0);
-        m1= reinterpret_cast<mpz_t*>(&a1);
-        m2= reinterpret_cast<mpz_t*>(&a2);
-        m3= reinterpret_cast<mpz_t*>(&a3);
-        mp_limb_t *m0_d,*m1_d,*m2_d,*m3_d;
-        m0_d = m0[0]->_mp_d;
-        m1_d = m1[0]->_mp_d;
-        m2_d = m2[0]->_mp_d;
-        m3_d = m3[0]->_mp_d;
-        m0[0]->_mp_alloc = m1[0]->_mp_alloc = m2[0]->_mp_alloc = m3[0]->_mp_alloc = (int32_t)(k4*8/sizeof(mp_limb_t)); // to ensure 32 bits portability
-        m0[0]->_mp_size  = m1[0]->_mp_size  = m2[0]->_mp_size  = m3[0]->_mp_size  = (int32_t)(k4*8/sizeof(mp_limb_t)); // to ensure 32 bits portability
-        for (size_t j=0;j<n;j++)
-            for(size_t i=0;i<m;i++){
+}
+  }
 
 
-                size_t idx=i+j*m;
-                for (size_t l=0;l<k;l++){
-                    uint64_t tmp=(uint64_t)A_beta[l+idx*k];
-                    uint16_t* tptr= reinterpret_cast<uint16_t*>(&tmp);
-#ifdef __FFLASFFPACK_HAVE_LITTLE_ENDIAN
-                    A0[l  ]= tptr[0];
-                    A1[l+1]= tptr[1];
-                    A2[l+2]= tptr[2];
-                    A3[l+3]= tptr[3];
-#else
-                    A0[l     ^ ((sizeof(mp_limb_t)/2U) - 1U)] = tptr[3];
-                    A1[(l+1) ^ ((sizeof(mp_limb_t)/2U) - 1U)] = tptr[2];
-                    A2[(l+2) ^ ((sizeof(mp_limb_t)/2U) - 1U)] = tptr[1];
-                    A3[(l+3) ^ ((sizeof(mp_limb_t)/2U) - 1U)] = tptr[0];
-#endif
-                }
-                // see A0,A1,A2,A3 as a the gmp integers a0,a1,a2,a3
-                m0[0]->_mp_d= reinterpret_cast<mp_limb_t*>(&A0[0]);
-                m1[0]->_mp_d= reinterpret_cast<mp_limb_t*>(&A1[0]);
-                m2[0]->_mp_d= reinterpret_cast<mp_limb_t*>(&A2[0]);
-                m3[0]->_mp_d= reinterpret_cast<mp_limb_t*>(&A3[0]);
-                res = a0;res+= a1;res+= a2;res+= a3;
-                res%=_M;
 
-                // get the correct result according to the expected sign of A
-                if (res>hM)
-                    res-=_M;
-                if (gamma==0)
-                    Aiter[j+i*lda]=res;
-                else
-                    if (gamma==integer(1))
-                        Aiter[j+i*lda]+=res;
-                    else
-                        if (gamma==integer(-1))
-                            Aiter[j+i*lda]=res-Aiter[j+i*lda];
-                        else{
-                            Aiter[j+i*lda]*=gamma;
-                            Aiter[j+i*lda]+=res;
-                        }
+  // reduce entries of Arns to be less than the rns basis elements
+  inline void rns_double::reduce(size_t n, double* Arns, size_t rda, bool RNS_MAJOR) const{
 
-            }
-        m0[0]->_mp_d = m0_d;
-        m1[0]->_mp_d = m1_d;
-        m2[0]->_mp_d = m2_d;
-        m3[0]->_mp_d = m3_d;
-        m0[0]->_mp_alloc = m1[0]->_mp_alloc = m2[0]->_mp_alloc= m3[0]->_mp_alloc = 1;
-        m0[0]->_mp_size  = m1[0]->_mp_size  = m2[0]->_mp_size = m3[0]->_mp_size  = 0;
-        FFLAS::fflas_delete( A_beta);
-#ifdef CHECK_RNS
-        bool ok=true;
-        for (size_t i=0;i<m;i++)
-            for(size_t j=0;j<n;j++)
-                for(size_t k=0;k<_size;k++){
-                    ok&= (((A[i*lda+j] % (int64_t) _basis[k])+(A[i*lda+j]% (int64_t) _basis[k]<0?(int64_t)_basis[k]:0)) == (int64_t) Arns[i+j*m+k*rda]);
-                    //std::cout<<A[i*lda+j]<<" mod "<<(int64_t) _basis[k]<<"="<<(int64_t) Arns[i*n+j+k*rda]<<";"<<std::endl;
-                }
-        std::cout<<"RNS convert ... "<<(ok?"OK":"ERROR")<<std::endl;
-#endif // CHECK_RNS
-        }
-    }
-
-    // reduce entries of Arns to be less than the rns basis elements
-    inline void rns_double::reduce(size_t n, double* Arns, size_t rda, bool RNS_MAJOR) const{
-
-        if (RNS_MAJOR) {
+    if (RNS_MAJOR) {
 #ifdef __FFLASFFPACK_HAVE_SSE4_1_INSTRUCTIONS
-            using simd = Simd<double>;
-            using vect_t = typename simd::vect_t;
-
-            if(_size % simd::vect_size == 0){
-                for(size_t i = 0 ; i < n ; i++){
-                    vect_t tmp1, tmp2, tmp3, v, max, basis, inv, neg;
-                    for(size_t j = 0 ; j < _size ; j+=simd::vect_size){
-                        basis = simd::load(_basis.data()+j);
-                        inv   = simd::load(_invbasis.data()+j);
-                        max   = simd::load(_basisMax.data()+j);
-                        neg   = simd::load(_negbasis.data()+j);
-                        v     = simd::load(Arns+i*_size+j);
-                        tmp1  = simd::floor(simd::mul(v, inv));
-                        tmp2  = simd::fnmadd(v, tmp1, basis);
-                        tmp1  = simd::greater(tmp2, max);
-                        tmp3  = simd::lesser(tmp2, simd::zero());
-                        tmp1  = simd::vand(tmp1, neg);
-                        tmp3  = simd::vand(tmp3, basis);
-                        tmp1  = simd::vor(tmp1, tmp3);
-                        tmp2  = simd::add(tmp2, tmp1);
-                        simd::store(Arns+i*_size+j, tmp2);
-                    }
-                }
-            }else{
-                for(size_t i = 0 ; i < n ; i++){
-                    vect_t tmp1, tmp2, tmp3, v, max, basis, inv, neg;
-                    size_t j = 0;
-                    for( ; j < ROUND_DOWN(_size, simd::vect_size) ; j+=simd::vect_size){
-                        basis = simd::load(_basis.data()+j);
-                        inv   = simd::load(_invbasis.data()+j);
-                        max   = simd::load(_basisMax.data()+j);
-                        neg   = simd::load(_negbasis.data()+j);
-                        v     = simd::loadu(Arns+i*_size+j);
-                        tmp1  = simd::floor(simd::mul(v, inv));
-                        tmp2  = simd::fnmadd(v, tmp1, basis);
-                        tmp1  = simd::greater(tmp2, max);
-                        tmp3  = simd::lesser(tmp2, simd::zero());
-                        tmp1  = simd::vand(tmp1, neg);
-                        tmp3  = simd::vand(tmp3, basis);
-                        tmp1  = simd::vor(tmp1, tmp3);
-                        tmp2  = simd::add(tmp2, tmp1);
-                        simd::storeu(Arns+i*_size+j, tmp2);
-                    }
-                    for( ; j < _size ; ++j){
-                        // std::cout << j << std::endl;
-                        // auto x = std::floor(Arns[i*_size+j] * _invbasis[j]);
-                        Arns[i*_size+j] -= std::floor(Arns[i*_size+j]*_invbasis[j])*_basis[j];
-                        // Arns[i*_size+j] = std::fma(Arns[i*_size+j], -x, _basis[j]);
-                        if(Arns[i*_size+j] >= _basis[j]){
-                            Arns[i*_size+j] -= _basis[j];
-                        }else if(Arns[i*_size+j] < 0){
-                            Arns[i*_size+j] += _basis[j];
-                        }
-                    }
-                }
-            }
-#else
-            for(size_t i = 0 ; i < n ; i+= _size){
-                for(size_t j = 0 ; j < _size ; ++j){
-                    //_field_rns.reduce(Arns+i*_size+j);
-                    _field_rns[i].reduce(Arns[i*_size+j]);
-                }
+      using simd = Simd<double>;
+      using vect_t = typename simd::vect_t;
+
+      if(_size % simd::vect_size == 0){
+        for(size_t i = 0 ; i < n ; i++){
+          vect_t tmp1, tmp2, tmp3, v, max, basis, inv, neg;
+          for(size_t j = 0 ; j < _size ; j+=simd::vect_size){
+            basis = simd::load(_basis.data()+j);
+            inv   = simd::load(_invbasis.data()+j);
+            max   = simd::load(_basisMax.data()+j);
+            neg   = simd::load(_negbasis.data()+j);
+            v     = simd::load(Arns+i*_size+j);
+            tmp1  = simd::floor(simd::mul(v, inv));
+            tmp2  = simd::fnmadd(v, tmp1, basis);
+            tmp1  = simd::greater(tmp2, max); 
+            tmp3  = simd::lesser(tmp2, simd::zero());
+            tmp1  = simd::vand(tmp1, neg);
+            tmp3  = simd::vand(tmp3, basis);
+            tmp1  = simd::vor(tmp1, tmp3);
+            tmp2  = simd::add(tmp2, tmp1);
+            simd::store(Arns+i*_size+j, tmp2);
+          }
+        }
+      }else{
+        for(size_t i = 0 ; i < n ; i++){
+          vect_t tmp1, tmp2, tmp3, v, max, basis, inv, neg;
+          size_t j = 0;
+          for( ; j < ROUND_DOWN(_size, simd::vect_size) ; j+=simd::vect_size){
+            basis = simd::load(_basis.data()+j);
+            inv   = simd::load(_invbasis.data()+j);
+            max   = simd::load(_basisMax.data()+j);
+            neg   = simd::load(_negbasis.data()+j);
+            v     = simd::loadu(Arns+i*_size+j);
+            tmp1  = simd::floor(simd::mul(v, inv));
+            tmp2  = simd::fnmadd(v, tmp1, basis);
+            tmp1  = simd::greater(tmp2, max);
+            tmp3  = simd::lesser(tmp2, simd::zero());
+            tmp1  = simd::vand(tmp1, neg);
+            tmp3  = simd::vand(tmp3, basis);
+            tmp1  = simd::vor(tmp1, tmp3);
+            tmp2  = simd::add(tmp2, tmp1);
+            simd::storeu(Arns+i*_size+j, tmp2);
+          }
+          for( ; j < _size ; ++j){
+            // std::cout << j << std::endl;
+            // auto x = std::floor(Arns[i*_size+j] * _invbasis[j]);
+            Arns[i*_size+j] -= std::floor(Arns[i*_size+j]*_invbasis[j])*_basis[j];
+            // Arns[i*_size+j] = std::fma(Arns[i*_size+j], -x, _basis[j]);
+            if(Arns[i*_size+j] >= _basis[j]){
+              Arns[i*_size+j] -= _basis[j];
+            }else if(Arns[i*_size+j] < 0){
+              Arns[i*_size+j] += _basis[j];
             }
-#endif
-        }
-        else { // NOT IN RNS MAJOR
-            // #ifndef __FFLASFFPACK_SEQUENTIAL
-            // 			auto sp=SPLITTER(MAX_THREADS);
-            // #else
-            // 			auto sp=SPLITTER(1);
-            // #endif
-            PARFOR1D(i,_size,SPLITTER(NUM_THREADS),
-                     //for(size_t i=0;i<_size;i++)
-                     FFLAS::freduce (_field_rns[i],n,Arns+i*rda,1);
-                    );
+          }
         }
-
-    }
-
-
-    // TODO: less naive implementation
-    inline void rns_double_extended::init(size_t m, double* Arns, const integer* A, size_t lda) const{
-        for(size_t i = 0 ; i < m ; ++i){
-            for(size_t j = 0 ; j < _size ; ++j){
-                Arns[i*_size+j] = (double)((A[i*lda]%integer(_basis[j]))[0]);
-            }
+      }
+#else
+      for(size_t i = 0 ; i < n ; i++){
+        for(size_t j = 0 ; j < _size ; ++j){
+          //_field_rns.reduce(Arns+i*_size+j);
+          _field_rns[j].reduce(Arns[i*_size+j]);
         }
+      }
+#endif
     }
-
-    // TODO: less naive implementation
-    inline void rns_double_extended::convert(size_t m, integer *A, const double *Arns) const{
-        integer hM= (_M-1)/2;
-        for(size_t i = 0 ; i < m ; ++i){
-            A[i] = 0;
-            integer tmp;
-            for(size_t j = 0 ; j < _size ; ++j){
-                A[i] += ((integer(Arns[i*_size+j])*integer(_MMi[j]))%integer(_basis[j]))*integer(_Mi[j]);
-            }
-            A[i] %= _M;
-            if(A[i] > hM)
-                A[i] -= _M;
-        }
+    else { // NOT IN RNS MAJOR
+      // #ifndef __FFLASFFPACK_SEQUENTIAL
+      // 			auto sp=SPLITTER(MAX_THREADS);
+      // #else
+      // 			auto sp=SPLITTER(1);
+      // #endif
+      PARFOR1D(i,_size,SPLITTER(NUM_THREADS),
+               //for(size_t i=0;i<_size;i++)
+               FFLAS::freduce (_field_rns[i],n,Arns+i*rda,1);
+               );
     }
 
-    // reduce entries of Arns to be less than the rns basis elements
-    inline void rns_double_extended::reduce(size_t n, double* Arns, size_t rda, bool RNS_MAJOR) const{
-
-#ifdef __FFLASFFPACK_HAVE_SSE4_1_INSTRUCTIONS
-        using simd = Simd<double>;
-        using vect_t = typename simd::vect_t;
+  }
 
-        if(_size % simd::vect_size == 0){
-            //#pragma omp parallel for schedule(static, 256)
-            for(size_t i = 0 ; i < n ; i++){
-                vect_t tmp1, tmp2, tmp3, v, max, basis, inv, neg;
-                for(size_t j = 0 ; j < _size ; j+=simd::vect_size){
-                    basis = simd::load(_basis.data()+j);
-                    inv   = simd::load(_invbasis.data()+j);
-                    max   = simd::load(_basisMax.data()+j);
-                    neg   = simd::load(_negbasis.data()+j);
-                    v     = simd::load(Arns+i*_size+j);
-                    tmp2 = modSimd(v, basis, inv, neg);
-                    tmp1  = simd::greater(tmp2, max);
-                    tmp3  = simd::lesser(tmp2, simd::zero());
-                    tmp1  = simd::vand(tmp1, neg);
-                    tmp3  = simd::vand(tmp3, basis);
-                    tmp1  = simd::vor(tmp1, tmp3);
-                    tmp2  = simd::add(tmp2, tmp1);
-                    simd::store(Arns+i*_size+j, tmp2);
-                }
-            }
-        }else{
-            //#pragma omp parallel for schedule(static, 256)
-            for(size_t i = 0 ; i < n ; i++){
-                vect_t tmp1, tmp2, tmp3, v, max, basis, inv, neg;
-                size_t j = 0;
-                for( ; j < ROUND_DOWN(_size, simd::vect_size) ; j+=simd::vect_size){
-                    basis = simd::load(_basis.data()+j);
-                    inv   = simd::load(_invbasis.data()+j);
-                    max   = simd::load(_basisMax.data()+j);
-                    neg   = simd::load(_negbasis.data()+j);
-                    v     = simd::loadu(Arns+i*_size+j);
-                    tmp2 = modSimd(v, basis, inv, neg);
-                    tmp1  = simd::greater(tmp2, max);
-                    tmp3  = simd::lesser(tmp2, simd::zero());
-                    tmp1  = simd::vand(tmp1, neg);
-                    tmp3  = simd::vand(tmp3, basis);
-                    tmp1  = simd::vor(tmp1, tmp3);
-                    tmp2  = simd::add(tmp2, tmp1);
-                    simd::storeu(Arns+i*_size+j, tmp2);
-                }
-                for( ; j < _size ; ++j){
-                    _field_rns[j].reduce(Arns[i*_size+j]);
-                }
-            }
-        }
-#else
-
-        // TODO : SIMD version
-        for(size_t i = 0 ; i < n ; i+= _size){
-            for(size_t j = 0 ; j < _size ; ++j){
-                //_field_rns.reduce(Arns+i*_size+j);
-                _field_rns[i].reduce(Arns[i*_size+j]);
-            }
-        }
-
-#endif
-
-    }
 
-    } // FFPACK
+} // FFPACK
 
 #endif // __FFLASFFPACK_field_rns_double_INL
-    /* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
-    // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
+ 
diff --git a/fflas-ffpack/field/rns-integer.h b/fflas-ffpack/field/rns-integer.h
index 27782c83f..2ed6bdbac 100644
--- a/fflas-ffpack/field/rns-integer.h
+++ b/fflas-ffpack/field/rns-integer.h
@@ -34,6 +34,7 @@
 #include <givaro/givinteger.h>
 
 #include "fflas-ffpack/field/rns-double.h"
+#include "fflas-ffpack/field/rns-double-extended.h"
 
 namespace FFPACK {
 
@@ -111,7 +112,7 @@ namespace FFPACK {
         }
         Element& init(Element& x, const Givaro::Integer& y) const{
             init(x);
-            size_t k =(y.bitsize())/16+((y.bitsize())%16?1:0);
+            uint64_t k =(y.bitsize())/16+((y.bitsize())%16?1:0);
             _rns->init(1,1,x._ptr,x._stride, &y,1,k);
             return x;
         }
@@ -153,19 +154,44 @@ namespace FFPACK {
 
 namespace FFLAS {
 
-    // specialization for the fflas alloc function
-    template<>
-    inline FFPACK::rns_double_elt_ptr
-    fflas_new(const FFPACK::RNSInteger<FFPACK::rns_double> &F, const size_t m, const Alignment align){
-        double *ptr=FFLAS::fflas_new<double>(m*F.size(), align);
-        return FFPACK::rns_double_elt_ptr(ptr,m);
-    }
-
-    template<>
-    inline FFPACK::rns_double_elt_ptr
-    fflas_new(const FFPACK::RNSInteger<FFPACK::rns_double> &F, const size_t m, const size_t n,  const Alignment align){
-        return fflas_new(F, m*n, align);
-    }
+  // specialization for the fflas alloc function
+  template<>
+  inline FFPACK::rns_double_elt_ptr
+  fflas_new(const FFPACK::RNSInteger<FFPACK::rns_double> &F, const size_t m, const Alignment align){
+    double *ptr=FFLAS::fflas_new<double>(m*F.size(), align);
+    return FFPACK::rns_double_elt_ptr(ptr,m);
+  }
+
+  template<>
+  inline FFPACK::rns_double_elt_ptr
+  fflas_new(const FFPACK::RNSInteger<FFPACK::rns_double> &F, const size_t m, const size_t n,  const Alignment align){
+    return fflas_new(F, m*n, align);
+  }
+  // specialization for the fflas alloc function
+  template<>
+  inline FFPACK::rns_double_elt_ptr
+  fflas_new(const FFPACK::RNSInteger<FFPACK::rns_double_extended_V1> &F, const size_t m, const Alignment align){
+    double *ptr=FFLAS::fflas_new<double>(m*F.size(), align);
+    return FFPACK::rns_double_elt_ptr(ptr,m);
+  }
+
+  template<>
+  inline FFPACK::rns_double_elt_ptr
+  fflas_new(const FFPACK::RNSInteger<FFPACK::rns_double_extended_V1> &F, const size_t m, const size_t n,  const Alignment align){
+    return fflas_new(F, m*n, align);
+  }
+  template<>
+  inline FFPACK::rns_double_elt_ptr
+  fflas_new(const FFPACK::RNSInteger<FFPACK::rns_double_extended_V2> &F, const size_t m, const Alignment align){
+    double *ptr=FFLAS::fflas_new<double>(m*F.size(), align);
+    return FFPACK::rns_double_elt_ptr(ptr,m);
+  }
+
+  template<>
+  inline FFPACK::rns_double_elt_ptr
+  fflas_new(const FFPACK::RNSInteger<FFPACK::rns_double_extended_V2> &F, const size_t m, const size_t n,  const Alignment align){
+    return fflas_new(F, m*n, align);
+  }
 
     // function to convert from integer to RNS (note: this is not the finit function from FFLAS, extra k)
     template<typename RNS>
diff --git a/tests/Makefile.am b/tests/Makefile.am
index c3c4338e7..ca18a1937 100755
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -77,6 +77,7 @@ BASIC_TESTS =               \
 		test-fgemv \
 		test-nullspace \
 		test-storage-transpose	\
+	  test-rns-double	\
 		regression-check
 
 if FFLASFFPACK_PRECOMPILED
@@ -162,6 +163,7 @@ test_interfaces_c_SOURCES = test-interfaces-c.c
 test_maxdelayeddim_SOURCES = test-maxdelayeddim.C
 #  test_fspmv_SOURCES = test-fspmv.C
 test_storage_transpose_SOURCES = test-storage-transpose.C
+test_rns_double_SOURCES = test-rns-double.C
 regression_check_SOURCES = regression-check.C
 
 test_solve_SOURCES = test-solve.C
diff --git a/tests/test-rns-double.C b/tests/test-rns-double.C
new file mode 100644
index 000000000..507345439
--- /dev/null
+++ b/tests/test-rns-double.C
@@ -0,0 +1,236 @@
+/*
+ * Copyright (C) 2023 the FFLAS-FFPACK group
+ *
+ * Written by Pascal Giorgi <pascal.giorgi@lirmm.fr>
+ *
+ * This file is Free Software and part of FFLAS-FFPACK.
+ *
+ * ========LICENCE========
+ * This file is part of the library FFLAS-FFPACK.
+ *
+ * FFLAS-FFPACK is free software: you can redistribute it and/or modify
+ * it under the terms of the  GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * This library 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
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+ * ========LICENCE========
+ *.
+ */
+
+//#define CHECK_RNS
+//#define RNS_DEBUG
+//#define BENCH_RNS
+ 
+#define __FFLASFFPACK_SEQUENTIAL
+
+
+#include "fflas-ffpack/fflas-ffpack-config.h"
+
+
+#include <iostream>
+#include <vector>
+#include <string>
+#include <chrono>
+#include <algorithm>
+using namespace std; 
+
+#include "fflas-ffpack/utils/timer.h"
+#include "fflas-ffpack/utils/fflas_io.h"
+#include "fflas-ffpack/fflas/fflas.h"
+#include "fflas-ffpack/field/rns-double.h"
+#include "fflas-ffpack/field/rns-double-extended.h"
+#include "fflas-ffpack/field/rns-integer.h"
+#include "fflas-ffpack/utils/args-parser.h"
+#include "givaro/givinteger.h" 
+#include "givaro/modular-integer.h"
+#include "givaro/zring.h"
+#include <givaro/givrns.h>
+
+template<typename RNS_t>
+void print_rns();
+
+template<>
+void print_rns<FFPACK::rns_double>(){ std::cout<<std::setw(22)<<"(rns_double)";}
+template<>
+void print_rns<FFPACK::rns_double_extended>(){ std::cout<<std::setw(22)<<"(rns_double_extended)";}
+
+
+
+template <typename RNS_t>
+bool run_check(size_t n, size_t primes_bits, size_t b, size_t iters, size_t seed){
+	Givaro::Integer p;
+
+	size_t bits=b;
+  bool ok=true;
+
+  std::cout<<"Checking "; print_rns<RNS_t>();
+  std::cout<<" int bits="<<std::left<<std::setw(10)<<bits<<" mi bits="<<std::left<<std::setw(10)<<primes_bits<<std::endl;
+  Givaro::ZRing<Givaro::Integer> ZZ;
+  Givaro::ZRing<double> ZD;
+
+#define SPC  20
+
+	for (size_t loop=0;(loop<iters);loop++){
+    using Field =  Givaro::Modular<Givaro::Integer>;
+    Givaro::Integer::seeding(seed);
+		Givaro::Integer::random_exact_2exp(p, bits);
+
+    Field F(p);
+		Field::RandIter Rand(F,seed);
+		Field::Element_ptr A,At,B;
+		A     = FFLAS::fflas_new(F,n,n);
+    At    = FFLAS::fflas_new(F,n,n);
+    B     = FFLAS::fflas_new(F,n,n);
+
+		FFLAS::frand(F,Rand,n,n,A,n);
+    Givaro::Integer shift= (p-1)/2;
+    for(size_t i=0;i<n*n;i++)
+      A[i]-=shift;
+
+
+    FFLAS::ftranspose(F,n,n,A,n,At,n);
+    FFLAS::fassign(F,n,n,A,n,B,n);
+
+
+
+    /* FFLAS::WriteMatrix(std::cout,F,n,n,A, n); */
+    /* FFLAS::WriteMatrix(std::cout,F,n,n,At, n); */
+
+    Givaro::Integer rns_bound = 3*p;
+
+    // construct an RNS structure and its associated Domain
+    RNS_t RNS(rns_bound, primes_bits,false,seed);
+    typedef FFPACK::RNSInteger<RNS_t> RnsDomain;
+    RnsDomain Zrns(RNS);
+    typename RnsDomain::Element_ptr mod_A  = FFLAS::fflas_new(Zrns,n,n);
+    typename RnsDomain::Element_ptr mod_At = FFLAS::fflas_new(Zrns,n,n);
+    size_t n2=n*n, s=RNS._size;
+
+    bool RNSMAJOR;
+    /**********************
+     ** check RNSMajor=false
+     **********************/
+#if 1
+    RNSMAJOR=false;
+    
+    RNS.init(n,n,mod_A._ptr,mod_A._stride,A,n,p, RNSMAJOR);
+    ok = ok and  FFPACK::RNS_COMMON::check_integers_to_rns(FFLAS::FflasNoTrans, n, n, mod_A._ptr, mod_A._stride, A, n, RNS._basis, RNSMAJOR, false);
+
+    RNS.init_transpose(n,n,mod_At._ptr,mod_At._stride,At,n,p,RNSMAJOR);
+    ok = ok and  FFPACK::RNS_COMMON::check_integers_to_rns(FFLAS::FflasTrans, n, n, mod_At._ptr, mod_At._stride, At, n, RNS._basis, RNSMAJOR, false);
+
+    ok = ok and (FFLAS::fequal(ZD,s,n2,mod_A._ptr,n2,mod_At._ptr,n2));
+
+      /* std::cout<<std::boolalpha<<ok<<std::endl; */
+      /* FFLAS::WriteMatrix(std::cout,ZD,n2, s, mod_A._ptr, n)<<std::endl;//FFLAS::FflasSageMath)<<std::endl; */
+      /* FFLAS::WriteMatrix(std::cout,ZD,n2, s, mod_At._ptr, n)<<std::endl;//FFLAS::FflasSageMath)<<std::endl; */
+
+    RNS.convert(n,n,Givaro::Integer(0), B, n, mod_A._ptr,mod_A._stride,RNSMAJOR);
+    ok= ok and (FFLAS::fequal(ZZ,n,n,A,n,B,n));
+
+      /* FFLAS::WriteMatrix(std::cout,ZZ,n, n, A, n)<<std::endl;//FFLAS::FflasSageMath)<<std::endl; */
+      /* FFLAS::WriteMatrix(std::cout,ZZ,n, n, B, n)<<std::endl;//FFLAS::FflasSageMath)<<std::endl; */
+
+    RNS.convert_transpose(n,n,Givaro::Integer(0), B, n, mod_A._ptr,mod_A._stride,RNSMAJOR);
+    ok= ok and (FFLAS::fequal(ZZ,n,n,At,n,B,n));
+
+      /* FFLAS::WriteMatrix(std::cout,ZZ,n, n, At, n)<<std::endl;//FFLAS::FflasSageMath)<<std::endl; */
+      /* FFLAS::WriteMatrix(std::cout,ZZ,n, n, B, n)<<std::endl;//FFLAS::FflasSageMath)<<std::endl; */
+
+      std::cout<<std::right<<std::setw(SPC)<<"RNSMajor=false ..."  <<  (ok? "PASSED\n": "FAILED\n");
+
+#endif
+
+#if 1
+      /**********************
+       ** check RNSMajor=true
+       **********************/
+      RNSMAJOR=true;
+
+      RNS.init(n,n,mod_A._ptr,mod_A._stride,A,n,p,RNSMAJOR);
+      ok = ok and  FFPACK::RNS_COMMON::check_integers_to_rns(FFLAS::FflasNoTrans, n, n, mod_A._ptr, mod_A._stride, A, n, RNS._basis, RNSMAJOR, false);
+
+      RNS.init_transpose(n,n,mod_At._ptr,mod_At._stride,At,n,p,RNSMAJOR);
+      ok = ok and  FFPACK::RNS_COMMON::check_integers_to_rns(FFLAS::FflasTrans, n, n, mod_At._ptr, mod_At._stride, At, n, RNS._basis, RNSMAJOR, false);
+
+      ok = ok and (FFLAS::fequal(ZD,n2,s,mod_A._ptr,s,mod_At._ptr,s));
+
+      /* std::cout<<std::boolalpha<<ok<<std::endl; */
+      /* FFLAS::WriteMatrix(std::cout,ZD,n2, s, mod_A._ptr, s)<<std::endl<<"--------------\n";;//FFLAS::FflasSageMath)<<std::endl; */
+      /* FFLAS::WriteMatrix(std::cout,ZD,n2, s, mod_At._ptr, s)<<std::endl;//FFLAS::FflasSageMath)<<std::endl; */
+
+
+			RNS.convert(n,n,Givaro::Integer(0), B, n, mod_A._ptr,mod_A._stride,RNSMAJOR);
+      ok= ok and (FFLAS::fequal(ZZ,n,n,A,n,B,n));
+
+      /* std::cout<<std::boolalpha<<ok<<std::endl; */
+      /* FFLAS::WriteMatrix(std::cout,ZZ,n, n, A, n)<<std::endl;//FFLAS::FflasSageMath)<<std::endl; */
+      /* FFLAS::WriteMatrix(std::cout,ZZ,n, n, B, n)<<std::endl;//FFLAS::FflasSageMath)<<std::endl; */ 
+      
+
+      RNS.convert_transpose(n,n,Givaro::Integer(0), B, n, mod_A._ptr,mod_A._stride,RNSMAJOR);
+      ok= ok and (FFLAS::fequal(ZZ,n,n,At,n,B,n));
+
+      /* std::cout<<std::boolalpha<<ok<<std::endl; */
+      /* FFLAS::WriteMatrix(std::cout,ZZ,n, n, At, n)<<std::endl;//FFLAS::FflasSageMath)<<std::endl; */
+      /* FFLAS::WriteMatrix(std::cout,ZZ,n, n, B, n)<<std::endl;//FFLAS::FflasSageMath)<<std::endl; */
+
+
+      std::cout<<std::right<<std::setw(SPC)<<"RNSMajor=true  ..."  <<  (ok? "PASSED\n": "FAILED\n");
+
+
+#endif
+
+			FFLAS::fflas_delete(mod_A);
+      FFLAS::fflas_delete(mod_At);
+      FFLAS::fflas_delete(A);
+      FFLAS::fflas_delete(At);
+      FFLAS::fflas_delete(B);
+		}
+
+  std::cout<<"Test suite for rns-double ...."<<(ok?"PASSED\n":"FAILED\n");
+  if (!ok) std::cout<<"seed was "<<seed<<std::endl;
+  return ok;
+}
+
+
+
+int main(int argc, char** argv){
+	static size_t iters = 1 ;
+	static unsigned long b = 512 ;
+	static unsigned long p = 20 ;
+	static size_t m = 16 ;
+	static size_t seed= time(NULL);	
+	static Argument as[] = {
+		{ 'b', "-b B", "Set the bitsize of the matrix entries.",         TYPE_INT , &b },
+		{ 'p', "-p B", "Set the bitsize of the RNS Prime.",              TYPE_INT , &p },
+		{ 'm', "-m M", "Set the dimension m of the matrix.",             TYPE_INT , &m },
+		{ 'i', "-i R", "Set minimal number of repetitions.",             TYPE_INT , &iters },
+		{ 's', "-s S", "Set the seeding value",                          TYPE_INT , &seed },
+		END_OF_ARGUMENTS
+	};
+	FFLAS::parseArguments(argc,argv,as);
+
+	/* int logb= log((double)n)/log(2.0); */
+	/* size_t fflas_primes_bits       = (primes_bits ? primes_bits : 41-logb); */
+	/* size_t fflas_primes_bits_ext   = 2* fflas_primes_bits; */
+
+  bool ok=true;
+  ok = ok and run_check<FFPACK::rns_double>         (m,  p,b,iters,seed);
+  ok = ok and run_check<FFPACK::rns_double_extended>(m,2*p,b,iters,seed);
+
+	return !ok;
+}
+
+
+/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
+