[libmath-prime-util-perl] 04/06: Try another 64-bit detection method
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:20 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.04
in repository libmath-prime-util-perl.
commit 8ebc5ae46fe4b7341665d5fe1a941bc20760143a
Author: Dana Jacobsen <dana at acm.org>
Date: Thu Jun 7 22:59:36 2012 -0600
Try another 64-bit detection method
---
TODO | 2 ++
examples/bench-random-prime.pl | 5 +++--
factor.c | 28 +++++++++++++---------------
factor.h | 2 +-
lib/Math/Prime/Util.pm | 15 +++++++++------
ptypes.h | 12 ++++++++++++
util.h | 6 +++---
7 files changed, 43 insertions(+), 27 deletions(-)
diff --git a/TODO b/TODO
index 840db11..e3aa794 100644
--- a/TODO
+++ b/TODO
@@ -20,3 +20,5 @@
- prime count bounds are <=, nth_prime bounds are <. Pick one.
- test file for random_prime
+
+- speed up random_prime for large numbers
diff --git a/examples/bench-random-prime.pl b/examples/bench-random-prime.pl
index 3e5c9e9..e60a221 100755
--- a/examples/bench-random-prime.pl
+++ b/examples/bench-random-prime.pl
@@ -6,15 +6,16 @@ use Math::Prime::Util qw/random_prime random_ndigit_prime/;
use Benchmark qw/:all/;
use List::Util qw/min max/;
my $count = shift || -3;
+my $maxdigits = (~0 <= 4294967295) ? 10 : 20;
srand(29);
-test_at_digits($_) for (2 .. 10);
+test_at_digits($_) for (2 .. $maxdigits);
sub test_at_digits {
my $digits = shift;
die "Digits must be > 0" unless $digits > 0;
cmpthese($count,{
- "$digits digits" => sub { random_ndigit_prime($digits) for (1..10000) },
+ "$digits digits" => sub { random_ndigit_prime($digits) for (1..1000) },
});
}
diff --git a/factor.c b/factor.c
index 76a7ba5..fef8b55 100644
--- a/factor.c
+++ b/factor.c
@@ -1,15 +1,13 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <limits.h>
-#define __STDC_LIMIT_MACROS
-#include <stdint.h>
-#include <math.h>
+//#include <stdio.h>
+//#include <stdlib.h>
+//#include <string.h>
+//#include <limits.h>
+//#include <math.h>
+#include "ptypes.h"
#include "factor.h"
#include "util.h"
#include "sieve.h"
-#include "ptypes.h"
/*
* You need to remember to use UV for unsigned and IV for signed types that
@@ -79,7 +77,7 @@ int trial_factor(UV n, UV *factors, UV maxtrial)
}
-#if UINT64_MAX > UV_MAX
+#if (BITS_PER_WORD == 32) && HAVE_STD_U64
/* We have 64-bit available, but UV is 32-bit. Do the math in 64-bit.
* Even if it is emulated, it should be as fast or faster than us doing it.
@@ -162,7 +160,7 @@ int trial_factor(UV n, UV *factors, UV maxtrial)
* Returns 1 if probably prime relative to the bases, 0 if composite.
* Bases must be between 2 and n-2
*/
-int miller_rabin(UV n, const UV *bases, UV nbases)
+int miller_rabin(UV n, const UV *bases, int nbases)
{
int b;
int s = 0;
@@ -313,7 +311,7 @@ int fermat_factor(UV n, UV *factors, UV rounds)
} while (r > 0);
}
r = (x-y)/2;
- if ( (r != 1) && (r != n) ) {
+ if ( (r != 1) && ((UV)r != n) ) {
factors[0] = r;
factors[1] = n/r;
MPUassert( factors[0] * factors[1] == n , "incorrect factoring");
@@ -477,7 +475,7 @@ static void enqu(IV q, IV *iter) {
int squfof_factor(UV n, UV *factors, UV rounds)
{
- UV rounds2 = rounds/16;
+ IV rounds2 = (IV) (rounds/16);
UV temp;
IV iq,ll,l2,p,pnext,q,qlast,r,s,t,i;
IV jter, iter;
@@ -507,7 +505,7 @@ int squfof_factor(UV n, UV *factors, UV rounds)
/* the end of the loop. Is there a faster way? The current way is */
/* EXPENSIVE! (many branches and double prec sqrt) */
- for (jter=0; jter < rounds; jter++) {
+ for (jter=0; (UV)jter < rounds; jter++) {
iq = (s + p)/q;
pnext = iq*q - p;
if (q <= ll) {
@@ -538,7 +536,7 @@ int squfof_factor(UV n, UV *factors, UV rounds)
break;
} /* end of main loop */
- if (jter >= rounds) {
+ if ((UV)jter >= rounds) {
factors[0] = n; return 1;
}
@@ -582,7 +580,7 @@ int squfof_factor(UV n, UV *factors, UV rounds)
if ((q & 1) == 0) q/=2; /* q was factor or 2*factor */
- if ( (q == 1) || (q == n) ) {
+ if ( (q == 1) || ((UV)q == n) ) {
factors[0] = n; return 1;
}
diff --git a/factor.h b/factor.h
index ad0a757..d474203 100644
--- a/factor.h
+++ b/factor.h
@@ -19,7 +19,7 @@ extern int prho_factor(UV n, UV *factors, UV maxrounds);
extern int pminus1_factor(UV n, UV *factors, UV maxrounds);
-extern int miller_rabin(UV n, const UV *bases, UV nbases);
+extern int miller_rabin(UV n, const UV *bases, int nbases);
extern int is_prob_prime(UV n);
static UV gcd_ui(UV x, UV y) {
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 357400e..01dfb85 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -109,6 +109,7 @@ sub random_prime {
$low = 2 if $low < 2;
# Make sure we have a valid range.
+ # TODO: this is is killing performance with large numbers
$low = next_prime($low-1);
$high = ($high < ~0) ? prev_prime($high+1) : prev_prime($high);
return $low if ($low == $high) && is_prime($low);
@@ -463,9 +464,10 @@ to a BPSW test, depending on speed.
my $small_prime = random_prime(1000); # random prime <= limit
my $rand_prime = random_prime(100, 10000); # random prime within a range
-Returns a randomly selected prime that will be greater than or equal to the
-lower limit and less than or equal to the upper limit. If no lower limit is
-given, 2 is implied. Returns undef if no primes exist within the range.
+Returns a psuedo-randomly selected prime that will be greater than or equal
+to the lower limit and less than or equal to the upper limit. If no lower
+limit is given, 2 is implied. Returns undef if no primes exist within the
+range. The L<rand> function is called one or more times for selection.
This will return a uniform distribution of the primes in the range, meaning
for each prime in the range, the chances are equally likely that it will be
@@ -474,9 +476,10 @@ seen.
The current algorithm does a random index selection for small numbers, which
is deterministic. For larger numbers, this can be very slow, so the
obvious Monte Carlo method is used, where random numbers in the range are
-selected until one is prime.
-
-The L<rand> function is called one or more times for selection.
+selected until one is prime. That also gets slow as the number of digits
+increases, but not something that impacts us in 64-bit. A GMP version would
+likely need to consider using C<next_prime>/C<prev_prime> -- giving up a
+small bit of uniformity for reasonable performance.
=head2 random_ndigit_prime
diff --git a/ptypes.h b/ptypes.h
index cd4a2e2..e7f604a 100644
--- a/ptypes.h
+++ b/ptypes.h
@@ -1,6 +1,9 @@
#ifndef MPU_PTYPES_H
#define MPU_PTYPES_H
+#define __STDC_LIMIT_MACROS
+#include <stdint.h>
+
#include "EXTERN.h"
#include "perl.h"
@@ -37,6 +40,15 @@
#define UVCONST(x) U32_CONST(x)
#endif
+/* Try to determine if we have 64-bit available via uint64_t */
+#define HAVE_STD_U64 0
+#if defined(UINT64_MAX) && defined(__UINT64_C)
+ #if (UINT64_MAX >= __UINT64_C(18446744073709551615))
+ #undef HAVE_STD_U64
+ #define HAVE_STD_U64 1
+ #endif
+#endif
+
#define MAXBIT (BITS_PER_WORD-1)
#define NWORDS(bits) ( ((bits)+BITS_PER_WORD-1) / BITS_PER_WORD )
#define NBYTES(bits) ( ((bits)+8-1) / 8 )
diff --git a/util.h b/util.h
index 0ecd27c..b14c175 100644
--- a/util.h
+++ b/util.h
@@ -25,10 +25,10 @@ extern void free_prime_segment(void);
/* Above this value, is_prime will do deterministic Miller-Rabin */
/* With 64-bit math, we can do much faster mulmods from 2^16-2^32 */
-#if (BITS_PER_WORD == 32) && (UINT64_MAX <= UV_MAX)
-#define MPU_PROB_PRIME_BEST UVCONST(100000000)
+#if (BITS_PER_WORD == 64) || HAVE_STD_U64
+ #define MPU_PROB_PRIME_BEST UVCONST(100000)
#else
-#define MPU_PROB_PRIME_BEST UVCONST(100000)
+ #define MPU_PROB_PRIME_BEST UVCONST(100000000)
#endif
#endif
--
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