/usr/lib/x86_64-linux-gnu/perl5/5.24/Math/Prime/Util/GMP.pm is in libmath-prime-util-gmp-perl 0.41-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 1466 1467 1468 1469 1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1535 1536 1537 1538 1539 1540 1541 1542 1543 1544 1545 1546 1547 1548 1549 1550 1551 1552 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562 1563 1564 1565 1566 1567 1568 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581 1582 1583 1584 1585 1586 1587 1588 1589 1590 1591 1592 1593 1594 1595 1596 1597 1598 1599 1600 1601 1602 1603 1604 1605 1606 1607 1608 1609 1610 1611 1612 1613 1614 1615 1616 1617 1618 1619 1620 1621 1622 1623 1624 1625 1626 1627 1628 1629 1630 1631 1632 1633 1634 1635 1636 1637 1638 1639 1640 1641 1642 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 1670 1671 1672 1673 1674 1675 1676 1677 1678 1679 1680 1681 1682 1683 1684 1685 1686 1687 1688 1689 1690 1691 1692 1693 1694 1695 1696 1697 1698 1699 1700 1701 1702 1703 1704 1705 1706 1707 1708 1709 1710 1711 1712 1713 1714 1715 1716 1717 1718 1719 1720 1721 1722 1723 1724 1725 1726 1727 1728 1729 1730 1731 1732 1733 1734 1735 1736 1737 1738 1739 1740 1741 1742 1743 1744 1745 1746 1747 1748 1749 | package Math::Prime::Util::GMP;
use strict;
use warnings;
use Carp qw/croak confess carp/;
BEGIN {
$Math::Prime::Util::GMP::AUTHORITY = 'cpan:DANAJ';
$Math::Prime::Util::GMP::VERSION = '0.41';
}
# parent is cleaner, and in the Perl 5.10.1 / 5.12.0 core, but not earlier.
# use parent qw( Exporter );
use base qw( Exporter );
our @EXPORT_OK = qw(
is_prime
is_prob_prime
is_bpsw_prime
is_provable_prime
is_provable_prime_with_cert
is_aks_prime
is_nminus1_prime
is_nplus1_prime
is_bls75_prime
is_ecpp_prime
is_pseudoprime
is_euler_pseudoprime
is_euler_plumb_pseudoprime
is_strong_pseudoprime
is_lucas_pseudoprime
is_strong_lucas_pseudoprime
is_extra_strong_lucas_pseudoprime
is_almost_extra_strong_lucas_pseudoprime
is_perrin_pseudoprime
is_frobenius_pseudoprime
is_frobenius_underwood_pseudoprime
is_frobenius_khashin_pseudoprime
is_mersenne_prime
is_llr_prime
is_proth_prime
is_miller_prime
miller_rabin_random
lucas_sequence lucasu lucasv
primes
sieve_primes
sieve_twin_primes
sieve_prime_cluster
sieve_range
next_prime
prev_prime
surround_primes
trial_factor
prho_factor
pbrent_factor
pminus1_factor
pplus1_factor
holf_factor
squfof_factor
ecm_factor
qs_factor
factor
divisors
sigma
chinese
moebius
prime_count
primorial
pn_primorial
factorial
consecutive_integer_lcm
partitions bernfrac bernreal harmfrac harmreal stirling
zeta riemannr
gcd lcm kronecker valuation binomial gcdext
invmod sqrtmod addmod mulmod divmod powmod
vecsum vecprod
exp_mangoldt
liouville
totient
jordan_totient
carmichael_lambda
sqrtint rootint
is_power is_prime_power
is_primitive_root
znorder
znprimroot
ramanujan_tau
Pi
todigits
);
# Should add:
# nth_prime
our %EXPORT_TAGS = (all => [ @EXPORT_OK ]);
BEGIN {
eval {
require XSLoader;
XSLoader::load(__PACKAGE__, $Math::Prime::Util::GMP::VERSION);
_GMP_init();
1;
} or do {
die $@;
}
}
END {
_GMP_destroy();
}
sub _validate_positive_integer {
my($n, $min, $max) = @_;
croak "Parameter must be defined" if !defined $n;
if (ref($n) eq 'Math::BigInt' && $n->can("sign")) {
croak "Parameter '$n' must be a positive integer" unless $n->sign() eq '+';
} else {
my $sn = "$n";
croak "Parameter '$sn' must be a positive integer"
if $sn eq '' || $sn =~ tr/0123456789//c;
}
croak "Parameter '$n' must be >= $min" if defined $min && $n < $min;
croak "Parameter '$n' must be <= $max" if defined $max && $n > $max;
1;
}
sub is_provable_prime {
my ($n) = @_;
return 0 if $n < 2;
return _is_provable_prime($n);
}
sub is_provable_prime_with_cert {
my ($n) = @_;
my @composite = (0, '');
return @composite if $n < 2;
my ($result, $text) = _is_provable_prime($n, 1);
return @composite if $result == 0;
return ($result, '') if $result != 2;
$text = "Type Small\nN $n\n" if !defined $text || $text eq '';
$text =~ s/\n$//;
$text = "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN $n\n\n$text";
return ($result, $text);
}
sub primes {
my($low,$high) = (scalar(@_) == 1) ? (2,$_[0]) : ($_[0], $_[1]);
_validate_positive_integer($low);
_validate_positive_integer($high);
[ sieve_primes($low, $high, 0) ];
}
1;
__END__
# ABSTRACT: Utilities related to prime numbers and factoring, using GMP
=pod
=encoding utf8
=for stopwords Möbius Deléglise Bézout gcdext vecsum vecprod moebius totient liouville znorder znprimroot bernfrac bernreal harmfrac harmreal stirling zeta riemannr lucasu lucasv OpenPFGW gmpy2 nonresidue chinese tuplets sqrtmod addmod mulmod powmod divmod superset sqrtint rootint todigits
=head1 NAME
Math::Prime::Util::GMP - Utilities related to prime numbers and factoring, using GMP
=head1 VERSION
Version 0.41
=head1 SYNOPSIS
use Math::Prime::Util::GMP ':all';
my $n = "115792089237316195423570985008687907853269984665640564039457584007913129639937";
# This doesn't impact the operation of the module at all, but does let you
# enter big number arguments directly as well as enter (e.g.): 2**2048 + 1.
use bigint;
# These return 0 for composite, 2 for prime, and 1 for probably prime
# Numbers under 2^64 will return 0 or 2.
# is_prob_prime does a BPSW primality test for numbers > 2^64
# is_prime adds some MR tests and a quick test to try to prove the result
# is_provable_prime will spend a lot of effort on proving primality
say "$n is probably prime" if is_prob_prime($n);
say "$n is ", qw(composite prob_prime def_prime)[is_prime($n)];
say "$n is definitely prime" if is_provable_prime($n) == 2;
# Miller-Rabin and strong Lucas-Selfridge pseudoprime tests
say "$n is a prime or spsp-2/7/61" if is_strong_pseudoprime($n, 2, 7, 61);
say "$n is a prime or slpsp" if is_strong_lucas_pseudoprime($n);
say "$n is a prime or eslpsp" if is_extra_strong_lucas_pseudoprime($n);
# Return array reference to primes in a range.
my $aref = primes( 10 ** 200, 10 ** 200 + 10000 );
$next = next_prime($n); # next prime > n
$prev = prev_prime($n); # previous prime < n
# Primorials and lcm
say "23# is ", primorial(23);
say "The product of the first 47 primes is ", pn_primorial(47);
say "lcm(1..1000) is ", consecutive_integer_lcm(1000);
# Find prime factors of big numbers
@factors = factor(5465610891074107968111136514192945634873647594456118359804135903459867604844945580205745718497);
# Finer control over factoring.
# These stop after finding one factor or exceeding their limit.
# # optional arguments o1, o2, ...
@factors = trial_factor($n); # test up to o1
@factors = prho_factor($n); # no more than o1 rounds
@factors = pbrent_factor($n); # no more than o1 rounds
@factors = holf_factor($n); # no more than o1 rounds
@factors = squfof_factor($n); # no more than o1 rounds
@factors = pminus1_factor($n); # o1 = smoothness limit, o2 = stage 2 limit
@factors = ecm_factor($n); # o1 = B1, o2 = # of curves
@factors = qs_factor($n); # (no arguments)
=head1 DESCRIPTION
A module for number theory in Perl using GMP. This includes primality tests,
getting primes in a range, factoring, and more.
While it certainly can be used directly, the main purpose of this
module is for L<Math::Prime::Util>. That module will automatically
load this one if it is installed, greatly speeding up many of its
operations on big numbers.
Inputs and outputs for big numbers are via strings, so you do not need
to use a bigint package in your program. However if you do use bigints,
inputs will be converted internally so there is no need to convert
before a call. Output results are returned as either Perl scalars
(for native-size) or strings (for bigints). L<Math::Prime::Util> tries
to reconvert all strings back into the callers bigint type if possible,
which makes it more convenient for calculations.
The various C<is_*_pseudoprime> tests are more appropriately called
C<is_*_probable_prime> or C<is_*_prp>. They return 1 if the input is a
probable prime based on their test. The naming convention is historical
and follows Pari, L<Math::Primality>, and some other math packages.
The modern definition of pseudoprime is a I<composite> that passes the
test, rather than any number.
=head1 FUNCTIONS
=head2 is_prob_prime
my $prob_prime = is_prob_prime($n);
# Returns 0 (composite), 2 (prime), or 1 (probably prime)
Takes a positive number as input and returns back either 0 (composite),
2 (definitely prime), or 1 (probably prime).
For inputs below C<2^64> the test is deterministic, so the possible
return values are 0 (composite) or 2 (definitely prime).
For inputs above C<2^64>, a probabilistic test is performed. Only 0
(composite) and 1 (probably prime) are returned. The current
implementation uses the Baillie-PSW (BPSW) test. There is a
possibility that composites may be returned marked prime, but since
the test was published in 1980, not a single BPSW pseudoprime has
been found, so it is extremely likely to be prime.
While we believe (Pomerance 1984) that an infinite number of
counterexamples exist, there is a weak conjecture (Martin) that
none exist under 10000 digits.
In more detail, we are using the extra-strong Lucas test
(Grantham 2000) using the Baillie parameter selection method
(see OEIS A217719). Previous versions of this module used the
strong Lucas test with Selfridge parameters, but the extra-strong
version produces fewer pseudoprimes while running 1.2 - 1.5x faster.
It is slightly stronger than the test used in
L<Pari|http://pari.math.u-bordeaux.fr/faq.html#primetest>.
=head2 is_prime
say "$n is prime!" if is_prime($n);
Takes a positive number as input and returns back either 0 (composite),
2 (definitely prime), or 1 (probably prime). Composites will act
exactly like C<is_prob_prime>, as will numbers less than C<2^64>.
For numbers larger than C<2^64>, some additional tests are performed
on probable primes to see if they can be proven by another means.
This call walks the line between the performance of L</is_prob_prime>
and the certainty of L</is_provable_prime>. Those calls may be more
appropriate in some cases. What this function does is give most of
the performance of the former, while adding more certainty. For finer
tuning of this tradeoff, you may instead use L</is_prob_prime> followed
by additional probable prime tests such as L</miller_rabin_random>
and/or L</is_frobenius_underwood_pseudoprime>.
As with L</is_prob_prime>, a BPSW test is first performed. This is
deterministic for all 64-bit numbers. Next, if the number is a
Proth or LLR form, then a proof is constructed. If the result is
still "probably prime" and the input is smaller than the
Sorenson/Webster (2015) deterministic Miller-Rabin limit
(approximately 82 bits) then the 11 or 12 Miller-Rabin tests are
performed and the result is confirmed. For larger inputs that are
still "probably prime" but under 200 bits, a quick
BLS75 C<n-1> primality proof is attempted. This is tuned to give up
if the result cannot be quickly determined, and results in success
rates of ~80% at 80 bits, ~30% at 128 bits, and ~13% at 160 bits.
Lastly, for results still "probably prime", an additional random-base
Miller-Rabin test is performed.
The result is that many numbers will return 2 (definitely prime),
and the numbers that return 1 (probably prime) have gone through
more tests than L</is_prob_prime> while not taking too long.
For cryptographic key generation, you may want even more testing for
probable primes (NIST recommends a few more additional M-R tests than
we perform). The function L</miller_rabin_random> is made for this.
Alternately, a different test such as
L</is_frobenius_underwood_pseudoprime> can be used.
Even better, use L</is_provable_prime> which should be reasonably
fast for sizes under 2048 bits.
Typically for key generation one wants random primes, and there are
many functions for that.
=head2 is_provable_prime
say "$n is definitely prime!" if is_provable_prime($n) == 2;
Takes a positive number as input and returns back either 0 (composite),
2 (definitely prime), or 1 (probably prime). A great deal of effort is
taken to return either 0 or 2 for all numbers.
The current method first uses BPSW to find composites and provide a
deterministic answer for tiny numbers (under C<2^64>). If no
certificate is required, LLR and Proth tests can be run, and small
numbers (under approximately C<2^82>) can be satisfied with a
deterministic Miller-Rabin test. If the result is still not determined,
a quick BLS75 C<n-1> test is attempted, followed by ECPP.
The time required for primes of different input sizes on a circa-2009
workstation averages about C<3ms> for 30-digits, C<5ms> for 40-digit,
C<20ms> for 60-digit, C<50ms> for 80-digit, C<100ms> for 100-digit,
C<2s> for 200-digit, and 400-digit inputs about a minute.
Expect a lot of time variation for larger inputs. You can see progress
indication if verbose is turned on (some at level 1, and a lot at level 2).
A certificate can be obtained along with the result using the
L</is_provable_prime_with_cert> method. There is no appreciable extra
performance cost for returning a certificate.
=head2 is_provable_prime_with_cert
Takes a positive number as input and returns back an array with two
elements. The result will be one of:
(0, '') The input is composite.
(1, '') The input is probably prime but we could not prove it.
This is a failure in our ability to factor some necessary
element in a reasonable time, not a significant proof
failure (in other words, it remains a probable prime).
(2, '...') The input is prime, and the certificate contains all the
information necessary to verify this.
The certificate is a text representation containing all the necessary
information to verify the primality of the input in a reasonable time.
The result can be used with L<Math::Prime::Util/verify_prime> for
verification. Proof types used include:
ECPP
BLS3
BLS15
BLS5
Small
=head2 is_pseudoprime
Takes a positive number C<n> and one or more non-zero positive bases as input.
Returns C<1> if the input is a probable prime to each base, C<0> if not.
This is the simple Fermat primality test.
Removing primes, given base 2 this produces the sequence L<OEIS A001567|http://oeis.org/A001567>.
=head2 is_euler_pseudoprime
Takes a positive number C<n> and one or more non-zero positive bases as input.
Returns C<1> if the input is an Euler probable prime to each base, C<0> if not.
This is the Euler test, sometimes called the Euler-Jacobi test.
Removing primes, given base 2 this produces the sequence L<OEIS A047713|http://oeis.org/A047713>.
=head2 is_strong_pseudoprime
my $maybe_prime = is_strong_pseudoprime($n, 2);
my $probably_prime = is_strong_pseudoprime($n, 2, 3, 5, 7, 11, 13, 17);
Takes a positive number C<n> and one or more non-zero positive bases as input.
Returns C<1> if the input is a strong probable prime to each base, C<0> if not.
This is often called the Miller-Rabin test.
If 0 is returned, then the number really is a composite. If 1 is
returned, then it is either a prime or a strong pseudoprime to all
the given bases. Given enough distinct bases, the chances become
very strong that the number is actually prime.
Both the input number and the bases may be big integers. If base
modulo n E<lt>= 1 or base modulo n = n-1, then the result will be 1.
This allows the bases to be larger than n if desired, while still
returning meaningful results. For example,
is_strong_pseudoprime(367, 1101)
would incorrectly return 0 if this was not done properly. A 0 result
should be returned only if n is composite, regardless of the base.
This is usually used in combination with other tests to make either stronger
tests (e.g. the strong BPSW test) or deterministic results for numbers less
than some verified limit (e.g. Jaeschke showed in 1993 that no more than three
selected bases are required to give correct primality test results for any
32-bit number). Given the small chances of passing multiple bases, there
are some math packages that just use multiple MR tests for primality testing,
though in the early 1990s almost all serious software switched to the
BPSW test.
Even numbers other than 2 will always return 0 (composite). While the
algorithm works with even input, most sources define it only on odd input.
Returning composite for all non-2 even input makes the function match most
other implementations including L<Math::Primality>'s C<is_strong_pseudoprime>
function.
=head2 miller_rabin_random
my $maybe_prime = miller_rabin_random($n, 10); # 10 random bases
Takes a positive number (C<n>) as input and a positive number (C<k>) of bases
to use. Performs C<k> Miller-Rabin tests using uniform random bases
between 2 and C<n-2>. This is the correct way to perform C<k> Miller-Rabin
tests, rather than the common but broken method of using the first C<k>
primes.
An optional third argument may be given, which is a seed to use. The seed
should be a number either in decimal, binary with a leading C<0b>, hex with
a leading C<0x>, or octal with a leading C<0>. It will be converted to a
GMP integer, so may be large. Typically this is not necessary, but
cryptographic applications may prefer the ability to use this, and it
allows repeatable test results.
There is no check for duplicate bases. Input sizes below 65-bits make
little sense for this function since L<is_prob_prime> is deterministic
at that size. For numbers of 65+ bits, the chance of duplicate bases
is quite small. The exponentiation approximation for the birthday
problem gives a probability of less than 2e-16 for 100 random bases to have
a duplicate with a 65-bit input, and less than 2e-35 with a 128-bit input.
=head2 is_lucas_pseudoprime
=head2 is_strong_lucas_pseudoprime
Takes a positive number as input, and returns 1 if the input is a standard
or strong Lucas probable prime. The Selfridge method of choosing D, P, and
Q are used (some sources call this a Lucas-Selfridge test). This is one
half of the BPSW primality test (the Miller-Rabin strong probable prime test
with base 2 being the other half). The canonical BPSW test (page 1401 of
Baillie and Wagstaff (1980)) uses the strong Lucas test with Selfridge
parameters, but in practice a variety of Lucas tests with different
parameters are used by tests calling themselves BPSW.
The standard Lucas test implemented here corresponds to the Lucas test
described in FIPS 186-4 section C.3.3, though uses a slightly more
efficient calculation. Since the standard Lucas-Selfridge test is a
subset of the strong Lucas-Selfridge test, I recommend using the strong
test rather than the standard test for cryptographic purposes. It is
often slightly faster, has over 4x fewer pseudoprimes, and is the method
recommended by Baillie and Wagstaff in their 1980 paper.
=head2 is_extra_strong_lucas_pseudoprime
Takes a positive number as input, and returns 1 if the input is an
extra-strong Lucas probable prime. This is defined in Grantham (2000),
and is a slightly more stringent test than the strong Lucas test, though
because different parameters are used the pseudoprimes are not a subset.
As expected by the extra conditions, the number of pseudoprimes is less
than 2/3 that of the strong Lucas-Selfridge test.
Runtime performance is 1.2 to 1.5x faster than the strong Lucas test.
The parameters are selected using the Baillie-OEIS method:
P = 3;
Q = 1;
while ( jacobi( P*P-4, n ) != -1 )
P += 1;
=head2 is_almost_extra_strong_lucas_pseudoprime
Takes a positive number as input and returns 1 if the input is an "almost"
extra-strong Lucas probable prime. This is the classic extra-strong Lucas
test but without calculating the U sequence. This makes it very fast,
although as the input increases in size the time converges to the conventional
extra-strong implementation: at 30 digits this routine is about 15% faster,
at 300 digits it is only 2% faster.
With the current implementations, there is little reason to prefer this unless
trying to reproduce specific results. The extra-strong implementation has been
optimized to use similar features, removing most of the performance advantage.
An optional second argument (must be between 1 and 256) indicates the
increment amount for P parameter selection. The default value of one yields
the method described in L</is_extra_strong_lucas_pseudoprime>. A value of
2 yields the method used in
L<Pari|http://pari.math.u-bordeaux.fr/faq.html#primetest>.
Because the C<U = 0> condition is ignored, this produces about 5% more
pseudoprimes than the extra-strong Lucas test. However this is still only
66% of the number produced by the strong Lucas-Selfridge test. No BPSW
counterexamples have been found with any of the Lucas tests described.
=head2 is_euler_plumb_pseudoprime
Takes a positive number C<n> as input and returns 1 if C<n> passes
Colin Plumb's Euler Criterion primality test. Pseudoprimes to this test
are a subset of the base 2 Fermat and Euler tests, but a superset
of the base 2 strong pseudoprime (Miller-Rabin) test.
The main reason for this test is that is a bit more efficient
than other probable prime tests.
=head2 is_perrin_pseudoprime
Takes a positive number C<n> as input and returns 1 if C<n> divides C<P(n)>
where C<P(n)> is the Perrin number of C<n>. The Perrin sequence is defined by
C<P(n) = P(n-2) + P(n-3)> with C<P(0) = 3, P(1) = 0, P(2) = 2>.
This is not a commonly used test, as it runs slower than most
of the other probable prime tests and offers little benefit, especially over
combined tests like L</is_bpsw_prime>,
L</is_frobenius_underwood_pseudoprime>, and
L</is_frobenius_khashin_pseudoprime>.
An optional second argument C<r> indicates whether to run additional tests.
With C<r=1>, C<P(-n) = -1 mod n> is also verified, creating the
"minimal restricted" test.
With C<r=2>, the full signature is also tested using the Adams and Shanks (1982)
rules (without the quadratic form test).
With C<r=3>, the full signature is tested using the Grantham (2000) test, which
additionally does not allow pseudoprimes to be divisible by 2 or 23.
The minimal restricted pseudoprime sequence is L<OEIS A018187|http://oeis.org/A018187>.
=head2 is_frobenius_pseudoprime
Takes a positive number C<n> as input, and two optional parameters C<a> and
C<b>, and returns 1 if the C<n> is a Frobenius probable prime with respect
to the polynomial C<x^2 - ax + b>. Without the parameters, C<b = 2> and
C<a> is the least positive odd number such that C<(a^2-4b|n) = -1>.
This selection has no pseudoprimes below C<2^64> and none known. In any
case, the discriminant C<a^2-4b> must not be a perfect square.
=head2 is_frobenius_underwood_pseudoprime
Takes a positive number as input, and returns 1 if the input passes the
efficient Frobenius test of Paul Underwood. This selects a parameter C<a>
as the least non-negative integer such that C<(a^2-4|n)=-1>, then verifies that
C<(x+2)^(n+1) = 2a + 5 mod (x^2-ax+1,n)>. This combines a Fermat and Lucas
test at a computational cost of about 2.5x a strong pseudoprime test. This
makes it similar to, but faster than, a standard Frobenius test.
This test is deterministic (no randomness is used). There are no known
pseudoprimes to this test. This test also has no overlap with the BPSW
test, making it a very effective method for adding additional certainty.
=head2 is_frobenius_khashin_pseudoprime
Takes a positive number as input, and returns 1 if the input passes the
Frobenius test of Sergey Khashin. This ensures C<n> is not a perfect square,
selects the parameter C<c> as the smallest odd prime such that C<(c|n)=-1>,
then verifies that C<(1+D)^n = (1-D) mod n> where C<D = sqrt(c) mod n>.
This test is deterministic (no randomness is used).
There are no known pseudoprimes to this test.
=head2 is_bpsw_prime
Given a positive number input, returns 0 (composite), 2 (definitely prime),
or 1 (probably prime), using the BPSW primality test (extra-strong variant).
This function does the extra-strong BPSW test and nothing more. That is,
it will skip all pretests and any extra work that the L</is_prob_prime>
test may add.
=head2 is_aks_prime
say "$n is definitely prime" if is_aks_prime($n);
Takes a positive number as input, and returns 1 if the input passes the
Agrawal-Kayal-Saxena (AKS) primality test. This is a deterministic
unconditional primality test which runs in polynomial time for general input.
The particular method used is theorem 4.1 from Bernstein (2003). This is
substantially faster than the original AKS publication, the later version
with improvements by Lenstra (sometimes called the V6 paper), or the later
improvements of Voloch and Bornemann. It is, by a large order, faster than
any other known implementation as of early 2016.
For theoretical analysis of the primality task, AKS is extremely important.
In practice, it is essentially useless. Estimated run time for a 150 digit
input is over 2 days, making the case that while the algorithmic complexity
I<growth> is polynomial, the constants are extremely high. It will take
years for for numbers that ECPP or APR-CL can prove in seconds.
Typically you should use L</is_provable_prime> and let it decide the method.
=head2 is_mersenne_prime
say "2^607-1 (M607) is a Mersenne prime" if is_mersenne_prime(607);
Takes a positive number C<p> as input and returns 1 if C<2^p-1> is prime.
After some pre-testing, the Lucas-Lehmer test is performed.
This is a deterministic unconditional test that runs very fast compared
to other primality methods for numbers of comparable size, and vastly
faster than any known general-form primality proof methods.
=head2 is_llr_prime
Takes a positive number C<n> as input and returns one of: 0 (definitely
composite), 2 (definitely prime), or -1 (test does not indicate anything).
This implements the Lucas-Lehmer-Riesel test for
fast deterministic primality testing on numbers of the form C<k * 2^n - 1>.
If C<k = 1> then this is a Mersenne number and the Lucas-Lehmer test is used.
If the number is not of this form, or if C<k E<lt>= 2^n>, then C<-1> will
be returned as the test does not apply.
Otherwise, the LLR
test is performed. While not as fast as the Lucas-Lehmer test for Mersenne
numbers, it is almost as fast as a single strong pseudoprime test (i.e.
Miller-Rabin test) while giving a certain answer.
=head2 is_proth_prime
Takes a positive number C<n> as input and returns one of: 0 (definitely
composite), 2 (definitely prime), or -1 (test does not indicate anything).
This applies Proth's theorem for fast Las Vegas
primality testing on numbers of the form C<k * 2^n + 1>.
If the number is not of this form, or if C<k E<lt>= 2^n>, then C<-1> will
be returned as the test does not apply.
Otherwise, a search is performed to find a quadratic nonresidue modulo C<n>.
If none can be found after a brief search, C<-1> is returned as no conclusion
can be reached. Otherwise, Proth's theorem is checked which conclusively
indicates primality.
While not as fast as the Lucas-Lehmer test for Mersenne
numbers, it is almost as fast as a single strong pseudoprime test (i.e.
Miller-Rabin test) while giving a certain answer.
=head2 is_miller_prime
say "$n is definitely prime" if is_miller_prime($n);
say "$n is definitely prime assuming the GRH" if is_miller_prime($n, 1);
Takes a positive number as input, and returns 1 if the input passes the
deterministic Miller test. An optional second argument indicates whether
the Generalized Riemann Hypothesis should be assumed, and defaults to 0.
Setting the verbose flag to 2 or higher will show how many bases are used.
The unconditional test is exponential time, while the conditional test
(assuming the GRH) is polynomial time.
This is a very slow method in practice, and generally should not be used.
The asymptotic complexity of the GRH version is good in theory, matching
ECPP, but in practice it is much slower. The number of bases used by
the unconditional test grows quite rapidly, impractically many past about
160 bits, and overflows a 64-bit integer at 456 bits -- sizes that are
trivial for the unconditional APR-CL and ECPP tests.
=head2 is_nminus1_prime
say "$n is definitely prime" if is_nminus1_prime($n);
Takes a positive number as input, and returns 1 if the input passes either
theorem 5 or theorem 7 of the Brillhart-Lehmer-Selfridge primality test.
This is a deterministic unconditional primality test which requires factoring
C<n-1> to a linear factor less than the cube root of the input. For small
inputs (under 40 digits) this is typically very easy, and some numbers will
naturally lead to this being very fast. As the input grows, this method
slows down rapidly.
This method is most appropriate for numbers of the form C<k+1> where C<k>
can be easily factored.
Typically you should use L</is_provable_prime> and let it decide the method.
=head2 is_nplus1_prime
Takes a positive number as input, and returns 1 if the input passes either
theorem 17 or theorem 19 of the Brillhart-Lehmer-Selfridge primality test.
This is a deterministic unconditional primality test which requires factoring
C<n+1> to a linear factor less than the cube root of the input. For small
inputs (under 40 digits) this is typically very easy, and some numbers will
naturally lead to this being very fast. As the input grows, this method
slows down rapidly.
Disregarding factoring, this is slightly slower than the C<n-1> methods.
It is most appropriate for numbers of the form C<k-1> where C<k> can be
easily factored.
=head2 is_bls75_prime
Takes a positive number as input, and returns 1 if the input passes one of
the tests from the Brillhart-Lehmer-Selfridge (1975) paper. These use
partial factoring of C<n-1> and C<n+1>. Currently the implementation will
use one of:
N-1 Corollary 1, Theorem 5, Theorem 7
N+1 Corollary 8, Theorem 17, Theorem 19
Comb Theorem 20
This is appropriate for cases where either C<n-1> or C<n+1> can be easily
factored, or when both of them have many small factors.
=head2 is_ecpp_prime
say "$n is definitely prime" if is_ecpp_prime($n);
Takes a positive number as input, and returns 1 if the input passes the
ECPP primality test. This is the Atkin-Morain Elliptic Curve Primality
Proving algorithm. It is the fastest primality proving method in
Math::Prime::Util.
This implementation uses a "factor all strategy" (FAS) with backtracking.
A limited set of about 500 precalculated discriminants are used, which works
well for inputs up to 300 digits, and for many inputs up to one thousand
digits. Having a larger set will help with large numbers (a set of 2650
is available on github in the C<xt/> directory). A future implementation
may include code to generate class polynomials as needed.
Typically you should use L</is_provable_prime> and let it decide the method.
=head2 primes
my $aref1 = primes( 1_000_000 );
my $aref2 = primes( 2 ** 448, 2 ** 448 + 10000 );
say join ",", @{primes( 2**2048, 2**2048 + 10000 )};
Returns all the primes between the lower and upper limits (inclusive), with
a lower limit of C<2> if none is given.
An array reference is returned, matching the signature of the function
of the same name in L<Math::Prime::Util>.
Values above 64-bit are extra-strong BPSW probable primes.
=head2 sieve_primes
my @primes = sieve_primes(2**100, 2**100 + 10000);
my @candidates = sieve_primes(2**1000, 2**1000 + 10000, 40000);
Given two arguments C<low> and C<high>, this returns the primes in the
interval (inclusive) as a list. It operates similar to L<primes>, though
must always have an lower and upper bound and returns a list.
With three arguments C<low>, C<high>, and C<limit>, this does a partial
sieve over the inclusive range and returns the list that pass the sieve.
If C<limit> is less than C<2> then it is identical to the two-argument
version, in that a primality test will be performed after sieving.
Otherwise, sieving is performed up to C<limit>.
The two-argument version is typically only used internally and adds little
functionality. The three-argument version is quite useful for applications
that want to apply their own primality or other tests, and wish to have a
list of values in the range with no small factors. This is quite common
for applications involving prime gaps.
Also see L</sieve_range>.
=head2 sieve_range
my @candidates = sieve_range(2**1000, 10000, 40000);
Given a start value C<n>, and native unsigned integers C<width> and C<depth>,
a sieve of maximum depth C<depth> is done for the C<width> consecutive
numbers beginning with C<n>. An array of offsets from the start is returned.
The returned list contains those offsets in the range C<n> to C<n+width-1>
where C<n + offset> has no prime factors less than C<depth>.
This function is very similar to the three argument form of L</sieve_primes>.
The differences are using C<(n,width)> instead of C<(low,high)>, and most
importantly returning small offsets from the start value rather than the
values themselves. This can substantially reduce overhead for
multi-thousand digit numbers.
=head2 sieve_twin_primes
my @primes = sieve_twin_primes(2**1000, 2**1000 + 500000);
Given two arguments C<low> and C<high>, this returns each lower twin prime
in the interval (inclusive). The result is a list, not a reference.
This does a partial sieve of the range, removes any non-twin candidates,
then checks that each pair are both BPSW probable primes. This is
substantially more efficient than sieving for all primes followed by
removing those that are not twin primes.
=head2 sieve_prime_cluster
# Find some prime septuplets
my @s = sieve_prime_cluster(2**100, 2**100+1e12, 2,6,8,12,18,20);
Efficiently finds prime clusters between the first two arguments C<low>
and C<high> (inclusive). The remaining arguments describe the cluster.
The cluster values must be even, less than 31 bits, and strictly increasing.
Given a cluster set C<C>, the returned values are all primes in the
range where C<p+c> is prime for all C<c> in the cluster set C<C>.
The cluster is described as offsets from 0, with the implicit prime
at 0. Hence an empty list is asking for all primes (the cluster
C<p+0>). A list with the single value C<2> will find all twin primes
(the cluster where C<p+0> and C<p+2> are prime). The list C<2,6,8>
will find prime quadruplets. Note that there is no requirement that
the list denote a constellation (a cluster with minimal distance) --
the list C<42,92,606> is just fine.
For long clusters, e.g. L<OEIS series A213601|http://oeis.org/A213601>
prime 12-tuplets, this will be immensely more efficient than filtering
out the cluster from a list of primes. For that example, a range of
C<10^13> takes less than a second to search -- thousands of times faster
than filtering results from primes or twin primes.
Shorter clusters are not quite this efficient, and the overhead for
returning large arrays should not be ignored.
=head2 next_prime
$n = next_prime($n);
Returns the prime following the input number (the smallest prime number
that is greater than the input number).
The function L</is_prob_prime> is used to determine when a prime is found,
hence the result is a probable prime (using BPSW).
For large inputs this function is quite a bit faster than GMP's
C<mpz_nextprime> or Pari's C<nextprime>.
=head2 prev_prime
$n = prev_prime($n);
Returns the prime preceding the input number (the largest prime number
that is less than the input number).
0 is returned if the input is C<2> or lower.
The function L</is_prob_prime> is used to determine when a prime is found,
hence the result is a probable prime (using BPSW).
=head2 surround_primes
($dprev, $dnext) = surround_primes($n);
Returns the distances to the previous and next primes of the input C<n>.
This is slightly more efficient than calling both L</prev_prime> and
L</next_prime>, and returning the native integer distances can be more
efficient with large inputs.
If an optional second argument C<d> is given, and the input C<n> is larger
than C<2^64>, then if a SPSP-2 is found in the range C<n-d> to C<n+d>
(inclusive) then it will be returned with the other argument set to C<0>.
Otherwise, the first SPSP-2 values found are returned.
This feature is especially useful for prime gap searches as well as
finding the nearest prime to a value.
Note that with a non-zero second argument, the values returned have not
undergone a full BPSW test; just sieving and a SPSP-2 test.
=head2 lucasu
say "Fibonacci($_) = ", lucasu(1,-1,$_) for 0..100;
Given integers C<P>, C<Q>, and the non-negative integer C<k>,
computes C<U_k> for the Lucas sequence defined by C<P>,C<Q>. These include
the Fibonacci numbers (C<1,-1>), the Pell numbers (C<2,-1>), the Jacobsthal
numbers (C<1,-2>), the Mersenne numbers (C<3,2>), and more.
This corresponds to OpenPFGW's C<lucasU> function and gmpy2's C<lucasu>
function.
=head2 lucasv
say "Lucas($_) = ", lucasv(1,-1,$_) for 0..100;
Given integers C<P>, C<Q>, and the non-negative integer C<k>,
computes C<V_k> for the Lucas sequence defined by C<P>,C<Q>. These include
the Lucas numbers (C<1,-1>).
This corresponds to OpenPFGW's C<lucasV> function and gmpy2's C<lucasv>
function.
=head2 lucas_sequence
my($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $k)
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<< P > 0 >>
- C<< P < n >>
- C<< Q < n >>
- C<< k >= 0 >>
- C<< n >= 2 >>
=head2 primorial
$p = primorial($n);
Given an unsigned integer argument, returns the product of the prime numbers
which are less than or equal to C<n>. This definition of C<n#> follows
L<OEIS series A034386|http://oeis.org/A034386> and
L<Wikipedia: Primorial definition for natural numbers|http://en.wikipedia.org/wiki/Primorial#Definition_for_natural_numbers>.
=head2 pn_primorial
$p = pn_primorial($n)
Given an unsigned integer argument, returns the product of the first C<n>
prime numbers. This definition of C<p_n#> follows
L<OEIS series A002110|http://oeis.org/A002110> and
L<Wikipedia: Primorial definition for prime numbers|http://en.wikipedia.org/wiki/Primorial#Definition_for_prime_numbers>.
The two are related with the relationships:
pn_primorial($n) == primorial( nth_prime($n) )
primorial($n) == pn_primorial( prime_count($n) )
=head2 factorial
Given positive integer argument C<n>, returns the factorial of C<n>,
defined as the product of the integers 1 to C<n> with the special case
of C<factorial(0) = 1>. This corresponds to Pari's C<factorial(n)>
and Mathematica's C<Factorial[n]> functions.
=head2 gcd
Given a list of integers, returns the greatest common divisor. This is
often used to test for L<coprimality|https://oeis.org/wiki/Coprimality>.
=head2 lcm
Given a list of integers, returns the least common multiple.
=head2 gcdext
Given two integers C<x> and C<y>, returns C<u,v,d> such that C<d = gcd(x,y)>
and C<u*x + v*y = d>. This uses the extended Euclidian algorithm to compute
the values satisfying Bézout's Identity.
This corresponds to Pari's C<gcdext> function, which was renamed from
C<bezout> out Pari 2.6. The results will hence match L<Math::Pari/bezout>.
=head2 chinese
say chinese( [14,643], [254,419], [87,733] ); # 87041638
Solves a system of simultaneous congruences using the Chinese Remainder
Theorem (with extension to non-coprime moduli). A list of C<[a,n]> pairs
are taken as input, each representing an equation C<x ≡ a mod n>. If no
solution exists, C<undef> is returned. If a solution is returned, the
modulus is equal to the lcm of all the given moduli (see L</lcm>. In
the standard case where all values of C<n> are coprime, this is just the
product. The C<n> values must be positive integers, while the C<a> values
are integers.
=head2 vecsum
Returns the sum of all arguments, each of which must be an integer.
=head2 vecprod
Returns the product of all arguments, each of which must be an integer.
=head2 kronecker
Returns the Kronecker symbol C<(a|n)> for two integers. The possible
return values with their meanings for odd positive C<n> are:
0 a = 0 mod n
1 a is a quadratic residue modulo n (a = x^2 mod n for some x)
-1 a is a quadratic non-residue modulo n
The Kronecker symbol is an extension of the Jacobi symbol to all integer
values of C<n> from the latter's domain of positive odd values of C<n>.
The Jacobi symbol is itself an extension of the Legendre symbol, which is
only defined for odd prime values of C<n>. This corresponds to Pari's
C<kronecker(a,n)> function and Mathematica's C<KroneckerSymbol[n,m]>
function.
=head2 binomial
Given integer arguments C<n> and C<k>, returns the binomial coefficient
C<n*(n-1)*...*(n-k+1)/k!>, also known as the choose function. Negative
arguments use the L<Kronenburg extensions|http://arxiv.org/abs/1105.3689/>.
This corresponds to Mathematica's C<Binomial[n,k]> function, Pari's
C<binomial(n,k)> function, and GMP's C<mpz_bin_ui> function.
For negative arguments, this matches Mathematica. Pari does not implement
the C<n E<lt> 0, k E<lt>= n> extension and instead returns C<0> for this
case. GMP's API does not allow negative C<k> but otherwise matches.
L<Math::BigInt> does not implement any extensions and the results for
C<n E<lt> 0, k > 0> are undefined.
=head2 bernfrac
Returns the Bernoulli number C<B_n> for an integer argument C<n>, as a
rational number. Two values are returned, the numerator and denominator.
B_1 = 1/2.
This corresponds to Pari's C<bernfrac(n)> and Mathematica's C<BernoulliB>
functions.
=head2 bernreal
Returns the Bernoulli number C<B_n> for an integer argument C<n>, as a
string floating point. An optional second argument indicates the number
of significant digits to be used, with the result rounded. The default
is 40 digits.
=head2 harmfrac
Returns the Harmonic number C<H_n> for an integer argument C<n>, as a
rational number. Two values are returned, the numerator and denominator.
numbers are the sum of reciprocals of the first C<n> natural numbers:
C<1 + 1/2 + 1/3 + ... + 1/n>.
=head2 harmreal
Returns the Harmonic number C<H_n> for an integer argument C<n>, as a
string floating point. An optional second argument indicates the number
of digits to be preserved past the decimal place, with a default of 40.
=head2 stirling
say "s(14,2) = ", stirling(14, 2);
say "S(14,2) = ", stirling(14, 2, 2);
Returns the Stirling numbers of either the first kind (default), the
second kind, or the third kind (the unsigned Lah numbers), with the kind
selected as an optional third argument. It takes two non-negative integer
arguments C<n> and C<k> plus the optional C<type>. This corresponds to Pari's
C<stirling(n,k,{type})> function and Mathematica's
C<StirlingS1> / C<StirlingS2> functions.
Stirling numbers of the first kind are C<-1^(n-k)> times the number of
permutations of C<n> symbols with exactly C<k> cycles. Stirling numbers
of the second kind are the number of ways to partition a set of C<n>
elements into C<k> non-empty subsets. The Lah numbers are the number of
ways to split a set of C<n> elements into C<k> non-empty lists.
=head2 zeta
Given a positive integer or float C<n>, returns the real Riemann Zeta
value as a string floating point. An optional second argument indicates
the number of digits past the decimal point (default 40).
The implementation is algorithm 2 of Borwein (1991). Performance with
integer inputs is good, but floating point arguments with high precision
will be much slower than methods using MPFR. L<Math::Prime::Util> will
try to use L<Math::MPFR> if possible.
=head2 riemannr
Given a positive integer or float C<n>, returns the real Riemann R function
as a string floating point. An optional second argument indicates the number
of significant digits (default 40) with the result rounded.
The implementation is the standard Gram series.
=head2 znorder
$order = znorder(17, "100000000000000000000000065");
Given two positive integers C<a> and C<n>, returns the multiplicative order
of C<a> modulo C<n>. This is the smallest positive integer C<k> such that
C<a^k ≡ 1 mod n>. Returns 1 if C<a = 1>. Returns undef if C<a = 0> or if
C<a> and C<n> are not coprime, since no value will result in 1 mod n.
This corresponds to Pari's C<znorder(Mod(a,n))> function and Mathematica's
C<MultiplicativeOrder[a,n]> function.
=head2 znprimroot
Given a positive integer C<n>, returns the smallest primitive root
of C<(Z/nZ)^*>, or C<undef> if no root exists. A root exists when
C<euler_phi($n) == carmichael_lambda($n)>, which will be true for
all prime C<n> and some composites.
L<OEIS A033948|http://oeis.org/A033948> is a sequence of integers where
the primitive root exists, while L<OEIS A046145|http://oeis.org/A046145>
is a list of the smallest primitive roots, which is what this function
produces.
=head2 is_primitive_root
Given two non-negative numbers C<a> and C<n>, returns C<1> if C<a> is a
primitive root modulo C<n>, and C<0> if not. If C<a> is a primitive root,
then C<euler_phi(n)> is the smallest C<e> for which C<a^e = 1 mod n>.
=head2 sigma
say "Sum of divisors of $n:", sigma( $n );
say "sigma_2($n) = ", sigma($n, 2);
say "Number of divisors: sigma_0($n) = ", sigma($n, 0);
This function takes a positive integer as input and returns the sum of
its divisors, including 1 and itself. An optional second argument C<k>
may be given, which will result in the sum of the C<k-th> powers of the
divisors to be returned.
This is known as the sigma function (see Hardy and Wright section 16.7,
or OEIS A000203). The API is identical to Pari/GP's C<sigma> function.
This function is useful for calculating things like aliquot sums, abundant
numbers, perfect numbers, etc.
=head2 ramanujan_tau
Takes a positive integer as input and returns the value of Ramanujan's tau
function. The result is a signed integer.
This corresponds to Mathematica's C<RamanujanTau> function.
=head2 valuation
say "$n is divisible by 2 ", valuation($n,2), " times.";
Given integers C<n> and C<k>, returns the numbers of times C<n> is divisible
by C<k>. This is a very limited version of the algebraic valuation meaning,
just applied to integers.
This corresponds to Pari's C<valuation> function.
C<0> is returned if C<n> or C<k> is one of the values C<-1>, C<0>, or C<1>.
=head2 moebius
say "$n is square free" if moebius($n) != 0;
$sum += moebius($_) for (1..200); say "Mertens(200) = $sum";
say "Mertens(2000) = ", vecsum(moebius(0,2000));
Returns μ(n), the Möbius function (also known as the Moebius, Mobius, or
MoebiusMu function) for an integer input. This function is 1 if
C<n = 1>, 0 if C<n> is not square free (i.e. C<n> has a repeated factor),
and C<-1^t> if C<n> is a product of C<t> distinct primes. This is an
important function in prime number theory. Like SAGE, we define
C<moebius(0) = 0> for convenience.
If called with two arguments, they define a range C<low> to C<high>, and the
function returns an array with the value of the Möbius function for every n
from low to high inclusive.
=head2 invmod
say "The inverse of 42 mod 2017 = ", invmod(42,2017);
Given two integers C<a> and C<n>, return the inverse of C<a> modulo C<n>.
If not defined, undef is returned. If defined, then the return value
multiplied by C<a> equals C<1> modulo C<n>.
=head2 sqrtmod
Given two integers C<a> and C<p>, return the square root of C<a> mod C<p>.
If no square root exists, undef is returned. If defined, the return value
C<s> will always satisfy C<mulmod(s,s,p) = a>.
If C<p> is not a prime, it is possible no result will be returned even
though a modular root exists.
Only one root is returned, even though there are at least two. In the
case of C<p> a prime and a return value C<s>, then both C<+s mod n> and
C<-s mod n> are roots. The least C<s> will be returned. In the case of
composites, many roots may exist, but only one will be returned.
=head2 addmod
Given three integers C<a>, C<b>, and C<n> where C<a> and C<n> are unsigned,
return C<(a+b) mod n>. This is particularly useful when dealing with
numbers that are larger than a half-word but still native size. No
bigint package is needed and this can be 10-200x faster than using one.
=head2 mulmod
Given three integers C<a>, C<b>, and C<n> where C<a> and C<n> are unsigned,
return C<(a*b) mod n>. This is particularly useful when C<n> fits in a
native integer. No bigint package is needed and this can be 10-200x
faster than using one.
=head2 powmod
Given three integers C<a>, C<b>, and C<n> where C<a> and C<n> are unsigned,
return C<(a ** b) mod n>. Typically binary exponentiation is used, so
the process is very efficient. With native size inputs, no bigint
library is needed.
=head2 divmod
Given three integers C<a>, C<b>, and C<n> where C<a> and C<n> are unsigned,
return C<(a/b) mod n>. This is done as C<(a * (1/b mod n)) mod n>. If
no inverse of C<b> mod C<n> exists then undef if returned.
=head2 consecutive_integer_lcm
$lcm = consecutive_integer_lcm($n);
Given an unsigned integer argument, returns the least common multiple of all
integers from 1 to C<n>. This can be done by manipulation of the primes up
to C<n>, resulting in much faster and memory-friendly results than using
factorials.
=head2 partitions
Calculates the partition function p(n) for a non-negative integer input.
This is the number of ways of writing the integer n as a sum of positive
integers, without restrictions. This corresponds to Pari's C<numbpart>
function and Mathematica's C<PartitionsP> function. The values produced
in order are L<OEIS series A000041|http://oeis.org/A000041>.
This uses a combinatorial calculation, which means it will not be very
fast compared to Pari, Mathematica, or FLINT which use the Rademacher
formula using multi-precision floating point. In 10 seconds, the pure
Perl version can produce C<partitions(10_000)> while with
L<Math::Prime::Util::GMP> it can do C<partitions(220_000)>. In contrast,
in about 10 seconds Pari can solve C<numbpart(22_000_000)>.
If you want the enumerated partitions, see L<Integer::Partition>. It is
very fast and uses an extremely memory efficient iterator. It is not,
however, practical for producing the partition I<number> for values
over 100 or so.
=head2 Pi
Takes a positive integer argument C<n> and returns the constant Pi with that
many digits (including the leading 3). Rounding is performed.
The implementation uses AGM and is only slightly slower than MPFR (which has
tighter bounds on the intermediate bits and exit conditions).
=head2 exp_mangoldt
say "exp(lambda($_)) = ", exp_mangoldt($_) for 1 .. 100;
Returns EXP(Λ(n)), the exponential of the Mangoldt function (also known
as von Mangoldt's function) for an integer value.
The Mangoldt function is equal to log p if n is prime or a power of a prime,
and 0 otherwise. We return the exponential so all results are integers.
Hence the return value for C<exp_mangoldt> is:
p if n = p^m for some prime p and integer m >= 1
1 otherwise.
=head2 totient
say "The Euler totient of $n is ", totient($n);
Returns φ(n), the Euler totient function (also called Euler's phi or phi
function) for an integer value. This is an arithmetic function which counts
the number of positive integers less than or equal to C<n> that are relatively
prime to C<n>. Given the definition used, C<totient> will return 0 for all
C<n E<lt> 1>. This follows the logic used by SAGE. Mathematica and Pari
return C<totient(-n)> for C<n E<lt> 0>. Mathematica returns 0 for C<n = 0>,
Pari pre-2.6.2 raises and exception, and Pari 2.6.2 and newer returns 2.
=head2 jordan_totient
say "Jordan's totient J_$k($n) is ", jordan_totient($k, $n);
Returns Jordan's totient function for a given integer value. Jordan's totient
is a generalization of Euler's totient, where
C<jordan_totient(1,$n) == euler_totient($n)>
This counts the number of k-tuples less than or equal to n that form a coprime
tuple with n. As with C<totient>, 0 is returned for all C<n E<lt> 1>.
This function can be used to generate some other useful functions, such as
the Dedekind psi function, where C<psi(n) = J(2,n) / J(1,n)>.
=head2 carmichael_lambda
Returns the Carmichael function (also called the reduced totient function,
or Carmichael λ(n)) of a positive integer argument. It is the smallest
positive integer C<m> such that C<a^m = 1 mod n> for every integer C<a>
coprime to C<n>. This is L<OEIS series A002322|http://oeis.org/A002322>.
=head2 liouville
Returns λ(n), the Liouville function for a non-negative integer input.
This is -1 raised to Ω(n) (the total number of prime factors).
=head2 is_power
say "$n is a perfect square" if is_power($n, 2);
say "$n is a perfect cube" if is_power($n, 3);
say "$n is a ", is_power($n), "-th power";
Given a single positive integer input C<n>, returns k if C<n = p^k> for
some integer C<p E<gt> 1, k E<gt> 1>, and 0 otherwise. The k returned is
the largest possible. This can be used in a boolean statement to
determine if C<n> is a perfect power.
If given two arguments C<n> and C<k>, returns 1 if C<n> is a C<k-th> power,
and 0 otherwise. For example, if C<k=2> then this detects perfect squares.
This corresponds to Pari/GP's C<ispower> function, with the limitations of
only integer arguments and no third argument may be given to return the root.
=head2 is_prime_power
Given an integer input C<n>, returns C<k> if C<n = p^k> for some prime p,
and zero otherwise.
This corresponds to Pari/GP's C<isprimepower> function.
=head2 sqrtint
Returns the truncated integer part of the square root of C<n>.
This corresponds to Pari/GP's C<sqrtint> function.
=head2 rootint
Given C<n> and C<k>, returns the truncated integer part of the C<k-th> root
of C<n>.
This corresponds to Pari/GP's C<sqrtnint> function.
=head2 factor
@factors = factor(640552686568398413516426919223357728279912327120302109778516984973296910867431808451611740398561987580967216226094312377767778241368426651540749005659);
# Returns an array of 11 factors
Returns a list of prime factors of a positive number, in numerical order. The
special cases of C<n = 0> and C<n = 1> will return C<n>.
Like most advanced factoring programs, a mix of methods is used. This
includes trial division for small factors, perfect power detection,
Pollard's Rho, Pollard's P-1 with various smoothness and stage settings,
Hart's OLF (a Fermat variant), ECM (elliptic curve method), and
QS (quadratic sieve).
Certainly improvements could be designed for this algorithm
(suggestions are welcome).
In practice, this factors 26-digit semiprimes in under C<100ms>, 36-digit
semiprimes in under one second. Arbitrary integers are factored faster.
It is many orders of magnitude faster than any other factoring module on
CPAN circa 2013. It is comparable in speed to Math::Pari's C<factorint>
for most inputs.
If you want better factoring in general, I recommend looking at the
standalone programs
L<yafu|http://sourceforge.net/projects/yafu/>,
L<msieve|http://sourceforge.net/projects/msieve/>,
L<gmp-ecm|http://ecm.gforge.inria.fr/>, and
L<GGNFS|http://sourceforge.net/projects/ggnfs/>.
=head2 divisors
my @divisors = divisors(30); # returns (1, 2, 3, 5, 6, 10, 15, 30)
Produces all the divisors of a positive number input, including 1 and
the input number. The divisors are a power set of multiplications of
the prime factors, returned as a sorted list with no duplications.
The result is identical to that of Pari's C<divisors> and
Mathematica's C<Divisors[n]> functions.
In scalar context this returns the sigma0 function (OEIS A000005), and
is the same result as evaluating the returned array in scalar context
(but much more efficient).
The result then corresponds to Pari's C<numdiv> and Mathematica's
C<DivisorSigma[0,n]> functions.
=head2 trial_factor
my @factors = trial_factor($n);
my @factors = trial_factor($n, 1000);
Given a positive number input, tries to discover a factor using trial division.
The resulting array will contain either two factors (it succeeded) or the
original number (no factor was found). In either case, multiplying @factors
yields the original input. An optional divisor limit may be given as the
second parameter. Factoring will stop when the input is a prime, one factor
is found, or the input has been tested for divisibility with all primes less
than or equal to the limit. If no limit is given, then C<2**31-1> will be used.
This is a good and fast initial test, and will be very fast for small numbers
(e.g. under 1 million). For larger numbers, faster methods for complete
factoring have been known since the 17th century.
For inputs larger than about 1000 digits, a dynamic product/remainder tree
is used, which is faster than GMP's native methods. This helps when pruning
composites or looking for very small factors.
=head2 prho_factor
my @factors = prho_factor($n);
my @factors = prho_factor($n, 100_000_000);
Given a positive number input, tries to discover a factor using Pollard's Rho
method. The resulting array will contain either two factors (it succeeded)
or the original number (no factor was found). In either case, multiplying
@factors yields the original input. An optional number of rounds may be
given as the second parameter. Factoring will stop when the input is a prime,
one factor has been found, or the number of rounds has been exceeded.
This is the Pollard Rho method with C<f = x^2 + 3> and default rounds 64M. It
is very good at finding small factors. Typically L</pbrent_factor> will be
preferred as it behaves similarly but runs quite a bit faster. They use
different parameters however, so are not completely identical.
=head2 pbrent_factor
my @factors = pbrent_factor($n);
my @factors = pbrent_factor($n, 100_000_000);
Given a positive number input, tries to discover a factor using Pollard's Rho
method with Brent's algorithm. The resulting array will contain either two
factors (it succeeded) or the original number (no factor was found). In
either case, multiplying @factors yields the original input. An optional
number of rounds may be given as the second parameter. Factoring will stop
when the input is a prime, one factor has been found, or the number of
rounds has been exceeded.
This is the Pollard Rho method using Brent's modified cycle detection,
delayed C<gcd> computations, and backtracking. It is essentially
Algorithm P''2 from Brent (1980). Parameters used are C<f = x^2 + 3>
and default rounds 64M. It is very good at finding small factors.
=head2 pminus1_factor
my @factors = pminus1_factor($n);
# Set B1 smoothness to 10M, second stage automatically set.
my @factors = pminus1_factor($n, 10_000_000);
# Run p-1 with B1 = 10M, B2 = 100M.
my @factors = pminus1_factor($n, 10_000_000, 100_000_000);
Given a positive number input, tries to discover a factor using Pollard's
C<p-1> method. The resulting array will contain either two factors (it
succeeded) or the original number (no factor was found). In either case,
multiplying @factors yields the original input. An optional first stage
smoothness factor (B1) may be given as the second parameter. This will be
the smoothness limit B1 for the first stage, and will use C<10*B1> for
the second stage limit B2. If a third parameter is given, it will be used
as the second stage limit B2.
Factoring will stop when the input is a prime, one factor has been found, or
the algorithm fails to find a factor with the given smoothness.
This is Pollard's C<p-1> method using a default smoothness of 5M and a
second stage of C<B2 = 10 * B1>. It can quickly find a factor C<p> of the input
C<n> if the number C<p-1> factors into small primes. For example
C<n = 22095311209999409685885162322219> has the factor C<p = 3916587618943361>,
where C<p-1 = 2^7 * 5 * 47 * 59 * 3137 * 703499>, so this method will find
a factor in the first stage if C<B1 E<gt>= 703499> or in the second stage if
C<B1 E<gt>= 3137> and C<B2 E<gt>= 703499>.
The implementation is written from scratch using the basic algorithm including
a second stage as described in Montgomery 1987. It is faster than most simple
implementations I have seen (many of which are written assuming native
precision inputs), but slower than Ben Buhrow's code used in earlier
versions of L<yafu|http://sourceforge.net/projects/yafu/>, and nowhere close
to the speed of the version included with modern GMP-ECM with large B values
(it is actually quite a bit faster than GMP-ECM with small smoothness values).
=head2 pplus1_factor
my @factors = pplus1_factor($n);
Given a positive number input, tries to discover a factor using Williams'
C<p+1> method. The resulting array will contain either two factors (it
succeeded) or the original number (no factor was found). In either case,
multiplying @factors yields the original input. An optional first stage
smoothness factor (B1) may be given as the second parameter. This will be
the smoothness limit B1 for the first stage.
Factoring will stop when the input is a prime, one factor has been found, or
the algorithm fails to find a factor with the given smoothness.
=head2 holf_factor
my @factors = holf_factor($n);
my @factors = holf_factor($n, 100_000_000);
Given a positive number input, tries to discover a factor using Hart's OLF
method. The resulting array will contain either two factors (it succeeded)
or the original number (no factor was found). In either case, multiplying
@factors yields the original input. An optional number of rounds may be
given as the second parameter. Factoring will stop when the input is a
prime, one factor has been found, or the number of rounds has been exceeded.
This is Hart's One Line Factorization method, which is a variant of Fermat's
algorithm. A premultiplier of 480 is used. It is very good at factoring
numbers that are close to perfect squares, or small numbers. Very naive
methods of picking RSA parameters sometimes yield numbers in this form, so
it can be useful to run a few rounds to check. For example, the number:
18548676741817250104151622545580576823736636896432849057 \
10984160646722888555430591384041316374473729421512365598 \
29709849969346650897776687202384767704706338162219624578 \
777915220190863619885201763980069247978050169295918863
was proposed by someone as an RSA key. It is indeed composed of two distinct
prime numbers of similar bit length. Most factoring methods will take a
B<very> long time to break this. However one factor is almost exactly 5x
larger than the other, allowing HOLF to factor this 222-digit semiprime in
only a few milliseconds.
=head2 squfof_factor
my @factors = squfof_factor($n);
my @factors = squfof_factor($n, 100_000_000);
Given a positive number input, tries to discover a factor using Shanks'
square forms factorization method (usually known as SQUFOF). The resulting
array will contain either two factors (it succeeded) or the original number
(no factor was found). In either case, multiplying @factors yields the
original input. An optional number of rounds may be given as the second
parameter. Factoring will stop when the input is a prime, one factor has
been found, or the number of rounds has been exceeded.
This is Daniel Shanks' SQUFOF (square forms factorization) algorithm. The
particular implementation is a non-racing multiple-multiplier version, based
on code ideas of Ben Buhrow and Jason Papadopoulos as well as many others.
SQUFOF is often the preferred method for small numbers, and L<Math::Prime::Util>
as well as many other packages use it was the default method for native size
(e.g. 32-bit or 64-bit) numbers after trial division. The GMP version used
in this module will work for larger values, but my testing indicates it is
generally slower than the C<prho> and C<pbrent> implementations.
=head2 ecm_factor
my @factors = ecm_factor($n);
my @factors = ecm_factor($n, 12500); # B1 = 12500
my @factors = ecm_factor($n, 12500, 10); # B1 = 12500, curves = 10
Given a positive number input, tries to discover a factor using ECM. The
resulting array will contain either two factors (it succeeded) or the original
number (no factor was found). In either case, multiplying @factors yields the
original input. An optional maximum smoothness may be given as the second
parameter, which relates to the size of factor to search for. An optional
third parameter indicates the number of random curves to use at each
smoothness value being searched.
This is an implementation of Hendrik Lenstra's elliptic curve factoring
method, usually referred to as ECM. The implementation is reasonable,
using projective coordinates, Montgomery's PRAC heuristic for EC
multiplication, and two stages.
It is much slower than the latest GMP-ECM, but still quite useful for
factoring reasonably sized inputs.
=head2 qs_factor
my @factors = qs_factor($n);
Given a positive number input, tries to discover factors using QS (the
quadratic sieve). The resulting array will contain one or more numbers such
that multiplying @factors yields the original input. Typically multiple
factors will be produced, unlike the other C<..._factor> routines.
The current implementation is a modified version of SIMPQS, a predecessor to
the QS in FLINT, and was written by William Hart in 2006. It will not operate
on input less than 30 digits. The memory use for large inputs is more than
desired, so other methods such as L</pbrent_factor>, L</pminus1_factor>, and
L</ecm_factor> are recommended to begin with to filter out small factors.
However, it is substantially faster than the other methods on large inputs
having large factors, and is the method of choice for 35+ digit semiprimes.
=head2 todigits
Given an integer C<n>, return an array of digits of C<|n|>. An optional
second integer argument specifies a base (default 10). For example,
given a base of 2, this returns an array of binary digits of C<n>.
An optional third argument specifies a length for the returned array.
The result will be either have upper digits truncated or have leading
zeros added.
C<todigits(0)> returns an empty array.
The base must be at least 2, and is limited to an int.
Length must be at least zero and is limited to an int.
=head1 SEE ALSO
=over 4
=item L<Math::Prime::Util>
Has many more functions, lots of fast code for dealing with native-precision
arguments (including much faster primes using sieves), and will use this
module when needed for big numbers. Using L<Math::Prime::Util> rather than
this module directly is recommended.
=item L<Math::Primality> (version 0.08)
A Perl module with support for the strong Miller-Rabin test, strong
Lucas-Selfridge test, the BPSW probable prime test, next_prime / prev_prime,
the AKS primality test, and prime_count. It uses L<Math::GMPz> to do all
the calculations, so is faster than pure Perl bignums, but a little slower
than XS+GMP. The prime_count function is only usable for very small inputs,
but the other functions are quite good for big numbers. Make sure to use
version 0.05 or newer.
=item L<Math::Pari>
Supports quite a bit of the same functionality (and much more). See
L<Math::Prime::Util/"SEE ALSO"> for more detailed information on how the
modules compare.
=item L<yafu|http://sourceforge.net/projects/yafu/>,
L<msieve|http://sourceforge.net/projects/msieve/>,
L<gmp-ecm|http://ecm.gforge.inria.fr/>,
L<GGNFS|http://sourceforge.net/projects/ggnfs/>
Good general purpose factoring utilities. These will be faster than this
module, and B<much> better as the factor increases in size.
=item L<Primo|http://www.ellipsa.eu/public/primo/primo.html>
is the state of the art in freely available (though not open source!)
primality proving programs. If you have 1000+ digit numbers to prove,
you want to use this.
=item L<mpz_aprcl|http://sourceforge.net/projects/mpzaprcl/>
Open source APR-CL primality proof implementation.
Fast primality proving, though without certificates.
=item L<GMP-ECPP|http://sourceforge.net/projects/gmp-ecpp/>.
An open source ECPP primality proving program. Slower than this
module's ECPP for all inputs when the large polynomial set from
github is used. Extremely slow once past 300 or so digits.
There are now better alternatives.
=back
=head1 REFERENCES
=over 4
=item Robert Baillie and Samuel S. Wagstaff, Jr., "Lucas Pseudoprimes", Mathematics of Computation, v35 n152, October 1980, pp 1391-1417. L<http://mpqs.free.fr/LucasPseudoprimes.pdf>
=item Daniel J. Bernstein, "Proving Primality After Agrawal-Kayal-Saxena",
preprint, Jan 2003. L<http://cr.yp.to/papers/aks.pdf>
=item Jon Grantham, "Frobenius Pseudoprimes", Mathematics of Computation, v70 n234, March 2000, pp 873-891. L<http://www.ams.org/journals/mcom/2001-70-234/S0025-5718-00-01197-2/>
=item John Brillhart, D. H. Lehmer, and J. L. Selfridge, "New Primality Criteria and Factorizations of 2^m +/- 1", Mathematics of Computation, v29, n130, Apr 1975, pp 620-647. L<http://www.ams.org/journals/mcom/1975-29-130/S0025-5718-1975-0384673-1/S0025-5718-1975-0384673-1.pdf>
=item Richard P. Brent, "An improved Monte Carlo factorization algorithm", BIT 20, 1980, pp. 176-184. L<http://www.cs.ox.ac.uk/people/richard.brent/pd/rpb051i.pdf>
=item Peter L. Montgomery, "Speeding the Pollard and Elliptic Curve Methods of Factorization", Mathematics of Computation, v48, n177, Jan 1987, pp 243-264. L<http://www.ams.org/journals/mcom/1987-48-177/S0025-5718-1987-0866113-7/>
=item Richard P. Brent, "Parallel Algorithms for Integer Factorisation", in Number Theory and Cryptography, Cambridge University Press, 1990, pp 26-37. L<http://www.cs.ox.ac.uk/people/richard.brent/pd/rpb115.pdf>
=item Richard P. Brent, "Some Parallel Algorithms for Integer Factorisation", in Proc. Third Australian Supercomputer Conference, 1999. (Note: there are multiple versions of this paper) L<http://www.cs.ox.ac.uk/people/richard.brent/pd/rpb193.pdf>
=item William B. Hart, "A One Line Factoring Algorithm", preprint. L<http://wstein.org/home/wstein/www/home/wbhart/onelinefactor.pdf>
=item Daniel Shanks, "SQUFOF notes", unpublished notes, transcribed by Stephen McMath. L<http://www.usna.edu/Users/math/wdj/mcmath/shanks_squfof.pdf>
=item Jason E. Gower and Samuel S. Wagstaff, Jr, "Square Form Factorization", Mathematics of Computation, v77, 2008, pages 551-588. L<http://homes.cerias.purdue.edu/~ssw/squfof.pdf>
=item A.O.L. Atkin and F. Morain, "Elliptic Curves and primality proving", Mathematics of Computation, v61, 1993, pages 29-68. L<http://www.ams.org/journals/mcom/1993-61-203/S0025-5718-1993-1199989-X/>
=item R.G.E. Pinch, "Some Primality Testing Algorithms", June 1993. Describes the primality testing methods used by many CAS systems and how most were compromised. Gives recommendations for primality testing APIs. L<http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.33.4409>
=back
=head1 AUTHORS
Dana Jacobsen E<lt>dana@acm.orgE<gt>
William Hart wrote the SIMPQS code which is the basis for the QS code.
=head1 ACKNOWLEDGEMENTS
Obviously none of this would be possible without the mathematicians who
created and published their work. Eratosthenes, Gauss, Euler, Riemann,
Fermat, Lucas, Baillie, Pollard, Brent, Montgomery, Shanks, Hart, Wagstaff,
Dixon, Pomerance, A.K. Lenstra, H. W. Lenstra Jr., Atkin, Knuth, etc.
The GNU GMP team, whose product allows me to concentrate on coding high-level
algorithms and not worry about any of the details of how modular exponentiation
and the like happen, and still get decent performance for my purposes.
Ben Buhrow and Jason Papadopoulos deserve special mention for their open
source factoring tools, which are both readable and fast. In particular I am
leveraging their SQUFOF work in the current implementation. They are a huge
resource to the community.
Jonathan Leto and Bob Kuo, who wrote and distributed the L<Math::Primality>
module on CPAN. Their implementation of BPSW provided the motivation I needed
to do it in this module and L<Math::Prime::Util>. I also used their
module quite a bit for testing against.
Paul Zimmermann's papers and GMP-ECM code were of great value for my projective
ECM implementation, as well as the papers by Brent and Montgomery.
=head1 COPYRIGHT
Copyright 2011-2016 by Dana Jacobsen E<lt>dana@acm.orgE<gt>
This program is free software; you can redistribute it and/or modify it under the same terms as Perl itself.
SIMPQS Copyright 2006, William Hart. SIMPQS is distributed under GPL v2+.
=cut
|