Skip to content

Commit

Permalink
- Re-implemented powmod(b, n, m) and complex_powmod(a, b, n, m) f…
Browse files Browse the repository at this point in the history 
…or better performance when `a` and/or `b` are rational numbers.
  • Loading branch information
@trizen
trizen committed Sep 13, 2020
1 parent e41428c commit e025a0a
Show file tree
Hide file tree
Showing 3 changed files with 127 additions and 116 deletions.
218 changes: 102 additions & 116 deletions lib/Sidef/Types/Number/Number.pm
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7180,47 +7180,55 @@ package Sidef::Types::Number::Number {
bless \$t;
}

sub powmod {
my ($x, $y, $z) = @_;
sub _modular_rational {
my ($n, $m) = @_;

_valid(\$y, \$z);
if (ref($n) ne 'Math::GMPq') {
$n = _any2mpq($n) // return;
}

$x = $$x;
$y = _any2mpz($$y) // goto &nan;
$z = _any2mpz($$z) // goto &nan;
state $z = Math::GMPz::Rmpz_init_nobless();

if (ref($x) eq 'Math::GMPz' || __is_int__($x)) {
my $t = Math::GMPz::Rmpz_init();
Math::GMPq::Rmpq_get_den($z, $n);
Math::GMPz::Rmpz_invert($t, $z, $m) or return;
Math::GMPq::Rmpq_get_num($z, $n);
Math::GMPz::Rmpz_mul($t, $t, $z);

$x = _any2mpz($x) // goto &nan;
Math::GMPz::Rmpz_sgn($z) || goto &nan;
return $t;
}

if (Math::GMPz::Rmpz_sgn($y) < 0) {
my $t = Math::GMPz::Rmpz_init();
Math::GMPz::Rmpz_gcd($t, $x, $z);
Math::GMPz::Rmpz_cmp_ui($t, 1) == 0 or goto &nan;
}
sub powmod {
my ($n, $k, $m) = @_;

my $r = Math::GMPz::Rmpz_init();
_valid(\$k, \$m);

Math::GMPz::Rmpz_fits_ulong_p($y)
? Math::GMPz::Rmpz_powm_ui($r, $x, Math::GMPz::Rmpz_get_ui($y), $z)
: Math::GMPz::Rmpz_powm($r, $x, $y, $z);
$n = $$n;
$k = _any2mpz($$k) // goto &nan;
$m = _any2mpz($$m) // goto &nan;

return bless \$r;
}
Math::GMPz::Rmpz_sgn($m) || goto &nan;

$x = _any2mpq($x) // goto &nan;
if (ref($n) ne 'Math::GMPz') {
if (__is_int__($n)) {
$n = _any2mpz($n) // goto &nan;
}
else {
$n = _modular_rational($n, $m) // goto &nan;
}
}

my $num = Math::GMPz::Rmpz_init();
my $den = Math::GMPz::Rmpz_init();
my $r = Math::GMPz::Rmpz_init();

Math::GMPq::Rmpq_get_num($num, $x);
Math::GMPq::Rmpq_get_den($den, $x);
if (Math::GMPz::Rmpz_sgn($k) < 0) {
Math::GMPz::Rmpz_invert($r, $n, $m) or goto &nan;
}

my $numpow = (bless \$num)->powmod((bless \$y), (bless \$z));
my $denpow = (bless \$den)->powmod((bless \$y), (bless \$z));
Math::GMPz::Rmpz_fits_ulong_p($k)
? Math::GMPz::Rmpz_powm_ui($r, $n, Math::GMPz::Rmpz_get_ui($k), $m)
: Math::GMPz::Rmpz_powm($r, $n, $k, $m);

$numpow->mul($denpow->invmod(bless \$z))->mod(bless \$z);
bless \$r;
}

*expmod = \&powmod;
Expand Down Expand Up @@ -7480,122 +7488,100 @@ package Sidef::Types::Number::Number {
$n = _any2mpz($$n) // return (nan(), nan());
$m = _any2mpz($$m) // return (nan(), nan());

if ( (ref($x) eq 'Math::GMPz' || __is_int__($x))
&& (ref($y) eq 'Math::GMPz' || __is_int__($y))) {

$x = _any2mpz($x) // return (nan(), nan());
$y = _any2mpz($y) // return (nan(), nan());

my $c0 = Math::GMPz::Rmpz_init_set_ui(1);
my $c1 = Math::GMPz::Rmpz_init_set_ui(0);

$x = Math::GMPz::Rmpz_init_set($x);
$y = Math::GMPz::Rmpz_init_set($y);

state $t = Math::GMPz::Rmpz_init_nobless();
Math::GMPz::Rmpz_sgn($m) || return (nan(), nan());

# Handle negative exponent
if (Math::GMPz::Rmpz_sgn($n) < 0) {
# Identities for fractional x = a/b, y = c/d:
# ((a/b) + (c/d)*i)^n = ((a*d + b*c*i) / (b*d))^n
# = (a*d + b*c*i)^n * (b*d)^(-n)

$n = Math::GMPz::Rmpz_init_set($n);
Math::GMPz::Rmpz_abs($n, $n);
# We use:
# ((a/b) + (c/d)*i)^n mod m = (a*invmod(b,m) + c*invmod(d, m)*i)^n mod m

my $t = Math::GMPz::Rmpz_init();
if (ref($x) ne 'Math::GMPz') {
if (__is_int__($x)) {
$x = _any2mpz($x) // return (nan(), nan());
}
else {
$x = _modular_rational($x, $m) // return (nan(), nan());
}
}

Math::GMPz::Rmpz_mul($t, $x, $x);
Math::GMPz::Rmpz_addmul($t, $y, $y);
if (ref($y) ne 'Math::GMPz') {
if (__is_int__($y)) {
$y = _any2mpz($y) // return (nan(), nan());
}
else {
$y = _modular_rational($y, $m) // return (nan(), nan());
}
}

if (Math::GMPz::Rmpz_invert($t, $t, $m)) {
$x = _any2mpz($x) // return (nan(), nan());
$y = _any2mpz($y) // return (nan(), nan());

Math::GMPz::Rmpz_mul($c0, $x, $t);
Math::GMPz::Rmpz_mul($c1, $y, $t);
Math::GMPz::Rmpz_neg($c1, $c1);
Math::GMPz::Rmpz_mod($c0, $c0, $m);
Math::GMPz::Rmpz_mod($c1, $c1, $m);
my $c0 = Math::GMPz::Rmpz_init_set_ui(1);
my $c1 = Math::GMPz::Rmpz_init_set_ui(0);

Math::GMPz::Rmpz_set($x, $c0);
Math::GMPz::Rmpz_set($y, $c1);
$x = Math::GMPz::Rmpz_init_set($x);
$y = Math::GMPz::Rmpz_init_set($y);

Math::GMPz::Rmpz_set_ui($c0, 1);
Math::GMPz::Rmpz_set_ui($c1, 0);
}
else { # no inverse exists
return (nan(), nan());
}
}
state $t = Math::GMPz::Rmpz_init_nobless();

foreach my $k (0 .. Math::GMPz::Rmpz_sizeinbase($n, 2) - 1) {
# Handle negative exponent
if (Math::GMPz::Rmpz_sgn($n) < 0) {

if (Math::GMPz::Rmpz_tstbit($n, $k)) {
Math::GMPz::Rmpz_set($t, $c0);
$n = Math::GMPz::Rmpz_init_set($n);
Math::GMPz::Rmpz_abs($n, $n);

Math::GMPz::Rmpz_mul($c0, $c0, $x);
Math::GMPz::Rmpz_submul($c0, $c1, $y);
my $t = Math::GMPz::Rmpz_init();

Math::GMPz::Rmpz_mul($c1, $c1, $x);
Math::GMPz::Rmpz_addmul($c1, $t, $y);
Math::GMPz::Rmpz_mul($t, $x, $x);
Math::GMPz::Rmpz_addmul($t, $y, $y);

Math::GMPz::Rmpz_mod($c0, $c0, $m);
Math::GMPz::Rmpz_mod($c1, $c1, $m);
}
if (Math::GMPz::Rmpz_invert($t, $t, $m)) {

Math::GMPz::Rmpz_mul($t, $x, $y);
Math::GMPz::Rmpz_mul_2exp($t, $t, 1);
Math::GMPz::Rmpz_mul($c0, $x, $t);
Math::GMPz::Rmpz_mul($c1, $y, $t);
Math::GMPz::Rmpz_neg($c1, $c1);
Math::GMPz::Rmpz_mod($c0, $c0, $m);
Math::GMPz::Rmpz_mod($c1, $c1, $m);

Math::GMPz::Rmpz_powm_ui($x, $x, 2, $m);
Math::GMPz::Rmpz_powm_ui($y, $y, 2, $m);
Math::GMPz::Rmpz_set($x, $c0);
Math::GMPz::Rmpz_set($y, $c1);

Math::GMPz::Rmpz_sub($x, $x, $y);
Math::GMPz::Rmpz_mod($y, $t, $m);
Math::GMPz::Rmpz_set_ui($c0, 1);
Math::GMPz::Rmpz_set_ui($c1, 0);
}
else { # no inverse exists
return (nan(), nan());
}

return ((bless \$c0), (bless \$c1));
}

# Handle fractional x = a/b, y = c/d:
# ((a/b) + (c/d)*i)^n = ((a*d + b*c*i) / (b*d))^n
# = (a*d + b*c*i)^n * (b*d)^(-n)
foreach my $k (0 .. Math::GMPz::Rmpz_sizeinbase($n, 2) - 1) {

$x = _any2mpq($x) // return (nan(), nan());
$y = _any2mpq($y) // return (nan(), nan());
if (Math::GMPz::Rmpz_tstbit($n, $k)) {
Math::GMPz::Rmpz_set($t, $c0);

my $A = Math::GMPz::Rmpz_init();
my $B = Math::GMPz::Rmpz_init();
my $C = Math::GMPz::Rmpz_init();
my $D = Math::GMPz::Rmpz_init();
Math::GMPz::Rmpz_mul($c0, $c0, $x);
Math::GMPz::Rmpz_submul($c0, $c1, $y);

Math::GMPq::Rmpq_get_num($A, $x);
Math::GMPq::Rmpq_get_den($B, $x);
Math::GMPz::Rmpz_mul($c1, $c1, $x);
Math::GMPz::Rmpz_addmul($c1, $t, $y);

Math::GMPq::Rmpq_get_num($C, $y);
Math::GMPq::Rmpq_get_den($D, $y);
Math::GMPz::Rmpz_mod($c0, $c0, $m);
Math::GMPz::Rmpz_mod($c1, $c1, $m);
}

state $t = Math::GMPz::Rmpz_init_nobless();
state $u = Math::GMPz::Rmpz_init_nobless();
Math::GMPz::Rmpz_mul($t, $x, $y);
Math::GMPz::Rmpz_mul_2exp($t, $t, 1);

Math::GMPz::Rmpz_mul($u, $B, $D);
Math::GMPz::Rmpz_powm_ui($x, $x, 2, $m);
Math::GMPz::Rmpz_powm_ui($y, $y, 2, $m);

if (Math::GMPz::Rmpz_invert($t, $u, $m)) {
Math::GMPz::Rmpz_powm($u, $u, __neg__($n), $m); # u = (b*d)^(-n) (mod m)
}
else { # no inverse exists
return (nan(), nan());
Math::GMPz::Rmpz_sub($x, $x, $y);
Math::GMPz::Rmpz_mod($y, $t, $m);
}

# Compute: (a*d + b*c*i)^n

Math::GMPz::Rmpz_mul($A, $A, $D); # A*D
Math::GMPz::Rmpz_mul($B, $B, $C); # B*C

my ($re, $im) = complex_powmod((bless \$A), (bless \$B), (bless \$n), (bless \$m));

Math::GMPz::Rmpz_mul($A, $$re, $u);
Math::GMPz::Rmpz_mul($B, $$im, $u);

Math::GMPz::Rmpz_mod($A, $A, $m);
Math::GMPz::Rmpz_mod($B, $B, $m);

return ((bless \$A), (bless \$B));
((bless \$c0), (bless \$c1));
}

*cpowmod = \&complex_powmod;
Expand Down
3 changes: 3 additions & 0 deletions scripts/Tests/gaussian_integers.sf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -153,6 +153,9 @@ assert_eq(Gauss(Mod(13, 97), Mod(43, 97)), Gauss(Mod(13, 97), Mod(43, 97)))
assert_eq(Mod(Gauss(3, 4), 97)*1234, Mod(Gauss(16, 86), 97))
assert_eq(Mod(Gauss(3/4, 5/6), 1234567)**10 * Mod(Gauss(3/4, 5/6), 1234567)**-10, Mod(Gauss(1,0), 1234567))

assert_eq((Mod(Gauss(43/3, 97/5), 127)**(-11) * Mod(Gauss(43/3, 97/5), 127))**+9, Mod(Gauss(52, 73), 127))
assert_eq((Mod(Gauss(43/3, 97/5), 127)**(-11) * Mod(Gauss(43/3, 97/5), 127))**-9, Mod(Gauss(81, 89), 127))

assert_eq((3 + Gauss(4, 5)), Gauss(7, 5))
assert_eq((3 - Gauss(4, 5)), Gauss(-1, -5))
assert_eq((3 * Gauss(4, 5)), Gauss(12, 15))
Expand Down
22 changes: 22 additions & 0 deletions scripts/Tests/mod_class.sf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,28 @@ do {
assert_eq(a**(-1), Mod(14, 19)) # Mod(14, 19)
assert_eq(sqrt(a+1), Mod(4, 19)) # Mod(4, 19)

assert_eq(Mod(3/4, 4171)**1234, Mod(2138, 4171))
assert_eq(Mod(3/4, 4171)**(-1234), Mod(2304, 4171))

assert_eq(powmod(3/4, 1234, 4171), 2138)
assert_eq(powmod(43/97, -129, 57 * 123), 3970)
assert_eq(powmod(3/4, -1234, 4171), 2304)

assert_eq(Mod(3/4, 4171)**(1234), Mod(2138, 4171))
assert_eq(Mod(43/97, 57 * 123)**(-129), Mod(3970, 57*123))
assert_eq(Mod(3/4, 4171)**(-1234), Mod(2304, 4171))

assert_eq(powmod(43/97, 127, 43), 0)
assert_eq(powmod(43/97, 127, 43 * 2), 43)
assert(powmod(43/97, 127, 97 * 3) -> is_nan)

assert_eq(Mod(43/97, 43)**127, 0)
assert_eq(Mod(43/97, 43 * 2)**127, 43)
assert(Mod(43/97, 97 * 3)**128 -> lift.is_nan)

assert_eq(Mod(182398124/123124124, 4171)**(-3), Mod(1970, 4171))
assert_eq(Mod(182398124/123124124, 4171)**(-5), Mod(3635, 4171))

assert_eq(chinese(Mod(43, 19), Mod(13, 41)), Mod(423, 779)) # Mod(423, 779)

with (1e10.irand, 1e10.irand) {|n, m|
Expand Down

0 comments on commit e025a0a

Please sign in to comment.