Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

Alternate mulmod (polynomial multiplication via convolution). Add Big…

…Poly tests.
  • Loading branch information...
commit 4b1f4917d02fbce8376ddd6088c2adac81906df9 1 parent f6a2bc3
@danaj danaj authored
Showing with 51 additions and 52 deletions.
  1. +18 −52 lib/Math/Primality/BigPolynomial.pm
  2. +33 −0 t/big_polynomial.t
View
70 lib/Math/Primality/BigPolynomial.pm
@@ -14,21 +14,22 @@ Math::Primality::BigPolynomials - Polynomials with BigInts
sub new {
- my $self = {};
- my $class = shift;
- my $construction_junk = shift;
+ my $self = {};
+ my $class = shift;
+ my $init = shift;
$self->{ZERO} = Math::GMPz->new(0);
- if ($construction_junk) {
- my $type = ref $construction_junk;
+ if ($init) {
+ my $type = ref $init;
if ( $type eq 'ARRAY' ) {
- $self->{COEF} = $construction_junk;
+ die "Initialization array must be non-empty" unless @$init > 1;
+ $self->{COEF} = $init;
} elsif ($type eq 'Math::Primality::BigPolynomial') {
- copy($self, $construction_junk);
+ $self->{COEF} = [ map { Math::GMPz->new($_) } $init->coef ];
} else {
my $a = [];
- for ( my $i = 0 ; $i < $construction_junk ; $i++ ) {
+ for ( my $i = 0 ; $i < $init ; $i++ ) {
push @$a, $self->{ZERO};
}
$self->{COEF} = $a;
@@ -43,10 +44,7 @@ sub new {
sub copy {
my $self = shift;
- my $other = shift;
- die "Copy called with a non-object" unless ref($other) eq 'Math::Primality::BigPolynomial';
- $self->{COEF} = [ map { Math::GMPz->new($_) } @{$other->{COEF}} ];
- 1;
+ return Math::Primality::BigPolynomial->new($self);
}
sub string {
@@ -135,49 +133,17 @@ sub mulmod {
my @xset = map { Rmpz_sgn($_) > 0 } @x;
my @yset = map { Rmpz_sgn($_) > 0 } @y;
- $self->clear();
-
- # We use these in the loop.
- my $sum = Math::GMPz->new(0);
- my $temp = Math::GMPz->new(0);
-
- for ( my $k = 0 ; $k < $polymod ; $k++ ) {
- Rmpz_set_ui($sum, 0);
- for ( my $i = 0 ; $i <= $k ; $i++ ) {
- # sum += x[i] * (y[k-i] + y[k+polymod-i])
- # it turns out one, two, or three of the terms above are typically
- # zero, so only do as much work as needed.
- next unless $xset[$i];
- if ($yset[$k-$i]) {
- if ($yset[$k+$polymod-$i]) {
- Rmpz_add( $temp, $y[$k-$i], $y[$k+$polymod-$i] );
- Rmpz_addmul( $sum, $x[$i], $temp );
- } else {
- Rmpz_addmul( $sum, $x[$i], $y[$k-$i] );
- }
- } else {
- if ($yset[$k+$polymod-$i]) {
- Rmpz_addmul( $sum, $x[$i], $y[$k+$polymod-$i] );
- } else {
- # Nothing
- }
- }
- }
-
- # Only do this loop if there are any non-zero components to multiply.
- if ( ($k+1) <= $xdeg ) {
- for ( my $i = $k+1 ; $i <= ( $k + $polymod ) ; $i++ ) {
- next unless $xset[$i] && $yset[$k+$polymod-$i];
- Rmpz_addmul( $sum, $x[$i], $y[$k+$polymod-$i] );
- }
- }
-
- Rmpz_mod( $temp, $sum, $mod );
- $self->setCoef( $temp, $k );
+ my @res = map { Math::GMPz->new(0) } 0 .. $polymod-1;
- last if $k > $maxdeg && $sum == 0;
+ for (my $ix = 0; $ix <= $xdeg; $ix++) {
+ next unless $xset[$ix];
+ for (my $iy = 0; $iy <= $ydeg; $iy++) {
+ next unless $yset[$iy];
+ Rmpz_addmul( $res[ ($ix + $iy) % $polymod ], $x[$ix], $y[$iy] );
+ }
}
+ $self->{COEF} = [ map { $_ % $mod } @res ];
$self->compact();
}
View
33 t/big_polynomial.t
@@ -11,6 +11,9 @@ BEGIN {
use_ok ('Math::Primality::BigPolynomial' );
}
+sub new_poly {
+ return Math::Primality::BigPolynomial->new([ map { Math::GMPz->new($_) } @_ ]);
+}
my $b = Math::Primality::BigPolynomial->new([1,3,7]);
isa_ok($b, 'Math::Primality::BigPolynomial');
@@ -23,4 +26,34 @@ cmp_ok($b->getCoef(3),'==', Math::GMPz->new(0), 'coef(3) is 0');
is($b->getCoef(-1), undef, 'coef(-1) is undef');
is($b->getCoef(-10), undef, 'coef(-10) is undef');
+# Now do some testing for multiplication
+{
+ my $p1 = new_poly(3,1);
+ my $p2 = new_poly(2,1);
+ is($p1->copy()->mulmod($p2,1024,1024)->string,
+ "x^2 + 5x + 6",
+ "(x+3) * (x+2) = x^2 + 5x + 6");
+}
+{
+ my $p1 = new_poly(4, 0, 2, 1);
+ my $p2 = new_poly(1, 1, 0, 2);
+ is($p1->copy()->mulmod($p2,1024,1024)->string,
+ "2x^6 + 4x^5 + x^4 + 11x^3 + 2x^2 + 4x + 4",
+ "(x^3 + 2x^2 + 4) * (2x^3 + x + 1)");
+ is($p1->copy()->powmod(3, 1024, 1024)->string,
+ "x^9 + 6x^8 + 12x^7 + 20x^6 + 48x^5 + 48x^4 + 48x^3 + 96x^2 + 64",
+ "(x^3 + 2x^2 + 4) ^ 3");
+}
+{
+ # http://reference.wolfram.com/mathematica/tutorial/PolynomialsModuloPrimes.html
+ my $p = new_poly(1,1);
+ is($p->copy->powmod(6, 1024, 1024)->string,
+ "x^6 + 6x^5 + 15x^4 + 20x^3 + 15x^2 + 6x + 1",
+ "(1+x)^6");
+ is($p->copy->powmod(6, 2, 1024)->string,
+ "x^6 + x^4 + x^2 + 1",
+ "(1+x)^6 modulo 2");
+}
+
done_testing;
+
Please sign in to comment.
Something went wrong with that request. Please try again.