diff --git a/src/apply.cc b/src/apply.cc
index 9adf1158c6e83adee1cae204efe638c8de5cb199..01f222eba286b8e08db35c95752b2f25050167fe 100644
--- a/src/apply.cc
+++ b/src/apply.cc
@@ -33,18 +33,10 @@ SEXP gmpMatToListZ(SEXP X, SEXP line)
bigvec oneLine ;
for(unsigned int j = 0; j < ncol; ++j)
{
- oneLine.value.push_back(matrix.value[i+j*nrow]);
-
- // modulus, if one by cell
- if(matrix.modulus.size() ==matrix.value.size() )
- oneLine.modulus.push_back(matrix.modulus[i+j*nrow]);
+ oneLine.push_back(matrix[i+j*nrow]);
}
- // modulus, if one by line
- if(((matrix.modulus.size() == nrow ) || (matrix.modulus.size() == 1) ) && (matrix.modulus.size() !=matrix.value.size()) )
- oneLine.modulus.push_back(matrix.modulus[i % matrix.modulus.size() ]);
-
SET_VECTOR_ELT(ans, i,bigintegerR::create_SEXP(oneLine));
@@ -60,19 +52,10 @@ SEXP gmpMatToListZ(SEXP X, SEXP line)
bigvec oneLine ;
for(unsigned int i = 0; i < nrow; ++i)
{
- oneLine.value.push_back(matrix.value[i+j*nrow]);
-
- // modulus, if one by cell
- if(matrix.modulus.size() ==matrix.value.size() )
- oneLine.modulus.push_back(matrix.modulus[i+j*nrow]);
+ oneLine.push_back(matrix[i+j*nrow]);
}
- // modulus, if one by line
- if( (matrix.modulus.size() == 1) && (matrix.modulus.size() !=matrix.value.size()) )
- oneLine.modulus.push_back(matrix.modulus[0 ]);
-
-
SET_VECTOR_ELT(ans, j,bigintegerR::create_SEXP(oneLine));
}
diff --git a/src/biginteger.cc b/src/biginteger.cc
index 0beb5885fea584f38cf8c6f7b7e74c3bff48e0c1..409fa50d9cadf3e8e2a94017abef38d303606b73 100644
--- a/src/biginteger.cc
+++ b/src/biginteger.cc
@@ -4,7 +4,7 @@
* \version 1
*
* \date Created: 27/10/04
- * \date Last modified: Time-stamp: <2023-01-16 18:48:41 (antoine)>
+ * \date Last modified: Time-stamp: <2023-01-21 21:01:28 (antoine)>
*
* \author Immanuel Scholz
*
@@ -17,10 +17,13 @@
#include "biginteger.h"
using std::string;
+static int count=0;
+static int countALL=0;
biginteger::biginteger()
: na(true) {
-
+ count++;
+ countALL++;
mpz_init(value);}
@@ -29,6 +32,8 @@ biginteger::biginteger()
* Construct a biginteger from a int value.
*/
biginteger::biginteger(const int value_) : na(false) {
+ count++;
+ countALL++;
if(value_ == NA_INTEGER)
{mpz_init(value); na = true ;}
else
@@ -38,6 +43,8 @@ biginteger::biginteger(const int value_) : na(false) {
* Construct a biginteger from a long value.
*/
biginteger::biginteger(const long int value_) : na(false) {
+ count++;
+ countALL++;
if(value_ == NA_INTEGER)
{mpz_init(value); na = true ;}
else
@@ -47,6 +54,8 @@ biginteger::biginteger(const long int value_) : na(false) {
* Construct a biginteger from a unsigned long value.
*/
biginteger:: biginteger(const unsigned long int value_) : na(false) {
+ count++;
+ countALL++;
if(value_ == (unsigned long int) NA_INTEGER)
{mpz_init(value); na = true ;}
else
@@ -56,6 +65,8 @@ biginteger:: biginteger(const unsigned long int value_) : na(false) {
* Construct a biginteger from a double value.
*/
biginteger::biginteger(const double value_) : na(false) {
+ count++;
+ countALL++;
if( R_FINITE(value_) )
mpz_init_set_d(value, value_);
else
@@ -67,6 +78,8 @@ biginteger::biginteger(const double value_) : na(false) {
*/
biginteger::biginteger(const std::string& value_) : na(false)
{
+ count++;
+ countALL++;
/* mpz_init.. return -1 when error, 0: ok */
if(mpz_init_set_str(value, value_.c_str(), 0))
{
@@ -82,6 +95,8 @@ biginteger::biginteger(const std::string& value_) : na(false)
*/
biginteger::biginteger(const biginteger& rhs) : na(rhs.na)
{
+ count++;
+ countALL++;
mpz_init_set(value, rhs.getValueTemp());
}
@@ -89,10 +104,16 @@ biginteger::biginteger(const biginteger& rhs) : na(rhs.na)
/**
* Free the owned mpz_t structs
*/
-biginteger::~biginteger() {mpz_clear(value);}
+biginteger::~biginteger() {
+ count--;
+ // printf("nb : %d %d \n",count,countALL);
+ mpz_clear(value);
+}
biginteger::biginteger(const char* raw)
{
+ count++;
+ countALL++;
mpz_init(value);
na = true;
const int* r = (int*)(raw);
@@ -117,6 +138,8 @@ biginteger::biginteger(const char* raw)
biginteger::biginteger(const mpz_t value_)
: na(false)
{
+ count++;
+ countALL++;
mpz_init_set(value, value_);
}
diff --git a/src/biginteger.h b/src/biginteger.h
index 3277ab97299986ea29165542b2721d7dda2d67dd..6c9fecd36876f8ad1a22f0b384e7dcfcb7349030 100644
--- a/src/biginteger.h
+++ b/src/biginteger.h
@@ -2,7 +2,7 @@
* \brief Description of class biginteger
*
* \date Created: 2004
- * \date Last modified: Time-stamp: <2023-01-16 19:17:13 (antoine)>
+ * \date Last modified: Time-stamp: <2023-01-17 18:20:27 (antoine)>
*
* \author Immanuel Scholz
*
diff --git a/src/bigintegerR.cc b/src/bigintegerR.cc
index 106e7b7f27eef1626a17a8c37c83da29a309003c..666510850002e4356917db17d805843604d27670 100644
--- a/src/bigintegerR.cc
+++ b/src/bigintegerR.cc
@@ -31,7 +31,6 @@ using namespace std;
static gmp_randstate_t seed_state;
static int seed_init=0;
-
namespace bigintegerR
{
// \brief create a vector of bigvecs, all without a modulus.
@@ -48,10 +47,10 @@ namespace bigintegerR
int pos = sizeof(int); // position in raw[]. Starting after header.
int sizevec = ((int*)raw)[0];
//std::cout << "nb element a lire " << sizevec << std::endl;
- v.value.resize(sizevec);
+ v.resize(sizevec);
for (int i = 0; i < sizevec; ++i) {
- v.value[i] = biginteger(&raw[pos]);
- pos += v.value[i].raw_size(); // increment number of bytes read.
+ v[i].setValue(make_shared<biginteger>(&raw[pos]));
+ pos += v[i].getValue().raw_size(); // increment number of bytes read.
}
return v;
}
@@ -60,30 +59,33 @@ namespace bigintegerR
double* d = REAL(param);
//bigvec v(d,d+LENGTH(param));
bigvec v;
- v.value.resize(LENGTH(param));
+ v.resize(LENGTH(param));
for (int j = 0; j < LENGTH(param); ++j) {
/// New: numeric '+- Inf' give +- "Large" instead of NA
double dj = d[j];
if(R_FINITE(dj) || ISNAN(dj))
- v.value[j] = dj;
+ v[j].setValue(make_shared<biginteger>( dj ));
else { // dj is +- Inf : use LARGE ( = +- 2 ^ 80000 -- arbitrarily )
- mpz_t LARGE;
- mpz_init(LARGE);
- // FIXME: Keep 'LARGE' a static const; initialized only once
+ mpz_t LARGE;
+ mpz_init(LARGE);
+ // FIXME: Keep 'LARGE' a static const; initialized only once
+
+ if(dj == R_PosInf){
mpz_ui_pow_ui (LARGE, (unsigned long int) 2, (unsigned long int) 8000);
- if(dj == R_PosInf)
- v.value[j] = LARGE;
- else if(dj == R_NegInf) {
- mpz_t neg_L;
- mpz_init(neg_L);
- mpz_neg(neg_L, LARGE);
- v.value[j] = neg_L;
- mpz_clear(neg_L);
- }
- else// should never happen
- v.value[j] = dj;
-
- mpz_clear(LARGE);
+ v[j].setValue(make_shared<biginteger>(LARGE));
+ }
+ else if(dj == R_NegInf) {
+ mpz_t neg_L;
+ mpz_init(neg_L);
+ mpz_neg(neg_L, LARGE);
+ v[j].setValue( make_shared<biginteger>(neg_L));
+ mpz_clear(neg_L);
+ }
+ else// should never happen
+ v[j].setValue( make_shared<biginteger>(dj));
+
+ mpz_clear(LARGE);
+
}
}
return v;
@@ -94,22 +96,20 @@ namespace bigintegerR
int* i = INTEGER(param);
//bigvec v(i,i+LENGTH(param));
bigvec v;
- v.value.resize(LENGTH(param));
+ v.resize(LENGTH(param));
for (int j = 0; j < LENGTH(param); ++j)
- v.value[j] = i[j];
+ v[j].setValue( make_shared<biginteger>(i[j]));
return v;
}
case STRSXP:
{
bigvec v;
- v.value.resize(LENGTH(param));
+ v.resize(LENGTH(param));
for (int i = 0; i < LENGTH(param); ++i) {
- if (STRING_ELT(param,i) == NA_STRING)
- v.value[i]= biginteger();
- else
- v.value[i]=biginteger(std::string(CHAR(STRING_ELT(param,i))));
+ if (STRING_ELT(param,i) != NA_STRING)
+ v[i].setValue(make_shared<biginteger>(std::string(CHAR(STRING_ELT(param,i)))));
}
return v;
}
@@ -141,7 +141,11 @@ namespace bigintegerR
if (TYPEOF(modAttr) != NILSXP)
{
//std::cout << "import value" << std::endl;
- v.modulus = bigintegerR::create_vector(modAttr).value;
+ bigvec vMod = bigintegerR::create_vector(modAttr);
+ for (int i = 0 ; i < v.size(); i++){
+ v[i].setModulus(vMod[i % vMod.size()].getValuePtr());
+ }
+ v.setType(vMod.size() == 1 ? TYPE_MODULUS::MODULUS_GLOBAL : TYPE_MODULUS::MODULUS_BY_CELL);
}
return v;
@@ -155,7 +159,6 @@ namespace bigintegerR
double* d = REAL(param);
// copy vector manually to avoid stupid conversion warning in STL-code :-/
vector<int> v;
- v.reserve(LENGTH(param));
for (int i = 0; i < LENGTH(param); ++i)
v.push_back(static_cast<int>(d[i]));
return v;
@@ -174,26 +177,35 @@ namespace bigintegerR
}
- SEXP create_SEXP(const std::vector<biginteger>& v)
+ SEXP create_SEXP(const bigvec & v , bigmod_accessor_fn fct, int size)
{
unsigned int i;
- int size = sizeof(int); // starting with vector-size-header
- for (i = 0; i < v.size(); ++i)
- size += v[i].raw_size(); // adding each bigint's needed size
- SEXP ans = PROTECT(Rf_allocVector(RAWSXP, size));
+ int sizeRaw = sizeof(int); // starting with vector-size-header
+ for (i = 0; i < size; ++i)
+ sizeRaw += fct(v[i]).raw_size(); // adding each bigint's needed size
+ SEXP ans = PROTECT(Rf_allocVector(RAWSXP, sizeRaw));
// Rprintf(" o create_SEXP(vect<biginteger>): size=%d, v.size()=%d\n", size, v.size());
char* r = (char*)(RAW(ans));
- ((int*)(r))[0] = v.size(); // first int is vector-size-header
+ ((int*)(r))[0] = size; // first int is vector-size-header
int pos = sizeof(int); // current position in r[] (starting after vector-size-header)
- for (i = 0; i < v.size(); ++i)
- pos += v[i].as_raw(&r[pos]);
+ for (i = 0; i < size; ++i)
+ pos += fct(v[i]).as_raw(&r[pos]);
UNPROTECT(1);
return(ans);
}
- SEXP create_SEXP(const bigvec& v) {
- SEXP ans = PROTECT(create_SEXP(v.value));
+ const biginteger & bigModToValue(const bigmod& v){
+ return v.getValue();
+ }
+
+ const biginteger & bigModToModulus(const bigmod& v){
+ return v.getModulus();
+ }
+
+ SEXP create_SEXP(const bigvec& v) {
+ int size = v.size();
+ SEXP ans = PROTECT(create_SEXP(v, bigModToValue, size));
// set the class attribute to "bigz"
Rf_setAttrib(ans, R_ClassSymbol, Rf_mkString("bigz"));
@@ -208,8 +220,9 @@ namespace bigintegerR
UNPROTECT(2);
}
// set the mod attribute
- if(v.modulus.size() > 0) {
- SEXP mod = PROTECT(create_SEXP(v.modulus)); // and set *its* class
+ if(v.getType() != TYPE_MODULUS::NO_MODULUS && v.size() > 0 ) {
+ int sizeM = v.getType() == TYPE_MODULUS::MODULUS_GLOBAL ? 1 : size;
+ SEXP mod = PROTECT(create_SEXP(v, bigModToModulus, sizeM )); // and set *its* class
Rf_setAttrib(mod, R_ClassSymbol, Rf_mkString("bigz"));
Rf_setAttrib(ans, Rf_mkString("mod"), mod);
UNPROTECT(1);
@@ -228,12 +241,21 @@ namespace bigintegerR
{
bigvec va = bigintegerR::create_bignum(a);
bigvec vb = bigintegerR::create_bignum(b), result;
- int size = (va.value.empty() || vb.value.empty()) ? 0 : max(va.value.size(), vb.value.size());
- result.value.reserve(size);
- for (int i = 0; i < size; ++i)
+ int size = (va.size() ==0 || vb.size()==0 ) ? 0 : max(va.size(), vb.size());
+
+
+ int nrow = matrixz::checkDims(va.nrow,vb.nrow);
+ if(result.nrow == -2){
+ va.clear();
+ vb.clear();
+ error(_("Matrix dimensions do not match"));
+ }
+
+ for (int i = 0; i < size; ++i)
result.push_back(f(va[i%va.size()], vb[i%vb.size()]));
- result.nrow = matrixz::checkDims(va.nrow,vb.nrow);
+ result.nrow = nrow;
+
// Rprintf(" o bigI_b_op(.); size=%d -> nrow=%d\n", size, result.nrow);
return bigintegerR::create_SEXP(result);
}
@@ -243,22 +265,28 @@ namespace bigintegerR
{
bigvec va = bigintegerR::create_bignum(a);
bigvec vb = bigintegerR::create_bignum(b), result;
- int size = (va.value.empty() || vb.value.empty()) ? 0 : max(va.value.size(), vb.value.size());
+
+ int nrow = matrixz::checkDims(va.nrow,vb.nrow) ;
+ if(nrow == -2){
+ va.clear();
+ vb.clear();
+ error(_("Matrix dimensions do not match"));
+ }
+
+ int size = (va.size()==0 || vb.size() == 0) ? 0 : max(va.size(), vb.size());
// int sizemod = max(va.modulus.size(), vb.modulus.size());
SEXP ans = PROTECT(Rf_allocVector(LGLSXP, size));
int *r = LOGICAL(ans);
/* TODO: this kind of situation 5 == (5 %% 17)*/
for (int i = 0; i < size; ++i) {
- biginteger am = va.value[i % va.value.size()];
- biginteger bm = vb.value[i % vb.value.size()];
+ biginteger & am = va[i % va.size()].getValue();
+ biginteger & bm = vb[i % vb.size()].getValue();
if (am.isNA() || bm.isNA())
r[i] = NA_LOGICAL;
else
r[i] = f(am, bm) ? 1 : 0;
}
- int nrow = matrixz::checkDims(va.nrow,vb.nrow) ;
-
// Add dimension parameter when available
if(nrow >= 0)
{
@@ -306,32 +334,43 @@ SEXP biginteger_div (SEXP a, SEXP b) { // called from "/.bigz" == div.bigz
bigvec A = bigintegerR::create_bignum(a),
B = bigintegerR::create_bignum(b);
// Note: a or b may be simple numbers (e.g. '1') => work with (A,B)
- int len_m_a = A.modulus.size(),
- len_m_b = B.modulus.size();
- if(len_m_a == 0 && len_m_b == 0) // deal with important case quickly:
+ TYPE_MODULUS type_a = A.getType();
+ TYPE_MODULUS type_b = B.getType();
+
+ if(type_a == TYPE_MODULUS::NO_MODULUS && type_b == TYPE_MODULUS::NO_MODULUS) // deal with important case quickly:
+ {
return bigrational_div(a, b);
- else if(len_m_a == 0) { // and len_m_b > 0:
+ }
+ else if(type_a == TYPE_MODULUS::NO_MODULUS) { // and len_m_b > 0:
// should work directly using b's "mod" --> compute a * b^(-1)
}
- else if(len_m_b == 0) { // and len_m_a > 0:
- // should use a's "mod" for b: div_via_inv() need's b's modulus
- B.modulus = A.modulus;
- return bigintegerR::biginteger_binary_operation(a,
- bigintegerR::create_SEXP(B),
- div_via_inv);
+ else if(type_b == TYPE_MODULUS::NO_MODULUS) { // and len_m_a > 0:
+ // should use a's "mod" for b: div_via_inv() need's b's modulus
+ if(type_a == TYPE_MODULUS::MODULUS_GLOBAL){
+ B.setGlobalModulus(A.getGlobalModulus());
+ } else{
+ for (unsigned int i = 0 ; i < B.size(); i++){
+ B[i].setModulus(A[i % A.size()].getModulusPtr());
+ }
+ }
+ return bigintegerR::biginteger_binary_operation(a,
+ bigintegerR::create_SEXP(B),
+ div_via_inv);
}
else { // len_m_a > 0 and len_m_b > 0:
bool same_mod = true;// are the two mods the "same" (after recycling)?
+ int len_m_a = A.getModulusSize();
+ int len_m_b = B.getModulusSize();
int m = (len_m_a < len_m_b) ? len_m_b : len_m_a; // = max(l..a, l..b)
- for(int i = 0; i < m; i++)
- if(A.modulus[i % len_m_a] != B.modulus[i % len_m_b]) {
- same_mod = false; break;
- }
+ for(unsigned int i = 0; i < m; i++)
+ if(A[i % len_m_a].getModulus() != B[i % len_m_b].getModulus()) {
+ same_mod = false; break;
+ }
if(same_mod) {
- // compute a * b^(-1) ... should work w/o more
+ // compute a * b^(-1) ... should work w/o more
} else {
- // use *rational* a/b (not considering 'mod' anymore):
- return bigrational_div(a, b);
+ // use *rational* a/b (not considering 'mod' anymore):
+ return bigrational_div(a, b);
}
}
return bigintegerR::biginteger_binary_operation(a,b, div_via_inv);
@@ -340,11 +379,11 @@ SEXP biginteger_div (SEXP a, SEXP b) { // called from "/.bigz" == div.bigz
SEXP biginteger_pow (SEXP a, SEXP b) {
bigvec v = bigintegerR::create_bignum(a),
exp = bigintegerR::create_bignum(b);
- if(v.modulus.size() == 0) { /* has no modulus: now, if any b < 0, the
+ if(v.getType() == TYPE_MODULUS::NO_MODULUS) { /* has no modulus: now, if any b < 0, the
result must be (non-integer) bigrational */
bool use_rat = FALSE;
- for (unsigned int i = 0; i < exp.value.size(); ++i) {
- if(mpz_sgn(exp.value[i].getValueTemp()) < 0) {
+ for (unsigned int i = 0; i < exp.size(); ++i) {
+ if(mpz_sgn(exp[i].getValue().getValueTemp()) < 0) {
use_rat = TRUE;
break;
}
@@ -384,9 +423,11 @@ SEXP biginteger_as_character(SEXP a, SEXP b)
bigvec v = bigintegerR::create_bignum(a);
SEXP ans;
int base = Rf_asInteger(b);
- if (base < 2 || base > 36)
+ if (base < 2 || base > 36){
+ v.clear();
error(_("select a base between 2 and 36"));
-
+ }
+
PROTECT(ans = Rf_allocVector(STRSXP, v.size()));
for (unsigned int i = 0; i < v.size(); ++i)
SET_STRING_ELT(ans, i, Rf_mkChar(v.str(i,base).c_str()));
@@ -395,7 +436,7 @@ SEXP biginteger_as_character(SEXP a, SEXP b)
{
SEXP dim = PROTECT(Rf_allocVector(INTSXP, 2));
INTEGER(dim)[0] = v.nrow;
- INTEGER(dim)[1] = v.value.size() / v.nrow;
+ INTEGER(dim)[1] = v.size() / v.nrow;
Rf_setAttrib(ans, Rf_mkString("dim"), dim);
UNPROTECT(1);
}
@@ -410,7 +451,7 @@ SEXP biginteger_as_numeric(SEXP a)
SEXP ans = PROTECT(Rf_allocVector(REALSXP,v.size()));
double *r = REAL(ans);
for (unsigned int i = 0; i < v.size(); ++i)
- r[i] = v.value[i].isNA() ? NA_REAL : v.value[i].as_double();
+ r[i] = v[i].isNA() ? NA_REAL : v[i].getValue().as_double();
UNPROTECT(1);
return ans;
}
@@ -421,14 +462,14 @@ SEXP biginteger_as_integer(SEXP a)
SEXP ans = PROTECT(Rf_allocVector(INTSXP,v.size()));
int *r = INTEGER(ans);
for (unsigned int i = 0; i < v.size(); ++i) {
- if(v.value[i].isNA()) {
+ if(v[i].isNA()) {
r[i] = NA_INTEGER;
}
- else if(!mpz_fits_slong_p(v.value[i].getValueTemp())) {
+ else if(!mpz_fits_slong_p(v[i].getValue().getValueTemp())) {
Rf_warning("NAs introduced by coercion from big integer");
r[i] = NA_INTEGER;
} else {
- r[i] = mpz_get_si(v.value[i].getValueTemp());
+ r[i] = mpz_get_si(v[i].getValue().getValueTemp());
}
}
UNPROTECT(1);
@@ -476,6 +517,8 @@ SEXP biginteger_set_at(SEXP src, SEXP idx, SEXP value)
}
if(vvalue.size() == 0) {
+ vvalue.clear();
+ result.clear();
error(_("replacement has length zero"));
}
@@ -560,8 +603,9 @@ SEXP biginteger_c(SEXP args)
bigvec result;
for(int i=0; i < LENGTH(args); i++) {
bigvec v = bigintegerR::create_bignum(VECTOR_ELT(args,i));
- for(unsigned int j=0; j < v.size(); j++)
+ for(unsigned int j=0; j < v.size(); j++){
result.push_back(v[j]);
+ }
v.clear();
}
return bigintegerR::create_SEXP(result);
@@ -574,7 +618,7 @@ SEXP biginteger_rep(SEXP x, SEXP times)
result;
int rep = Rf_asInteger(times);
- result.value.reserve(v.size()*rep);
+
for(int i = 0 ; i < rep ; i++)
for(unsigned int j = 0 ; j < v.size() ; j++)
result.push_back(v[j]);
@@ -605,7 +649,7 @@ SEXP biginteger_nextprime(SEXP a)
{
bigvec v = bigintegerR::create_bignum(a),
result;
- result.value.reserve(v.size());
+
mpz_t val;
mpz_init(val);
@@ -613,7 +657,7 @@ SEXP biginteger_nextprime(SEXP a)
for (unsigned int i = 0; i < v.size(); ++i) {
mpz_nextprime(val,v[i].getValue().getValueTemp());
- result.push_back(DefaultBigMod(val));
+ result.push_back(bigmod(val));
}
return bigintegerR::create_SEXP(result);
}
@@ -622,7 +666,7 @@ SEXP biginteger_abs(SEXP a)
{
bigvec v = bigintegerR::create_bignum(a),
result;
- result.value.reserve(v.size());
+
mpz_t val;
mpz_init(val);
@@ -631,14 +675,9 @@ SEXP biginteger_abs(SEXP a)
for (unsigned int i = 0; i < v.size(); ++i)
{
mpz_abs(val,v[i].getValue().getValueTemp());
- result.push_back(DefaultBigMod(val));
-
- // TODO: understand why following lines don't work.
- //result.push_back(bigmod());
- //result[i].value.setValue(val);
+ result.push_back(bigmod(make_shared<biginteger>(val),v[i].getModulusPtr()));
}
- result.modulus = v.modulus;
return bigintegerR::create_SEXP(result);
}
@@ -657,7 +696,6 @@ SEXP biginteger_gcdex(SEXP a, SEXP b)
if(va.size() != vb.size())
return bigintegerR::create_SEXP(bigvec());
- result.value.reserve(3*va.size());
mpz_t g;
mpz_t s;
@@ -668,17 +706,15 @@ SEXP biginteger_gcdex(SEXP a, SEXP b)
mpz_t_sentry val_g(g);
mpz_t_sentry val_s(s);
mpz_t_sentry val_t(t);
+ shared_ptr<biginteger> mod = make_shared<biginteger>(); // NA
for(unsigned int i=0; i < va.size(); i++)
{
mpz_gcdext (g,s,t,va[i].getValue().getValueTemp(),vb[i].getValue().getValueTemp());
- result.value.push_back(biginteger(g)); // Hem... not very elegant !
- result.value.push_back(biginteger(s));
- result.value.push_back(biginteger(t));
- /* result[i*3].value.setValue(g);
- result[i*3+1].value.setValue(s);
- result[i*3+2].value.setValue(t);*/
-
+ result.push_back(bigmod(make_shared<biginteger>(g)));
+ result.push_back(bigmod(make_shared<biginteger>(s)));
+ result.push_back(bigmod(make_shared<biginteger>(t)));
+
}
return bigintegerR::create_SEXP(result);
}
@@ -710,7 +746,6 @@ SEXP biginteger_rand_u (SEXP nb, SEXP length, SEXP newseed, SEXP ok)
size = Rf_asInteger(nb);
UNPROTECT(3);
- result.value.reserve(size);
/* Random seed initialisation */
@@ -734,7 +769,7 @@ SEXP biginteger_rand_u (SEXP nb, SEXP length, SEXP newseed, SEXP ok)
{
/* Random number generation */
mpz_urandomb(bz,seed_state,len);
- result.push_back(DefaultBigMod(bz));
+ result.push_back(bigmod(bz));
}
return bigintegerR::create_SEXP(result);
}
@@ -764,11 +799,11 @@ SEXP bigI_factorial(SEXP n)
{
bigvec result;
int *nn = INTEGER(AS_INTEGER(n)), size = Rf_length(n);
- result.value.resize(size);
+ result.resize(size);
for (int i = 0; i < size; ++i) {
- result.value[i].NA(false);
+ result[i].getValue().NA(false);
if(nn[i] != NA_INTEGER && nn[i] >= 0) {
- mpz_fac_ui(result.value[i].getValue(), (unsigned long int)nn[i]);
+ mpz_fac_ui(result[i].getValue().getValue(), (unsigned long int)nn[i]);
}
}
return bigintegerR::create_SEXP(result);
@@ -782,20 +817,23 @@ SEXP bigI_choose(SEXP n, SEXP k)
{
bigvec result, n_ = bigintegerR::create_bignum(n);
int *kk = INTEGER(AS_INTEGER(k)), n_k = Rf_length(k);
- int size = (n_.value.empty() || n_k == 0) ? 0 :
+ int size = (n_.size() == 0 || n_k == 0) ? 0 :
// else: max(n_.value.size(), n_k)
- (((int)n_.value.size() <= n_k) ? n_k : n_.value.size());
+ (((int)n_.size() <= n_k) ? n_k : n_.size());
- result.value.resize(size);
+ result.resize(size);
+ // if(n_.getType() == MODULUS_GLOBAL){
+ // result.setGlobalModulus(n_.getGlobalModulus());
+ //}
for (int i = 0; i < size; ++i) {
- result.value[i].NA(false);
+ result[i].getValue().NA(false);
int ik_i = kk[i % n_k];
// check if k in range:
if(ik_i != NA_INTEGER && ik_i >= 0) {
unsigned long int k_i = (unsigned long int)ik_i;
/* void mpz_bin_ui (mpz_t ROP, mpz_t N, unsigned long int K) */
- mpz_bin_ui(result.value[i].getValue(),
- n_.value[i % n_.value.size()].getValueTemp(), k_i);
+ mpz_bin_ui(result[i].getMpValue(),
+ n_[i % n_.size()].getValueTemp(), k_i);
}
}
return bigintegerR::create_SEXP(result);
@@ -818,7 +856,7 @@ SEXP bigI_fibnum(SEXP n)
mpz_t_sentry val_s(val);
mpz_fib_ui(val,num);
- result.push_back(DefaultBigMod(val));
+ result.push_back(bigmod(val));
// result[0].value.setValue(val);
}
// else
@@ -840,7 +878,7 @@ SEXP bigI_fibnum2(SEXP n)
unsigned long int num = nn;
if(nn < 0 || nn == NA_INTEGER)
error(_("argument must be non-negative"));
- result.value.reserve(1);
+
mpz_t val;
mpz_init(val);
mpz_t_sentry val_s(val);
@@ -849,8 +887,8 @@ SEXP bigI_fibnum2(SEXP n)
mpz_t_sentry val_s2(val2);
mpz_fib2_ui(val,val2, num);
- result.push_back(DefaultBigMod(val2));
- result.push_back(DefaultBigMod(val));
+ result.push_back(bigmod(val2));
+ result.push_back(bigmod(val));
}
else
error(_("argument must not be an empty list"));
@@ -877,7 +915,7 @@ SEXP bigI_lucnum(SEXP n)
mpz_t_sentry val_s(val);
mpz_lucnum_ui(val,num);
- result.push_back(DefaultBigMod(val));
+ result.push_back(bigmod(val));
}
// else
// error(_("argument must not be an empty list"));
@@ -906,8 +944,8 @@ SEXP bigI_lucnum2(SEXP n)
mpz_t_sentry val_s2(val2);
mpz_lucnum2_ui(val,val2,num);
- result.push_back(DefaultBigMod(val2));
- result.push_back(DefaultBigMod(val));
+ result.push_back(bigmod(val2));
+ result.push_back(bigmod(val));
}
else
error(_("argument must not be an empty list"));
@@ -932,31 +970,16 @@ SEXP biginteger_max(SEXP a, SEXP narm)
for(unsigned int i = 1 ; i < va.size(); ++i)
{
- if(va.value[i].isNA() && !na_remove)
+ if(va[i].isNA() && !na_remove)
return(bigintegerR::create_SEXP(result));
else
- if(!(va.value[i] < va.value[maximum] ))
+ if(va[i].getValue() > va[maximum].getValue() )
maximum = i; // if va.value[maximum = 0] is NA => false for the "<" => maximum changed = good
}
- result.push_back(va.value[maximum]);
+ result.push_back(va[maximum]);
-
- // now the modulus !
- if(va.modulus.size() == 1)
- result.modulus.push_back(va.modulus[0]);
-
- if(va.modulus.size()>1)
- {
- biginteger modulus ;
- modulus.setValue(va.modulus[0].getValueTemp());
- for(unsigned int i = 1 ; i < va.modulus.size() ; ++i)
- {
- if(modulus != va.modulus[i]) // if one is different: no modulus
- return(bigintegerR::create_SEXP(result));
- }
- result.modulus.push_back(modulus);
- }
+ if(va.getType() == TYPE_MODULUS::MODULUS_BY_CELL) result[0].getModulus().NA(TRUE);
return(bigintegerR::create_SEXP(result));
}
@@ -976,30 +999,17 @@ SEXP biginteger_min(SEXP a, SEXP narm)
for(unsigned int i = 1 ; i < va.size(); ++i)
{
- if(va.value[i].isNA() && !na_remove)
+ if(va[i].isNA() && !na_remove)
return bigintegerR::create_SEXP(result);
else
- if(!(va.value[i] > va.value[minimum] ))
+ if(va[i].getValue() < va[minimum].getValue() )
minimum = i;
}
- result.push_back(va.value[minimum]);
+ result.push_back(va[minimum]);
// now the modulus !
- if(va.modulus.size() == 1)
- result.modulus.push_back(va.modulus[0]);
-
- if(va.modulus.size()>1)
- {
- biginteger modulus ;
- modulus.setValue(va.modulus[0].getValueTemp());
- for(unsigned int i = 1 ; i < va.modulus.size() ; ++i)
- {
- if(modulus != va.modulus[i]) // if one is different: no modulus
- return bigintegerR::create_SEXP(result);
- }
- result.modulus.push_back(modulus);
- }
+ if(va.getType() == TYPE_MODULUS::MODULUS_BY_CELL) result[0].getModulus().NA(TRUE);
return bigintegerR::create_SEXP(result);
}
@@ -1009,48 +1019,31 @@ SEXP biginteger_cumsum(SEXP a)
{
bigvec result, va = bigintegerR::create_bignum(a);
- result.value.resize(va.value.size());
+ result.resize(va.size());
mpz_t val;
mpz_init(val);
mpz_t_sentry val_s(val);
- bool hasmodulus = true;
-
- // first the modulus !
- if(va.modulus.size() > 1) {
- biginteger modulus ;
- modulus.setValue(va.modulus[0].getValueTemp());
-
- for(unsigned int i = 1 ; i < va.modulus.size() ; ++i) {
- if(modulus != va.modulus[i]) { // if one is different: no modulus
- hasmodulus = false;
- break;
- }
- }
- if(hasmodulus)
- result.modulus.push_back(modulus);
- }
- else if(va.modulus.size() == 1) {
- result.modulus.push_back(va.modulus[0]);
- hasmodulus = true;
- }
- else hasmodulus = false;
+ bool hasmodulus = va.getType() == TYPE_MODULUS::MODULUS_GLOBAL;
+
for(unsigned int i = 0 ; i < va.size(); ++i)
{
{
- if(va.value[i].isNA() )
+ if(va[i].isNA() )
{
break; // all last values are NA.
}
- mpz_add(val,val,va.value[i].getValueTemp());
-
- if(hasmodulus)
- mpz_mod(val,val,va.modulus[0].getValueTemp() );
+ mpz_add(val,val,va[i].getValue().getValueTemp());
- result.value[i].setValue(val);
+ if(hasmodulus){
+ mpz_mod(val,val,va.getGlobalModulus()->getValueTemp() );
+ result[i].setModulus(va.getGlobalModulus());
+ }
+ result[i].getValue().setValue(val);
+
}
}
@@ -1063,55 +1056,33 @@ SEXP biginteger_sum(SEXP a)
{
bigvec result, va = bigintegerR::create_bignum(a);
- result.value.resize(1);
+ result.resize(1);
mpz_t val;
mpz_init(val);
mpz_t_sentry val_s(val);
- bool hasmodulus = true;
-
- // first the modulus !
- if(va.modulus.size() > 1) {
- biginteger modulus ;
- modulus.setValue(va.modulus[0].getValueTemp());
-
- for(unsigned int i = 1 ; i < va.modulus.size() ; ++i) {
- if(modulus != va.modulus[i]) { // if one is different: no modulus
- hasmodulus = false;
- break;
- }
- }
- if(hasmodulus)
- result.modulus.push_back(modulus);
- }
- else if(va.modulus.size() == 1) {
- result.modulus.push_back(va.modulus[0]);
- hasmodulus = true;
- }
- else
- hasmodulus = false;
-
-
+ bool hasmodulus = va.getType() == TYPE_MODULUS::MODULUS_GLOBAL;
for(unsigned int i = 0 ; i < va.size(); ++i)
{
{
- if(va.value[i].isNA() )
+ if(va[i].isNA() )
{
break; // all last values are NA.
}
- mpz_add(val,val,va.value[i].getValueTemp());
+ mpz_add(val,val,va[i].getValueTemp());
if(hasmodulus)
- mpz_mod(val,val,va.modulus[0].getValueTemp() );
+ mpz_mod(val,val,va.getGlobalModulus()->getValueTemp() );
}
}
- result.value[0].setValue(val);
+ result[0].setValue(val);
+ if(hasmodulus) result[0].setModulus(va.getGlobalModulus());
return(bigintegerR::create_SEXP(result));
}
@@ -1122,54 +1093,31 @@ SEXP biginteger_prod(SEXP a)
bigvec result;
bigvec va = bigintegerR::create_bignum(a);
- result.value.resize(1);
+ result.resize(1);
mpz_t val;
mpz_init(val);
mpz_set_ui(val,1);
mpz_t_sentry val_s(val);
- bool hasmodulus = true;
-
- // first the modulus !
- if(va.modulus.size()>1)
- {
- biginteger modulus ;
- modulus.setValue(va.modulus[0].getValueTemp());
-
- for(unsigned int i = 1 ; i < va.modulus.size() ; ++i) {
- if(modulus != va.modulus[i]) { // if one is different: no modulus
- hasmodulus = false;
- break;
- }
- }
- if(hasmodulus)
- result.modulus.push_back(modulus);
-
- }
- else
- hasmodulus = false;
-
- if(va.modulus.size() == 1)
- {
- result.modulus.push_back(va.modulus[0]);
- hasmodulus = true;
- }
+ bool hasmodulus = va.getType() == TYPE_MODULUS::MODULUS_GLOBAL;
+
for(unsigned int i = 0 ; i < va.size(); ++i)
{
- if(va.value[i].isNA() )
+ if(va[i].isNA() )
{
return (bigintegerR::create_SEXP(result));
}
- mpz_mul(val,val,va.value[i].getValueTemp());
+ mpz_mul(val,val,va[i].getValue().getValueTemp());
if(hasmodulus)
- mpz_mod(val,val,va.modulus[0].getValueTemp() );
+ mpz_mod(val,val,va.getGlobalModulus()->getValueTemp() );
}
- result.value[0].setValue(val);
+ result[0].setValue(val);
+ if(hasmodulus) result[0].setModulus(va.getGlobalModulus());
return(bigintegerR::create_SEXP(result));
}
@@ -1183,18 +1131,18 @@ SEXP biginteger_powm(SEXP x, SEXP y, SEXP n)
bigvec vy = bigintegerR::create_bignum(y);
bigvec vn = bigintegerR::create_bignum(n);
- result.value.resize(vx.value.size());
-
- for (unsigned int i = 0 ; i < vx.value.size(); i++)
+ result.resize(vx.size());
+ result.setGlobalModulus(vn[0].getValuePtr());
+ for (unsigned int i = 0 ; i < vx.size(); i++)
{
- result.value[i].NA(false);
+ result[i].getValue().NA(false);
// check if n != 0
- if(mpz_sgn(vn.value[i % vn.value.size()].getValueTemp()) != 0)
- mpz_powm(result.value[i].getValue(),
- vx.value[i].getValueTemp(),
- vy.value[i % vy.value.size()].getValueTemp(),
- vn.value[i % vn.value.size()].getValueTemp());
+ if(mpz_sgn(vn[i % vn.size()].getValueTemp()) != 0)
+ mpz_powm(result[i].getMpValue(),
+ vx[i].getValue().getValueTemp(),
+ vy[i % vy.size()].getValueTemp(),
+ vn[i % vn.size()].getValueTemp());
}
@@ -1220,7 +1168,7 @@ SEXP bigI_frexp(SEXP x)
const char *nms[] = {"d", "exp", ""};
SEXP ans, d_R, exp_R;
bigvec vx = bigintegerR::create_bignum(x);
- int n = vx.value.size();
+ int n = vx.size();
PROTECT(ans = Rf_mkNamed(VECSXP, nms)); // =~= R list(d = . , exp= .)
d_R = Rf_allocVector(REALSXP, n); SET_VECTOR_ELT(ans, 0, d_R);
@@ -1229,11 +1177,14 @@ SEXP bigI_frexp(SEXP x)
int *exp_ = INTEGER(exp_R);
for (int i = 0 ; i < n; i++) {
signed long int ex;
- d_[i] = mpz_get_d_2exp(&ex, vx.value[i].getValueTemp());
- if(abs(ex) < INT_MAX)
+ d_[i] = mpz_get_d_2exp(&ex, vx[i].getValueTemp());
+ if(abs(ex) < INT_MAX){
exp_[i] = (int) ex;
- else
+ }
+ else{
+ vx.clear();
error(_("exponent too large to fit into an integer"));
+ }
}
UNPROTECT(1);
return ans;
@@ -1248,7 +1199,7 @@ SEXP biginteger_log2(SEXP x)
for (unsigned int i = 0; i < v.size(); ++i) {
signed long int ex;
- double di = mpz_get_d_2exp(&ex, v.value[i].getValueTemp());
+ double di = mpz_get_d_2exp(&ex, v[i].getValueTemp());
// xi = di * 2 ^ ex ==> log2(xi) = log2(di) + ex :
r[i] = log(di) / M_LN2 + (double)ex;
}
@@ -1264,7 +1215,7 @@ SEXP biginteger_log(SEXP x)
for (unsigned int i = 0; i < v.size(); ++i) {
signed long int ex;
- double di = mpz_get_d_2exp(&ex, v.value[i].getValueTemp());
+ double di = mpz_get_d_2exp(&ex, v[i].getValueTemp());
// xi = di * 2 ^ ex ==> log(xi) = log(di) + ex*log(2) :
r[i] = log(di) + M_LN2*(double)ex;
}
diff --git a/src/bigintegerR.h b/src/bigintegerR.h
index f5a8b073b34f2768df7b4edc742e0b7f84f9330d..35c2f495919a3af3cf71741c476c07008bb63f83 100644
--- a/src/bigintegerR.h
+++ b/src/bigintegerR.h
@@ -4,7 +4,7 @@
* \version 1
*
* \date Created: 2006
- * \date Last modified: Time-stamp: <2023-01-16 19:11:24 (antoine)>
+ * \date Last modified: Time-stamp: <2023-01-19 19:43:37 (antoine)>
*
*
* \note Licence: GPL (>= 2)
@@ -18,7 +18,7 @@
#ifndef T_BIGMOD_BINARY_OPERATION
#define T_BIGMOD_BINARY_OPERATION 1
/// A pointer to a binary operator for bigintegers
-typedef DefaultBigMod (*biginteger_binary_fn)(const bigmod&, const bigmod&);
+typedef bigmod (*biginteger_binary_fn)(const bigmod&, const bigmod&);
#endif
#ifndef T_BIGMOD_BINARY_OPERATION_LOGICAL
@@ -26,6 +26,7 @@ typedef DefaultBigMod (*biginteger_binary_fn)(const bigmod&, const bigmod&);
typedef bool (*biginteger_logical_binary_fn)(const biginteger&, const biginteger&);
#endif
+typedef const biginteger & (*bigmod_accessor_fn)(const bigmod&);
struct lockSexp {
lockSexp(const SEXP & value) {
@@ -61,7 +62,7 @@ namespace bigintegerR{
/**
* \brief export vector of biginteger to R value
*/
- SEXP create_SEXP(const std::vector<biginteger>& v);
+ SEXP create_SEXP(const bigvec& v,bigmod_accessor_fn fct, int size);
/**
* \brief export bigvec to R value
diff --git a/src/bigmod.cc b/src/bigmod.cc
index 99ae1dc606ec5193e42a0dffbdbd04ed8c4ed683..69a8aba39ef1abf5f06ef192cdf0168f04d5389e 100644
--- a/src/bigmod.cc
+++ b/src/bigmod.cc
@@ -12,21 +12,21 @@
// for GetOption etc:
#include <R.h>
#include <Rinternals.h>
-
+using namespace std;
// string representation of (.) wrt base 'b' :
std::string bigmod::str(int b) const
{
- if (value.isNA())
+ if (value->isNA())
return "NA";
std::string s; // sstream seems to collide with libgmp :-(
- if (!modulus.isNA())
+ if (!modulus->isNA())
s = "(";
- s += value.str(b);
- if (!modulus.isNA()) {
+ s += value->str(b);
+ if (!modulus->isNA()) {
s += " %% ";
- s += modulus.str(b);
+ s += modulus->str(b);
s += ")";
}
return s;
@@ -36,15 +36,15 @@ bigmod & bigmod::operator= (const bigmod& rhs)
{
if(this != &rhs)
{
- modulus.setValue( rhs.getModulus() );
- value.setValue(rhs.value );
+ modulus = make_shared<biginteger>(rhs.getModulus());
+ value = make_shared<biginteger>(rhs.getValue() );
}
return(*this);
}
bigmod bigmod::inv () const
{
- if(value.isNA() || modulus.isNA()) {
+ if(value->isNA() || modulus->isNA()) {
return bigmod();
}
@@ -59,7 +59,7 @@ bigmod bigmod::inv () const
return bigmod(); // return NA; was
}
- return bigmod(val, modulus);
+ return bigmod(std::make_shared<biginteger>(val), std::make_shared<biginteger>(modulus->getValueTemp()));
}
@@ -72,23 +72,25 @@ bool operator!=(const bigmod& rhs, const bigmod& lhs)
bool operator==(const bigmod& rhs, const bigmod& lhs)
{
+
if(rhs.getValue() != lhs.getValue())
return(false);
+ // if(rhs.getModulusPtr().get() == lhs.getModulusPtr().get()) return true; // same pointer
return(!(rhs.getModulus() != lhs.getModulus()));
}
-DefaultBigMod operator+(const bigmod& lhs, const bigmod& rhs)
+bigmod operator+(const bigmod& lhs, const bigmod& rhs)
{
return create_bigmod(lhs, rhs, mpz_add);
}
-DefaultBigMod operator-(const bigmod& lhs, const bigmod& rhs)
+bigmod operator-(const bigmod& lhs, const bigmod& rhs)
{
return create_bigmod(lhs, rhs, mpz_sub);
}
-DefaultBigMod operator*(const bigmod& lhs, const bigmod& rhs)
+bigmod operator*(const bigmod& lhs, const bigmod& rhs)
{
return create_bigmod(lhs, rhs, mpz_mul);
}
@@ -98,9 +100,9 @@ DefaultBigMod operator*(const bigmod& lhs, const bigmod& rhs)
* ~~~~~~~~~~~~~~
* itself called from "/.bigz" = div.bigz()
*/
-DefaultBigMod div_via_inv(const bigmod& a, const bigmod& b) {
+bigmod div_via_inv(const bigmod& a, const bigmod& b) {
// compute a/b as a * b^(-1)
- return operator*(a, pow(b, DefaultBigMod(-1)));
+ return operator*(a, pow(b, bigmod(-1)));
}
@@ -125,17 +127,17 @@ void integer_div(mpz_t result,const mpz_t a, const mpz_t b) {
/* called via biginteger_binary_operation(.) from R's
* .Call(biginteger_divq, a, b) , itself called from '%/%.bigz' = divq.bigz()
*/
-DefaultBigMod operator/(const bigmod& lhs, const bigmod& rhs) {
+bigmod operator/(const bigmod& lhs, const bigmod& rhs) {
return create_bigmod(lhs, rhs, integer_div, false);
}
-DefaultBigMod operator%(const bigmod& lhs, const bigmod& rhs)
+bigmod operator%(const bigmod& lhs, const bigmod& rhs)
{
if (lhs.getValue().isNA() || rhs.getValue().isNA())
- return DefaultBigMod();
+ return bigmod();
if (mpz_sgn(rhs.getValue().getValueTemp()) == 0) {
warning(_("biginteger division by zero: returning NA"));
- return DefaultBigMod();
+ return bigmod();
}
biginteger mod;
if (!lhs.getModulus().isNA() || !rhs.getModulus().isNA())
@@ -145,13 +147,13 @@ DefaultBigMod operator%(const bigmod& lhs, const bigmod& rhs)
mpz_init(val);
mpz_t_sentry val_s(val);
mpz_mod(val, lhs.getValue().getValueTemp(), rhs.getValue().getValueTemp());
- return DefaultBigMod(val, mod);
+ return bigmod(val, mod);
}
// Either 'base' has a modulus, or it has not *and* exp >= 0 :
-DefaultBigMod pow(const bigmod& base, const bigmod& exp)
+bigmod pow(const bigmod& base, const bigmod& exp)
{
biginteger mod = get_modulus(base, exp);
#ifdef DEBUG_bigmod
@@ -165,9 +167,9 @@ DefaultBigMod pow(const bigmod& base, const bigmod& exp)
if(mod.isNA() &&
((!base.getValue().isNA() && !mpz_cmp_si(base.getValue().getValueTemp(), 1)) ||
(! exp.getValue().isNA() && !mpz_cmp_si( exp.getValue().getValueTemp(), 0))))
- return DefaultBigMod(biginteger(1));
+ return bigmod(biginteger(1));
if (base.getValue().isNA() || exp.getValue().isNA())
- return DefaultBigMod();
+ return bigmod();
int sgn_exp = mpz_sgn(exp.getValue().getValueTemp());
bool neg_exp = (sgn_exp < 0); // b ^ -|e| = 1 / b^|e|
mpz_t val; mpz_init(val); mpz_t_sentry val_s(val);
@@ -177,7 +179,9 @@ DefaultBigMod pow(const bigmod& base, const bigmod& exp)
mod.str(10).c_str());
#endif
if (mod.isNA()) { // <==> (both have no mod || both have mod. but differing)
- if(neg_exp) error(_("** internal error (negative powers for Z/nZ), please report!"));
+ if(neg_exp){
+ error(_("** internal error (negative powers for Z/nZ), please report!"));
+ }
if (!mpz_fits_ulong_p(exp.getValue().getValueTemp()))
error(_("exponent e too large for pow(z,e) = z^e"));// FIXME? return( "Inf" )
// else :
@@ -190,7 +194,7 @@ DefaultBigMod pow(const bigmod& base, const bigmod& exp)
SEXP wOpt = Rf_GetOption1(Rf_install("gmp:warnNoInv"));
if(wOpt != R_NilValue && Rf_asInteger(wOpt))
warning(_("pow(x, -|n|) returning NA as x has no inverse wrt modulus"));
- return(DefaultBigMod()); // return NA; was
+ return(bigmod()); // return NA; was
} // else: val = x^(-1) already: ==> result = val ^ |exp| = val ^ (-exp) :
// nExp := - exp
mpz_t nExp; mpz_init(nExp); mpz_neg(nExp, exp.getValue().getValueTemp());
@@ -199,18 +203,18 @@ DefaultBigMod pow(const bigmod& base, const bigmod& exp)
mpz_powm(val, base.getValue().getValueTemp(), exp.getValue().getValueTemp(), mod.getValueTemp());
}
}
- return DefaultBigMod(val, mod);
+ return bigmod(val, mod);
}
-DefaultBigMod inv(const bigmod& x, const bigmod& m)
+bigmod inv(const bigmod& x, const bigmod& m)
{
if (x.getValue().isNA() || m.getValue().isNA())
- return DefaultBigMod();
+ return bigmod();
SEXP wOpt = Rf_GetOption1(Rf_install("gmp:warnNoInv"));
bool warnI = (wOpt != R_NilValue && Rf_asInteger(wOpt));
if (mpz_sgn(m.getValue().getValueTemp()) == 0) {
if(warnI) warning(_("inv(0) returning NA"));
- return DefaultBigMod();
+ return bigmod();
}
biginteger mod = get_modulus(x, m);
mpz_t val;
@@ -218,30 +222,30 @@ DefaultBigMod inv(const bigmod& x, const bigmod& m)
mpz_t_sentry val_s(val);
if (mpz_invert(val, x.getValue().getValueTemp(), m.getValue().getValueTemp()) == 0) {
if(warnI) warning(_("inv(x,m) returning NA as x has no inverse modulo m"));
- return(DefaultBigMod()); // return NA; was
+ return(bigmod()); // return NA; was
}
- return DefaultBigMod(val, mod);
+ return bigmod(val, mod);
}
// R as.bigz() :
-DefaultBigMod set_modulus(const bigmod& x, const bigmod& m)
+bigmod set_modulus(const bigmod& x, const bigmod& m)
{
if (!m.getValue().isNA() && mpz_sgn(m.getValue().getValueTemp()) == 0)
error(_("modulus 0 is invalid"));
// if (!m.getValue().isNA() && mpz_cmp(x.getValue().getValueTemp(),m.getValue().getValueTemp())>=0) {
if (!m.getValue().isNA() ) {
- DefaultBigMod t(x%m);
- return DefaultBigMod(t.getValue(), m.getValue());
+ bigmod t(x%m);
+ return bigmod(t.getValue(), m.getValue());
} else
- return DefaultBigMod(x.getValue(), m.getValue());
+ return bigmod(x.getValue(), m.getValue());
}
-DefaultBigMod gcd(const bigmod& lhs, const bigmod& rhs)
+bigmod gcd(const bigmod& lhs, const bigmod& rhs)
{
return create_bigmod(lhs, rhs, mpz_gcd);
}
-DefaultBigMod lcm(const bigmod& lhs, const bigmod& rhs)
+bigmod lcm(const bigmod& lhs, const bigmod& rhs)
{
return create_bigmod(lhs, rhs, mpz_lcm);
}
@@ -271,13 +275,13 @@ biginteger get_modulus(const bigmod& b1, const bigmod& b2)
// Create a bigmod from a binary combination of two other bigmods
-DefaultBigMod create_bigmod(const bigmod& lhs, const bigmod& rhs, gmp_binary f,
+bigmod create_bigmod(const bigmod& lhs, const bigmod& rhs, gmp_binary f,
bool zeroRhsAllowed) {
if (lhs.getValue().isNA() || rhs.getValue().isNA())
- return DefaultBigMod();
+ return bigmod();
if (!zeroRhsAllowed && mpz_sgn(rhs.getValue().getValueTemp()) == 0) {
warning(_("returning NA for (modulus) 0 in RHS"));
- return DefaultBigMod();
+ return bigmod();
}
biginteger mod = get_modulus(lhs, rhs);
mpz_t val;
@@ -313,5 +317,5 @@ DefaultBigMod create_bigmod(const bigmod& lhs, const bigmod& rhs, gmp_binary f,
}
#endif
}
- return DefaultBigMod(val, mod);
+ return bigmod(val, mod);
}
diff --git a/src/bigmod.h b/src/bigmod.h
index efbd2645048a4be8bed94f8f27b220d2145cdedb..3df7808917b51a80348840b8274478a6a66a5ce1 100644
--- a/src/bigmod.h
+++ b/src/bigmod.h
@@ -2,7 +2,7 @@
* \brief Description of class bigmod
*
* \date Created: 22/05/06
- * \date Last modified: Time-stamp: <2019-03-10 10:30:30 (antoine)>
+ * \date Last modified: Time-stamp: <2023-01-21 11:41:06 (antoine)>
*
* \author Immanuel Scholz
*
@@ -13,6 +13,7 @@
#ifndef BIGMOD_HEADER_
#define BIGMOD_HEADER_ 1
+#include <memory>
#include "biginteger.h"
typedef void (*gmp_binary)(mpz_t, const mpz_t, const mpz_t);
@@ -36,32 +37,25 @@ extern "C" {
*/
class bigmod {
private:
- /** optional source */
- biginteger * value_ptr;
- biginteger * modulus_ptr;
-
- protected:
- /** \brief Value of our bigmod -- only references*/
- biginteger & value;
+ /** \brief Value of our bigmod */
+ std::shared_ptr<biginteger> value;
/** \brief modulus of our bigmod representation: value %% modulus */
- biginteger & modulus;
+ std::shared_ptr<biginteger> modulus;
public:
- /** keep both references value / modulus
+ /** keep both pointers value / modulus
*/
- bigmod(biginteger& value_,
- biginteger& modulus_) :
- value_ptr(NULL),
- modulus_ptr(NULL),
- value(value_),modulus(modulus_) {};
+ bigmod(biginteger * value_,
+ biginteger * modulus_) :
+ value(value_),
+ modulus(modulus_) {};
/** keep references value / modulus is new object.
*/
- bigmod(biginteger& value_) :
- value_ptr(NULL),
- modulus_ptr(new biginteger()),
- value(value_),modulus(*modulus_ptr) {};
+ bigmod(biginteger* value_) :
+ value(value_),
+ modulus(std::make_shared<biginteger>()) {};
/**
@@ -69,33 +63,48 @@ class bigmod {
*/
bigmod(const biginteger& value_,
const biginteger& modulus_) :
- value_ptr(new biginteger(value_)),
- modulus_ptr(new biginteger(modulus_)),
- value(*value_ptr),modulus(*modulus_ptr) {};
+ value(std::make_shared<biginteger>(value_)),
+ modulus(std::make_shared<biginteger>(modulus_)) {};
+
bigmod(const biginteger& value_) :
- value_ptr(new biginteger(value_)),
- modulus_ptr(new biginteger()),
- value(*value_ptr),modulus(*modulus_ptr) {};
+ value(std::make_shared<biginteger>(value_)),
+ modulus(std::make_shared<biginteger>()) {};
bigmod() :
- value_ptr(new biginteger()),
- modulus_ptr(new biginteger()),
- value(*value_ptr),modulus(*modulus_ptr) {};
-
+ value(std::make_shared<biginteger>()),
+ modulus(std::make_shared<biginteger>()) {};
+
+ bigmod(std::shared_ptr<biginteger> value_, std::shared_ptr<biginteger> modulus_ ) :
+ value(),
+ modulus()
+ {
+ value = value_;
+ modulus = modulus_;
+ };
+
+ bigmod(std::shared_ptr<biginteger> value_ ) :
+ value(),
+ modulus(std::make_shared<biginteger>())
+ {
+ value = value_;
+ };
+
/** \brief copy operator */
- bigmod(const bigmod & rhs) :
- value_ptr(new biginteger(rhs.getValue())),
- modulus_ptr(new biginteger(rhs.getModulus())),
- value(*value_ptr),modulus(*modulus_ptr) {
+ bigmod(const bigmod & rhs) :
+ value(),
+ modulus()
+ {
+ value = rhs.value;
+ modulus = rhs.modulus;
};
+
virtual ~bigmod(){
- if(value_ptr != NULL) delete value_ptr;
- if(modulus_ptr != NULL) delete modulus_ptr;
+
};
/**
@@ -116,86 +125,78 @@ class bigmod {
bigmod inv () const;
- biginteger & getValue() {
- return value;
+ inline biginteger & getValue() {
+ return *value;
+ }
+
+ inline mpz_t & getMpValue()
+ {
+ return value->getValue();
+ }
+
+ inline const mpz_t& getValueTemp() const
+ {
+ return value->getValueTemp();
}
biginteger & getModulus() {
- return modulus;
+ return *modulus;
}
- const biginteger & getValue() const{
+
+ inline bool isNA(){
+ return value->isNA();
+ }
+
+ std::shared_ptr<biginteger> & getValuePtr() {
return value;
}
-
- const biginteger & getModulus() const {
+
+ const std::shared_ptr<biginteger> & getModulusPtr() const {
return modulus;
}
-};
+ inline void setValue(const std::shared_ptr<biginteger> & v){
+ value = v;
+ }
+
+ inline void setValue(const bigmod & v){
+ value = std::make_shared<biginteger>(v.getValue());
+ modulus = std::make_shared<biginteger>(v.getModulus());
+ }
+ inline void setValue(const mpz_t & v){
+ value->setValue(v);
+ }
-class DefaultBigMod : public bigmod {
- private:
- /** \brief Value of our bigmod */
- biginteger valueLocal;
- /** \brief modulus of our bigmod representation: value %% modulus */
- biginteger modulusLocal;
-
- public:
-/** \brief creator
- */
- DefaultBigMod(const biginteger& value_ = biginteger(),
- const biginteger& modulus_ = biginteger()) :
- bigmod(valueLocal,modulusLocal),
- valueLocal(value_),modulusLocal(modulus_) {
- value = valueLocal;
- modulus = modulusLocal;
-}
+ inline void setValue(int v){
+ value->setValue(v);
+ }
- /** \brief copy operator */
- DefaultBigMod(const bigmod & rhs) :
- bigmod(valueLocal,modulusLocal),
- valueLocal(rhs.getValue()),modulusLocal(rhs.getModulus()) {
- value = valueLocal;
- modulus = modulusLocal;
-}
+ inline void setValue(double v){
+ value->setValue(v);
+ }
- /** \brief copy operator */
- DefaultBigMod(const DefaultBigMod & rhs) :
- bigmod(valueLocal,modulusLocal),
- valueLocal(rhs.getValue()),modulusLocal(rhs.getModulus()) {
- value = valueLocal;
- modulus = modulusLocal;
-}
- ~DefaultBigMod(){};
-
+ inline void setValue(){
+ value->setValue();
+ modulus->setValue();
+ }
+
+ inline void setModulus(const std::shared_ptr<biginteger> & v){
+ modulus = v;
+ }
+ const biginteger & getValue() const{
+ return *value;
+ }
+
+ const biginteger & getModulus() const {
+ return *modulus;
+ }
};
-/**
- * a bigmod that has only integer.
- */
-class BigModInt : public bigmod {
- private:
- /** \brief modulus of our bigmod representation */
- biginteger modulusLocal;
-
- public:
-/** \brief creator
- */
- BigModInt(biginteger& value_) :
- bigmod(value_,modulusLocal),
- modulusLocal() {
- modulus = modulusLocal;
-}
-
- ~BigModInt(){};
-
-
-};
@@ -217,33 +218,33 @@ bool operator== (const bigmod& rhs, const bigmod& lhs);
* a modulus set. If none modulus for either bigmod is set, the result will not
* have a modulus as well.
*/
-DefaultBigMod operator+(const bigmod& rhs, const bigmod& lhs);
+bigmod operator+(const bigmod& rhs, const bigmod& lhs);
/**
* \brief Subtract two bigmods.
*
* For modulus description, see operator+(bigmod, bigmod)
*/
-DefaultBigMod operator-(const bigmod& rhs, const bigmod& lhs);
+bigmod operator-(const bigmod& rhs, const bigmod& lhs);
/**
* \brief Multiply two bigmods.
*
* For modulus description, see operator+(bigmod, bigmod)
*/
-DefaultBigMod operator*(const bigmod& rhs, const bigmod& lhs);
+bigmod operator*(const bigmod& rhs, const bigmod& lhs);
/**
* \brief Divide two bigmods a / b := a * b^(-1)
*/
-DefaultBigMod div_via_inv(const bigmod& a, const bigmod& b);
+bigmod div_via_inv(const bigmod& a, const bigmod& b);
/**
* \brief Divide two bigmods.
*
* For modulus description, see operator+(bigmod, bigmod)
*/
-DefaultBigMod operator/(const bigmod& rhs, const bigmod& lhs);
+bigmod operator/(const bigmod& rhs, const bigmod& lhs);
/**
* \brief Calculate the modulus (remainder) of two bigmods.
@@ -253,7 +254,7 @@ DefaultBigMod operator/(const bigmod& rhs, const bigmod& lhs);
* was before, except if rhs and lhs has both no modulus set,
* in which case the resulting modulus will be unset too.
*/
-DefaultBigMod operator%(const bigmod& rhs, const bigmod& lhs);
+bigmod operator%(const bigmod& rhs, const bigmod& lhs);
/**
* \brief Return the power of "exp" to the base of "base" (return = base^exp).
@@ -265,14 +266,14 @@ DefaultBigMod operator%(const bigmod& rhs, const bigmod& lhs);
*
* For other modulus description, see operator+(bigmod, bigmod)
*/
-DefaultBigMod pow(const bigmod& base, const bigmod& exp);
+bigmod pow(const bigmod& base, const bigmod& exp);
/**
* \brief Return the modulo inverse to x mod m. (return = x^-1 % m)
*
* For modulus description, see operator+(bigmod, bigmod)
*/
-DefaultBigMod inv(const bigmod& x, const bigmod& m);
+bigmod inv(const bigmod& x, const bigmod& m);
/**
* \brief Return a bigmod with value (x % m) and the intern modulus set to m.
@@ -280,7 +281,7 @@ DefaultBigMod inv(const bigmod& x, const bigmod& m);
*
* Do not confuse this with operator%(bigmod, bigmod).
*/
-DefaultBigMod set_modulus(const bigmod& x, const bigmod& m);
+bigmod set_modulus(const bigmod& x, const bigmod& m);
biginteger get_modulus(const bigmod& b1, const bigmod& b2);
@@ -289,20 +290,20 @@ biginteger get_modulus(const bigmod& b1, const bigmod& b2);
*
* For modulus description, see operator+(bigmod, bigmod)
*/
-DefaultBigMod gcd(const bigmod& rhs, const bigmod& lhs);
+bigmod gcd(const bigmod& rhs, const bigmod& lhs);
/**
* \brief Return the least common multiply of both parameter.
*
* For modulus description, see operator+(bigmod, bigmod)
*/
-DefaultBigMod lcm(const bigmod& rhs, const bigmod& lhs);
+bigmod lcm(const bigmod& rhs, const bigmod& lhs);
/**
* \brief function used to make any binary operation between
* two bigmod that return a bigmod (addition substraction... )
*/
-DefaultBigMod create_bigmod(const bigmod& lhs, const bigmod& rhs, gmp_binary f,
+bigmod create_bigmod(const bigmod& lhs, const bigmod& rhs, gmp_binary f,
bool zeroRhsAllowed = true) ;
#endif
diff --git a/src/bigrational.h b/src/bigrational.h
index 2264f006e7e9ef1614f55288df4608f540bb0c66..0575afc964b0e8487a35ca1d8f5db48963dc28ef 100644
--- a/src/bigrational.h
+++ b/src/bigrational.h
@@ -2,7 +2,7 @@
* \brief Description of class bigrational
*
* \date Created: 22/05/06
- * \date Last modified: Time-stamp: <2009-10-27 17:28:06 antoine>
+ * \date Last modified: Time-stamp: <2023-01-21 13:13:10 (antoine)>
*
* \author Antoine Lucas (adapted from biginteger class made by
* Immanuel Scholz)
@@ -134,6 +134,17 @@ class bigrational
Rf_error("Not a valid number"); */
}
+ /**
+ * Copy constructor (mpz_t aren't standard-copyable)
+ */
+ bigrational(const biginteger & rhs) :
+ value(),
+ na(rhs.isNA())
+ {
+ mpq_init(value);
+ mpq_set_z(value, rhs.getValueTemp());
+ }
+
/**
* Copy constructor (mpz_t aren't standard-copyable)
*/
diff --git a/src/bigrationalR.cc b/src/bigrationalR.cc
index d776403cd5a0b364dc1f5ca0a151a599e6f08269..0f36ef0ce9ea5d06d1ee33c480a79f7bb7d7f6a7 100644
--- a/src/bigrationalR.cc
+++ b/src/bigrationalR.cc
@@ -127,7 +127,7 @@ namespace bigrationalR
if (TYPEOF(denAttr) != NILSXP)
{
bigvec_q attrib = bigrationalR::create_vector(denAttr);
- if (!attrib.value.empty()) // sanity check
+ if (!attrib.size() == 0) // sanity check
for (unsigned int i = 0; i < v.size(); ++i)
{
if( attrib[i%attrib.size()].sgn() != 0)
@@ -212,14 +212,23 @@ namespace bigrationalR
{
bigvec_q va = bigrationalR::create_bignum(a);
bigvec_q vb = bigrationalR::create_bignum(b), result;
+
+ int nrow = matrixz::checkDims(va.nrow,vb.nrow) ;
+ if(nrow == -2){
+ va.clear();
+ vb.clear();
+ error(_("Matrix dimensions do not match"));
+ }
+
// MM: I think, 0 length vector operations should work too!
// if (va.value.empty() || vb.value.empty())
// error(_("argument must not be an empty list"));
- int size = (va.value.empty() || vb.value.empty()) ? 0 : max(va.size(), vb.size());
+ int size = (va.size() ==0 || vb.size() == 0) ? 0 : max(va.size(), vb.size());
result.value.reserve(size);
for (int i = 0; i < size; ++i)
result.push_back(f(va.value[i%va.size()], vb.value[i%vb.size()]));
- result.nrow = matrixz::checkDims(va.nrow,vb.nrow) ;
+ result.nrow = nrow;
+
return bigrationalR::create_SEXP(result);
}
@@ -227,12 +236,20 @@ namespace bigrationalR
{
bigvec_q va = bigrationalR::create_bignum(a), result;
bigvec vb = bigintegerR::create_bignum(b);
- int size = (va.value.empty() || vb.value.empty()) ? 0 : max(va.size(), vb.size());
- result.value.reserve(size);
- for (int i = 0; i < size; ++i)
- result.push_back(f(va.value[i%va.size()], vb.value[i%vb.size()]));
- result.nrow = matrixz::checkDims(va.nrow,vb.nrow) ;
- return bigrationalR::create_SEXP(result);
+ int size = (va.size()==0 || vb.size()==0) ? 0 : max(va.size(), vb.size());
+
+ int nrow = matrixz::checkDims(va.nrow,vb.nrow) ;
+ if(nrow == -2){
+ va.clear();
+ vb.clear();
+ error(_("Matrix dimensions do not match"));
+ }
+
+ for (int i = 0; i < size; ++i)
+ result.push_back(f(va.value[i%va.size()], vb[i%vb.size()].getValue()));
+ result.nrow = nrow;
+
+ return bigrationalR::create_SEXP(result);
}
SEXP bigrational_logical_binary_operation(SEXP a, SEXP b, bigrational_logical_binary_fn f)
@@ -240,9 +257,17 @@ namespace bigrationalR
bigvec_q va = bigrationalR::create_bignum(a);
bigvec_q vb = bigrationalR::create_bignum(b), result;
// MM: I think, 0 length vector operations should work too!
- // if (va.value.empty() || vb.value.empty())
+ // if (va.size() == 0 || vb.size() == 0)
// error(_("argument must not be an empty list"));
- int size = (va.value.empty() || vb.value.empty()) ? 0 : max(va.size(), vb.size());
+
+ int nrow = matrixz::checkDims(va.nrow,vb.nrow) ;
+ if(nrow == -2){
+ va.clear();
+ vb.clear();
+ error(_("Matrix dimensions do not match"));
+ }
+
+ int size = (va.size() == 0 || vb.size() == 0) ? 0 : max(va.size(), vb.size());
SEXP ans = PROTECT(Rf_allocVector(LGLSXP, size));
for (int i = 0; i < size; ++i) {
bigrational am = va.value[i % va.size()];
@@ -253,8 +278,6 @@ namespace bigrationalR
LOGICAL(ans)[i] = f(va[i%va.size()], vb[i%vb.size()]) ? 1 : 0;
}
- int nrow = matrixz::checkDims(va.nrow,vb.nrow) ;
-
// Add dimension parameter when available
if(nrow >= 0)
{
@@ -484,11 +507,11 @@ SEXP bigrational_den(SEXP a)
mpz_init(z_tmp);
bigvec_q v =bigrationalR::create_bignum(a);
bigvec result;
- result.value.resize(v.size());
+ result.resize(v.size());
for (unsigned int i = 0; i < v.size(); ++i) {
- mpq_get_den(z_tmp,v.value[i].getValueTemp());
- result.value[i].setValue(z_tmp);
+ mpq_get_den(z_tmp,v[i].getValueTemp());
+ result[i].setValue(z_tmp);
}
mpz_clear(z_tmp);
return bigintegerR::create_SEXP(result);
@@ -503,9 +526,9 @@ SEXP bigrational_num(SEXP a)
result.resize(v.size());
for (unsigned int i = 0; i < v.size(); ++i) {
- if(!v.value[i].isNA()) {
- mpq_get_num(z_tmp,v.value[i].getValueTemp());
- result.value[i].setValue(z_tmp);
+ if(!v[i].isNA()) {
+ mpq_get_num(z_tmp,v[i].getValueTemp());
+ result[i].setValue(z_tmp);
} // else: uninitialized, i.e., NA
}
mpz_clear(z_tmp);
@@ -543,7 +566,7 @@ SEXP bigrational_setlength(SEXP vec, SEXP value)
error(_("invalid second argument"));
}
bigvec_q v =bigrationalR::create_bignum(vec);
- v.value.resize(len);
+ v.resize(len);
return bigrationalR::create_SEXP(v);
}
@@ -553,7 +576,7 @@ SEXP bigrational_is_na(SEXP a)
SEXP ans = PROTECT(Rf_allocVector(LGLSXP, v.size()));
int *a_ = LOGICAL(ans);
for (unsigned int i = 0; i < v.size(); ++i)
- a_[i] = v.value[i].isNA();
+ a_[i] = v[i].isNA();
UNPROTECT(1);
return ans;
}
@@ -568,7 +591,7 @@ SEXP bigrational_is_int(SEXP a)
mpz_init(z_tmp);
for (unsigned int i = 0; i < v.size(); ++i) {
- mpq_get_den(z_tmp,v.value[i].getValueTemp());
+ mpq_get_den(z_tmp,v[i].getValueTemp());
a_[i] = mpz_cmp_ui (z_tmp, 1) == 0; // <==> numerator == 1
}
mpz_clear(z_tmp);
@@ -584,7 +607,7 @@ SEXP bigrational_c(SEXP args)
for(int i = 0; i < Rf_length(args); i++) {
bigvec_q v = bigrationalR::create_bignum(VECTOR_ELT(args,i));
for(unsigned int j=0; j < v.size(); j++)
- result.push_back(v.value[j]);
+ result.push_back(v[j]);
v.value.clear();
}
@@ -600,7 +623,7 @@ SEXP bigrational_rep(SEXP x, SEXP times)
result.value.reserve(v.size()*rep);
for(i = 0 ; i< rep ; i++)
for(j = 0 ; j < v.size() ; j++)
- result.push_back(v.value[j]);
+ result.push_back(v[j]);
return bigrationalR::create_SEXP(result);
}
@@ -619,14 +642,14 @@ SEXP bigrational_max(SEXP a, SEXP narm)
for(unsigned int i = 1 ; i < va.size(); ++i)
{
- if(va.value[i].isNA() && ! na_remove)
+ if(va[i].isNA() && ! na_remove)
return(bigrationalR::create_SEXP(result));
else
- if(!(va.value[i] < va.value[maximum] ))
+ if(!(va[i] < va[maximum] ))
maximum = i; // if va.value[maximum = 0] is NA => false for the "<" => maximum changed = good
}
- result.push_back(va.value[maximum]);
+ result.push_back(va[maximum]);
return bigrationalR::create_SEXP(result);
}
@@ -645,14 +668,14 @@ SEXP bigrational_min(SEXP a, SEXP narm)
for(unsigned int i = 1 ; i < va.size(); ++i)
{
- if(va.value[i].isNA() && !na_remove)
+ if(va[i].isNA() && !na_remove)
return(bigrationalR::create_SEXP(result));
else
- if(!(va.value[i] > va.value[minimum] ))
+ if(!(va[i] > va[minimum] ))
minimum = i; // if va.value[maximum = 0] is NA => false for the "<" => maximum changed = good
}
- result.push_back(va.value[minimum]);
+ result.push_back(va[minimum]);
return bigrationalR::create_SEXP(result);
}
@@ -661,19 +684,19 @@ SEXP bigrational_min(SEXP a, SEXP narm)
SEXP bigrational_cumsum(SEXP a)
{
bigvec_q result, va = bigrationalR::create_bignum(a);
- result.value.resize(va.value.size());
+ result.resize(va.size());
mpq_t val;
mpq_init(val);
mpq_t_sentry val_s(val);
for(unsigned int i = 0 ; i < va.size(); ++i) {
- if(va.value[i].isNA() ) {
+ if(va[i].isNA() ) {
break; // all last values are NA.
}
- mpq_add(val,val,va.value[i].getValueTemp());
+ mpq_add(val,val,va[i].getValueTemp());
- result.value[i].setValue(val);
+ result[i].setValue(val);
}
return(bigrationalR::create_SEXP(result));
}
@@ -683,19 +706,19 @@ SEXP bigrational_cumsum(SEXP a)
SEXP bigrational_sum(SEXP a)
{
bigvec_q result, va = bigrationalR::create_bignum(a);
- result.value.resize(1);
+ result.resize(1);
mpq_t val;
mpq_init(val);
mpq_t_sentry val_s(val);
for(unsigned int i = 0 ; i < va.size(); ++i) {
- if(va.value[i].isNA()) {
+ if(va[i].isNA()) {
break; // all last values are NA.
}
- mpq_add(val,val,va.value[i].getValueTemp());
+ mpq_add(val,val,va[i].getValueTemp());
}
- result.value[0].setValue(val);
+ result[0].setValue(val);
return(bigrationalR::create_SEXP(result));
}
@@ -705,7 +728,7 @@ SEXP bigrational_sum(SEXP a)
SEXP bigrational_prod(SEXP a)
{
bigvec_q result, va = bigrationalR::create_bignum(a);
- result.value.resize(1);
+ result.resize(1);
mpq_t val;
mpq_init(val);
@@ -713,12 +736,12 @@ SEXP bigrational_prod(SEXP a)
mpq_t_sentry val_s(val);
for(unsigned int i = 0 ; i < va.size(); ++i) {
- if(va.value[i].isNA() ) {
+ if(va[i].isNA() ) {
return (bigrationalR::create_SEXP(result));
}
- mpq_mul(val,val,va.value[i].getValueTemp());
+ mpq_mul(val,val,va[i].getValueTemp());
}
- result.value[0].setValue(val);
+ result[0].setValue(val);
return(bigrationalR::create_SEXP(result));
}
@@ -728,33 +751,33 @@ SEXP bigrational_R_pow(SEXP x, SEXP y)
{
bigvec_q result, vx = bigrationalR::create_bignum(x);
bigvec vy = bigintegerR::create_bignum(y);
- int size = (vx.value.empty() || vy.value.empty()) ? 0 : max(vx.size(), vy.size());
+ int size = (vx.size() ==0 || vy.size() == 0) ? 0 : max(vx.size(), vy.size());
mpq_t val; mpq_init(val); mpq_t_sentry val_s(val);
mpz_t num, den; mpz_init(num); mpz_init(den);
mpz_t_sentry val_n(num); mpz_t_sentry val_d(den);
- result.value.resize(size);
+ result.resize(size);
for (int i = 0 ; i < size; i++) {
- int i_x = i % vx.value.size(),
- i_y = i % vy.value.size();
+ int i_x = i % vx.size(),
+ i_y = i % vy.size();
//fails: val.NA(false);
- if(vx.value[i_x].isNA() ||
- vy.value[i_y].isNA())
+ if(vx[i_x].isNA() ||
+ vy[i_y].isNA())
break;
// else -- res[i] := x[i] ^ y[i] = (num ^ y_i) / (den ^ y_i) ----
- if(mpz_sgn(vy.value[i_y].getValueTemp()) < 0)
+ if(mpz_sgn(vy[i_y].getValueTemp()) < 0)
error(_("Negative powers not yet implemented [i = %d]"), i_y +1);
- if (!mpz_fits_ulong_p(vy.value[i_y].getValueTemp()))
+ if (!mpz_fits_ulong_p(vy[i_y].getValueTemp()))
error(_("exponent too large for pow [i = %d]"), i_y +1);
- int y_i = mpz_get_ui(vy.value[i_y].getValueTemp());
- // *num = mpq_numref(vx.value[i].getValueTemp());
- // *den = mpq_denref(vx.value[i].getValueTemp());
- mpq_get_num(num, vx.value[i_x].getValueTemp());
- mpq_get_den(den, vx.value[i_x].getValueTemp());
+ int y_i = mpz_get_ui(vy[i_y].getValueTemp());
+ // *num = mpq_numref(vx[i].getValueTemp());
+ // *den = mpq_denref(vx[i].getValueTemp());
+ mpq_get_num(num, vx[i_x].getValueTemp());
+ mpq_get_den(den, vx[i_x].getValueTemp());
mpz_pow_ui(num, num, y_i); // num := num ^ y_i
mpz_pow_ui(den, den, y_i); // den := den ^ y_i
@@ -764,7 +787,7 @@ SEXP bigrational_R_pow(SEXP x, SEXP y)
/* Simplify numerator and denominator */
mpq_canonicalize(val);
- result.value[i].setValue(val);
+ result[i].setValue(val);
}
return bigrationalR::create_SEXP(result);
diff --git a/src/bigvec.cc b/src/bigvec.cc
index 46930e1f85f738c7c8a20b4efc2c6bc11f9ee7fa..2a4b8544369b5be7c4e51fc38bdda4647240b474 100644
--- a/src/bigvec.cc
+++ b/src/bigvec.cc
@@ -1,69 +1,57 @@
#include "bigvec.h"
+static int count = 0;
+static int countAll = 0;
+
/** \brief constructor
*
*/
bigvec::bigvec(unsigned int size) :
math::Matrix<bigmod>(),
- value(),
- modulus(),
- valuesMod(),
- nrow(-1)
-{
-
+ values(),
+ type(NO_MODULUS),
+ nrow(-1),
+ modulus()
+{
+ count++;
+ countAll++;
+
for (int i=0 ; i < size; i++){
- value.push_back(biginteger());
+ values.push_back(bigmod());
}
}
bigvec::bigvec(const bigvec & vecteur) :
math::Matrix<bigmod>(),
- value(),
- modulus(),
- valuesMod(),
- nrow(vecteur.nrow)
-{
- // *this = vecteur;
- std::vector<biginteger>::const_iterator it;
-
- for(it = vecteur.modulus.begin();it != vecteur.modulus.end(); ++it)
+ values(),
+ type(vecteur.type),
+ nrow(vecteur.nrow),
+ modulus(vecteur.modulus)
+{
+ count++;
+ countAll++;
+ // *this = vecteur;
+ values.reserve(vecteur.size());
+ for(std::vector<bigmod>::const_iterator it= vecteur.values.begin(); it != vecteur.values.end(); ++it)
{
- modulus.push_back(*it);
- }
-
- for(it= vecteur.value.begin(); it != vecteur.value.end(); ++it)
- {
- value.push_back(*it);
+ values.push_back(*it);
}
}
bigvec::~bigvec(){
- clearValuesMod();
+ count--;
+ // printf("bv %d %d \n",count, countAll);
+
+
}
//
std::string bigvec::str(int i,int b) const
{
- if (value[i].isNA())
- return ("NA");
-
- std::string s; // sstream seems to collide with libgmp :-(
- bool onemodulus = modulus.size()>0;
- if(onemodulus)
- onemodulus = !modulus[i%modulus.size()].isNA() ;
- if (onemodulus)
- s = "(";
- s += value[i].str(b);
- if (onemodulus) {
- s += " %% ";
- s += modulus[i%modulus.size()].str(b);
- s += ")";
- }
-
- return s;
+ return values[i].str(b);
}
bigmod & bigvec::get(unsigned int row, unsigned int col) {
@@ -73,116 +61,47 @@ bigmod & bigvec::get(unsigned int row, unsigned int col) {
bigmod & bigvec::operator[] (unsigned int i)
{
- checkValuesMod();
- return *valuesMod[i];
+ return values[i];
}
const bigmod & bigvec::operator[] (unsigned int i) const
{
- bigvec * nonConst = const_cast<bigvec*>(this);
- nonConst->checkValuesMod();
- return *valuesMod[i];
+ return values[i];
}
void bigvec::set(unsigned int row, unsigned int col, const bigmod & val) {
set( row + col*nRows(),val);
}
-void bigvec::checkValuesMod() {
- if (valuesMod.size() != value.size()){
- // reconstruct bigmod that are references to values and modulus:
- clearValuesMod();
- if(modulus.size()>0) {
- for (unsigned int i = 0 ; i < value.size(); i++)
- valuesMod.push_back(new bigmod(value[i], modulus[i%modulus.size()]));
- } else {
- for (unsigned int i= 0 ; i < value.size(); i++)
- valuesMod.push_back(new BigModInt(value[i]));
- }
-
- }
-}
-
-void bigvec::clearValuesMod(){
- for (unsigned int i = 0 ; i < valuesMod.size(); i++){
- delete valuesMod[i];
- }
- valuesMod.clear();
-}
void bigvec::set(unsigned int i,const bigmod & val)
{
- value[i] = val.getValue();
+ values[i] = val;
- if(!val.getModulus().isNA())
- {
- if(modulus.size()> i)
- {
- modulus[i] = val.getModulus() ;
- return;
- }
-
- unsigned int nrow_mod = nrow;
- if(nrow<1)
- nrow_mod = 1;
- if( (modulus.size() == (unsigned int)nrow_mod ) || (modulus.size() == 1) )
- {
- // check "row-mod" or "global" modulus
- if(!(val.getModulus() != modulus[i % modulus.size()])) {
- return;
- }
- }
-
- // we add "previous"
- nrow_mod = modulus.size();
- for(unsigned int j = nrow_mod; j< i;++j)
- modulus.push_back(modulus[j % nrow_mod]);
- modulus.push_back(val.getModulus());
- clearValuesMod();
+
+ if(type == NO_MODULUS){
+ if(val.getModulus().isNA()) return;
+ if(i == 0 && values.size() == 1) {
+ type = MODULUS_GLOBAL;
+ modulus = val.getModulusPtr();
+ } else {
+ type = MODULUS_BY_CELL;
}
+ }
+ if(type == MODULUS_GLOBAL){
+ if(values.size() == 1){
+ modulus = val.getModulusPtr();
+ } else if(val.getModulus() != *modulus ){
+ type = MODULUS_BY_CELL;
+ }
+ }
}
void bigvec::push_back(const bigmod & number)
{
- unsigned int nrow_mod = (nrow < 0) ? 1 : nrow;
- clearValuesMod();
-
- value.push_back(number.getValue());
-
- if((!number.getModulus().isNA()) || (modulus.size()>0) )
- {
- // pathologic case: value.size >0 ; modulus.size =0
- // we assume previous mod where NA
- if((modulus.size() ==0) && (value.size()>0))
- {
- modulus.resize(value.size()-1);
- modulus.push_back( number.getModulus());
- return;
- }
-
- // standard cas
- if((modulus.size() != 1 ) && (modulus.size() != nrow_mod))
- {
- modulus.push_back(number.getModulus());
- return;
- }
-
- // pathologic case:
- // value modulus [nrow=4]
- // 2 2 push_back: ok
- // 2 2 push_back: nothing
- // 2 1 push_back: shoud add previous modulus then 1
- // note nrow_mod is either 1 ither nrow (when nrow >1)
- nrow_mod = modulus.size();
- if(modulus[(value.size() -1)% nrow_mod ] != number.getModulus())
- {
- // we add "previous"
- for(unsigned int i = nrow_mod; i< value.size()-1;i++)
- modulus.push_back(modulus[i % nrow_mod]);
- modulus.push_back(number.getModulus());
- }
- }
+ values.push_back(bigmod());
+ set(values.size()-1, number);
}
/**
@@ -190,8 +109,7 @@ void bigvec::push_back(const bigmod & number)
*/
void bigvec::push_back(int value_p)
{
- clearValuesMod();
- value.push_back(biginteger(value_p));
+ push_back(biginteger(value_p));
}
/**
@@ -199,8 +117,7 @@ void bigvec::push_back(int value_p)
*/
void bigvec::push_back(biginteger & value_p)
{
- clearValuesMod();
- value.push_back(value_p);
+ push_back(bigmod(value_p));
}
/**
@@ -208,15 +125,14 @@ void bigvec::push_back(biginteger & value_p)
*/
void bigvec::push_back(const __mpz_struct * value_p)
{
- clearValuesMod();
- value.push_back(biginteger(value_p));
+ push_back(biginteger(value_p));
}
// return size of value
unsigned int bigvec::size() const
{
- return(value.size());
+ return(values.size());
}
@@ -228,20 +144,16 @@ unsigned int bigvec::nRows() const {
// hummm. to avoid !
void bigvec::resize(unsigned int i)
{
- clearValuesMod();
-
- value.resize(i);
- if(i < modulus.size())
- modulus.resize(i);
+ values.resize(i);
}
// clear all
void bigvec::clear()
{
- clearValuesMod();
- value.clear();
- modulus.clear();
+ values.clear();
+ type=NO_MODULUS;
+ modulus=nullptr;
nrow = -1;
}
@@ -251,22 +163,16 @@ bigvec & bigvec::operator= (const bigvec & rhs)
{
if(this != &rhs)
{
- clearValuesMod();
- value.resize(0);
- modulus.resize(0);
- std::vector<biginteger>::const_iterator it;
-
- for(it = rhs.modulus.begin();it != rhs.modulus.end(); ++it)
- {
- modulus.push_back(*it);
- }
+ values.resize(0);
+ modulus = rhs.modulus;
+ type=rhs.type;
- for(it= rhs.value.begin(); it != rhs.value.end(); ++it)
- {
- value.push_back(*it);
- }
+ for (unsigned int i = 0 ; i < rhs.size(); i++){
+ values.push_back(rhs[i]);
+ }
nrow = rhs.nrow;
+
}
return(*this);
}
@@ -276,24 +182,21 @@ bigvec & bigvec::operator= (const bigvec & rhs)
// Comparison operator
bool operator!=(const bigvec & rhs, const bigvec& lhs)
{
- std::vector<biginteger>::const_iterator it = rhs.value.begin();
- std::vector<biginteger>::const_iterator jt = lhs.value.begin();
- if( (rhs.value.size() != lhs.value.size()) || \
+
+ if( (rhs.size() != lhs.size()) || \
(rhs.nrow != lhs.nrow ) )
return(false);
+
+
// check value
- while(it != rhs.value.end())
+ for (unsigned int i = 0 ; i < rhs.size(); i++)
{
- if(*it != *jt)
+ if(rhs[i] != lhs[i])
return(false);
- it++; jt++;
}
- for(unsigned int i=0; i < std::max( rhs.modulus.size() ,lhs.modulus.size() ); ++i)
- if(rhs.modulus[i % rhs.modulus.size()] != \
- lhs.modulus[i % lhs.modulus.size()] )
- return(false);
+
return(true);
}
@@ -305,14 +208,34 @@ void bigvec::print()
if(nrow > 0) {
for(int i=0; i < nrow; ++i)
{
- for(unsigned int j=0; j < (value.size() / nrow); ++j)
- Rprintf("%s\t", value[i+j* nrow].str(10).c_str() );
+ for(unsigned int j=0; j < (values.size() / nrow); ++j)
+ Rprintf("%s\t", values[i+j* nrow].str(10).c_str() );
Rprintf("\n");
}
}
else {
- for(unsigned int i=0; i < value.size(); ++i)
- Rprintf("%s\t", value[i].str(10).c_str() );
+ for(unsigned int i=0; i < values.size(); ++i)
+ Rprintf("%s\t", values[i].str(10).c_str() );
Rprintf("\n");
}
}
+
+
+void bigvec::setGlobalModulus(std::shared_ptr<biginteger> & val){
+ modulus = val;
+ type = MODULUS_GLOBAL;
+ for (int i = 0 ; i < values.size(); i++){
+ values[i].setModulus(val);
+ }
+}
+
+
+std::shared_ptr<biginteger> bigvec::getGlobalModulus(bigvec & first, bigvec & second){
+ std::shared_ptr<biginteger> empty(nullptr);
+ if(first.getType() == MODULUS_GLOBAL && second.getType() == NO_MODULUS) return first.getGlobalModulus();
+ if(first.getType() == NO_MODULUS && second.getType() == MODULUS_GLOBAL) return second.getGlobalModulus();
+ if(first.getType() == MODULUS_GLOBAL && second.getType() == MODULUS_GLOBAL) {
+ return (*first.getGlobalModulus()) == (*second.getGlobalModulus()) ? first.getGlobalModulus() : empty;
+ }
+ return empty;
+}
diff --git a/src/bigvec.h b/src/bigvec.h
index b0f16e7b0a3f598a83f575b8ad72d87c478e2212..8aab56a77db50fad2d9a21d56421e476108db827 100644
--- a/src/bigvec.h
+++ b/src/bigvec.h
@@ -4,7 +4,7 @@
* \version 1
*
* \date Created: 2005
- * \date Last modified: Time-stamp: <2023-01-16 19:09:32 (antoine)>
+ * \date Last modified: Time-stamp: <2023-01-21 11:22:09 (antoine)>
*
*
* \note Licence: GPL (>= 2)
@@ -19,6 +19,12 @@
#include "bigmod.h"
#include "templateMatrix.h"
+enum TYPE_MODULUS{
+ NO_MODULUS,
+ MODULUS_GLOBAL,
+ MODULUS_BY_CELL
+};
+
/** \brief class bigvec
*
* It a class composed of 2 vectors, (value & modulus) that
@@ -26,15 +32,14 @@
* parameter (for matrix support)
*/
class bigvec : public math::Matrix<bigmod> {
- public:
- /** \brief value */
- std::vector<biginteger> value;
- /** \brief modulus */
- std::vector<biginteger> modulus;
+protected:
/** array with all bigmod, that are references to values in vector. */
- std::vector<bigmod *> valuesMod;
+ std::vector<bigmod > values;
+ TYPE_MODULUS type;
+ std::shared_ptr<biginteger> modulus;
+public:
/** \brief optional parameter used with matrix -- set to -1 for non-matrix */
int nrow ;
@@ -89,6 +94,10 @@ class bigvec : public math::Matrix<bigmod> {
*/
void push_back(int value_p);
+ inline void erase(unsigned int index){
+ values.erase(values.begin() + index);
+ }
+
/**
* Insert Big Integer value
*/
@@ -128,8 +137,33 @@ class bigvec : public math::Matrix<bigmod> {
*/
void print();
+ inline unsigned int getModulusSize(){
+ if(type == NO_MODULUS) return 0;
+ if(type == MODULUS_GLOBAL) return 1;
+ return size();
+ }
+ void setGlobalModulus(std::shared_ptr<biginteger> & modulus);
+ inline std::shared_ptr<biginteger> & getGlobalModulus(){
+ return modulus;
+ }
+
+ const TYPE_MODULUS getType() const {
+ return type;
+ }
+
+ inline void setType(TYPE_MODULUS val){
+ type = val;
+ if(val == MODULUS_GLOBAL && size() > 0){
+ modulus = (*this)[0].getModulusPtr();
+ }
+ }
+
+
+ /** return a global modulus if possible, null if not */
+ static std::shared_ptr<biginteger> getGlobalModulus(bigvec & first, bigvec & second);
+
private:
void checkValuesMod() ;
void clearValuesMod() ;
diff --git a/src/bigvec_q.cc b/src/bigvec_q.cc
index 51c4480db00f41e18472da84fb0b5e32dd9e5769..551b639620c24815c71cba2ff75d731a11a7a5b0 100644
--- a/src/bigvec_q.cc
+++ b/src/bigvec_q.cc
@@ -11,12 +11,12 @@ bigvec_q::bigvec_q(const bigvec_q & rhs):
bigvec_q::bigvec_q(const bigvec & rhs):
- value(rhs.value.size()),
+ value(rhs.size()),
nrow(rhs.nrow)
{
for(unsigned int i=0; i< rhs.size(); ++i)
{
- value[i].setValue(rhs.value[i]);
+ value[i].setValue(rhs[i].getValue());
}
}
diff --git a/src/bigvec_q.h b/src/bigvec_q.h
index d33df00210b8326dd37c101cbcbb3eb613c3f8ed..2feacb77ceb20185e662a6efe199eb952911fb3f 100644
--- a/src/bigvec_q.h
+++ b/src/bigvec_q.h
@@ -4,7 +4,7 @@
* \version 1
*
* \date Created: 2005
- * \date Last modified: Time-stamp: <2006-06-16 20:19:01 antoine>
+ * \date Last modified: Time-stamp: <2023-01-20 19:13:13 (antoine)>
*
*
* \note Licence: GPL (>= 2)
@@ -90,6 +90,9 @@ class bigvec_q : public math::Matrix<bigrational> {
*/
void set(unsigned int i,const mpq_t & val);
+ inline void erase(unsigned int index){
+ value.erase(value.begin() + index);
+ }
/**
* \brief add a value
diff --git a/src/extract_matrix.cc b/src/extract_matrix.cc
index 5524dc1c73f356ee7db4a5a1748931524af131df..d2dc3cf72e874a6a4944cf1f139edbfd41bbf397 100644
--- a/src/extract_matrix.cc
+++ b/src/extract_matrix.cc
@@ -13,12 +13,13 @@ SEXP matrix_get_at_q(SEXP A,SEXP INDI, SEXP INDJ)
// for something like x = A[indi, indj], but also simply A[ind]
SEXP matrix_get_at_z(SEXP A,SEXP INDI, SEXP INDJ)
{
+ // printf("ici\n");
bigvec mat = bigintegerR::create_bignum(A);
bigvec mat2 = extract_gmp_R::get_at( mat,INDI,INDJ);
-
+ /*
// now modulus !
// cell based modulus
- if(mat.modulus.size() == mat.value.size())
+ if(mat.getType() == TYPE_MODULUS:: MODULUS_BY_CELL)
{
for(unsigned int i = 0; i< mat.size(); ++i)
mat.value[i] = mat.modulus[i];
@@ -49,8 +50,7 @@ SEXP matrix_get_at_z(SEXP A,SEXP INDI, SEXP INDJ)
{
mat2.modulus.resize(1);
mat2.modulus[0] = mat.modulus[0];
- }
-
+ }*/
return(bigintegerR::create_SEXP(mat2) );
}
diff --git a/src/extract_matrix.h b/src/extract_matrix.h
index daf5763be1a0d453a8f878fcdf7edada9569f255..792a14bf31429f345b0c565c2e572fef34d6068a 100644
--- a/src/extract_matrix.h
+++ b/src/extract_matrix.h
@@ -75,19 +75,23 @@ namespace extract_gmp_R
if(A.nrow < 0)
A.nrow = A.size();
else // check that size is a multiple of row
- if((A.size() / A.nrow) != static_cast<float>(A.size()) / static_cast<float>(A.nrow))
- Rf_error("malformed matrix");
-
+ {
+ if((A.size() / A.nrow) != static_cast<float>(A.size()) / static_cast<float>(A.nrow)){
+ A.clear();
+ Rf_error("malformed matrix");
+ }
+ }
+
retour.resize(A.size() / A.nrow);
for(i = 0; i < retour.size(); ++i)
{
retour[i] = new T();
- retour[i]->value.resize(A.nrow);
+ retour[i]->resize(A.nrow);
}
// go !
- for(i= 0 ; i < A.value.size(); ++i)
+ for(i= 0 ; i < A.size(); ++i)
// retour[col ] [row ]
- (retour[i / A.nrow ])->value[ i % A.nrow].setValue(A.value[i]);
+ (*retour[i / A.nrow ])[ i % A.nrow].setValue(A[i]);
//return(retour);
@@ -116,14 +120,15 @@ namespace extract_gmp_R
*/
template <class T> T get_at (T & A, SEXP& INDI, SEXP& INDJ)
{
-
+ //DEBUG
+ //return A;
// result = A[indi,indj]
int oldnrow = A.nrow;
std::vector<int> indJ;
typename std::vector<T*> matricetmp ;
typename std::vector<T*> matricetmp2;
-
+
toVecVec(A,matricetmp);
typename std::vector<T*> copyAdress(matricetmp);
@@ -189,9 +194,11 @@ namespace extract_gmp_R
it = std::unique(indJ.begin(),indJ.end());
//indJ.erase(it,indJ.end());
- if ( indJ.back() > 0)
+ if ( indJ.back() > 0){
+ clearVec<T>(copyAdress);
+ A.clear();
Rf_error("only 0's may mix with negative subscripts");
-
+ }
it=indJ.begin();
@@ -225,13 +232,18 @@ namespace extract_gmp_R
// for all [new] rows
for( it=indJ.begin(); it != indJ.end(); it++)
{
- if (*it < 0)
+ if (*it < 0){
+ clearVec<T>(copyAdress);
+ A.clear();
Rf_error("only 0's may mix with negative subscripts");
+ }
if( static_cast<unsigned int>(*it-1) < matrice->size() ) {
//printf("on sort %s",(*matrice)[(*it)-1][0].str(10).c_str());
matricework->push_back( (*matrice)[*it-1] );
} else {
- Rf_error("column subscript out of bounds");
+ clearVec<T>(copyAdress);
+ A.clear();
+ Rf_error("column subscript out of bounds");
}
}
@@ -265,7 +277,7 @@ namespace extract_gmp_R
{
// for all columns j delete row i
for (j = 0; j < matrice->size(); ++j)
- (*matrice)[j]->value.erase(i+(*matrice)[j]->value.begin());
+ (*matrice)[j]->erase(i);
//++newnrow;
--i; // i: new indice in modified matrix
@@ -295,9 +307,11 @@ namespace extract_gmp_R
std::unique(indI.begin(),indI.end());
//indI.erase(it,indI.end());
- if ( indI.back() > 0)
+ if ( indI.back() > 0){
+ clearVec<T>(copyAdress);
+ A.clear();
Rf_error("only 0's may mix with negative subscripts");
-
+ }
//newnrow = A.nrow;
@@ -313,7 +327,7 @@ namespace extract_gmp_R
// for all columns j remove row i
for (j = 0; j < matrice->size(); ++j)
{
- (*matrice)[j]->value.erase(i+(*matrice)[j]->value.begin());
+ (*matrice)[j]->erase(i);
}
//--newnrow;
--i; // i: new indice in modified matrix
@@ -337,6 +351,8 @@ namespace extract_gmp_R
/* it = indI.erase(it); */
/* --it; */
} else { // matrix case:
+ clearVec<T>(copyAdress);
+ A.clear();
Rf_error("subscript out of bounds");
}
}
@@ -365,18 +381,21 @@ namespace extract_gmp_R
// for all [new] rows
for( i=0; i < indI.size(); ++i)
{
- if (indI[i] < 0)
+ if (indI[i] < 0){
+ clearVec<T>(copyAdress);
+ A.clear();
Rf_error("only 0's may mix with negative subscripts");
+ }
if( static_cast<unsigned int>(indI[i]-1) < (*matrice)[0]->size() )
{
// for all columns
for (j = 0; j < matricework->size(); ++j)
//newmat[i,j] = oldmat[ind[i],j]
- ( (*matricework)[j])->value[i] = ((*matrice)[j])->value[indI[i]-1];
+ ( *(*matricework)[j])[i] = (*(*matrice)[j])[indI[i]-1];
}
else
for (j = 0; j < matricework->size(); ++j)
- ( (*matricework)[j])->value[i].setValue();
+ ( *(*matricework)[j])[i].setValue();
}
matrice = matricework; // change addresses
@@ -393,7 +412,7 @@ namespace extract_gmp_R
retour.resize(matrice->size() * newnrow);
for(i=0; i< newnrow ; ++i)
for(j=0; j < matrice->size() ; ++j)
- retour.value[i + j* newnrow ] = ((*matrice)[j])->value[i] ;
+ retour[i + j* newnrow ] = (*(*matrice)[j])[i] ;
retour.nrow = (oldnrow < 0) ? -1 : newnrow;
@@ -413,9 +432,12 @@ namespace extract_gmp_R
src.nrow = src.size();
// check that size is a multiple of row
- if((src.size() / src.nrow) != static_cast<float>(src.size()) / static_cast<float>(src.nrow))
+ if((src.size() / src.nrow) != static_cast<float>(src.size()) / static_cast<float>(src.nrow)){
+ src.clear();
+ T & vval = const_cast<T&>(value);
+ vval.clear();
Rf_error("malformed matrix");
-
+ }
unsigned int ncol = src.size() / src.nrow; // number of col
std::vector<int> vidx = indice_get_at ( src.nrow, IDX);
std::vector<int> vjdx = indice_get_at ( ncol, JDX);
@@ -428,7 +450,12 @@ namespace extract_gmp_R
for(unsigned int i = 0 ; i < vidx.size(); ++i)
{
unsigned int index = vidx[i] + vjdx[j] * src.nrow;
- if(index >= src.size()) Rf_error("indice out of bounds");
+ if(index >= src.size()) {
+ src.clear();
+ T & vval = const_cast<T&>(value);
+ vval.clear();
+ Rf_error("indice out of bounds");
+ }
src.set(index, value[k % value.size()] );
++k;
diff --git a/src/factor.cc b/src/factor.cc
index cfb6d2b2a684717ac23169be0964b6862288dc21..a4487133eebd0d34a0a1f326ad573db0dfd1cc40 100644
--- a/src/factor.cc
+++ b/src/factor.cc
@@ -4,7 +4,7 @@
* \version 1
*
* \date Created: 04/12/04
- * \date Last modified: Time-stamp: <2013-06-08 22:47:18 antoine>
+ * \date Last modified: Time-stamp: <2023-01-21 16:28:55 (antoine)>
*
* \author Antoine Lucas (help from Immanuel Scholz) (R adaptation)
* Original C code from libgmp.
@@ -33,12 +33,14 @@ SEXP factorR (SEXP n)
mpz_set(val,v[0].getValue().getValueTemp());
int sgn = mpz_sgn(val);
- if(sgn == 0)
+ if(sgn == 0){
+ v.clear();
error(_("Cannot factorize 0"));
+ }
if(sgn<0)
{
mpz_abs(val,val);
- result.value.push_back(biginteger(-1));
+ result.push_back(bigmod(biginteger(-1)));
}
//
// function from gmplib, in demo/factorize.c
diff --git a/src/factorize.cc b/src/factorize.cc
index de6445b351b375d4bbc28c0f63cff8bc831098a9..bff20c14e8e16f096971ccc66897f5060d02ed00 100644
--- a/src/factorize.cc
+++ b/src/factorize.cc
@@ -180,7 +180,7 @@ mp_prime_p (mpz_t n)
}
}
-
+ factors.clear();
error( "Lucas prime test failure. This should not happen\n");
ret1:
if (flag_prove_primality)
diff --git a/src/matrix.cc b/src/matrix.cc
index 4a8e1c3b40710fc5534217172b69b73b2596f033..3b9f349ebd30a79d7744e1337919d4d93cd9d425 100644
--- a/src/matrix.cc
+++ b/src/matrix.cc
@@ -5,7 +5,7 @@
* \version 1
*
* \date Created: 19/02/06
- * \date Last modified: Time-stamp: <2023-01-16 18:58:59 (antoine)>
+ * \date Last modified: Time-stamp: <2023-01-21 21:18:14 (antoine)>
*
* \author A. Lucas
*
@@ -47,29 +47,34 @@ SEXP as_matrixz (SEXP x, SEXP nrR, SEXP ncR, SEXP byrowR, SEXP mod)
// get "bigz" vector, this makes all conversion int to bigz etc...
bigvec mat = bigintegerR::create_bignum(x);
- int lendat = mat.value.size();
+ int lendat = mat.size();
// int sizemod = mat.modulus.size();
// when modulus specified
bigvec modulus = bigintegerR::create_bignum(mod);
- if(modulus.value.size()>0) // should be allways the case
- if(!modulus.value[0].isNA()) {
- mat.modulus.resize(modulus.size());
- for (unsigned int i = 0; i < modulus.size(); i++)
- mat.modulus[i] = modulus.value[i];
- // sizemod = modulus.size();
- }
// A copy from R project to get correct dimension
// all warnings...
//
- if (nr == NA_INTEGER) // This is < 0
+ if (nr == NA_INTEGER){ // This is < 0
+ mat.clear();
+ modulus.clear();
error(_("matrix: invalid 'nrow' value (too large or NA)"));
- if (nr < 0)
+ }
+ if (nr < 0){
+ mat.clear();
+ modulus.clear();
error(_("matrix: invalid 'nrow' value (< 0)"));
- if (nc < 0)
+ }
+ if (nc < 0){
+ mat.clear();
+ modulus.clear();
error(_("matrix: invalid 'ncol' value (< 0)"));
- if (nc == NA_INTEGER)
+ }
+ if (nc == NA_INTEGER){
+ mat.clear();
+ modulus.clear();
error(_("matrix: invalid 'ncol' value (too large or NA)"));
+ }
if(lendat > 0 ) {
if (lendat > 1 && (nr * nc) % lendat != 0) {
@@ -96,11 +101,29 @@ SEXP as_matrixz (SEXP x, SEXP nrR, SEXP ncR, SEXP byrowR, SEXP mod)
// when we extend "x"
if(nc*nr > lendat)
{
- mat.value.resize(nr*nc);
+ mat.resize(nr*nc);
for(int i = lendat; i < nr*nc; i++)
- mat.value[i] = mat.value[i % lendat];
+ mat[i] = mat[i % lendat];
}
mat.nrow = nr;
+
+
+
+ if(modulus.size()>0) // should be allways the case
+ if(modulus[0].isNA()) {
+ // nothing but mat could have already a modulus
+ } else {
+ for (unsigned int i = 0; i < mat.size(); i++){
+ mat[i].setModulus(modulus[i % modulus.size()].getValuePtr());
+ }
+ if(modulus.size() == 1){
+ mat.setGlobalModulus(modulus[0].getValuePtr());
+ }
+ mat.setType( modulus.size() == 1 ? TYPE_MODULUS::MODULUS_GLOBAL : TYPE_MODULUS::MODULUS_BY_CELL);
+ // sizemod = modulus.size();
+ }
+
+
if(byrow)
{
bigvec mat2 = matrixz::bigint_transpose (mat);
@@ -128,6 +151,7 @@ SEXP bigint_transposeR(SEXP x)
} else if (TYPEOF(dimAttr) == INTSXP) {
nr = INTEGER(dimAttr)[0];
} else { nr = -1;// -Wall
+ mat.clear();
error(_("argument must be a matrix of class \"bigz\""));
}
UNPROTECT(2);
@@ -150,7 +174,7 @@ SEXP matrix_crossp_z (SEXP a, SEXP trans)
bool useMod = FALSE,
tr = (bool)Rf_asLogical(trans);
bigvec mat_a = bigintegerR::create_bignum(a);
- int sizemod = mat_a.modulus.size(),
+ int sizemod = mat_a.getModulusSize(),
a_nrow = mat_a.nrow,
a_len = mat_a.size();
@@ -175,13 +199,11 @@ SEXP matrix_crossp_z (SEXP a, SEXP trans)
mpz_init(tt);
mpz_t common_modulus; mpz_init(common_modulus);
- if(sizemod <= 1) { // maybe 'useMod' i.e., can use common modulus:
- if(sizemod == 1) {
- mpz_set(common_modulus, mat_a.modulus[0].getValueTemp());
- if(!mat_a.modulus[0].isNA())
- useMod = TRUE;
- }
+ if(sizemod == 1) {
+ mpz_set(common_modulus, mat_a.getGlobalModulus()->getValueTemp());
+ useMod = TRUE;
}
+
// here the computation:
for(int i=0; i < m; i++)
@@ -202,16 +224,16 @@ SEXP matrix_crossp_z (SEXP a, SEXP trans)
if(tr) {//------------- tcrossprod ---------------------------
-#define A_I_K mat_a.value [ i + k *a_nrow]
-#define B_K_J mat_a.value [ j + k *a_nrow]
+#define A_I_K mat_a [ i + k *a_nrow]
+#define B_K_J mat_a [ j + k *a_nrow]
K_LOOP
#undef A_I_K
#undef B_K_J
} else {//------------- crossprod ---------------------------
-#define A_I_K mat_a.value [ k + i *a_nrow]
-#define B_K_J mat_a.value [ k + j *a_nrow]
+#define A_I_K mat_a [ k + i *a_nrow]
+#define B_K_J mat_a [ k + j *a_nrow]
K_LOOP
#undef A_I_K
#undef B_K_J
@@ -219,17 +241,18 @@ SEXP matrix_crossp_z (SEXP a, SEXP trans)
}
if(isna) {
- res.value[i + j*m].setValue(0);
- res.value[i + j*m].NA(true);
+ res[i + j*m].setValue(0);
+ res[i + j*m].getValue().NA(true);
}
else
- res.value[i + j*m].setValue(R_ij);
+ res[i + j*m].setValue(R_ij);
} // for(i ..) for(j ..)
- if(useMod)
- res.modulus.push_back(biginteger(common_modulus));
-
+ if(useMod) {
+ std::shared_ptr<biginteger> ptr = std::make_shared<biginteger>(common_modulus);
+ res.setGlobalModulus(ptr);
+ }
mpz_clear(R_ij);
mpz_clear(tt);
mpz_clear(common_modulus);
@@ -248,6 +271,7 @@ SEXP matrix_crossp_z (SEXP a, SEXP trans)
*/
SEXP matrix_mul_z (SEXP a, SEXP b, SEXP op)
{
+
if(!strcmp(class_P(b), "bigq")) { // b "bigq", --> use q arithm:
return(matrix_mul_q(a, b, op));
}
@@ -258,8 +282,9 @@ SEXP matrix_mul_z (SEXP a, SEXP b, SEXP op)
bigvec mat_a = bigintegerR::create_bignum(a),
mat_b = bigintegerR::create_bignum(b);
- int sizemod_a = mat_a.modulus.size(),
- sizemod_b = mat_b.modulus.size();
+
+ int sizemod_a = mat_a.getModulusSize(),
+ sizemod_b = mat_b.getModulusSize();
int
a_nrow = mat_a.nrow, a_len = mat_a.size(),
@@ -340,9 +365,11 @@ SEXP matrix_mul_z (SEXP a, SEXP b, SEXP op)
if(((o_ == 0) && (a_ncol != b_nrow)) ||
((o_ == 1) && (a_nrow != b_nrow)) || // crossprod()
((o_ == 2) && (a_ncol != b_ncol)) // tcrossprod()
- )
+ ){
+ mat_a.clear();
+ mat_b.clear();
error(_("Matrix dimensions do not match"));
-
+ }
// Result R is R[1..n, 1..m] -- and R_{ij} = sum_{k=1} ^ p A.. B..
int n,m, p;
if(o_ == 0) {
@@ -352,6 +379,8 @@ SEXP matrix_mul_z (SEXP a, SEXP b, SEXP op)
}else if (o_ == 2) {
n= a_nrow; m= b_nrow; p= a_ncol;// = b_ncol
}else {
+ mat_a.clear();
+ mat_b.clear();
error(_("invalid 'op' code in matrix_mul_z()"));
n = m = p = -1;// -Wall
}
@@ -366,22 +395,10 @@ SEXP matrix_mul_z (SEXP a, SEXP b, SEXP op)
/* modulus when modulus are "global" (i.e. of size 1) and
* either are the same, or only one of a or b is specified
*/
- if( !(sizemod_a > 1 || sizemod_b > 1)) {
- if((sizemod_a == 1) && (sizemod_b == 1)) {
- mpz_set(common_modulus, mat_a.modulus[0].getValueTemp());
- if(mpz_cmp(common_modulus, mat_b.modulus[0].getValueTemp()) == 0
- && !mat_a.modulus[0].isNA())
- useMod = TRUE;
- }
- else { // at least one of the sizemod_* is > 1 :
- if ((sizemod_a == 1) && !mat_a.modulus[0].isNA()) {
- mpz_set(common_modulus, mat_a[0].getModulus().getValueTemp());
- useMod = TRUE;
- } else if ((sizemod_b == 1) && !mat_b.modulus[0].isNA()) {
- mpz_set(common_modulus, mat_b[0].getModulus().getValueTemp());
- useMod = TRUE;
- }
- }
+ std::shared_ptr<biginteger> globalMod = bigvec::getGlobalModulus(mat_a,mat_b);
+ useMod = globalMod.get() != nullptr;
+ if(useMod){
+ mpz_init_set(common_modulus,globalMod->getValueTemp());
}
// bigmod tmp;
@@ -389,12 +406,12 @@ SEXP matrix_mul_z (SEXP a, SEXP b, SEXP op)
for(int i=0; i < n; i++)
for(int j=0; j < m; j++)
{
-#define R_IJ res.value[ i + j*n]
+#define R_IJ res[ i + j*n]
#define K_LOOP \
for(int k=0; k < p; k++) \
{ \
if(A_I_K.isNA() || B_K_J.isNA()) { \
- R_IJ.setValue(0); R_IJ.NA(true); \
+ R_IJ.setValue(0); R_IJ.getValue().NA(true); \
break; \
} \
/* Z = A_I_K * B_K_J */ \
@@ -410,33 +427,34 @@ SEXP matrix_mul_z (SEXP a, SEXP b, SEXP op)
if(o_ == 0) { //------------- %*% ---------------------------
-#define A_I_K mat_a.value [ i + k *a_nrow]
-#define B_K_J mat_b.value [ k + j *b_nrow]
+#define A_I_K mat_a [ i + k *a_nrow]
+#define B_K_J mat_b [ k + j *b_nrow]
K_LOOP
#undef A_I_K
#undef B_K_J
} else if(o_ == 1){//------------- crossprod ---------------------------
-#define A_I_K mat_a.value [ k + i *a_nrow]
-#define B_K_J mat_b.value [ k + j *b_nrow]
+#define A_I_K mat_a [ k + i *a_nrow]
+#define B_K_J mat_b [ k + j *b_nrow]
K_LOOP
#undef A_I_K
#undef B_K_J
} else {//(o_ == 2) ------------- tcrossprod ---------------------------
-#define A_I_K mat_a.value [ i + k *a_nrow]
-#define B_K_J mat_b.value [ j + k *b_nrow]
+#define A_I_K mat_a [ i + k *a_nrow]
+#define B_K_J mat_b [ j + k *b_nrow]
K_LOOP
#undef A_I_K
#undef B_K_J
}
}
- if(useMod)
- res.modulus.push_back(biginteger(common_modulus));
-
+ if(useMod){
+ std::shared_ptr<biginteger> ptr = std::make_shared<biginteger>(common_modulus);
+ res.setGlobalModulus(ptr);
+ }
mpz_clear(tt);
mpz_clear(common_modulus);
@@ -542,6 +560,10 @@ namespace matrixz
bigvec matbis (mat.size());
matbis.nrow = mat.nCols();
+ if(mat.getType() == MODULUS_GLOBAL){
+ matbis.setGlobalModulus(mat.getGlobalModulus());
+ }
+
/* we compute transpose */
for(unsigned int i =0; i<mat.nRows(); i++)
for(unsigned int j =0; j<mat.nCols(); j++)
@@ -550,12 +572,12 @@ namespace matrixz
return matbis;
}
- /* return dimension in dima */
+ /* return dimension in dima return -2 if dimension does not match */
int checkDims(int dima, int dimb)
{
if(dima > 0 && dimb > 0) {
- if (dimb != dima)
- error(_("Matrix dimensions do not match"));
+ if (dimb != dima) return -2;
+ //error(_("Matrix dimensions do not match"));
}
else { /* either a or b is a matrix */
if(dima == -1)
diff --git a/src/matrixq.cc b/src/matrixq.cc
index 8563cc95d4a76fc763282254a942753093bf6519..ef8d6e7f827ad0e809d62e0bf050f25e86da9511 100644
--- a/src/matrixq.cc
+++ b/src/matrixq.cc
@@ -5,7 +5,7 @@
* \version 1
*
* \date Created: 19/02/06
- * \date Last modified: Time-stamp: <2022-12-09 15:55:39 (antoine)>
+ * \date Last modified: Time-stamp: <2023-01-21 16:26:18 (antoine)>
*
* \author A. Lucas
*
@@ -52,15 +52,27 @@ SEXP as_matrixq (SEXP x, SEXP nrR, SEXP ncR, SEXP byrowR, SEXP den)
/* A copy from R project to get correct dimension
* all warnings...
*/
- if (nr == NA_INTEGER) /* This is < 0 */
+ if (nr == NA_INTEGER){ /* This is < 0 */
+ mat.clear();
+ denominator.clear();
error(_("matrix: invalid 'nrow' value (too large or NA)"));
- if (nr < 0)
+ }
+ if (nr < 0){
+ mat.clear();
+ denominator.clear();
error(_("matrix: invalid 'nrow' value (< 0)"));
- if (nc < 0)
+ }
+ if (nc < 0){
+ mat.clear();
+ denominator.clear();
error(_("matrix: invalid 'ncol' value (< 0)"));
- if (nc == NA_INTEGER)
+ }
+ if (nc == NA_INTEGER){
+ mat.clear();
+ denominator.clear();
error(_("matrix: invalid 'ncol' value (too large or NA)"));
-
+ }
+
if(lendat > 0 ) {
if (lendat > 1 && (nr * nc) % lendat != 0) {
if (((lendat > nr) && (lendat / nr) * nr != lendat) ||
@@ -115,6 +127,7 @@ SEXP bigq_transposeR(SEXP x)
} else if (TYPEOF(dimAttr) == INTSXP) {
nr = INTEGER(dimAttr)[0];
} else {
+ mat.clear();
error(_("argument must be a matrix of class \"bigq\""));
nr = -1;// -Wall
}
@@ -306,9 +319,11 @@ SEXP matrix_mul_q (SEXP a, SEXP b, SEXP op)
if(((o_ == 0) && (a_ncol != b_nrow)) ||
((o_ == 1) && (a_nrow != b_nrow)) || // crossprod()
((o_ == 2) && (a_ncol != b_ncol)) // tcrossprod()
- )
+ ){
+ mat_a.clear();
+ mat_b.clear();
error(_("Matrix dimensions do not match"));
-
+ }
// Result R is R[1..n, 1..m] -- and R_{ij} = sum_{k=1} ^ p A.. B..
int n,m, p;
if(o_ == 0) {
@@ -318,6 +333,8 @@ SEXP matrix_mul_q (SEXP a, SEXP b, SEXP op)
}else if (o_ == 2) {
n= a_nrow; m= b_nrow; p= a_ncol;// = b_ncol
}else {
+ mat_a.clear();
+ mat_b.clear();
error(_("invalid 'op' code in matrix_mul_z()"));
n = m = p = -1;// -Wall
}
diff --git a/src/solve.cc b/src/solve.cc
index 7f46c70e3c84edd8afb311f154a9008512397d36..b4d87b734e8fca1e06664285ed8151ed8d46040c 100644
--- a/src/solve.cc
+++ b/src/solve.cc
@@ -4,7 +4,7 @@
* \version 1
*
* \date Created: 25/05/06
- * \date Last modified: Time-stamp: <2006-05-25 23:05:20 antoine>
+ * \date Last modified: Time-stamp: <2023-01-21 13:18:55 (antoine)>
*
* \author A. Lucas
*
@@ -28,15 +28,17 @@ SEXP inverse_q(SEXP A)
SEXP solve_gmp_R::inverse_q(bigvec_q a)
{
- if(a.nrow * a.nrow != (int) a.size())
+ if(a.nrow * a.nrow != (int) a.size()){
+ a.clear();
error(_("Argument 1 must be a square matrix"));
+ }
bigvec_q b (a.size());
b.nrow = a.nrow;
// initialize b to identity
for(int i=0; i<b.nrow ; ++i)
for(int j=0; j<b.nrow ; ++j)
- b.value[i+j*b.nrow].setValue((i == j) ? 1 : 0);
+ b[i+j*b.nrow].setValue((i == j) ? 1 : 0);
solve_gmp_R::solve(a,b);
@@ -47,18 +49,22 @@ SEXP solve_gmp_R::inverse_q(bigvec_q a)
SEXP inverse_z (SEXP A)
{
bigvec a = bigintegerR::create_bignum(A);
- if(a.modulus.size() == 1 && !a.modulus[0].isNA()) {
+
+ if(a.nrow * a.nrow != (int) a.size()){
+ a.clear();
+ error(_("Argument 1 must be a square matrix"));
+ }
+
+ if(a.getType() == TYPE_MODULUS::MODULUS_GLOBAL) {
bigvec b (a.size() );
b.nrow = a.nrow;
- if(a.nrow * a.nrow != (int) a.size())
- error(_("Argument 1 must be a square matrix"));
-
- b.modulus.push_back(a.modulus[0]);
+
// initialize b to identity
for(int i=0; i<b.nrow ; ++i)
for(int j=0; j<b.nrow ; ++j)
- b.value[i+j*b.nrow].setValue((i == j) ? 1 : 0);
+ b[i+j*b.nrow].setValue((i == j) ? 1 : 0);
+ b.setGlobalModulus(a.getGlobalModulus());
solve_gmp_R::solve(a,b);
return(bigintegerR::create_SEXP(b));
@@ -75,37 +81,41 @@ SEXP solve_z(SEXP A,SEXP B)
bool common_modulus=true;
bigvec a = bigintegerR::create_bignum(A);
bigvec b = bigintegerR::create_bignum(B);
- if(a.modulus.size() == 1 )
- if(!a.modulus[0].isNA())
+
+ // case: b a vector
+ if(b.nrow<1)
+ b.nrow = b.size();
+
+ if(a.nrow * a.nrow != (int) a.size()){
+ a.clear();
+ b.clear();
+ error(_("Argument 1 must be a square matrix"));
+ }
+
+ if(a.nrow != b.nrow){
+ a.clear();
+ b.clear();
+ error(_("Dimensions do not match"));
+ }
+
+ if(a.getType() == TYPE_MODULUS::MODULUS_GLOBAL && b.getType() != TYPE_MODULUS::MODULUS_BY_CELL ) {
+
+ if(b.getType() == TYPE_MODULUS::NO_MODULUS){
+ b.setGlobalModulus(a.getGlobalModulus());
+ }
+
+ if(*a.getGlobalModulus() == *b.getGlobalModulus())
{
- if(b.modulus.size()>1)
- common_modulus = false; // -> solve with rational
- if(b.modulus.size() == 1)
- {
- if(b.modulus[0] == a.modulus[0])
- common_modulus = false; // -> solve with rational
- }
- else
- b.modulus.push_back(a.modulus[0]);
// solve in Z/nZ
- if(common_modulus)
- {
- // case: b a vector
- if(b.nrow<1)
- b.nrow = b.size();
-
- if(a.nrow * a.nrow != (int) a.size())
- error(_("Argument 1 must be a square matrix"));
-
- if(a.nrow != b.nrow)
- error(_("Dimensions do not match"));
-
- solve_gmp_R::solve(a,b);
-
- return(bigintegerR::create_SEXP(b));
- }
+
+
+ solve_gmp_R::solve(a,b);
+
+ return(bigintegerR::create_SEXP(b));
+
}
-
+ }
+ // Solve as rational numbers.
bigvec_q aq (a);
bigvec_q bq (b);
return(solve_gmp_R::solve_q(aq,bq));
@@ -124,16 +134,21 @@ SEXP solve_q(SEXP A,SEXP B)
// solve AX = B
SEXP solve_gmp_R::solve_q(bigvec_q a, bigvec_q b)
{
- if(a.nrow * a.nrow != (int) a.size())
+ if(a.nrow * a.nrow != (int) a.size()){
+ a.clear();
+ b.clear();
error(_("Argument 1 must be a square matrix"));
-
+ }
// case: b a vector
if(b.nrow<0)
b.nrow = b.size();
- if(a.nrow != b.nrow)
+ if(a.nrow != b.nrow){
+ a.clear();
+ b.clear();
error(_("Dimensions do not match"));
-
+ }
+
solve_gmp_R::solve(a,b);
return(bigrationalR::create_SEXP(b));
diff --git a/src/solve.h b/src/solve.h
index 200f0e658cf66ff0c13c7dd5bbc4dfb1ba597039..0d42b6c6ecf0d6354eb8cab60360952bb09e74e2 100644
--- a/src/solve.h
+++ b/src/solve.h
@@ -58,9 +58,12 @@ namespace solve_gmp_R
// A [ i ,j] = A[ i + j * A.nrow]
for(unsigned int k = 0 ; k < A.nRows(); ++k)
{
- if(A.get(k, k).sgn() == 0 )
+ if(A.get(k, k).sgn() == 0 ){
+ A.clear();
+ B.clear();
Rf_error("System is singular");
-
+ }
+
// l_k <- (1/akk) l_k
T tmpValeur =A.get(k , k).inv() ;
A.mulLine(k,tmpValeur);
diff --git a/src/templateMatrix.h b/src/templateMatrix.h
index c7fceed371cc681032501cb27e610384145e24c8..89a952ccf000358c2b303f787f304b8882fc9d1b 100644
--- a/src/templateMatrix.h
+++ b/src/templateMatrix.h
@@ -16,6 +16,8 @@ namespace math{
virtual const T & operator[](unsigned int i) const=0;
virtual T & operator[](unsigned int i) =0;
+
+ virtual ~Vector(){};
};
template< class T>
@@ -44,6 +46,7 @@ namespace math{
/** return true if matrix is supposed to be a vector and not a matrix n x 1 */
virtual bool isVector() const = 0;
+ virtual void clear() = 0;
/**
* \brief assign a value at indice i
*/