Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Maintenance - Mise à jour mensuelle Lundi 6 Février entre 7h00 et 9h00
Open sidebar
Sylvain Jasson
gmp
Commits
c75a94c4
Commit
c75a94c4
authored
Mar 11, 2019
by
antoine lucas
Browse files
version 0.5-13.5
parent
ac08a0e2
Changes
8
Hide whitespace changes
Inline
Side-by-side
src/bigintegerR.cc
View file @
c75a94c4
...
...
@@ -201,7 +201,10 @@ namespace bigintegerR
if
(
v
.
nrow
>=
0
)
{
SEXP
nrowAttr
=
Rf_mkString
(
"nrow"
);
SEXP
nrowValue
=
Rf_ScalarInteger
((
int
)
v
.
nrow
);
PROTECT
(
nrowAttr
);
PROTECT
(
nrowValue
);
Rf_setAttrib
(
ans
,
nrowAttr
,
nrowValue
);
UNPROTECT
(
2
);
}
// set the mod attribute
if
(
v
.
modulus
.
size
()
>
0
)
{
...
...
@@ -348,9 +351,13 @@ 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
one
=
Rf_ScalarInteger
(
1
);
PROTECT
(
one
);
SEXP
aq
=
bigrational_as
(
a
,
one
);
PROTECT
(
aq
);
// 2) result = <bigq a> ^ b:
return
bigrational_pow
(
aq
,
b
);
SEXP
ans
=
bigrational_pow
(
aq
,
b
);
UNPROTECT
(
2
);
return
ans
;
}
}
// else, either, a has a modulus, or (no modulus *and* exp >= 0) :
...
...
@@ -375,7 +382,7 @@ SEXP biginteger_as_character(SEXP a, SEXP b)
{
bigvec
v
=
bigintegerR
::
create_bignum
(
a
);
SEXP
ans
;
int
base
=
INTEGER
(
AS_INTEGER
(
b
))[
0
]
;
int
base
=
Rf_asInteger
(
b
)
;
if
(
base
<
2
||
base
>
36
)
error
(
_
(
"select a base between 2 and 36"
));
...
...
@@ -486,7 +493,7 @@ bigvec bigintegerR::biginteger_get_at_C(bigvec va, SEXP ind)
result
.
push_back
(
va
[(
*
it
)
-
1
]);
}
else
result
.
push_back
(
b
ig
m
od
());
// NA for out of range's
result
.
push_back
(
DefaultB
ig
M
od
());
// NA for out of range's
}
}
}
...
...
@@ -563,7 +570,7 @@ SEXP biginteger_setlength(SEXP vec, SEXP value)
case
LGLSXP
:
if
(
LENGTH
(
value
)
!=
1
)
error
(
_
(
"invalid second argument"
));
len
=
*
INTEGER
(
value
);
len
=
Rf_asInteger
(
value
);
if
(
len
<
0
)
error
(
_
(
"vector size cannot be negative"
));
else
if
(
len
==
NA_INTEGER
)
...
...
@@ -595,7 +602,7 @@ SEXP biginteger_is_na(SEXP a)
bigvec
v
=
bigintegerR
::
create_bignum
(
a
);
SEXP
ans
=
PROTECT
(
Rf_allocVector
(
LGLSXP
,
v
.
size
()));
for
(
unsigned
int
i
=
0
;
i
<
v
.
size
();
++
i
)
LOGICAL
(
ans
)[
i
]
=
v
[
i
].
v
alue
.
isNA
();
LOGICAL
(
ans
)[
i
]
=
v
[
i
].
getV
alue
()
.
isNA
();
UNPROTECT
(
1
);
return
ans
;
}
...
...
@@ -607,7 +614,7 @@ SEXP biginteger_sgn(SEXP a)
SEXP
ans
=
PROTECT
(
Rf_allocVector
(
INTSXP
,
v
.
size
()));
int
*
r
=
INTEGER
(
ans
);
for
(
unsigned
int
i
=
0
;
i
<
v
.
size
();
++
i
)
r
[
i
]
=
mpz_sgn
(
v
[
i
].
v
alue
.
getValueTemp
());
r
[
i
]
=
mpz_sgn
(
v
[
i
].
getV
alue
()
.
getValueTemp
());
UNPROTECT
(
1
);
return
ans
;
}
...
...
@@ -647,7 +654,7 @@ SEXP biginteger_rep(SEXP x, SEXP times)
{
bigvec
v
=
bigintegerR
::
create_bignum
(
x
),
result
;
int
rep
=
INTEGER
(
AS_INTEGER
(
times
)
)[
0
]
;
int
rep
=
Rf_asInteger
(
times
);
result
.
value
.
reserve
(
v
.
size
()
*
rep
);
for
(
int
i
=
0
;
i
<
rep
;
i
++
)
...
...
@@ -667,10 +674,10 @@ SEXP biginteger_is_prime(SEXP a, SEXP reps)
int
*
r
=
INTEGER
(
ans
);
if
(
v
.
size
()
==
vb
.
size
())
for
(
i
=
0
;
i
<
v
.
size
();
++
i
)
r
[
i
]
=
v
[
i
].
v
alue
.
isprime
(
vb
[
i
]);
r
[
i
]
=
v
[
i
].
getV
alue
()
.
isprime
(
vb
[
i
]);
else
for
(
i
=
0
;
i
<
v
.
size
();
++
i
)
r
[
i
]
=
v
[
i
].
v
alue
.
isprime
(
vb
[
0
]);
r
[
i
]
=
v
[
i
].
getV
alue
()
.
isprime
(
vb
[
0
]);
UNPROTECT
(
1
);
return
ans
;
}
...
...
@@ -687,8 +694,8 @@ SEXP biginteger_nextprime(SEXP a)
mpz_t_sentry
val_s
(
val
);
for
(
unsigned
int
i
=
0
;
i
<
v
.
size
();
++
i
)
{
mpz_nextprime
(
val
,
v
[
i
].
v
alue
.
getValueTemp
());
result
.
push_back
(
b
ig
m
od
(
val
));
mpz_nextprime
(
val
,
v
[
i
].
getV
alue
()
.
getValueTemp
());
result
.
push_back
(
DefaultB
ig
M
od
(
val
));
}
return
bigintegerR
::
create_SEXP
(
result
);
}
...
...
@@ -705,8 +712,8 @@ SEXP biginteger_abs(SEXP a)
for
(
unsigned
int
i
=
0
;
i
<
v
.
size
();
++
i
)
{
mpz_abs
(
val
,
v
[
i
].
v
alue
.
getValueTemp
());
result
.
push_back
(
b
ig
m
od
(
val
));
mpz_abs
(
val
,
v
[
i
].
getV
alue
()
.
getValueTemp
());
result
.
push_back
(
DefaultB
ig
M
od
(
val
));
// TODO: understand why following lines don't work.
//result.push_back(bigmod());
...
...
@@ -746,7 +753,7 @@ SEXP biginteger_gcdex(SEXP a, SEXP b)
for
(
unsigned
int
i
=
0
;
i
<
va
.
size
();
i
++
)
{
mpz_gcdext
(
g
,
s
,
t
,
va
[
i
].
v
alue
.
getValueTemp
(),
vb
[
i
].
v
alue
.
getValueTemp
());
mpz_gcdext
(
g
,
s
,
t
,
va
[
i
].
getV
alue
()
.
getValueTemp
(),
vb
[
i
].
getV
alue
()
.
getValueTemp
());
result
.
value
.
push_back
(
biginteger
(
g
));
// Hem... not very elegant !
result
.
value
.
push_back
(
biginteger
(
s
));
result
.
value
.
push_back
(
biginteger
(
t
));
...
...
@@ -780,9 +787,9 @@ SEXP biginteger_rand_u (SEXP nb, SEXP length, SEXP newseed, SEXP ok)
PROTECT
(
ok
=
AS_INTEGER
(
ok
));
PROTECT
(
length
=
AS_INTEGER
(
length
));
PROTECT
(
nb
=
AS_INTEGER
(
nb
));
flag
=
INTEGER
(
ok
)
[
0
]
;
len
=
INTEGER
(
length
)
[
0
]
;
size
=
INTEGER
(
nb
)
[
0
]
;
flag
=
Rf_asInteger
(
ok
);
len
=
Rf_asInteger
(
length
);
size
=
Rf_asInteger
(
nb
);
UNPROTECT
(
3
);
result
.
value
.
reserve
(
size
);
...
...
@@ -796,7 +803,7 @@ SEXP biginteger_rand_u (SEXP nb, SEXP length, SEXP newseed, SEXP ok)
}
if
(
flag
==
1
)
{
gmp_randseed
(
seed_state
,
va
[
0
].
v
alue
.
getValueTemp
());
gmp_randseed
(
seed_state
,
va
[
0
].
getV
alue
()
.
getValueTemp
());
Rprintf
(
"Seed initialisation
\n
"
);
}
...
...
@@ -809,7 +816,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
(
b
ig
m
od
(
bz
));
result
.
push_back
(
DefaultB
ig
M
od
(
bz
));
}
return
bigintegerR
::
create_SEXP
(
result
);
}
...
...
@@ -822,11 +829,11 @@ SEXP biginteger_rand_u (SEXP nb, SEXP length, SEXP newseed, SEXP ok)
SEXP
biginteger_sizeinbase
(
SEXP
x
,
SEXP
base
)
{
bigvec
vx
=
bigintegerR
::
create_bignum
(
x
);
int
basesize
=
INTEGER
(
AS_INTEGER
(
base
)
)[
0
]
;
int
basesize
=
Rf_asInteger
(
base
);
SEXP
ans
=
PROTECT
(
Rf_allocVector
(
INTSXP
,
vx
.
size
()));
int
*
r
=
INTEGER
(
ans
);
for
(
unsigned
int
i
=
0
;
i
<
vx
.
size
();
i
++
)
r
[
i
]
=
mpz_sizeinbase
(
vx
[
i
].
v
alue
.
getValueTemp
(),
basesize
);
r
[
i
]
=
mpz_sizeinbase
(
vx
[
i
].
getV
alue
()
.
getValueTemp
(),
basesize
);
UNPROTECT
(
1
);
return
ans
;
}
...
...
@@ -838,7 +845,7 @@ SEXP biginteger_sizeinbase(SEXP x, SEXP base)
SEXP
bigI_factorial
(
SEXP
n
)
{
bigvec
result
;
int
*
nn
=
INTEGER
(
AS_INTEGER
(
n
)),
size
=
L
ength
(
n
);
int
*
nn
=
INTEGER
(
AS_INTEGER
(
n
)),
size
=
Rf_l
ength
(
n
);
result
.
value
.
resize
(
size
);
for
(
int
i
=
0
;
i
<
size
;
++
i
)
{
result
.
value
[
i
].
NA
(
false
);
...
...
@@ -856,7 +863,7 @@ SEXP bigI_factorial(SEXP n)
SEXP
bigI_choose
(
SEXP
n
,
SEXP
k
)
{
bigvec
result
,
n_
=
bigintegerR
::
create_bignum
(
n
);
int
*
kk
=
INTEGER
(
AS_INTEGER
(
k
)),
n_k
=
L
ength
(
k
);
int
*
kk
=
INTEGER
(
AS_INTEGER
(
k
)),
n_k
=
Rf_l
ength
(
k
);
int
size
=
(
n_
.
value
.
empty
()
||
n_k
==
0
)
?
0
:
// else: max(n_.value.size(), n_k)
(((
int
)
n_
.
value
.
size
()
<=
n_k
)
?
n_k
:
n_
.
value
.
size
());
...
...
@@ -882,9 +889,9 @@ SEXP bigI_choose(SEXP n, SEXP k)
SEXP
bigI_fibnum
(
SEXP
n
)
{
bigvec
result
;
if
(
L
ength
(
n
)
>
0
)
if
(
Rf_l
ength
(
n
)
>
0
)
{
int
nn
=
INTEGER
(
AS_INTEGER
(
n
))[
0
]
;
int
nn
=
Rf_asInteger
(
n
)
;
unsigned
long
int
num
=
nn
;
if
(
nn
<
0
||
nn
==
NA_INTEGER
)
error
(
_
(
"argument must be non-negative"
));
...
...
@@ -893,7 +900,7 @@ SEXP bigI_fibnum(SEXP n)
mpz_t_sentry
val_s
(
val
);
mpz_fib_ui
(
val
,
num
);
result
.
push_back
(
b
ig
m
od
(
val
));
result
.
push_back
(
DefaultB
ig
M
od
(
val
));
// result[0].value.setValue(val);
}
// else
...
...
@@ -909,9 +916,9 @@ SEXP bigI_fibnum(SEXP n)
SEXP
bigI_fibnum2
(
SEXP
n
)
{
bigvec
result
;
if
(
L
ength
(
n
)
>
0
)
if
(
Rf_l
ength
(
n
)
>
0
)
{
int
nn
=
INTEGER
(
AS_INTEGER
(
n
))[
0
]
;
int
nn
=
Rf_asInteger
(
n
)
;
unsigned
long
int
num
=
nn
;
if
(
nn
<
0
||
nn
==
NA_INTEGER
)
error
(
_
(
"argument must be non-negative"
));
...
...
@@ -924,8 +931,8 @@ SEXP bigI_fibnum2(SEXP n)
mpz_t_sentry
val_s2
(
val2
);
mpz_fib2_ui
(
val
,
val2
,
num
);
result
.
push_back
(
b
ig
m
od
(
val2
));
result
.
push_back
(
b
ig
m
od
(
val
));
result
.
push_back
(
DefaultB
ig
M
od
(
val2
));
result
.
push_back
(
DefaultB
ig
M
od
(
val
));
}
else
error
(
_
(
"argument must not be an empty list"
));
...
...
@@ -940,9 +947,9 @@ SEXP bigI_fibnum2(SEXP n)
SEXP
bigI_lucnum
(
SEXP
n
)
{
bigvec
result
;
if
(
L
ength
(
n
)
>
0
)
if
(
Rf_l
ength
(
n
)
>
0
)
{
int
nn
=
INTEGER
(
AS_INTEGER
(
n
))[
0
]
;
int
nn
=
Rf_asInteger
(
n
)
;
unsigned
long
int
num
=
nn
;
if
(
nn
<
0
||
nn
==
NA_INTEGER
)
error
(
_
(
"argument must be non-negative"
));
...
...
@@ -952,7 +959,7 @@ SEXP bigI_lucnum(SEXP n)
mpz_t_sentry
val_s
(
val
);
mpz_lucnum_ui
(
val
,
num
);
result
.
push_back
(
b
ig
m
od
(
val
));
result
.
push_back
(
DefaultB
ig
M
od
(
val
));
}
// else
// error(_("argument must not be an empty list"));
...
...
@@ -968,8 +975,8 @@ SEXP bigI_lucnum2(SEXP n)
{
bigvec
result
;
if
(
L
ength
(
n
)
>
0
)
{
int
nn
=
INTEGER
(
AS_INTEGER
(
n
))[
0
]
;
if
(
Rf_l
ength
(
n
)
>
0
)
{
int
nn
=
Rf_asInteger
(
n
)
;
unsigned
long
int
num
=
nn
;
if
(
nn
<
0
||
nn
==
NA_INTEGER
)
error
(
_
(
"argument must be non-negative"
));
...
...
@@ -981,8 +988,8 @@ SEXP bigI_lucnum2(SEXP n)
mpz_t_sentry
val_s2
(
val2
);
mpz_lucnum2_ui
(
val
,
val2
,
num
);
result
.
push_back
(
b
ig
m
od
(
val2
));
result
.
push_back
(
b
ig
m
od
(
val
));
result
.
push_back
(
DefaultB
ig
M
od
(
val2
));
result
.
push_back
(
DefaultB
ig
M
od
(
val
));
}
else
error
(
_
(
"argument must not be an empty list"
));
...
...
src/bigmod.cc
View file @
c75a94c4
...
...
@@ -36,57 +36,65 @@ bigmod & bigmod::operator= (const bigmod& rhs)
{
if
(
this
!=
&
rhs
)
{
modulus
.
setValue
(
rhs
.
m
odulus
);
modulus
.
setValue
(
rhs
.
getM
odulus
()
);
value
.
setValue
(
rhs
.
value
);
}
return
(
*
this
);
}
bigmod
bigmod
::
inv
()
const
bigmod
&
bigmod
::
inv
()
{
if
(
value
.
isNA
()
||
modulus
.
isNA
())
return
(
bigmod
());
if
(
inverse
!=
NULL
){
inverse
=
NULL
;
delete
inverse
;
}
if
(
value
.
isNA
()
||
modulus
.
isNA
())
{
inverse
=
new
DefaultBigMod
();
return
*
inverse
;
}
mpz_t
val
;
mpz_init
(
val
);
mpz_t_sentry
val_s
(
val
);
if
(
mpz_invert
(
val
,
v
alue
.
getValueTemp
(),
m
odulus
.
getValueTemp
())
==
0
)
{
if
(
mpz_invert
(
val
,
getV
alue
()
.
getValueTemp
(),
getM
odulus
()
.
getValueTemp
())
==
0
)
{
SEXP
wOpt
=
Rf_GetOption1
(
Rf_install
(
"gmp:warnNoInv"
));
if
(
wOpt
!=
R_NilValue
&&
Rf_asInteger
(
wOpt
))
warning
(
_
(
"inv(x) returning NA as x has no inverse"
));
return
(
bigmod
());
// return NA; was
inverse
=
new
DefaultBigMod
();
return
*
inverse
;
// return NA; was
}
return
bigmod
(
val
,
modulus
)
;
inverse
=
new
DefaultBigMod
(
val
,
modulus
);
return
*
inverse
;
}
bool
operator
!=
(
const
bigmod
&
rhs
,
const
bigmod
&
lhs
)
{
if
(
rhs
.
v
alue
!=
lhs
.
v
alue
)
if
(
rhs
.
getV
alue
()
!=
lhs
.
getV
alue
()
)
return
(
true
);
return
(
rhs
.
m
odulus
!=
lhs
.
m
odulus
);
return
(
rhs
.
getM
odulus
()
!=
lhs
.
getM
odulus
()
);
}
bool
operator
==
(
const
bigmod
&
rhs
,
const
bigmod
&
lhs
)
{
if
(
rhs
.
v
alue
!=
lhs
.
v
alue
)
if
(
rhs
.
getV
alue
()
!=
lhs
.
getV
alue
()
)
return
(
false
);
return
(
!
(
rhs
.
m
odulus
!=
lhs
.
m
odulus
));
return
(
!
(
rhs
.
getM
odulus
()
!=
lhs
.
getM
odulus
()
));
}
b
ig
m
od
operator
+
(
const
bigmod
&
lhs
,
const
bigmod
&
rhs
)
DefaultB
ig
M
od
operator
+
(
const
bigmod
&
lhs
,
const
bigmod
&
rhs
)
{
return
create_bigmod
(
lhs
,
rhs
,
mpz_add
);
}
b
ig
m
od
operator
-
(
const
bigmod
&
lhs
,
const
bigmod
&
rhs
)
DefaultB
ig
M
od
operator
-
(
const
bigmod
&
lhs
,
const
bigmod
&
rhs
)
{
return
create_bigmod
(
lhs
,
rhs
,
mpz_sub
);
}
b
ig
m
od
operator
*
(
const
bigmod
&
lhs
,
const
bigmod
&
rhs
)
DefaultB
ig
M
od
operator
*
(
const
bigmod
&
lhs
,
const
bigmod
&
rhs
)
{
return
create_bigmod
(
lhs
,
rhs
,
mpz_mul
);
}
...
...
@@ -96,9 +104,9 @@ bigmod operator*(const bigmod& lhs, const bigmod& rhs)
* ~~~~~~~~~~~~~~
* itself called from "/.bigz" = div.bigz()
*/
b
ig
m
od
div_via_inv
(
const
bigmod
&
a
,
const
bigmod
&
b
)
{
DefaultB
ig
M
od
div_via_inv
(
const
bigmod
&
a
,
const
bigmod
&
b
)
{
// compute a/b as a * b^(-1)
return
operator
*
(
a
,
pow
(
b
,
b
ig
m
od
(
-
1
)));
return
operator
*
(
a
,
pow
(
b
,
DefaultB
ig
M
od
(
-
1
)));
}
...
...
@@ -123,123 +131,123 @@ 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()
*/
b
ig
m
od
operator
/
(
const
bigmod
&
lhs
,
const
bigmod
&
rhs
)
{
DefaultB
ig
M
od
operator
/
(
const
bigmod
&
lhs
,
const
bigmod
&
rhs
)
{
return
create_bigmod
(
lhs
,
rhs
,
integer_div
,
false
);
}
b
ig
m
od
operator
%
(
const
bigmod
&
lhs
,
const
bigmod
&
rhs
)
DefaultB
ig
M
od
operator
%
(
const
bigmod
&
lhs
,
const
bigmod
&
rhs
)
{
if
(
lhs
.
v
alue
.
isNA
()
||
rhs
.
v
alue
.
isNA
())
return
b
ig
m
od
();
if
(
mpz_sgn
(
rhs
.
v
alue
.
getValueTemp
())
==
0
)
{
if
(
lhs
.
getV
alue
()
.
isNA
()
||
rhs
.
getV
alue
()
.
isNA
())
return
DefaultB
ig
M
od
();
if
(
mpz_sgn
(
rhs
.
getV
alue
()
.
getValueTemp
())
==
0
)
{
warning
(
_
(
"biginteger division by zero: returning NA"
));
return
b
ig
m
od
();
return
DefaultB
ig
M
od
();
}
biginteger
mod
;
if
(
!
lhs
.
m
odulus
.
isNA
()
||
!
rhs
.
m
odulus
.
isNA
())
mod
=
rhs
.
v
alue
;
if
(
!
lhs
.
getM
odulus
()
.
isNA
()
||
!
rhs
.
getM
odulus
()
.
isNA
())
mod
=
rhs
.
getV
alue
()
;
mpz_t
val
;
mpz_init
(
val
);
mpz_t_sentry
val_s
(
val
);
mpz_mod
(
val
,
lhs
.
v
alue
.
getValueTemp
(),
rhs
.
v
alue
.
getValueTemp
());
return
b
ig
m
od
(
val
,
mod
);
mpz_mod
(
val
,
lhs
.
getV
alue
()
.
getValueTemp
(),
rhs
.
getV
alue
()
.
getValueTemp
());
return
DefaultB
ig
M
od
(
val
,
mod
);
}
// Either 'base' has a modulus, or it has not *and* exp >= 0 :
b
ig
m
od
pow
(
const
bigmod
&
base
,
const
bigmod
&
exp
)
DefaultB
ig
M
od
pow
(
const
bigmod
&
base
,
const
bigmod
&
exp
)
{
biginteger
mod
=
get_modulus
(
base
,
exp
);
#ifdef DEBUG_bigmod
if
(
mod
.
isNA
()
&&
!
mpz_cmp_si
(
base
.
v
alue
.
getValueTemp
(),
1
))
Rprintf
(
"bigmod pow(1, exp=%d)
\n
"
,
mpz_get_si
(
exp
.
v
alue
.
getValueTemp
()));
else
if
(
mod
.
isNA
()
&&
!
mpz_cmp_si
(
exp
.
v
alue
.
getValueTemp
(),
0
))
Rprintf
(
"bigmod pow(base=%d, 0)
\n
"
,
mpz_get_si
(
base
.
v
alue
.
getValueTemp
()));
if
(
mod
.
isNA
()
&&
!
mpz_cmp_si
(
base
.
getV
alue
()
.
getValueTemp
(),
1
))
Rprintf
(
"bigmod pow(1, exp=%d)
\n
"
,
mpz_get_si
(
exp
.
getV
alue
()
.
getValueTemp
()));
else
if
(
mod
.
isNA
()
&&
!
mpz_cmp_si
(
exp
.
getV
alue
()
.
getValueTemp
(),
0
))
Rprintf
(
"bigmod pow(base=%d, 0)
\n
"
,
mpz_get_si
(
base
.
getV
alue
()
.
getValueTemp
()));
#endif
// if (base == 1 or exp == 0) return 1
if
(
mod
.
isNA
()
&&
((
!
base
.
v
alue
.
isNA
()
&&
!
mpz_cmp_si
(
base
.
v
alue
.
getValueTemp
(),
1
))
||
(
!
exp
.
v
alue
.
isNA
()
&&
!
mpz_cmp_si
(
exp
.
v
alue
.
getValueTemp
(),
0
))))
return
b
ig
m
od
(
biginteger
(
1
));
if
(
base
.
v
alue
.
isNA
()
||
exp
.
v
alue
.
isNA
())
return
b
ig
m
od
();
int
sgn_exp
=
mpz_sgn
(
exp
.
v
alue
.
getValueTemp
());
((
!
base
.
getV
alue
()
.
isNA
()
&&
!
mpz_cmp_si
(
base
.
getV
alue
()
.
getValueTemp
(),
1
))
||
(
!
exp
.
getV
alue
()
.
isNA
()
&&
!
mpz_cmp_si
(
exp
.
getV
alue
()
.
getValueTemp
(),
0
))))
return
DefaultB
ig
M
od
(
biginteger
(
1
));
if
(
base
.
getV
alue
()
.
isNA
()
||
exp
.
getV
alue
()
.
isNA
())
return
DefaultB
ig
M
od
();
int
sgn_exp
=
mpz_sgn
(
exp
.
getV
alue
()
.
getValueTemp
());
bool
neg_exp
=
(
sgn_exp
<
0
);
// b ^ -|e| = 1 / b^|e|
mpz_t
val
;
mpz_init
(
val
);
mpz_t_sentry
val_s
(
val
);
#ifdef DEBUG_bigmod
Rprintf
(
"bigmod pow(base=%3s, exp=%3s [mod=%3s]) ..
\n
"
,
base
.
v
alue
.
str
(
10
).
c_str
(),
exp
.
v
alue
.
str
(
10
).
c_str
(),
base
.
getV
alue
()
.
str
(
10
).
c_str
(),
exp
.
getV
alue
()
.
str
(
10
).
c_str
(),
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
(
!
mpz_fits_ulong_p
(
exp
.
v
alue
.
getValueTemp
()))
if
(
!
mpz_fits_ulong_p
(
exp
.
getV
alue
()
.
getValueTemp
()))
error
(
_
(
"exponent e too large for pow(z,e) = z^e"
));
// FIXME? return( "Inf" )
// else :
mpz_pow_ui
(
val
,
base
.
v
alue
.
getValueTemp
(),
mpz_get_ui
(
exp
.
v
alue
.
getValueTemp
()));
mpz_pow_ui
(
val
,
base
.
getV
alue
()
.
getValueTemp
(),
mpz_get_ui
(
exp
.
getV
alue
()
.
getValueTemp
()));
}
else
if
(
mpz_sgn
(
mod
.
getValueTemp
())
!=
0
)
{
// check modulus non-zero
if
(
neg_exp
)
{
// negative exponent -- only ok if inverse exists
if
(
mpz_invert
(
val
,
base
.
v
alue
.
getValueTemp
(),
mod
.
getValueTemp
())
==
0
)
{
if
(
mpz_invert
(
val
,
base
.
getV
alue
()
.
getValueTemp
(),
mod
.
getValueTemp
())
==
0
)
{
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
(
b
ig
m
od
());
// return NA; was
return
(
DefaultB
ig
M
od
());
// 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
.
v
alue
.
getValueTemp
());
mpz_t
nExp
;
mpz_init
(
nExp
);
mpz_neg
(
nExp
,
exp
.
getV
alue
()
.
getValueTemp
());
mpz_powm
(
val
,
val
,
nExp
,
mod
.
getValueTemp
());
}
else
{
// non-negative exponent
mpz_powm
(
val
,
base
.
v
alue
.
getValueTemp
(),
exp
.
v
alue
.
getValueTemp
(),
mod
.
getValueTemp
());
mpz_powm
(
val
,
base
.
getV
alue
()
.
getValueTemp
(),
exp
.
getV
alue
()
.
getValueTemp
(),
mod
.
getValueTemp
());
}
}
return
b
ig
m
od
(
val
,
mod
);
return
DefaultB
ig
M
od
(
val
,
mod
);
}
b
ig
m
od
inv
(
const
bigmod
&
x
,
const
bigmod
&
m
)
DefaultB
ig
M
od
inv
(
const
bigmod
&
x
,
const
bigmod
&
m
)
{
if
(
x
.
v
alue
.
isNA
()
||
m
.
v
alue
.
isNA
())
return
b
ig
m
od
();
if
(
x
.
getV
alue
()
.
isNA
()
||
m
.
getV
alue
()
.
isNA
())
return
DefaultB
ig
M
od
();
SEXP
wOpt
=
Rf_GetOption1
(
Rf_install
(
"gmp:warnNoInv"
));
bool
warnI
=
(
wOpt
!=
R_NilValue
&&
Rf_asInteger
(
wOpt
));
if
(
mpz_sgn
(
m
.
v
alue
.
getValueTemp
())
==
0
)
{
if
(
mpz_sgn
(
m
.
getV
alue
()
.
getValueTemp
())
==
0
)
{
if
(
warnI
)
warning
(
_
(
"inv(0) returning NA"
));
return
b
ig
m
od
();
return
DefaultB
ig
M
od
();
}
biginteger
mod
=
get_modulus
(
x
,
m
);
mpz_t
val
;
mpz_init
(
val
);
mpz_t_sentry
val_s
(
val
);
if
(
mpz_invert
(
val
,
x
.
v
alue
.
getValueTemp
(),
m
.
v
alue
.
getValueTemp
())
==
0
)
{
if
(
mpz_invert
(
val
,
x
.
getV
alue
()
.
getValueTemp
(),
m
.
getV
alue
()
.
getValueTemp
())
==
0
)
{
if
(
warnI
)
warning
(
_
(
"inv(x,m) returning NA as x has no inverse modulo m"
));
return
(
b
ig
m
od
());
// return NA; was
return
(
DefaultB
ig
M
od
());
// return NA; was
}
return
b
ig
m
od
(
val
,
mod
);
return
DefaultB
ig
M
od
(
val
,
mod
);
}
// R as.bigz() :
b
ig
m
od
set_modulus
(
const
bigmod
&
x
,
const
bigmod
&
m
)
DefaultB
ig
M
od
set_modulus
(
const
bigmod
&
x
,
const
bigmod
&
m
)
{
if
(
!
m
.
v
alue
.
isNA
()
&&
mpz_sgn
(
m
.
v
alue
.
getValueTemp
())
==
0
)
if
(
!
m
.
getV
alue
()
.
isNA
()
&&
mpz_sgn
(
m
.
getV
alue
()
.
getValueTemp
())
==
0
)
error
(
_
(
"modulus 0 is invalid"
));
// if (!m.
v
alue.isNA() && mpz_cmp(x.
v
alue.getValueTemp(),m.
v
alue.getValueTemp())>=0) {
if
(
!
m
.
v
alue
.
isNA
()
)
{
b
ig
m
od
t
(
x
%
m
);
return
b
ig
m
od
(
t
.
v
alue
,
m
.
v
alue
);
// if (!m.
getV
alue
()
.isNA() && mpz_cmp(x.
getV
alue
()
.getValueTemp(),m.
getV
alue
()
.getValueTemp())>=0) {
if
(
!
m
.
getV
alue
()
.
isNA
()
)
{
DefaultB
ig
M
od
t
(
x
%
m
);
return
DefaultB
ig
M
od
(
t
.
getV
alue
()
,
m
.
getV
alue
()
);
}
else
return
b
ig
m
od
(
x
.
v
alue
,
m
.
v
alue
);
return
DefaultB
ig
M
od
(
x
.
getV
alue
()
,
m
.
getV
alue
()
);
}
b
ig
m
od
gcd
(
const
bigmod
&
lhs
,
const
bigmod
&
rhs
)
DefaultB
ig
M
od
gcd
(
const
bigmod
&
lhs
,
const
bigmod
&
rhs
)
{
return
create_bigmod
(
lhs
,
rhs
,
mpz_gcd
);
}
b
ig
m
od
lcm
(
const
bigmod
&
lhs
,
const
bigmod
&
rhs
)
DefaultB
ig
M
od
lcm
(
const
bigmod
&
lhs
,
const
bigmod
&
rhs
)
{
return
create_bigmod
(
lhs
,
rhs
,
mpz_lcm
);
}
...
...
@@ -249,17 +257,17 @@ bigmod lcm(const bigmod& lhs, const bigmod& rhs)
// NA if incompatible.
biginteger
get_modulus
(
const
bigmod
&
b1
,
const
bigmod
&
b2
)
{
if
(
b1
.
m
odulus
.
isNA
())
// NA: means "no modulus" <==> R's is.null(modulus(.))
return
b2
.
m
odulus
;
// if b2 is NA too, the return is correct: NA
else
if
(
b2
.
m
odulus
.
isNA
())
return
b1
.
m
odulus
;
else
if
(
mpz_cmp
(
b1
.
m
odulus
.
getValueTemp
(),
b2
.
m
odulus
.
getValueTemp
()))
{
if
(
b1
.
getM
odulus
()
.
isNA
())
// NA: means "no modulus" <==> R's is.null(modulus(.))
return
b2
.
getM
odulus
()
;
// if b2 is NA too, the return is correct: NA
else
if
(
b2
.
getM
odulus
()
.
isNA
())
return
b1
.
getM
odulus
()
;
else
if
(
mpz_cmp
(
b1
.
getM
odulus
()
.
getValueTemp
(),
b2
.
getM
odulus
()
.
getValueTemp
()))
{
SEXP
wOpt
=
Rf_GetOption1
(
Rf_install
(
"gmp:warnModMismatch"
));
if
(
wOpt
!=
R_NilValue
&&
Rf_asInteger
(
wOpt
))
warning
(
_
(
"modulus mismatch in bigz.* arithmetic"
));
return
biginteger
();
// i.e. NA
}
else
// equal
return
b1
.
m
odulus
;
return
b1
.
getM
odulus
()
;
}
...
...
@@ -269,19 +277,19 @@ biginteger get_modulus(const bigmod& b1, const bigmod& b2)
// Create a bigmod from a binary combination of two other bigmods
b
ig
m
od
create_bigmod
(
const
bigmod
&
lhs
,
const
bigmod
&
rhs
,
gmp_binary
f
,
DefaultB
ig
M
od
create_bigmod
(
const
bigmod
&
lhs
,
const
bigmod
&
rhs
,
gmp_binary
f
,
bool
zeroRhsAllowed
)
{
if
(
lhs
.
v
alue
.
isNA
()
||
rhs
.
v
alue
.
isNA
())
return
b
ig
m
od
();
if
(
!
zeroRhsAllowed
&&
mpz_sgn
(
rhs
.
v
alue
.
getValueTemp
())
==
0
)
{
if
(
lhs
.
getV
alue
()
.
isNA
()
||
rhs
.
getV
alue
()
.
isNA
())
return
DefaultB
ig
M
od
();
if
(
!
zeroRhsAllowed
&&
mpz_sgn
(
rhs
.
getV
alue
()
.
getValueTemp
())
==
0
)
{
warning
(
_
(
"returning NA for (modulus) 0 in RHS"
));
return
b
ig
m
od
();
return
DefaultB
ig
M
od
();
}
biginteger
mod
=
get_modulus
(
lhs
,
rhs
);
mpz_t
val
;
mpz_init
(
val
);
mpz_t_sentry
val_s
(
val
);
f
(
val
,
lhs
.
v
alue
.
getValueTemp
(),
rhs
.
v
alue
.
getValueTemp
());
f
(
val
,
lhs
.
getV
alue
()
.
getValueTemp
(),
rhs
.
getV
alue
()
.
getValueTemp
());
//--- val := f(lhs, rhs)
#ifdef DEBUG_bigmod
bool
iNA
=
biginteger
(
val
).
isNA
();
...
...
@@ -294,8 +302,8 @@ bigmod create_bigmod(const bigmod& lhs, const bigmod& rhs, gmp_binary f,
mpz_get_str
(
buf
,
10
,
val
);
}
Rprintf
(
"create_bigmod(lhs=%3s, rhs=%3s [mod=%3s]) = %s%s"
,
lhs
.
v
alue
.
str
(
10
).
c_str
(),
rhs
.
v
alue
.
str
(
10
).
c_str
(),
lhs
.
getV
alue
()
.
str
(
10
).
c_str
(),
rhs
.
getV
alue
()
.
str
(
10
).
c_str
(),
mod
.
str
(
10
).
c_str
(),
(
iNA
)
?
"NA"
:
buf
,
(
mod
.
isNA
())
?
"
\n
"
:
" {before 'mod'}"
);
...
...
@@ -311,5 +319,5 @@ bigmod create_bigmod(const bigmod& lhs, const bigmod& rhs, gmp_binary f,
}
#endif
}
return
b
ig
m
od
(
val
,
mod
);
return
DefaultB
ig
M
od
(
val
,
mod
);
}
src/bigrationalR.cc
View file @
c75a94c4
...
...
@@ -106,13 +106,17 @@ namespace bigrationalR
if( CHAR(STRING_ELT(param,0)) == "bigz")
return(bigvec_q(bigintegerR::create_bignum(param)) );
*/
lockSexp
lock
(
param
);
PROTECT
(
param
);