1 /************************************************************/
2 /*! \file bigintegerR.cc
3 * \brief C function to interface R and libgmp with big integer values
4 *
5 * \version 1
6 *
7 * \date Created: 27/10/04
8 * \date Last modified: $Id: bigintegerR.cc,v 1.37 2014-09-16 07:33:43 mmaechler Exp $
9 *
10 * \author Immanuel Scholz (help from A. Lucas)
11 *
12 * \note Licence: GPL (>= 2)
13 */
14
15 #include "Rgmp.h"
16
17 #include "bigintegerR.h"
18 // only for <bigz> ^ -|n| :
19 #include "bigrationalR.h"
20
21 #include "matrix.h"
22
23 #include <vector>
24 #include <algorithm>
25
26 using namespace std;
27
28 #include <Rmath.h>
29
30 /* Globals variables */
31
32 static gmp_randstate_t seed_state;
33 static int seed_init=0;
34
35 namespace bigintegerR
36 {
37 // \brief create a vector of bigvecs, all without a modulus.
create_vector(const SEXP & param)38 bigvec create_vector(const SEXP & param) {
39 lockSexp lock (param);
40 switch (TYPEOF(param)) {
41 case NILSXP:
42 return bigvec(); // = bigz(0)
43 case RAWSXP:
44 {
45 // deserialise the vector. first int is the size.
46 bigvec v;
47 const char* raw = (char*)RAW(param);
48 int pos = sizeof(int); // position in raw[]. Starting after header.
49 int sizevec = ((int*)raw)[0];
50 //std::cout << "nb element a lire " << sizevec << std::endl;
51 v.value.resize(sizevec);
52 for (int i = 0; i < sizevec; ++i) {
53 v.value[i] = biginteger(&raw[pos]);
54 pos += v.value[i].raw_size(); // increment number of bytes read.
55 }
56 return v;
57 }
58 case REALSXP:
59 {
60 double* d = REAL(param);
61 //bigvec v(d,d+LENGTH(param));
62 bigvec v;
63 v.value.resize(LENGTH(param));
64 for (int j = 0; j < LENGTH(param); ++j) {
65 /// New: numeric '+- Inf' give +- "Large" instead of NA
66 double dj = d[j];
67 if(R_FINITE(dj) || ISNAN(dj))
68 v.value[j] = dj;
69 else { // dj is +- Inf : use LARGE ( = +- 2 ^ 80000 -- arbitrarily )
70 mpz_t LARGE;
71 mpz_init(LARGE);
72 // FIXME: Keep 'LARGE' a static const; initialized only once
73 mpz_ui_pow_ui (LARGE, (unsigned long int) 2, (unsigned long int) 8000);
74 if(dj == R_PosInf)
75 v.value[j] = LARGE;
76 else if(dj == R_NegInf) {
77 mpz_t neg_L;
78 mpz_init(neg_L);
79 mpz_neg(neg_L, LARGE);
80 v.value[j] = neg_L;
81 mpz_clear(neg_L);
82 }
83 else// should never happen
84 v.value[j] = dj;
85
86 mpz_clear(LARGE);
87 }
88 }
89 return v;
90 }
91 case INTSXP:
92 case LGLSXP:
93 {
94 int* i = INTEGER(param);
95 //bigvec v(i,i+LENGTH(param));
96 bigvec v;
97 v.value.resize(LENGTH(param));
98
99 for (int j = 0; j < LENGTH(param); ++j)
100 v.value[j] = i[j];
101
102 return v;
103 }
104 case STRSXP:
105 {
106 bigvec v;
107 v.value.resize(LENGTH(param));
108 for (int i = 0; i < LENGTH(param); ++i) {
109 if (STRING_ELT(param,i) == NA_STRING)
110 v.value[i]= biginteger();
111 else
112 v.value[i]=biginteger(std::string(CHAR(STRING_ELT(param,i))));
113 }
114 return v;
115 }
116 default:
117 // no longer: can be fatal later! return bigvec();
118 error(_("only logical, numeric or character (atomic) vectors can be coerced to 'bigz'"));
119 }
120 }
121
create_bignum(const SEXP & param)122 bigvec create_bignum(const SEXP & param) {
123 lockSexp lock (param);
124 SEXP
125 modAttr = Rf_getAttrib(param, Rf_mkString("mod")),
126 dimAttr = Rf_getAttrib(param, Rf_mkString("nrow"));
127
128 // try to catch biz-nrow dimension value
129 //std::cout << "import value" << std::endl;
130 bigvec v = bigintegerR::create_vector(param);
131
132 if (TYPEOF(dimAttr) == INTSXP)
133 v.nrow = INTEGER(dimAttr)[0];
134 else {
135 // catch to get std matrix dimensions value
136 dimAttr = Rf_getAttrib(param, Rf_mkString("dim"));
137 v.nrow = (TYPEOF(dimAttr) == INTSXP) ? INTEGER(dimAttr)[0] : -1;// -1: want support 0-row
138 }
139
140 if (TYPEOF(modAttr) != NILSXP)
141 {
142 //std::cout << "import value" << std::endl;
143 v.modulus = bigintegerR::create_vector(modAttr).value;
144 }
145 return v;
146
147 }
148
create_int(const SEXP & param)149 std::vector<int> create_int(const SEXP & param) {
150 lockSexp lock (param);
151 switch (TYPEOF(param)) {
152 case REALSXP:
153 {
154 double* d = REAL(param);
155 // copy vector manually to avoid stupid conversion warning in STL-code :-/
156 vector<int> v;
157 v.reserve(LENGTH(param));
158 for (int i = 0; i < LENGTH(param); ++i)
159 v.push_back(static_cast<int>(d[i]));
160 return v;
161 //return vector<int>(d, d+LENGTH(param));
162 }
163 case INTSXP:
164 case LGLSXP:
165 {
166 int* i = INTEGER(param);
167
168 return std::vector<int>(i, i+LENGTH(param));
169 }
170 default:
171 return std::vector<int>();
172 }
173 }
174
175
create_SEXP(const std::vector<biginteger> & v)176 SEXP create_SEXP(const std::vector<biginteger>& v)
177 {
178 unsigned int i;
179 int size = sizeof(int); // starting with vector-size-header
180 for (i = 0; i < v.size(); ++i)
181 size += v[i].raw_size(); // adding each bigint's needed size
182 SEXP ans = PROTECT(Rf_allocVector(RAWSXP, size));
183 // Rprintf(" o create_SEXP(vect<biginteger>): size=%d, v.size()=%d\n", size, v.size());
184 char* r = (char*)(RAW(ans));
185 ((int*)(r))[0] = v.size(); // first int is vector-size-header
186 int pos = sizeof(int); // current position in r[] (starting after vector-size-header)
187 for (i = 0; i < v.size(); ++i)
188 pos += v[i].as_raw(&r[pos]);
189 UNPROTECT(1);
190 return(ans);
191 }
192
create_SEXP(const bigvec & v)193 SEXP create_SEXP(const bigvec& v) {
194
195 SEXP ans = PROTECT(create_SEXP(v.value));
196 // set the class attribute to "bigz"
197 Rf_setAttrib(ans, R_ClassSymbol, Rf_mkString("bigz"));
198
199 // Rprintf(" o create_SEXP(<bigvec>): v.nrow=%d", v.nrow);
200 // set the dim attribute
201 if(v.nrow >= 0) {
202 SEXP nrowAttr = Rf_mkString("nrow");
203 PROTECT(nrowAttr);
204 SEXP nrowValue = Rf_ScalarInteger((int) v.nrow);
205 PROTECT(nrowValue);
206 Rf_setAttrib(ans, nrowAttr,nrowValue);
207 UNPROTECT(2);
208 }
209 // set the mod attribute
210 if(v.modulus.size() > 0) {
211 SEXP mod = PROTECT(create_SEXP(v.modulus)); // and set *its* class
212 Rf_setAttrib(mod, R_ClassSymbol, Rf_mkString("bigz"));
213 Rf_setAttrib(ans, Rf_mkString("mod"), mod);
214 UNPROTECT(1);
215 }
216 UNPROTECT(1);
217 return ans;
218 }
219
220 /**
221 * \brief Main function of doing a binary operation on bigintegers.
222 * It calls a function argument for doing the correct thing.
223 * This could also be written as a class functor (template)
224 * to save one function call, but then code bloat will happen.
225 */
biginteger_binary_operation(const SEXP & a,const SEXP & b,biginteger_binary_fn f)226 SEXP biginteger_binary_operation(const SEXP& a,const SEXP& b, biginteger_binary_fn f)
227 {
228 bigvec va = bigintegerR::create_bignum(a);
229 bigvec vb = bigintegerR::create_bignum(b), result;
230 int size = (va.value.empty() || vb.value.empty()) ? 0 : max(va.value.size(), vb.value.size());
231 result.value.reserve(size);
232 for (int i = 0; i < size; ++i)
233 result.push_back(f(va[i%va.size()], vb[i%vb.size()]));
234
235 result.nrow = matrixz::checkDims(va.nrow,vb.nrow);
236 // Rprintf(" o bigI_b_op(.); size=%d -> nrow=%d\n", size, result.nrow);
237 return bigintegerR::create_SEXP(result);
238 }
239
240
biginteger_logical_binary_operation(const SEXP & a,const SEXP & b,biginteger_logical_binary_fn f)241 SEXP biginteger_logical_binary_operation(const SEXP & a,const SEXP & b, biginteger_logical_binary_fn f)
242 {
243 bigvec va = bigintegerR::create_bignum(a);
244 bigvec vb = bigintegerR::create_bignum(b), result;
245 int size = (va.value.empty() || vb.value.empty()) ? 0 : max(va.value.size(), vb.value.size());
246 // int sizemod = max(va.modulus.size(), vb.modulus.size());
247 SEXP ans = PROTECT(Rf_allocVector(LGLSXP, size));
248 int *r = LOGICAL(ans);
249 /* TODO: this kind of situation 5 == (5 %% 17)*/
250 for (int i = 0; i < size; ++i) {
251 biginteger am = va.value[i % va.value.size()];
252 biginteger bm = vb.value[i % vb.value.size()];
253 if (am.isNA() || bm.isNA())
254 r[i] = NA_LOGICAL;
255 else
256 r[i] = f(am, bm) ? 1 : 0;
257 }
258
259 int nrow = matrixz::checkDims(va.nrow,vb.nrow) ;
260
261 // Add dimension parameter when available
262 if(nrow >= 0)
263 {
264 SEXP dimVal;
265 PROTECT(dimVal = Rf_allocVector(INTSXP, 2));
266 INTEGER(dimVal)[0] = (int) nrow;
267 INTEGER(dimVal)[1] = (int) size / nrow;
268 Rf_setAttrib(ans, Rf_mkString("dim"), dimVal);
269 UNPROTECT(1);
270 }
271
272 UNPROTECT(1);
273 return ans;
274 }
275
lt(const biginteger & lhs,const biginteger & rhs)276 bool lt(const biginteger& lhs, const biginteger& rhs)
277 {return mpz_cmp(lhs.getValueTemp(), rhs.getValueTemp()) < 0;}
gt(const biginteger & lhs,const biginteger & rhs)278 bool gt(const biginteger& lhs, const biginteger& rhs)
279 {return mpz_cmp(lhs.getValueTemp(), rhs.getValueTemp()) > 0;}
lte(const biginteger & lhs,const biginteger & rhs)280 bool lte(const biginteger& lhs, const biginteger& rhs)
281 {return mpz_cmp(lhs.getValueTemp(), rhs.getValueTemp()) <= 0;}
gte(const biginteger & lhs,const biginteger & rhs)282 bool gte(const biginteger& lhs, const biginteger& rhs)
283 {return mpz_cmp(lhs.getValueTemp(), rhs.getValueTemp()) >= 0;}
eq(const biginteger & lhs,const biginteger & rhs)284 bool eq(const biginteger& lhs, const biginteger& rhs)
285 {return mpz_cmp(lhs.getValueTemp(), rhs.getValueTemp()) == 0;}
neq(const biginteger & lhs,const biginteger & rhs)286 bool neq(const biginteger& lhs, const biginteger& rhs)
287 {return mpz_cmp(lhs.getValueTemp(), rhs.getValueTemp()) != 0;}
288
289 }
290
291 /* End of namespace bigintegerR*/
292
293
R_gmp_get_version()294 SEXP R_gmp_get_version() {
295 return Rf_mkString(gmp_version);
296 }
297
biginteger_add(SEXP a,SEXP b)298 SEXP biginteger_add (SEXP a, SEXP b) {return bigintegerR::biginteger_binary_operation(a,b,operator+);}
biginteger_sub(SEXP a,SEXP b)299 SEXP biginteger_sub (SEXP a, SEXP b) {return bigintegerR::biginteger_binary_operation(a,b,operator-);}
biginteger_mul(SEXP a,SEXP b)300 SEXP biginteger_mul (SEXP a, SEXP b) {return bigintegerR::biginteger_binary_operation(a,b,operator*);}
biginteger_divq(SEXP a,SEXP b)301 SEXP biginteger_divq(SEXP a, SEXP b) {return bigintegerR::biginteger_binary_operation(a,b,operator/);}
biginteger_mod(SEXP a,SEXP b)302 SEXP biginteger_mod (SEXP a, SEXP b) {return bigintegerR::biginteger_binary_operation(a,b,operator%);}
303
biginteger_div(SEXP a,SEXP b)304 SEXP biginteger_div (SEXP a, SEXP b) { // called from "/.bigz" == div.bigz
305 bigvec A = bigintegerR::create_bignum(a),
306 B = bigintegerR::create_bignum(b);
307 // Note: a or b may be simple numbers (e.g. '1') => work with (A,B)
308 int len_m_a = A.modulus.size(),
309 len_m_b = B.modulus.size();
310 if(len_m_a == 0 && len_m_b == 0) // deal with important case quickly:
311 return bigrational_div(a, b);
312 else if(len_m_a == 0) { // and len_m_b > 0:
313 // should work directly using b's "mod" --> compute a * b^(-1)
314 }
315 else if(len_m_b == 0) { // and len_m_a > 0:
316 // should use a's "mod" for b: div_via_inv() need's b's modulus
317 B.modulus = A.modulus;
318 return bigintegerR::biginteger_binary_operation(a,
319 bigintegerR::create_SEXP(B),
320 div_via_inv);
321 }
322 else { // len_m_a > 0 and len_m_b > 0:
323 bool same_mod = true;// are the two mods the "same" (after recycling)?
324 int m = (len_m_a < len_m_b) ? len_m_b : len_m_a; // = max(l..a, l..b)
325 for(int i = 0; i < m; i++)
326 if(A.modulus[i % len_m_a] != B.modulus[i % len_m_b]) {
327 same_mod = false; break;
328 }
329 if(same_mod) {
330 // compute a * b^(-1) ... should work w/o more
331 } else {
332 // use *rational* a/b (not considering 'mod' anymore):
333 return bigrational_div(a, b);
334 }
335 }
336 return bigintegerR::biginteger_binary_operation(a,b, div_via_inv);
337 }
338
biginteger_pow(SEXP a,SEXP b)339 SEXP biginteger_pow (SEXP a, SEXP b) {
340 bigvec v = bigintegerR::create_bignum(a),
341 exp = bigintegerR::create_bignum(b);
342 if(v.modulus.size() == 0) { /* has no modulus: now, if any b < 0, the
343 result must be (non-integer) bigrational */
344 bool use_rat = FALSE;
345 for (unsigned int i = 0; i < exp.value.size(); ++i) {
346 if(mpz_sgn(exp.value[i].getValueTemp()) < 0) {
347 use_rat = TRUE;
348 break;
349 }
350 }
351 if (use_rat) { // a ^ b with some b negative --> rational result
352 // 1) a := as.bigq(a, 1)
353 SEXP one = Rf_ScalarInteger(1);
354 PROTECT(one);
355 SEXP aq = bigrational_as(a, one);
356 PROTECT(aq);
357 // 2) result = <bigq a> ^ b:
358 SEXP ans = bigrational_pow(aq, b);
359 UNPROTECT(2);
360 return ans;
361 }
362 }
363 // else, either, a has a modulus, or (no modulus *and* exp >= 0) :
364 return bigintegerR::biginteger_binary_operation(a,b, pow); // -> pow() in ./bigmod.cc
365 }
biginteger_inv(SEXP a,SEXP b)366 SEXP biginteger_inv (SEXP a, SEXP b) {return bigintegerR::biginteger_binary_operation(a,b,inv);}
biginteger_gcd(SEXP a,SEXP b)367 SEXP biginteger_gcd (SEXP a, SEXP b) {return bigintegerR::biginteger_binary_operation(a,b,gcd);}
biginteger_lcm(SEXP a,SEXP b)368 SEXP biginteger_lcm (SEXP a, SEXP b) {return bigintegerR::biginteger_binary_operation(a,b,lcm);}
biginteger_as(SEXP a,SEXP mod)369 SEXP biginteger_as (SEXP a, SEXP mod){return bigintegerR::biginteger_binary_operation(a,mod,set_modulus);}
370 // set_modulus : -> ./bigmod.cc
371
biginteger_lt(SEXP a,SEXP b)372 SEXP biginteger_lt (SEXP a, SEXP b) {return bigintegerR::biginteger_logical_binary_operation(a,b,bigintegerR::lt);}
biginteger_gt(SEXP a,SEXP b)373 SEXP biginteger_gt (SEXP a, SEXP b) {return bigintegerR::biginteger_logical_binary_operation(a,b,bigintegerR::gt);}
biginteger_lte(SEXP a,SEXP b)374 SEXP biginteger_lte (SEXP a, SEXP b) {return bigintegerR::biginteger_logical_binary_operation(a,b,bigintegerR::lte);}
biginteger_gte(SEXP a,SEXP b)375 SEXP biginteger_gte (SEXP a, SEXP b) {return bigintegerR::biginteger_logical_binary_operation(a,b,bigintegerR::gte);}
biginteger_eq(SEXP a,SEXP b)376 SEXP biginteger_eq (SEXP a, SEXP b) {return bigintegerR::biginteger_logical_binary_operation(a,b,bigintegerR::eq);}
biginteger_neq(SEXP a,SEXP b)377 SEXP biginteger_neq (SEXP a, SEXP b) {return bigintegerR::biginteger_logical_binary_operation(a,b,bigintegerR::neq);}
378
379
380
biginteger_as_character(SEXP a,SEXP b)381 SEXP biginteger_as_character(SEXP a, SEXP b)
382 {
383 bigvec v = bigintegerR::create_bignum(a);
384 SEXP ans;
385 int base = Rf_asInteger(b);
386 if (base < 2 || base > 36)
387 error(_("select a base between 2 and 36"));
388
389 PROTECT(ans = Rf_allocVector(STRSXP, v.size()));
390 for (unsigned int i = 0; i < v.size(); ++i)
391 SET_STRING_ELT(ans, i, Rf_mkChar(v.str(i,base).c_str()));
392 // matrix part
393 if(v.nrow >= 0)
394 {
395 SEXP dim = PROTECT(Rf_allocVector(INTSXP, 2));
396 INTEGER(dim)[0] = v.nrow;
397 INTEGER(dim)[1] = v.value.size() / v.nrow;
398 Rf_setAttrib(ans, Rf_mkString("dim"), dim);
399 UNPROTECT(1);
400 }
401
402 UNPROTECT(1);
403 return ans;
404 }
405
biginteger_as_numeric(SEXP a)406 SEXP biginteger_as_numeric(SEXP a)
407 {
408 bigvec v = bigintegerR::create_bignum(a);
409 SEXP ans = PROTECT(Rf_allocVector(REALSXP,v.size()));
410 double *r = REAL(ans);
411 for (unsigned int i = 0; i < v.size(); ++i)
412 r[i] = v.value[i].isNA() ? NA_REAL : v.value[i].as_double();
413 UNPROTECT(1);
414 return ans;
415 }
416
biginteger_as_integer(SEXP a)417 SEXP biginteger_as_integer(SEXP a)
418 {
419 bigvec v = bigintegerR::create_bignum(a);
420 SEXP ans = PROTECT(Rf_allocVector(INTSXP,v.size()));
421 int *r = INTEGER(ans);
422 for (unsigned int i = 0; i < v.size(); ++i) {
423 if(v.value[i].isNA()) {
424 r[i] = NA_INTEGER;
425 }
426 else if(!mpz_fits_slong_p(v.value[i].getValueTemp())) {
427 Rf_warning("NAs introduced by coercion from big integer");
428 r[i] = NA_INTEGER;
429 } else {
430 r[i] = mpz_get_si(v.value[i].getValueTemp());
431 }
432 }
433 UNPROTECT(1);
434 return ans;
435 }
436
biginteger_get_at(SEXP a,SEXP i)437 SEXP biginteger_get_at(SEXP a, SEXP i)
438 {
439 //result = a [i]
440
441 bigvec va = bigintegerR::create_bignum(a);
442 return(bigintegerR::create_SEXP(bigintegerR::biginteger_get_at_C(va,i)));
443
444 }
445
446 // also called from matrix_get_at_z(.) in ./extract_matrix.cc :
biginteger_get_at_C(bigvec va,SEXP ind)447 bigvec bigintegerR::biginteger_get_at_C(bigvec va, SEXP ind)
448 {
449 vector<int> v_ind = bigintegerR::create_int(ind);
450 bigvec result;
451 // logical: ind = true/false
452 if (TYPEOF(ind) == LGLSXP) {
453 for (unsigned int i = 0; i < va.size(); ++i)
454 if (v_ind[i%v_ind.size()])
455 {
456 //std::cout << "cas LOGIC "<< std::endl;
457 result.push_back(va[i]);
458 }
459 return result;
460 }
461 else {
462 std::remove(v_ind.begin(), v_ind.end(), 0); // remove all zeroes from ind
463 if (v_ind.empty())
464 return bigvec();
465
466 // case: a[-ind]
467 if (v_ind[0] < 0) {
468 //std::cout << "cas ngatif" << std::cout;
469 for (vector<int>::iterator it = v_ind.begin(); it != v_ind.end(); ++it)
470 if (*it > 0)
471 error(_("only 0's may mix with negative subscripts"));
472 else if (-(*it)-1 >= (int)va.size())
473 error(_("subscript out of bounds"));
474
475 // TODO: This is optimized for large va.size and small v_ind.size.
476 // Maybe add a condition to use a different approach for large v_ind's
477 result.value.reserve(va.size()-v_ind.size());
478 for (int i = 0; i < (int)va.size(); ++i)
479 if (find(v_ind.begin(), v_ind.end(), -i-1) == v_ind.end())
480 {
481 result.push_back(va[i]);
482 }
483 }
484 else {
485 // standard case: a[ind] with ind: integers
486 result.value.reserve(v_ind.size());
487 for (vector<int>::iterator it = v_ind.begin(); it != v_ind.end(); ++it) {
488 if (*it <= 0)
489 error(_("only 0's may mix with negative subscripts"));
490 if (*it <= (int)va.size())
491 {
492 //std::cout << "on sort " << va.value[(*it)-1].str(10) << std::endl;
493 result.push_back(va[(*it)-1]);
494 }
495 else
496 result.push_back(DefaultBigMod()); // NA for out of range's
497 }
498 }
499 }
500 return (result);
501 }
502
biginteger_set_at(SEXP src,SEXP idx,SEXP value)503 SEXP biginteger_set_at(SEXP src, SEXP idx, SEXP value)
504 {
505 // return = ( src[idx] <- value )
506
507 bigvec result = bigintegerR::create_bignum(src);
508 bigvec vvalue = bigintegerR::create_bignum(value);
509 vector<int> vidx = bigintegerR::create_int(idx);
510
511 if(vvalue.size() == 0) {
512 if(result.size() == 0)
513 return bigintegerR::create_SEXP(result);
514 else
515 error(_("replacement has length zero"));
516 }
517 //case: logicals
518 if (TYPEOF(idx) == LGLSXP) {
519 int pos = 0;
520 for (unsigned int i = 0; i < result.size(); ++i)
521 if (vidx[i%vidx.size()])
522 result.set(i, vvalue[pos++ % vvalue.size()]);
523 return bigintegerR::create_SEXP(result);
524 }
525 else {
526 std::remove(vidx.begin(), vidx.end(), 0); // remove all zeroes
527 if (vidx.empty())
528 return bigintegerR::create_SEXP(result);
529 // return = (src[-idx] <- value)
530 if (vidx[0] < 0) {
531 for (vector<int>::iterator it = vidx.begin(); it != vidx.end(); ++it)
532 if (*it > 0)
533 error(_("only 0's may mix with negative subscripts"));
534 else if (-(*it)-1 >= (int)result.size())
535 error(_("subscript out of bounds"));
536 int pos = 0;
537 for (int i = 0; i < (int)result.size(); ++i)
538 if (find(vidx.begin(), vidx.end(), -i-1) == vidx.end())
539 result.set(i, vvalue[pos++%vvalue.size()]);
540 }
541 //standard case: return = (src[idx] <- value) with idx: positive integer
542 else {
543 // finding maximum to resize vector if needed
544 int maximum = INT_MIN;
545 for (vector<int>::iterator it = vidx.begin(); it != vidx.end(); ++it)
546 maximum = max(maximum, *it);
547 if (maximum > (int)result.size())
548 result.resize(maximum);
549 int pos = 0;
550 for (vector<int>::iterator it = vidx.begin(); it != vidx.end(); ++it) {
551 if (*it < 0)
552 error(_("only 0's may mix with negative subscripts"));
553 result.set((*it)-1,vvalue[pos++%vvalue.size()]);
554 }
555 }
556 }
557 return bigintegerR::create_SEXP(result);
558 }
559
biginteger_length(SEXP a)560 SEXP biginteger_length(SEXP a)
561 {
562 return Rf_ScalarInteger(bigintegerR::create_bignum(a).size());
563 }
564
biginteger_setlength(SEXP vec,SEXP value)565 SEXP biginteger_setlength(SEXP vec, SEXP value)
566 {
567 int len = 0;
568 switch (TYPEOF(value)) {
569 case INTSXP:
570 case LGLSXP:
571 if (LENGTH(value) != 1)
572 error(_("invalid second argument"));
573 len = Rf_asInteger(value);
574 if (len < 0)
575 error(_("vector size cannot be negative"));
576 else if (len == NA_INTEGER)
577 error(_("vector size cannot be NA"));
578 break;
579 case REALSXP:
580 if (LENGTH(value) != 1)
581 error(_("invalid second argument"));
582 len = (int)*REAL(value);
583 if (len < 0)
584 error(_("vector size cannot be negative"));
585 else if (! (R_FINITE (len ) ))
586 error(_("vector size cannot be NA, NaN of Inf"));
587 break;
588 case STRSXP:
589 // dunno why R spits out this strange error on "length(foo) <- -1"
590 // but I always follow the holy standard ;-)
591 error(_("negative length vectors are not allowed"));
592 default:
593 error(_("invalid second argument"));
594 }
595 bigvec v =bigintegerR::create_bignum(vec);
596 v.resize(len);
597 return bigintegerR::create_SEXP(v);
598 }
599
biginteger_is_na(SEXP a)600 SEXP biginteger_is_na(SEXP a)
601 {
602 bigvec v = bigintegerR::create_bignum(a);
603 SEXP ans = PROTECT(Rf_allocVector(LGLSXP, v.size()));
604 for (unsigned int i = 0; i < v.size(); ++i)
605 LOGICAL(ans)[i] = v[i].getValue().isNA();
606 UNPROTECT(1);
607 return ans;
608 }
609
610
biginteger_sgn(SEXP a)611 SEXP biginteger_sgn(SEXP a)
612 {
613 bigvec v = bigintegerR::create_bignum(a);
614 SEXP ans = PROTECT(Rf_allocVector(INTSXP, v.size()));
615 int *r = INTEGER(ans);
616 for (unsigned int i = 0; i < v.size(); ++i)
617 r[i] = mpz_sgn(v[i].getValue().getValueTemp());
618 UNPROTECT(1);
619 return ans;
620 }
621
biginteger_c(SEXP args)622 SEXP biginteger_c(SEXP args)
623 {
624 // if(TYPEOF(args) != VECSXP) error(_("should be a list"));
625 bigvec result;
626 for(int i=0; i < LENGTH(args); i++) {
627 bigvec v = bigintegerR::create_bignum(VECTOR_ELT(args,i));
628 for(unsigned int j=0; j < v.size(); j++)
629 result.push_back(v[j]);
630 v.clear();
631 }
632 return bigintegerR::create_SEXP(result);
633 }
634
biginteger_cbind(SEXP args)635 SEXP biginteger_cbind(SEXP args)
636 {
637 // if(TYPEOF(args) != VECSXP) error(_("should be a list"));
638 bigvec result = bigintegerR::create_bignum(VECTOR_ELT(args,0));
639 if(result.nrow <= 0)
640 result.nrow = result.size();
641
642 for(int i = 1; i < LENGTH(args);i++)
643 {
644 bigvec v = bigintegerR::create_bignum(VECTOR_ELT(args,i));
645 for(unsigned int j=0; j< v.size() ; j++)
646 result.push_back(v[j]);
647 v.clear();
648 }
649
650 return bigintegerR::create_SEXP(result);
651 }
652
biginteger_rep(SEXP x,SEXP times)653 SEXP biginteger_rep(SEXP x, SEXP times)
654 {
655 bigvec v = bigintegerR::create_bignum(x),
656 result;
657 int rep = Rf_asInteger(times);
658
659 result.value.reserve(v.size()*rep);
660 for(int i = 0 ; i < rep ; i++)
661 for(unsigned int j = 0 ; j < v.size() ; j++)
662 result.push_back(v[j]);
663
664 return bigintegerR::create_SEXP(result);
665 }
666
667
biginteger_is_prime(SEXP a,SEXP reps)668 SEXP biginteger_is_prime(SEXP a, SEXP reps)
669 {
670 bigvec v = bigintegerR::create_bignum(a);
671 vector<int> vb = bigintegerR::create_int(reps);
672 unsigned int i;
673 SEXP ans = PROTECT(Rf_allocVector(INTSXP, v.size()));
674 int *r = INTEGER(ans);
675 if(v.size() == vb.size())
676 for (i = 0; i < v.size(); ++i)
677 r[i] = v[i].getValue().isprime(vb[i]);
678 else
679 for (i = 0; i < v.size(); ++i)
680 r[i] = v[i].getValue().isprime(vb[0]);
681 UNPROTECT(1);
682 return ans;
683 }
684
685
biginteger_nextprime(SEXP a)686 SEXP biginteger_nextprime(SEXP a)
687 {
688 bigvec v = bigintegerR::create_bignum(a),
689 result;
690 result.value.reserve(v.size());
691
692 mpz_t val;
693 mpz_init(val);
694 mpz_t_sentry val_s(val);
695
696 for (unsigned int i = 0; i < v.size(); ++i) {
697 mpz_nextprime(val,v[i].getValue().getValueTemp());
698 result.push_back(DefaultBigMod(val));
699 }
700 return bigintegerR::create_SEXP(result);
701 }
702
biginteger_abs(SEXP a)703 SEXP biginteger_abs(SEXP a)
704 {
705 bigvec v = bigintegerR::create_bignum(a),
706 result;
707 result.value.reserve(v.size());
708
709 mpz_t val;
710 mpz_init(val);
711 mpz_t_sentry val_s(val);
712
713 for (unsigned int i = 0; i < v.size(); ++i)
714 {
715 mpz_abs(val,v[i].getValue().getValueTemp());
716 result.push_back(DefaultBigMod(val));
717
718 // TODO: understand why following lines don't work.
719 //result.push_back(bigmod());
720 //result[i].value.setValue(val);
721 }
722
723 result.modulus = v.modulus;
724 return bigintegerR::create_SEXP(result);
725 }
726
727
728 /** @brief Bezoult coefficients: compute g,s and t as as + bt = g
729 * @param a BigInteger
730 * @param b BigInteger
731 */
biginteger_gcdex(SEXP a,SEXP b)732 SEXP biginteger_gcdex(SEXP a, SEXP b)
733 {
734 bigvec
735 va = bigintegerR::create_bignum(a),
736 vb = bigintegerR::create_bignum(b),
737 result;
738
739 if(va.size() != vb.size())
740 return bigintegerR::create_SEXP(bigvec());
741
742 result.value.reserve(3*va.size());
743
744 mpz_t g;
745 mpz_t s;
746 mpz_t t;
747 mpz_init(g);
748 mpz_init(s);
749 mpz_init(t);
750 mpz_t_sentry val_g(g);
751 mpz_t_sentry val_s(s);
752 mpz_t_sentry val_t(t);
753
754 for(unsigned int i=0; i < va.size(); i++)
755 {
756 mpz_gcdext (g,s,t,va[i].getValue().getValueTemp(),vb[i].getValue().getValueTemp());
757 result.value.push_back(biginteger(g)); // Hem... not very elegant !
758 result.value.push_back(biginteger(s));
759 result.value.push_back(biginteger(t));
760 /* result[i*3].value.setValue(g);
761 result[i*3+1].value.setValue(s);
762 result[i*3+2].value.setValue(t);*/
763
764 }
765 return bigintegerR::create_SEXP(result);
766 }
767
768 /** @brief Random number generation
769 \note If seed is not initialised: generation of a new seed
770 @param nb Integer: number of number to generate
771 @param length Integer number will be of length 2^length
772 @param newseed Integer, seed initialisation (if exists)
773 @param ok Integer 1: seed generation 0 not
774 */
biginteger_rand_u(SEXP nb,SEXP length,SEXP newseed,SEXP ok)775 SEXP biginteger_rand_u (SEXP nb, SEXP length, SEXP newseed, SEXP ok)
776 {
777 mpz_t bz;
778 int i,flag,len,size;
779 bigvec result;
780
781 //extern int seed_init;
782 //extern gmp_randstate_t seed_state;
783
784
785 /* store input data into appropriate mode */
786 bigvec va = bigintegerR::create_bignum(newseed);
787 PROTECT (ok = AS_INTEGER(ok));
788 PROTECT (length = AS_INTEGER(length));
789 PROTECT (nb = AS_INTEGER(nb));
790 flag = Rf_asInteger(ok);
791 len = Rf_asInteger(length);
792 size = Rf_asInteger(nb);
793 UNPROTECT(3);
794
795 result.value.reserve(size);
796
797 /* Random seed initialisation */
798
799 if(seed_init==0)
800 {
801 gmp_randinit_default(seed_state);
802 Rprintf("Seed default initialisation\n");
803 }
804 if(flag == 1)
805 {
806 gmp_randseed(seed_state,va[0].getValue().getValueTemp());
807 Rprintf("Seed initialisation\n");
808 }
809
810 seed_init = 1;
811
812 mpz_init (bz);
813 mpz_t_sentry val_s(bz);
814
815 for(i= 0; i<size; i++)
816 {
817 /* Random number generation */
818 mpz_urandomb(bz,seed_state,len);
819 result.push_back(DefaultBigMod(bz));
820 }
821 return bigintegerR::create_SEXP(result);
822 }
823
824
825 /** @brief biginteger_sizeinbase return
826 * @param x BigInteger
827 * @param base BigInteger
828 */
biginteger_sizeinbase(SEXP x,SEXP base)829 SEXP biginteger_sizeinbase(SEXP x, SEXP base)
830 {
831 bigvec vx = bigintegerR::create_bignum(x);
832 int basesize= Rf_asInteger(base);
833 SEXP ans = PROTECT(Rf_allocVector(INTSXP,vx.size()));
834 int *r = INTEGER(ans);
835 for(unsigned int i=0; i < vx.size(); i++)
836 r[i] = mpz_sizeinbase(vx[i].getValue().getValueTemp(), basesize);
837 UNPROTECT(1);
838 return ans;
839 }
840
841
842 /** @brief bigI_factorial returns n!
843 * @param n non-negative integer vector
844 */
bigI_factorial(SEXP n)845 SEXP bigI_factorial(SEXP n)
846 {
847 bigvec result;
848 int *nn = INTEGER(AS_INTEGER(n)), size = Rf_length(n);
849 result.value.resize(size);
850 for (int i = 0; i < size; ++i) {
851 result.value[i].NA(false);
852 if(nn[i] != NA_INTEGER && nn[i] >= 0) {
853 mpz_fac_ui(result.value[i].getValue(), (unsigned long int)nn[i]);
854 }
855 }
856 return bigintegerR::create_SEXP(result);
857 } // bigI_factorial
858
859 /** @brief bigI_choose(n, k) returns binomial coefficient (n \choose k)
860 * @param n integer, either R "integer" (non-negative), or a "bigz"
861 * @param k non-negative integer
862 */
bigI_choose(SEXP n,SEXP k)863 SEXP bigI_choose(SEXP n, SEXP k)
864 {
865 bigvec result, n_ = bigintegerR::create_bignum(n);
866 int *kk = INTEGER(AS_INTEGER(k)), n_k = Rf_length(k);
867 int size = (n_.value.empty() || n_k == 0) ? 0 :
868 // else: max(n_.value.size(), n_k)
869 (((int)n_.value.size() <= n_k) ? n_k : n_.value.size());
870
871 result.value.resize(size);
872 for (int i = 0; i < size; ++i) {
873 result.value[i].NA(false);
874 int ik_i = kk[i % n_k];
875 // check if k in range:
876 if(ik_i != NA_INTEGER && ik_i >= 0) {
877 unsigned long int k_i = (unsigned long int)ik_i;
878 /* void mpz_bin_ui (mpz_t ROP, mpz_t N, unsigned long int K) */
879 mpz_bin_ui(result.value[i].getValue(),
880 n_.value[i % n_.value.size()].getValueTemp(), k_i);
881 }
882 }
883 return bigintegerR::create_SEXP(result);
884 }
885
886 /** @brief fibnum return nth Fibonacci number
887 * @param n integer
888 */
bigI_fibnum(SEXP n)889 SEXP bigI_fibnum(SEXP n)
890 {
891 bigvec result;
892 if(Rf_length(n) > 0)
893 {
894 int nn = Rf_asInteger(n);
895 unsigned long int num = nn;
896 if(nn < 0 || nn == NA_INTEGER)
897 error(_("argument must be non-negative"));
898 mpz_t val;
899 mpz_init(val);
900 mpz_t_sentry val_s(val);
901
902 mpz_fib_ui(val,num);
903 result.push_back(DefaultBigMod(val));
904 // result[0].value.setValue(val);
905 }
906 // else
907 // error(_("argument must not be an empty list"));
908
909 return bigintegerR::create_SEXP(result);
910
911 }
912
913 /** @brief fibnum2 return nth and n-1th Fibonacci number
914 * @param n integer
915 */
bigI_fibnum2(SEXP n)916 SEXP bigI_fibnum2(SEXP n)
917 {
918 bigvec result;
919 if(Rf_length(n) > 0)
920 {
921 int nn = Rf_asInteger(n);
922 unsigned long int num = nn;
923 if(nn < 0 || nn == NA_INTEGER)
924 error(_("argument must be non-negative"));
925 result.value.reserve(1);
926 mpz_t val;
927 mpz_init(val);
928 mpz_t_sentry val_s(val);
929 mpz_t val2;
930 mpz_init(val2);
931 mpz_t_sentry val_s2(val2);
932
933 mpz_fib2_ui(val,val2, num);
934 result.push_back(DefaultBigMod(val2));
935 result.push_back(DefaultBigMod(val));
936 }
937 else
938 error(_("argument must not be an empty list"));
939
940 return bigintegerR::create_SEXP(result);
941
942 }
943
944 /** @brief lucnum return nth lucas number
945 * @param n integer
946 */
bigI_lucnum(SEXP n)947 SEXP bigI_lucnum(SEXP n)
948 {
949 bigvec result;
950 if(Rf_length(n) > 0)
951 {
952 int nn = Rf_asInteger(n);
953 unsigned long int num = nn;
954 if(nn < 0 || nn == NA_INTEGER)
955 error(_("argument must be non-negative"));
956
957 mpz_t val;
958 mpz_init(val);
959 mpz_t_sentry val_s(val);
960
961 mpz_lucnum_ui(val,num);
962 result.push_back(DefaultBigMod(val));
963 }
964 // else
965 // error(_("argument must not be an empty list"));
966
967 return bigintegerR::create_SEXP(result);
968
969 }
970
971 /** @brief lucnum2 return nth and n-1th lucas number
972 * @param n integer
973 */
bigI_lucnum2(SEXP n)974 SEXP bigI_lucnum2(SEXP n)
975 {
976 bigvec result;
977
978 if(Rf_length(n) > 0) {
979 int nn = Rf_asInteger(n);
980 unsigned long int num = nn;
981 if(nn < 0 || nn == NA_INTEGER)
982 error(_("argument must be non-negative"));
983 mpz_t val;
984 mpz_init(val);
985 mpz_t_sentry val_s(val);
986 mpz_t val2;
987 mpz_init(val2);
988 mpz_t_sentry val_s2(val2);
989
990 mpz_lucnum2_ui(val,val2,num);
991 result.push_back(DefaultBigMod(val2));
992 result.push_back(DefaultBigMod(val));
993 }
994 else
995 error(_("argument must not be an empty list"));
996
997 return bigintegerR::create_SEXP(result);
998 }
999
1000
1001
1002 //Return max
1003
biginteger_max(SEXP a,SEXP narm)1004 SEXP biginteger_max(SEXP a, SEXP narm)
1005 {
1006 bigvec result;
1007 bigvec va = bigintegerR::create_bignum(a);
1008
1009 if( ! va.size())
1010 return bigintegerR::create_SEXP(result);
1011
1012 unsigned int maximum = 0;
1013 int na_remove = Rf_asInteger(narm);
1014
1015 for(unsigned int i = 1 ; i < va.size(); ++i)
1016 {
1017 if(va.value[i].isNA() && !na_remove)
1018 return(bigintegerR::create_SEXP(result));
1019 else
1020 if(!(va.value[i] < va.value[maximum] ))
1021 maximum = i; // if va.value[maximum = 0] is NA => false for the "<" => maximum changed = good
1022 }
1023
1024 result.push_back(va.value[maximum]);
1025
1026
1027 // now the modulus !
1028 if(va.modulus.size() == 1)
1029 result.modulus.push_back(va.modulus[0]);
1030
1031 if(va.modulus.size()>1)
1032 {
1033 biginteger modulus ;
1034 modulus.setValue(va.modulus[0].getValueTemp());
1035 for(unsigned int i = 1 ; i < va.modulus.size() ; ++i)
1036 {
1037 if(modulus != va.modulus[i]) // if one is different: no modulus
1038 return(bigintegerR::create_SEXP(result));
1039 }
1040 result.modulus.push_back(modulus);
1041 }
1042
1043 return(bigintegerR::create_SEXP(result));
1044 }
1045
1046
1047 // Return min
biginteger_min(SEXP a,SEXP narm)1048 SEXP biginteger_min(SEXP a, SEXP narm)
1049 {
1050 bigvec result;
1051 bigvec va = bigintegerR::create_bignum(a);
1052
1053 if( ! va.size())
1054 return bigintegerR::create_SEXP(result);
1055
1056 unsigned int minimum = 0;
1057 int na_remove = Rf_asInteger(narm);
1058
1059 for(unsigned int i = 1 ; i < va.size(); ++i)
1060 {
1061 if(va.value[i].isNA() && !na_remove)
1062 return bigintegerR::create_SEXP(result);
1063 else
1064 if(!(va.value[i] > va.value[minimum] ))
1065 minimum = i;
1066 }
1067
1068 result.push_back(va.value[minimum]);
1069
1070 // now the modulus !
1071 if(va.modulus.size() == 1)
1072 result.modulus.push_back(va.modulus[0]);
1073
1074 if(va.modulus.size()>1)
1075 {
1076 biginteger modulus ;
1077 modulus.setValue(va.modulus[0].getValueTemp());
1078 for(unsigned int i = 1 ; i < va.modulus.size() ; ++i)
1079 {
1080 if(modulus != va.modulus[i]) // if one is different: no modulus
1081 return bigintegerR::create_SEXP(result);
1082 }
1083 result.modulus.push_back(modulus);
1084 }
1085
1086 return bigintegerR::create_SEXP(result);
1087 }
1088
1089
biginteger_cumsum(SEXP a)1090 SEXP biginteger_cumsum(SEXP a)
1091 {
1092 bigvec result, va = bigintegerR::create_bignum(a);
1093
1094 result.value.resize(va.value.size());
1095
1096 mpz_t val;
1097 mpz_init(val);
1098 mpz_t_sentry val_s(val);
1099
1100 bool hasmodulus = true;
1101
1102 // first the modulus !
1103 if(va.modulus.size() > 1) {
1104 biginteger modulus ;
1105 modulus.setValue(va.modulus[0].getValueTemp());
1106
1107 for(unsigned int i = 1 ; i < va.modulus.size() ; ++i) {
1108 if(modulus != va.modulus[i]) { // if one is different: no modulus
1109 hasmodulus = false;
1110 break;
1111 }
1112 }
1113 if(hasmodulus)
1114 result.modulus.push_back(modulus);
1115 }
1116 else if(va.modulus.size() == 1) {
1117 result.modulus.push_back(va.modulus[0]);
1118 hasmodulus = true;
1119 }
1120 else hasmodulus = false;
1121
1122 for(unsigned int i = 0 ; i < va.size(); ++i)
1123 {
1124 {
1125 if(va.value[i].isNA() )
1126 {
1127 break; // all last values are NA.
1128 }
1129
1130 mpz_add(val,val,va.value[i].getValueTemp());
1131
1132 if(hasmodulus)
1133 mpz_mod(val,val,va.modulus[0].getValueTemp() );
1134
1135 result.value[i].setValue(val);
1136 }
1137 }
1138
1139 return(bigintegerR::create_SEXP(result));
1140 }
1141
1142
1143
biginteger_sum(SEXP a)1144 SEXP biginteger_sum(SEXP a)
1145 {
1146 bigvec result, va = bigintegerR::create_bignum(a);
1147
1148 result.value.resize(1);
1149
1150 mpz_t val;
1151 mpz_init(val);
1152 mpz_t_sentry val_s(val);
1153
1154 bool hasmodulus = true;
1155
1156 // first the modulus !
1157 if(va.modulus.size() > 1) {
1158 biginteger modulus ;
1159 modulus.setValue(va.modulus[0].getValueTemp());
1160
1161 for(unsigned int i = 1 ; i < va.modulus.size() ; ++i) {
1162 if(modulus != va.modulus[i]) { // if one is different: no modulus
1163 hasmodulus = false;
1164 break;
1165 }
1166 }
1167 if(hasmodulus)
1168 result.modulus.push_back(modulus);
1169 }
1170 else if(va.modulus.size() == 1) {
1171 result.modulus.push_back(va.modulus[0]);
1172 hasmodulus = true;
1173 }
1174 else
1175 hasmodulus = false;
1176
1177
1178
1179
1180
1181 for(unsigned int i = 0 ; i < va.size(); ++i)
1182 {
1183 {
1184 if(va.value[i].isNA() )
1185 {
1186 break; // all last values are NA.
1187 }
1188
1189 mpz_add(val,val,va.value[i].getValueTemp());
1190
1191 if(hasmodulus)
1192 mpz_mod(val,val,va.modulus[0].getValueTemp() );
1193 }
1194 }
1195
1196 result.value[0].setValue(val);
1197
1198 return(bigintegerR::create_SEXP(result));
1199 }
1200
1201
biginteger_prod(SEXP a)1202 SEXP biginteger_prod(SEXP a)
1203 {
1204 bigvec result;
1205 bigvec va = bigintegerR::create_bignum(a);
1206
1207 result.value.resize(1);
1208
1209 mpz_t val;
1210 mpz_init(val);
1211 mpz_set_ui(val,1);
1212 mpz_t_sentry val_s(val);
1213
1214 bool hasmodulus = true;
1215
1216 // first the modulus !
1217 if(va.modulus.size()>1)
1218 {
1219 biginteger modulus ;
1220 modulus.setValue(va.modulus[0].getValueTemp());
1221
1222 for(unsigned int i = 1 ; i < va.modulus.size() ; ++i) {
1223 if(modulus != va.modulus[i]) { // if one is different: no modulus
1224 hasmodulus = false;
1225 break;
1226 }
1227 }
1228 if(hasmodulus)
1229 result.modulus.push_back(modulus);
1230
1231 }
1232 else
1233 hasmodulus = false;
1234
1235 if(va.modulus.size() == 1)
1236 {
1237 result.modulus.push_back(va.modulus[0]);
1238 hasmodulus = true;
1239 }
1240
1241 for(unsigned int i = 0 ; i < va.size(); ++i)
1242 {
1243 if(va.value[i].isNA() )
1244 {
1245 return (bigintegerR::create_SEXP(result));
1246 }
1247
1248 mpz_mul(val,val,va.value[i].getValueTemp());
1249
1250 if(hasmodulus)
1251 mpz_mod(val,val,va.modulus[0].getValueTemp() );
1252 }
1253
1254 result.value[0].setValue(val);
1255
1256 return(bigintegerR::create_SEXP(result));
1257 }
1258
1259 // return x ^ y [n]
biginteger_powm(SEXP x,SEXP y,SEXP n)1260 SEXP biginteger_powm(SEXP x, SEXP y, SEXP n)
1261 {
1262 bigvec result;
1263
1264 bigvec vx = bigintegerR::create_bignum(x);
1265 bigvec vy = bigintegerR::create_bignum(y);
1266 bigvec vn = bigintegerR::create_bignum(n);
1267
1268 result.value.resize(vx.value.size());
1269
1270 for (unsigned int i = 0 ; i < vx.value.size(); i++)
1271 {
1272
1273 result.value[i].NA(false);
1274 // check if n != 0
1275 if(mpz_sgn(vn.value[i % vn.value.size()].getValueTemp()) != 0)
1276 mpz_powm(result.value[i].getValue(),
1277 vx.value[i].getValueTemp(),
1278 vy.value[i % vy.value.size()].getValueTemp(),
1279 vn.value[i % vn.value.size()].getValueTemp());
1280
1281 }
1282
1283 return bigintegerR::create_SEXP(result);
1284 } // ..._powm()
1285
1286
1287 // TODO: A version that only returns 'ex' {i.e. the binary precision}
1288 // ---- GMP manual suggests that size_t mpz_sizeinbase (mpz_t OP, int BASE)
1289 // i.e., mpz_sizeinbase (OP, 2) would give that
1290
1291 // (d, ex) where x = d * 2^ex, and 0.5 <= |d| < 1
bigI_frexp(SEXP x)1292 SEXP bigI_frexp(SEXP x)
1293 {
1294 // double mpz_get_d_2exp (signed long int *exp, mpz t op )
1295
1296 // Convert op to a double, truncating if necessary (ie. rounding towards zero), and returning
1297 // the exponent separately.
1298 // The return value is in the range 0.5 ≤ |d| < 1 and the exponent is stored to *exp . d ∗ 2exp is
1299 // the (truncated) op value. If op is zero, the return is 0.0 and 0 is stored to *exp .
1300 // This is similar to the standard C frexp function.
1301
1302 const char *nms[] = {"d", "exp", ""};
1303 SEXP ans, d_R, exp_R;
1304 bigvec vx = bigintegerR::create_bignum(x);
1305 int n = vx.value.size();
1306
1307 PROTECT(ans = Rf_mkNamed(VECSXP, nms)); // =~= R list(d = . , exp= .)
1308 d_R = Rf_allocVector(REALSXP, n); SET_VECTOR_ELT(ans, 0, d_R);
1309 exp_R = Rf_allocVector(INTSXP, n); SET_VECTOR_ELT(ans, 1, exp_R);
1310 double *d_ = REAL(d_R);
1311 int *exp_ = INTEGER(exp_R);
1312 for (int i = 0 ; i < n; i++) {
1313 signed long int ex;
1314 d_[i] = mpz_get_d_2exp(&ex, vx.value[i].getValueTemp());
1315 if(abs(ex) < INT_MAX)
1316 exp_[i] = (int) ex;
1317 else
1318 error(_("exponent too large to fit into an integer"));
1319 }
1320 UNPROTECT(1);
1321 return ans;
1322 } // bigI_frexp()
1323
1324
biginteger_log2(SEXP x)1325 SEXP biginteger_log2(SEXP x)
1326 {
1327 bigvec v = bigintegerR::create_bignum(x);
1328 SEXP ans = PROTECT(Rf_allocVector(REALSXP,v.size()));
1329 double *r = REAL(ans);
1330
1331 for (unsigned int i = 0; i < v.size(); ++i) {
1332 signed long int ex;
1333 double di = mpz_get_d_2exp(&ex, v.value[i].getValueTemp());
1334 // xi = di * 2 ^ ex ==> log2(xi) = log2(di) + ex :
1335 r[i] = log(di) / M_LN2 + (double)ex;
1336 }
1337 UNPROTECT(1);
1338 return ans;
1339 }
1340
biginteger_log(SEXP x)1341 SEXP biginteger_log(SEXP x)
1342 {
1343 bigvec v = bigintegerR::create_bignum(x);
1344 SEXP ans = PROTECT(Rf_allocVector(REALSXP,v.size()));
1345 double *r = REAL(ans);
1346
1347 for (unsigned int i = 0; i < v.size(); ++i) {
1348 signed long int ex;
1349 double di = mpz_get_d_2exp(&ex, v.value[i].getValueTemp());
1350 // xi = di * 2 ^ ex ==> log(xi) = log(di) + ex*log(2) :
1351 r[i] = log(di) + M_LN2*(double)ex;
1352 }
1353 UNPROTECT(1);
1354 return ans;
1355 }
1356