Math-Prime-Util
view release on metacpan or search on metacpan
Faster and smaller. Most primality and factoring code 10% faster.
- Speedup for factoring by running more Pollard-Rho-Brent, revising
SQUFOF, updating HOLF, updating recipe.
0.67 2017-09-23
[ADDED]
- lastfor stops forprimes (etc.) iterations
- is_square(n) returns 1 if n is a perfect square
- is_polygonal(n,k) returns 1 if n is a k-gonal number
[FUNCTIONALITY AND PERFORMANCE]
- shuffle prototype is @ instead of ;@, so matches List::Util.
- On Perl 5.8 and earlier we will call PP instead of trying
direct-to-GMP. Works around a bug in XS trying to turn the
result into an object where 5.8.7 and earlier gets lost.
if (items == 1)
XSRETURN_EMPTY;
randcxt = MY_CXT.randcxt;
/*
* Fisher-Yates shuffle with first 'k' selections returned.
*
* There is only one algorithm here, no shortcuts other than
* detecting an empty list.
*
* With a list input, the input is on the stack ST(1),ST(2),...
* We move the last item to ST(0) then shuffle 'k' iterations.
*
* With an array reference input, we cannot modify the input at all.
* We create an index array and shuffle using that. Remembering to
* act like the last item is at the front so we match the list results.
* We optimize by pushing each selection onto the return stack as
* we find it rather than pushing them all at the end with another loop.
*/
if (items > 2 || !SvROK(ST(1)) || SvTYPE(SvRV(ST(1))) != SVt_PVAV) {
/* Standard form, where we are given an array of items */
nitems = items-1;
UV rleft = (r > rounds) ? rounds : r;
UV saveXi = Xi;
/* Do rleft rounds, inner at a time */
while (rleft > 0) {
UV dorounds = (rleft > inner) ? inner : rleft;
saveXi = Xi;
rleft -= dorounds;
rounds -= dorounds;
Xi = sqraddmod(Xi, a, n); /* First iteration, no mulmod needed */
m = ABSDIFF(Xi,Xm);
while (--dorounds > 0) { /* Now do inner-1=63 more iterations */
Xi = sqraddmod(Xi, a, n);
f = ABSDIFF(Xi,Xm);
m = mulmod(m, f, n);
}
f = gcd_ui(m, n);
if (f != 1)
break;
}
/* If f == 1, then we didn't find a factor. Move on. */
if (f == 1) {
inverse_interpolate.c view on Meta::CPAN
for (iloopc = 1; (hi-lo) > 1 && rhi > n; iloopc++) {
MPUassert(lo < mid && mid < hi, "interpolation: assume 3 unique points");
rmid = MPU_CALLBACK(mid);
if (rmid >= n) { hi = mid; rhi = rmid; }
else { lo = mid; rlo = rmid; }
if (rhi == n) break;
mid += (IV)(((double)n-(double)rmid)*(double)(hi-lo) / (double)(rhi-rlo));
/* Sometimes we get stuck getting closer and closer but not bracketing.
* We could do Ridder's method of alternating bisection, or using a
* multiplier on mid on alternate iterations to reflect about n.
* What we're going to do instead is, every few loops, check if we're
* very close to one of the edges and try to pull in the other edge.
*/
if ((iloopc % 6) == 0) {
UV close = .003*(hi-lo) + 1.0;
if (lo+close > mid) mid = lo+close;
else if (hi-close < mid) mid = hi-close;
}
/* Alternately:
if (mid == lo) { mid = lo + .01*(hi-lo); }
lib/Math/Prime/Util.pm view on Meta::CPAN
=head2 nth_twin_prime_approx
Returns an approximation to the Nth twin prime. A curve fit is used for
small inputs (under 1200), while for larger inputs a binary search is done
on the approximate twin prime count.
=head2 nth_semiprime
Returns the Nth semiprime, similar to where a C<forsemiprimes> loop would
end after C<N> iterations, but much more efficiently.
=head2 nth_semiprime_approx
Returns an approximation to the Nth semiprime. The approximation is
orders of magnitude better than the simple C<n log n / log log n>
approximation for large C<n>. E.g. for C<n=10^12> the simple estimate
is within 3.6%, but this function is within 0.000012%.
=head2 nth_almost_prime
lib/Math/Prime/Util.pm view on Meta::CPAN
Given a single non-negative integer C<n>, returns a list containing each C<p>
for all prime pairs C<p> and C<q> where C<< p <= q >> and C<p + q = n>.
The number of elements returned is the same as L</goldbach_pair_count>.
If no such pairs exist, an empty list is returned.
=head2 is_happy
Given a single non-negative integer C<n>, returns the number of iterations
required for the map of sum of squared base-10 digits to converge to C<1>,
or C<0> if it does not converge to the value C<1>.
This returns the height using the OEIS A090425 definition of height, which is
zero for non-happy numbers, 1 for C<n=1>, 2 for numbers that produce 1 after
a single iteration, etc.
This is one more than the definitions used in many papers
(e.g. Cai and Zhou 2008) where C<n=1> is considered to have height 0.
An optional base and exponent may be given (default base 10 exponent 2).
lib/Math/Prime/Util/RandomPrimes.pm view on Meta::CPAN
my ($isp, $cert) = is_provable_prime_with_cert($c);
return (1,$c,$seed,$prime_gen_counter,$cert) if $isp;
return (0,0,0,0) if $prime_gen_counter > 10000 + 16*$k;
}
}
my($status,$c0,$seed,$prime_gen_counter,$cert)
= _ST_Random_prime( (($k+1)>>1)+1, $input_seed);
return (0,0,0,0) unless $status;
$cert = cmpint($c0,"18446744073709551615") <= 0
? "" : _strip_proof_header($cert);
my $iterations = int(($k + 255) / 256) - 1; # SHA256 generates 256 bits
my $old_counter = $prime_gen_counter;
my $c02 = lshiftint($c0); # $c02 = 2*$c0
my $xstr = '';
for my $i (0 .. $iterations) {
$xstr = Digest::SHA::sha256_hex($seed) . $xstr;
$seed = _seed_plus_one($seed);
}
my $x = fromdigits($xstr,16);
$x = $k2 + ($x % $k2);
my $t = cdivint($x, $c02);
_make_big_gcds() if $_big_gcd_use < 0;
while (1) {
my $c = add1int(mulint($t,$c02));
if ($c > 2*$k2) {
lib/Math/Prime/Util/RandomPrimes.pm view on Meta::CPAN
gcd($c, $_big_gcd[1]) == 1 &&
gcd($c, $_big_gcd[2]) == 1 &&
gcd($c, $_big_gcd[3]) == 1;
}
$looks_prime = 0 if $looks_prime && !is_strong_pseudoprime($c, 3);
}
if ($looks_prime) {
# We could use a in (2,3,5,7,11,13), but pedantically use FIPS 186-4.
my $astr = '';
for my $i (0 .. $iterations) {
$astr = Digest::SHA::sha256_hex($seed) . $astr;
$seed = _seed_plus_one($seed);
}
my $a = fromdigits($astr,16);
$a = addint(modint($a,$c-3),2);
my $z = powmod($a, lshiftint($t), $c);
if (gcd($z-1,$c) == 1 && powmod($z, $c0, $c) == 1) {
croak "Shawe-Taylor random prime failure at ($k): $c not prime"
unless is_prob_prime($c);
$cert = "[MPU - Primality Certificate]\nVersion 1.0\n\n" .
"Proof for:\nN $c\n\n" .
"Type Pocklington\nN $c\nQ $c0\nA $a\n" .
$cert;
return (1, $c, $seed, $prime_gen_counter, $cert);
}
} else {
# Update seed "as if" we performed the Pocklington check from FIPS 186-4
for my $i (0 .. $iterations) {
$seed = _seed_plus_one($seed);
}
}
return (0,0,0,0) if $prime_gen_counter > 10000 + 16*$k + $old_counter;
$t++;
}
}
sub random_safe_prime {
my $bits = int("$_[0]");
long double lastw = w;
for (i = 0; i < 100; i++) {
long double ew = expl(w);
long double wew = w * ew;
long double wewx = wew - x;
long double w1 = w + 1;
w = w - wewx / (ew * w1 - (w+2) * wewx/(2*w1));
if (w != 0.0L && fabsl((w-lastw)/w) <= 8*LDBL_EPSILON) break;
lastw = w;
}
#else /* Fritsch, see Veberic 2009. 1-2 iterations are enough. */
for (i = 0; i < 6 && w != 0.0L; i++) {
long double w1 = 1 + w;
long double zn = logl((long double)x/w) - w;
long double qn = 2 * w1 * (w1+(2.0L/3.0L)*zn);
long double en = (zn/w1) * (qn-zn)/(qn-2.0L*zn);
/* w *= 1.0L + en; if (fabsl(en) <= 16*LDBL_EPSILON) break; */
long double wen = w * en;
if (isnan(wen)) return 0;
w += wen;
if (fabsl(wen) <= 64*LDBL_EPSILON) break;
t/34-random.t view on Meta::CPAN
my $bits_on = 0;
my $bits_off = 0;
my $iter = 0;
for (1 .. 6400) {
$iter++;
my $v = irand64;
$bits_on |= $v;
$bits_off |= (~$v);
last if ~$bits_on == 0 && ~$bits_off == 0;
}
is( ~$bits_on, 0, "irand64 all bits on in $iter iterations" );
is( ~$bits_off, 0, "irand64 all bits off in $iter iterations" );
}
########
# This is really brute force, but it doesn't take too long.
{
my $mask = 0;
my $v;
for (1..1024) {
$v = drand;
xt/allrootmod.pl view on Meta::CPAN
die "xs rootmod fail for $a,$k,$n" unless $roots{$rxs};
die "pp rootmod fail for $a,$k,$n" unless $roots{$rpp};
if ($k == 2) {
my $sxs = sqrtmod($a,$n);
my $spp = Math::Prime::Util::PP::sqrtmod($a,$n);
die "xs sqrtmod fail for $a,$n" unless $roots{$sxs};
die "pp sqrtmod fail for $a,$n" unless $roots{$spp};
}
}
}
say "pass $iter iterations with max n $size for k $k";
}
( run in 0.469 second using v1.01-cache-2.11-cpan-71847e10f99 )