Commit 3d911e17 authored by Antoine Lucas's avatar Antoine Lucas
Browse files

version 0.5-13.4

parent ab119a96
Package: gmp
Version: 0.5-13.3
Version: 0.5-13.4
Date: 2019-02-21
Title: Multiple Precision Arithmetic
Author: Antoine Lucas, Immanuel Scholz, Rainer Boehme <rb-gmp@reflex-studio.de>,
......@@ -17,6 +17,6 @@ License: GPL
BuildResaveData: no
LazyDataNote: not available, as we use data/*.R *and* our classes
NeedsCompilation: yes
Packaged: 2019-02-21 15:31:59 UTC; l0413674
Packaged: 2019-03-01 10:51:01 UTC; l0413674
Repository: CRAN
Date/Publication: 2018-07-14 14:19:11 UTC
Date/Publication: 2019-03-01 22:00:10 UTC
......@@ -47,6 +47,7 @@ struct mpz_t_sentry {
};
/** \brief Class biginteger
*
* A big integer. Actually a wrapper for mpz_t to work with plus
......
......@@ -36,6 +36,7 @@ namespace bigintegerR
{
// \brief create a vector of bigvecs, all without a modulus.
bigvec create_vector(const SEXP & param) {
lockSexp lock (param);
switch (TYPEOF(param)) {
case NILSXP:
return bigvec(); // = bigz(0)
......@@ -119,7 +120,8 @@ namespace bigintegerR
}
bigvec create_bignum(const SEXP & param) {
SEXP
lockSexp lock (param);
SEXP
modAttr = Rf_getAttrib(param, Rf_mkString("mod")),
dimAttr = Rf_getAttrib(param, Rf_mkString("nrow"));
......@@ -145,6 +147,7 @@ namespace bigintegerR
}
std::vector<int> create_int(const SEXP & param) {
lockSexp lock (param);
switch (TYPEOF(param)) {
case REALSXP:
{
......@@ -195,10 +198,11 @@ namespace bigintegerR
// Rprintf(" o create_SEXP(<bigvec>): v.nrow=%d", v.nrow);
// set the dim attribute
if(v.nrow >= 0) // {
Rf_setAttrib(ans, Rf_mkString("nrow"), Rf_ScalarInteger((int) v.nrow));
// Rprintf(" *SET*\n");
// } else Rprintf(" no set\n");
if(v.nrow >= 0) {
SEXP nrowAttr = Rf_mkString("nrow");
SEXP nrowValue = Rf_ScalarInteger((int) v.nrow);
Rf_setAttrib(ans, nrowAttr,nrowValue);
}
// set the mod attribute
if(v.modulus.size() > 0) {
SEXP mod = PROTECT(create_SEXP(v.modulus)); // and set *its* class
......@@ -343,7 +347,8 @@ SEXP biginteger_pow (SEXP a, SEXP b) {
}
if (use_rat) { // a ^ b with some b negative --> rational result
// 1) a := as.bigq(a, 1)
SEXP aq = bigrational_as(a, Rf_ScalarInteger(1));
SEXP one = Rf_ScalarInteger(1);
SEXP aq = bigrational_as(a, one);
// 2) result = <bigq a> ^ b:
return bigrational_pow(aq, b);
}
......
......@@ -26,6 +26,18 @@ typedef bigmod (*biginteger_binary_fn)(const bigmod&, const bigmod&);
typedef bool (*biginteger_logical_binary_fn)(const biginteger&, const biginteger&);
#endif
struct lockSexp {
lockSexp(const SEXP & value) {
PROTECT(value);
}
~lockSexp(){
UNPROTECT(1);
}
};
/**
* \brief set of function useful for manipulation of SEXP and bigvec
*/
......
......@@ -51,6 +51,8 @@ class bigmod {
value(rhs.value),
modulus(rhs.modulus){}
~bigmod(){};
/**
* \brief Return as a human readible string
*/
......
......@@ -31,6 +31,7 @@ namespace bigrationalR
/** \brief create a vector of bigrationals, all without a denominator.
*/
bigvec_q create_vector(SEXP param) {
lockSexp lock (param);
switch (TYPEOF(param)) {
case NILSXP:
return bigvec_q(); // = bigq(0)
......@@ -105,15 +106,19 @@ namespace bigrationalR
if( CHAR(STRING_ELT(param,0)) == "bigz")
return(bigvec_q(bigintegerR::create_bignum(param)) );
*/
lockSexp lock (param);
bigvec_q v = bigrationalR::create_vector(param);
SEXP denAttr = Rf_getAttrib(param, Rf_mkString("denominator"));
SEXP dimAttr = Rf_getAttrib(param, Rf_mkString("nrow"));
SEXP denKey = Rf_mkString("denominator");
SEXP denAttr = Rf_getAttrib(param, denKey);
SEXP dimKey = Rf_mkString("nrow");
SEXP dimAttr = Rf_getAttrib(param,dimKey );
if (TYPEOF(dimAttr) == INTSXP)
v.nrow = INTEGER(dimAttr)[0];
else {
// catch to get std matrix dimensions value
dimAttr = Rf_getAttrib(param, Rf_mkString("dim"));
dimKey = Rf_mkString("dim");
dimAttr = Rf_getAttrib(param,dimKey );
v.nrow = (TYPEOF(dimAttr) == INTSXP) ? INTEGER(dimAttr)[0] : -1;
}
if (TYPEOF(denAttr) != NILSXP)
......
......@@ -288,7 +288,7 @@ void bigvec::print()
if(nrow > 0) {
for(int i=0; i < nrow; ++i)
{
for(int j=0; j < (value.size() / nrow); ++j)
for(unsigned int j=0; j < (value.size() / nrow); ++j)
Rprintf("%s\t", value[i+j* nrow].str(10).c_str() );
Rprintf("\n");
}
......
......@@ -41,6 +41,8 @@ class bigvec {
*/
bigvec(const bigvec & vecteur);
~bigvec(){};
/**
* \brief construct a bigmod at indice i
*
......
......@@ -55,6 +55,8 @@ class bigvec_q {
/** \brief create from bigintegers */
bigvec_q(const bigvec & rhs);
~bigvec_q(){};
/** \brief assignemt operator */
bigvec_q & operator= (const bigvec_q& rhs);
......
......@@ -325,18 +325,22 @@ namespace extract_gmp_R
}
else
// case: positive values: 1;5;7;7;9;10...
{
// case: positive values: 1;5;7;7;9;10...
// delete too large values or give error iff INDJ is non-NULL
for(it = indI.begin(); it != indI.end(); ++it)
if(*it > static_cast<int>((*matrice)[0]->size()))
if(oldnrow < 0) { // vector case: out-of-bound --> NA
for(it = indI.begin(); it != indI.end(); ++it)
{
if(*it > static_cast<int>((*matrice)[0]->size()))
{
if(oldnrow < 0) { // vector case: out-of-bound --> NA
/* it = indI.erase(it); */
/* --it; */
} else { // matrix case:
} else { // matrix case:
Rf_error("subscript out of bounds");
}
}
}
// note : can have duplicates and indices are not sorted
//newnrow = indI.size();
......
......@@ -46,16 +46,13 @@ factor_using_division (mpz_t t, bigvec & factors)
{
mpz_t q;
unsigned long int p;
int i;
if (flag_verbose > 0)
{
//printf ("[trial division] ");
}
mpz_init (q);
p = mpz_scan1 (t, 0);
mpz_div_2exp (t, t, p);
while (p)
{
......@@ -64,7 +61,7 @@ factor_using_division (mpz_t t, bigvec & factors)
}
p = 3;
for (i = 1; i <= PRIMES_PTAB_ENTRIES;)
for (unsigned int i = 1; i <= PRIMES_PTAB_ENTRIES;)
{
if (! mpz_divisible_ui_p (t, p))
{
......@@ -107,7 +104,7 @@ mp_millerrabin (mpz_srcptr n, mpz_srcptr nm1, mpz_ptr x, mpz_ptr y,
int
mp_prime_p (mpz_t n)
{
int k, r, is_prime;
int k, is_prime;
mpz_t q, a, nm1, tmp;
bigvec factors;
......@@ -115,10 +112,12 @@ mp_prime_p (mpz_t n)
if (mpz_cmp_ui (n, 1) <= 0)
return 0;
/* We have already casted out small primes. */
if (mpz_cmp_ui (n, (long) FIRST_OMITTED_PRIME * FIRST_OMITTED_PRIME) < 0)
return 1;
mpz_init (q);
mpz_init(a);
mpz_init( nm1);
......@@ -149,14 +148,12 @@ mp_prime_p (mpz_t n)
/* Loop until Lucas proves our number prime, or Miller-Rabin proves our
number composite. */
for (r = 0; r < PRIMES_PTAB_ENTRIES; r++)
for (unsigned int r = 0; r < PRIMES_PTAB_ENTRIES; r++)
{
int i;
if (flag_prove_primality)
{
is_prime = 1;
for (i = 0; i < factors.size() && is_prime; i++)
for (unsigned int i = 0; i < factors.size() && is_prime; i++)
{
mpz_divexact (tmp, nm1, factors[i].value.getValue());
mpz_powm (tmp, a, tmp, n);
......@@ -181,8 +178,7 @@ mp_prime_p (mpz_t n)
}
}
//fprintf (stderr, "Lucas prime test failure. This should not happen\n");
//abort ();
error( "Lucas prime test failure. This should not happen\n");
ret1:
if (flag_prove_primality)
......@@ -194,6 +190,7 @@ mp_prime_p (mpz_t n)
mpz_clear(nm1);
mpz_clear(tmp);
return is_prime;
}
......@@ -204,11 +201,6 @@ factor_using_pollard_rho (mpz_t n, unsigned long a, bigvec & factors)
mpz_t t, t2;
unsigned long k, l, i;
if (flag_verbose > 0)
{
//printf ("[pollard-rho (%lu)] ", a);
}
mpz_init (t);
mpz_init(t2);
mpz_init_set_si (y, 2);
......@@ -270,10 +262,6 @@ factor_using_pollard_rho (mpz_t n, unsigned long a, bigvec & factors)
if (!mp_prime_p (t))
{
if (flag_verbose > 0)
{
//printf ("[composite factor--restarting pollard-rho] ");
}
factor_using_pollard_rho (t, a + 1, factors);
}
else
......@@ -304,23 +292,20 @@ void
factor (mpz_t t, bigvec & factors)
{
if (mpz_sgn (t) != 0)
{
factor_using_division (t, factors);
if (mpz_cmp_ui (t, 1) != 0)
{
if (flag_verbose > 0)
{
//printf ("[is number prime?] ");
}
if (mp_prime_p (t))
factors.push_back( t);
else
factor_using_pollard_rho (t, 1, factors);
}
}
}
......@@ -31,7 +31,9 @@ using namespace std;
// given that x is "bigz" or "bigq",
// return TRUE if x is a bigz/q *matrix*: R's is.matrixZQ(.)
SEXP is_matrix_zq(SEXP x) {
return Rf_ScalarLogical(Rf_getAttrib(x, Rf_mkString("nrow")) != R_NilValue);
SEXP nrowSexp = Rf_mkString("nrow");
SEXP attributeRow = Rf_getAttrib(x,nrowSexp );
return Rf_ScalarLogical(attributeRow != R_NilValue);
}
// C++ side of R function matrix.bigz()
......@@ -114,7 +116,9 @@ SEXP as_matrixz (SEXP x, SEXP nrR, SEXP ncR, SEXP byrowR, SEXP mod)
*/
SEXP bigint_transposeR(SEXP x)
{
SEXP dimAttr = Rf_getAttrib(x, Rf_mkString("nrow"));
SEXP dimKey =Rf_mkString("nrow");
SEXP dimAttr = Rf_getAttrib(x,dimKey );
PROTECT(dimAttr);
bigvec mat = bigintegerR::create_bignum(x);
int nr, n = mat.size();
......@@ -125,6 +129,7 @@ SEXP bigint_transposeR(SEXP x)
} else { nr = -1;// -Wall
error(_("argument must be a matrix of class \"bigz\""));
}
UNPROTECT(1);
int nc = (int) n / nr;
// Rprintf(" o bigI_tr(<%d x %d>) ..\n", nr,nc);
return( bigintegerR::create_SEXP(matrixz::bigint_transpose(mat, nr,nc)));
......
......@@ -28,14 +28,14 @@ SEXP inverse_q(SEXP A)
SEXP solve_gmp_R::inverse_q(bigvec_q a)
{
if(a.nrow * a.nrow != a.size())
if(a.nrow * a.nrow != (int) a.size())
error(_("Argument 1 must be a square matrix"));
bigvec_q b (a.size());
b.nrow = a.nrow;
// initialize b to identity
for(unsigned int i=0; i<b.nrow ; ++i)
for(unsigned int j=0; j<b.nrow ; ++j)
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);
solve_gmp_R::solve(a,b);
......@@ -50,13 +50,13 @@ SEXP inverse_z (SEXP A)
if(a.modulus.size() == 1 && !a.modulus[0].isNA()) {
bigvec b (a.size() );
b.nrow = a.nrow;
if(a.nrow * a.nrow != a.size())
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(unsigned int i=0; i<b.nrow ; ++i)
for(unsigned int j=0; j<b.nrow ; ++j)
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);
solve_gmp_R::solve(a,b);
......@@ -94,7 +94,7 @@ SEXP solve_z(SEXP A,SEXP B)
if(b.nrow<1)
b.nrow = b.size();
if(a.nrow * a.nrow != a.size())
if(a.nrow * a.nrow != (int) a.size())
error(_("Argument 1 must be a square matrix"));
if(a.nrow != b.nrow)
......@@ -124,7 +124,7 @@ 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 != a.size())
if(a.nrow * a.nrow != (int) a.size())
error(_("Argument 1 must be a square matrix"));
// case: b a vector
......
......@@ -55,13 +55,13 @@ namespace solve_gmp_R
template< class T> void solve (T & A , T & B)
{
// n = A.nrow
unsigned int i;
int i;
// a temp value
T tmpValeur(1);
// A [ i ,j] = A[ i + j * A.nrow]
for(unsigned int k = 0 ; k < A.nrow; ++k)
for(int k = 0 ; k < A.nrow; ++k)
{
if(A[k + k* A.nrow ].sgn() == 0 )
Rf_error("System is singular");
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment