[libmath-prime-util-perl] 11/55: Add lucas_sequence tests
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:53:39 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.41
in repository libmath-prime-util-perl.
commit 3e76756552e7477ac06ab010330f71b4cb284d42
Author: Dana Jacobsen <dana at acm.org>
Date: Tue Apr 29 16:16:55 2014 -0700
Add lucas_sequence tests
---
Changes | 3 ++
MANIFEST | 1 +
lib/Math/Prime/Util.pm | 5 ++-
lib/Math/Prime/Util/PP.pm | 12 +++++--
primality.c | 8 ++++-
t/25-lucas_sequences.t | 79 +++++++++++++++++++++++++++++++++++++++++++++++
6 files changed, 101 insertions(+), 7 deletions(-)
diff --git a/Changes b/Changes
index 18ded9b..57a5773 100644
--- a/Changes
+++ b/Changes
@@ -12,6 +12,9 @@ Revision history for Perl module Math::Prime::Util
- Better AKS testing in xt/primality-aks.pl.
+ - Loosen requirements of lucas_sequence. More useful for general seqs.
+ Add tests for some common sequences.
+
0.40 2014-04-21
diff --git a/MANIFEST b/MANIFEST
index a5ed401..3d83c15 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -107,6 +107,7 @@ t/22-aks-prime.t
t/23-primality-proofs.t
t/23-random-certs.t
t/24-partitions.t
+t/25-lucas_sequences.t
t/30-relations.t
t/31-threading.t
t/32-iterators.t
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 40ea2f2..793fd32 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -1859,9 +1859,8 @@ Computes C<U_k>, C<V_k>, and C<Q_k> for the Lucas sequence defined by
C<P>,C<Q>, modulo C<n>. The modular Lucas sequence is used in a
number of primality tests and proofs.
The following conditions must hold:
-C< D = P*P - 4*Q != 0> ;
-C< 0 E<lt> P E<lt> n> ;
-C< Q E<lt> n> ;
+C< |P| E<lt> n> ;
+C< |Q| E<lt> n> ;
C< k E<gt>= 0> ;
C< n E<gt>= 2>.
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 0ebbadd..0ad7120 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1879,8 +1879,8 @@ sub lucas_sequence {
croak "lucas_sequence: n must be >= 2" if $n < 2;
croak "lucas_sequence: k must be >= 0" if $k < 0;
- croak "lucas_sequence: P out of range" if $P < 0 || $P >= $n;
- croak "lucas_sequence: Q out of range" if $Q >= $n;
+ croak "lucas_sequence: P out of range" if abs($P) >= $n;
+ croak "lucas_sequence: Q out of range" if abs($Q) >= $n;
$n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
@@ -1888,7 +1888,13 @@ sub lucas_sequence {
$P = $ZERO+$P unless ref($P) eq 'Math::BigInt';
$Q = $ZERO+$Q unless ref($Q) eq 'Math::BigInt';
my $D = $P*$P - BTWO*BTWO*$Q;
- croak "lucas_sequence: D is zero" if $D->is_zero;
+ if ($D->is_zero) {
+ my $S = ($ZERO+$P) >> 1;
+ my $U = $S->copy->bmodpow($k-1,$n)->bmul($k)->bmod($n);
+ my $V = $S->copy->bmodpow($k,$n)->bmul(BTWO)->bmod($n);
+ my $Qk = ($ZERO+$Q)->bmodpow($k, $n);
+ return ($U, $V, $Qk);
+ }
my $U = BONE->copy;
my $V = $P->copy;
my $Qk = $Q->copy;
diff --git a/primality.c b/primality.c
index 0f9b221..38d3e7f 100644
--- a/primality.c
+++ b/primality.c
@@ -331,7 +331,13 @@ void lucas_seq(UV* Uret, UV* Vret, UV* Qkret, UV n, IV P, IV Q, UV k)
Qmod = (Q < 0) ? (UV)(Q + (IV)n) : (UV)Q;
Pmod = (P < 0) ? (UV)(P + (IV)n) : (UV)P;
Dmod = submod( mulmod(Pmod, Pmod, n), mulmod(4, Qmod, n), n);
- MPUassert(Dmod != 0, "lucas_seq: D is 0");
+ if (Dmod == 0) {
+ b = Pmod >> 1;
+ *Uret = mulmod(k, powmod(b, k-1, n), n);
+ *Vret = mulmod(2, powmod(b, k, n), n);
+ *Qkret = powmod(Qmod, k, n);
+ return;
+ }
U = 1;
V = Pmod;
Qk = Qmod;
diff --git a/t/25-lucas_sequences.t b/t/25-lucas_sequences.t
new file mode 100644
index 0000000..41efed7
--- /dev/null
+++ b/t/25-lucas_sequences.t
@@ -0,0 +1,79 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+use Test::More;
+use Math::Prime::Util qw/lucas_sequence/;
+
+#my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
+my $usexs = Math::Prime::Util::prime_get_config->{'xs'};
+my $usegmp = Math::Prime::Util::prime_get_config->{'gmp'};
+
+# Values taken from the OEIS pages.
+my @lucas_seqs = (
+ [ [1, -1], 0, "U", "Fibonacci numbers",
+ [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610] ],
+ [ [1, -1], 0, "V", "Lucas numbers",
+ [2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, 199, 322, 521, 843] ],
+ [ [2, -1], 0, "U", "Pell numbers",
+ [0, 1, 2, 5, 12, 29, 70, 169, 408, 985, 2378, 5741, 13860, 33461, 80782] ],
+ [ [2, -1], 0, "V", "Pell-Lucas numbers",
+ [2, 2, 6, 14, 34, 82, 198, 478, 1154, 2786, 6726, 16238, 39202, 94642] ],
+ [ [1, -2], 0, "U", "Jacobsthal numbers",
+ [0, 1, 1, 3, 5, 11, 21, 43, 85, 171, 341, 683, 1365, 2731, 5461, 10923] ],
+ [ [1, -2], 0, "V", "Jacobsthal-Lucas numbers",
+ [2, 1, 5, 7, 17, 31, 65, 127, 257, 511, 1025, 2047, 4097, 8191, 16385] ],
+ [ [2, 2], 1, "U", "sin(x)*exp(x)",
+ [0, 1, 2, 2, 0, -4, -8, -8, 0, 16, 32, 32, 0, -64, -128, -128, 0, 256] ],
+ [ [2, 2], 1, "V", "offset sin(x)*exp(x)",
+ [2, 2, 0, -4, -8, -8, 0, 16, 32, 32, 0, -64, -128, -128, 0, 256, 512,512] ],
+ [ [2, 5], 1, "U", "A045873",
+ [0, 1, 2, -1, -12, -19, 22, 139, 168, -359, -1558, -1321, 5148, 16901] ],
+ [ [3,-5], 0, "U", "3*a(n-1)+5*a(n-2) [0,1]",
+ [0, 1, 3, 14, 57, 241, 1008, 4229, 17727, 74326, 311613, 1306469] ],
+ [ [3,-5], 0, "V", "3*a(n-1)+5*a(n-2) [2,3]",
+ [2, 3, 19, 72, 311, 1293, 5434, 22767, 95471, 400248, 1678099, 7035537] ],
+ [ [3,-4], 0, "U", "3*a(n-1)+4*a(n-2) [0,1]",
+ [0, 1, 3, 13, 51, 205, 819, 3277, 13107, 52429, 209715, 838861, 3355443] ],
+ [ [3,-4], 0, "V", "3*a(n-1)+4*a(n-2) [2,3]",
+ [2, 3, 17, 63, 257, 1023, 4097, 16383, 65537, 262143, 1048577, 4194303] ],
+ [ [3,-1], 0, "U", "A006190",
+ [0, 1, 3, 10, 33, 109, 360, 1189, 3927, 12970, 42837, 141481, 467280] ],
+ [ [3,-1], 0, "V", "A006497",
+ [2, 3, 11, 36, 119, 393, 1298, 4287, 14159, 46764, 154451, 510117,1684802]],
+ [ [3, 1], 0, "U", "Fibonacci(2n)",
+ [0, 1, 3, 8, 21, 55, 144, 377, 987, 2584, 6765, 17711, 46368, 121393] ],
+ [ [3, 1], 0, "V", "Lucas(2n)",
+ [2, 3, 7, 18, 47, 123, 322, 843, 2207, 5778, 15127, 39603, 103682, 271443]],
+ [ [3, 2], 0, "U", "2^n-1 Mersenne numbers (prime and composite)",
+ [0, 1, 3, 7, 15, 31, 63, 127, 255, 511, 1023, 2047, 4095, 8191, 16383] ],
+ [ [3, 2], 0, "V", "2^n+1",
+ [2, 3, 5, 9, 17, 33, 65, 129, 257, 513, 1025, 2049, 4097, 8193, 16385] ],
+ [ [4,-1], 0, "U", "Denominators of continued fraction convergents to sqrt(5)",
+ [0, 1, 4, 17, 72, 305, 1292, 5473, 23184, 98209, 416020, 1762289, 7465176]],
+ [ [4,-1], 0, "V", "Even Lucas numbers Lucas(3n)",
+ [2, 4, 18, 76, 322, 1364, 5778, 24476, 103682, 439204, 1860498, 7881196] ],
+ [ [4, 1], 0, "U", "A001353",
+ [0, 1, 4, 15, 56, 209, 780, 2911, 10864, 40545, 151316, 564719, 2107560] ],
+ [ [4, 1], 0, "V", "A003500",
+ [2, 4, 14, 52, 194, 724, 2702, 10084, 37634, 140452, 524174, 1956244] ],
+ [ [5, 4], 0, "U", "(4^n-1)/3",
+ [0, 1, 5, 21, 85, 341, 1365, 5461, 21845, 87381, 349525, 1398101, 5592405]],
+);
+
+# 4,4 has D=0. Old GMP won't handle that.
+if ($usexs || !$usegmp) {
+ push @lucas_seqs,
+ [ [4, 4], 0, "U", "n*2^(n-1)",
+ [0, 1, 4, 12, 32, 80, 192, 448, 1024, 2304, 5120, 11264, 24576, 53248] ],
+}
+
+plan tests => 0 + scalar(@lucas_seqs);
+
+foreach my $seqs (@lucas_seqs) {
+ my($apq, $isneg, $uorv, $name, $exp) = @$seqs;
+ my $idx = ($uorv eq 'U') ? 0 : 1;
+ my @seq = map { (lucas_sequence(2**32-1, @$apq, $_))[$idx] } 0 .. $#$exp;
+ do { for (@seq) { $_ -= (2**32-1) if $_ > 2**31; } } if $isneg;
+ is_deeply( [@seq], $exp, "${uorv}_n(@$apq) -- $name" );
+}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git
More information about the Pkg-perl-cvs-commits
mailing list