2016-06-23 23:51:08 +02:00
|
|
|
<?php
|
2016-07-14 15:15:50 +02:00
|
|
|
|
|
|
|
set_include_path(get_include_path().PATH_SEPARATOR.dirname(__FILE__).DIRECTORY_SEPARATOR.'libpy2php');
|
|
|
|
require_once 'libpy2php.php';
|
|
|
|
function everynth($array, $n)
|
|
|
|
{
|
|
|
|
$result = [];
|
2016-06-24 23:47:10 +02:00
|
|
|
$i = -1;
|
2016-07-14 15:15:50 +02:00
|
|
|
foreach ($array as $key => $value) {
|
|
|
|
if ($i++ == $n) {
|
2016-06-24 23:47:10 +02:00
|
|
|
$i = 0;
|
2016-07-14 15:15:50 +02:00
|
|
|
}
|
|
|
|
if ($i == 0) {
|
|
|
|
$result[$key] = $value;
|
2016-06-24 23:47:10 +02:00
|
|
|
}
|
|
|
|
}
|
2016-07-14 15:15:50 +02:00
|
|
|
|
2016-06-24 23:47:10 +02:00
|
|
|
return $result;
|
|
|
|
}
|
2016-07-14 15:15:50 +02:00
|
|
|
function array_merge_ignore_keys($array1, $array2)
|
|
|
|
{
|
|
|
|
if (count($array1) == count($array2)) {
|
|
|
|
$i = -1;
|
|
|
|
foreach ($array1 as $key => $val) {
|
|
|
|
$array1[$key] = $array2[$i++];
|
|
|
|
}
|
|
|
|
} else {
|
|
|
|
return;
|
|
|
|
}
|
|
|
|
|
|
|
|
return $array1;
|
2016-06-25 20:39:49 +02:00
|
|
|
}
|
2016-07-14 15:15:50 +02:00
|
|
|
function primesbelow($N)
|
|
|
|
{
|
2016-06-23 23:51:08 +02:00
|
|
|
$correction = (($N % 6) > 1);
|
2016-07-14 15:15:50 +02:00
|
|
|
$N = [0 => $N, 1 => ($N - 1), 2 => ($N + 4), 3 => ($N + 3), 4 => ($N + 2), 5 => ($N + 1)][($N % 6)];
|
2016-06-26 12:59:15 +02:00
|
|
|
$sieve = array_fill(0, floor($N / 3), true);
|
2016-06-23 23:51:08 +02:00
|
|
|
$sieve[0] = false;
|
2016-07-18 17:11:37 +02:00
|
|
|
foreach (pyjslib_range(floor((int)(pow($N, 0.5)) / 3) + 1) as $i) {
|
2016-06-23 23:51:08 +02:00
|
|
|
if ($sieve[$i]) {
|
|
|
|
$k = ((3 * $i) + 1) | 1;
|
2016-07-18 17:11:37 +02:00
|
|
|
foreach (pyjslib_range(floor(($k*$k) / 3), 2*$k) as $key) {
|
|
|
|
|
|
|
|
}
|
|
|
|
array_fill(0, floor((floor($N / 6) - floor(($k * $k) / 6) - 1) / $k) + 1, false)
|
2016-06-26 12:59:15 +02:00
|
|
|
$sieve = array_merge(
|
|
|
|
$sieve, array_merge(
|
|
|
|
everynth(
|
|
|
|
array_splice(
|
|
|
|
$sieve, floor(
|
|
|
|
(
|
2016-07-14 15:15:50 +02:00
|
|
|
$k * $k + 4 * $k - 2 * $k * ($i % 2)
|
2016-06-26 12:59:15 +02:00
|
|
|
) / 3
|
|
|
|
)
|
2016-07-14 15:15:50 +02:00
|
|
|
), 2 * $k
|
|
|
|
),
|
2016-06-26 12:59:15 +02:00
|
|
|
(
|
2016-07-18 17:11:37 +02:00
|
|
|
|
2016-06-26 12:59:15 +02:00
|
|
|
)
|
|
|
|
)
|
|
|
|
);
|
2016-06-23 23:51:08 +02:00
|
|
|
}
|
|
|
|
}
|
2016-07-18 17:11:37 +02:00
|
|
|
var_dump($sieve);
|
2016-07-14 15:15:50 +02:00
|
|
|
|
|
|
|
return [2, 3] + array_map(function ($i, $sieve) {
|
|
|
|
if ($sieve[$i]) {
|
|
|
|
return (3 * $i + 1) | 1;
|
|
|
|
}
|
|
|
|
}, pyjslib_range(1, (($N / 3) - $correction)), $sieve);
|
2016-06-23 23:51:08 +02:00
|
|
|
}
|
2016-06-26 12:59:15 +02:00
|
|
|
$smallprimeset = array_unique(primesbelow(100000));
|
2016-06-23 23:51:08 +02:00
|
|
|
$_smallprimeset = 100000;
|
2016-07-14 15:15:50 +02:00
|
|
|
function isprime($n, $precision = 7)
|
|
|
|
{
|
2016-06-23 23:51:08 +02:00
|
|
|
if (($n == 1) || (($n % 2) == 0)) {
|
|
|
|
return false;
|
2016-07-14 15:15:50 +02:00
|
|
|
} elseif (($n < 1)) {
|
2016-06-23 23:51:08 +02:00
|
|
|
throw new $ValueError('Out of bounds, first argument must be > 0');
|
2016-07-14 15:15:50 +02:00
|
|
|
} elseif (($n < $_smallprimeset)) {
|
2016-06-23 23:51:08 +02:00
|
|
|
return in_array($n, $smallprimeset);
|
|
|
|
}
|
|
|
|
$d = ($n - 1);
|
|
|
|
$s = 0;
|
|
|
|
while ((($d % 2) == 0)) {
|
2016-07-14 15:15:50 +02:00
|
|
|
$d = floor($d / 2);
|
|
|
|
$s += 1;
|
2016-06-23 23:51:08 +02:00
|
|
|
}
|
|
|
|
foreach (pyjslib_range($precision) as $repeat) {
|
|
|
|
$a = random::randrange(2, ($n - 2));
|
|
|
|
$x = pow($a, $d, $n);
|
|
|
|
if (($x == 1) || ($x == ($n - 1))) {
|
|
|
|
continue;
|
|
|
|
}
|
|
|
|
foreach (pyjslib_range(($s - 1)) as $r) {
|
|
|
|
$x = pow($x, 2, $n);
|
|
|
|
if (($x == 1)) {
|
|
|
|
return false;
|
|
|
|
}
|
|
|
|
if (($x == ($n - 1))) {
|
|
|
|
break;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2016-07-14 15:15:50 +02:00
|
|
|
|
2016-06-23 23:51:08 +02:00
|
|
|
return true;
|
|
|
|
}
|
2016-07-14 15:15:50 +02:00
|
|
|
function pollard_brent($n)
|
|
|
|
{
|
2016-06-23 23:51:08 +02:00
|
|
|
if ((($n % 2) == 0)) {
|
|
|
|
return 2;
|
|
|
|
}
|
|
|
|
if ((($n % 3) == 0)) {
|
|
|
|
return 3;
|
|
|
|
}
|
2016-07-14 15:15:50 +02:00
|
|
|
list($y, $c, $m) = [random::randint(1, ($n - 1)), random::randint(1, ($n - 1)), random::randint(1, ($n - 1))];
|
2016-06-23 23:51:08 +02:00
|
|
|
list($g, $r, $q) = [1, 1, 1];
|
|
|
|
while (($g == 1)) {
|
|
|
|
$x = $y;
|
|
|
|
foreach (pyjslib_range($r) as $i) {
|
|
|
|
$y = ((pow($y, 2, $n) + $c) % $n);
|
|
|
|
}
|
|
|
|
$k = 0;
|
|
|
|
while (($k < $r) && ($g == 1)) {
|
|
|
|
$ys = $y;
|
|
|
|
foreach (pyjslib_range(min($m, ($r - $k))) as $i) {
|
|
|
|
$y = ((pow($y, 2, $n) + $c) % $n);
|
|
|
|
$q = (($q * abs(($x - $y))) % $n);
|
|
|
|
}
|
|
|
|
$g = gcd($q, $n);
|
2016-07-14 15:15:50 +02:00
|
|
|
$k += $m;
|
2016-06-23 23:51:08 +02:00
|
|
|
}
|
2016-07-14 15:15:50 +02:00
|
|
|
$r *= 2;
|
2016-06-23 23:51:08 +02:00
|
|
|
}
|
|
|
|
if (($g == $n)) {
|
|
|
|
while (true) {
|
|
|
|
$ys = ((pow($ys, 2, $n) + $c) % $n);
|
|
|
|
$g = gcd(abs(($x - $ys)), $n);
|
|
|
|
if (($g > 1)) {
|
|
|
|
break;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2016-07-14 15:15:50 +02:00
|
|
|
|
2016-06-23 23:51:08 +02:00
|
|
|
return $g;
|
|
|
|
}
|
|
|
|
$smallprimes = primesbelow(10000);
|
2016-07-14 15:15:50 +02:00
|
|
|
function primefactors($n, $sort = false)
|
|
|
|
{
|
2016-06-26 15:28:39 +02:00
|
|
|
global $smallprimes;
|
2016-06-23 23:51:08 +02:00
|
|
|
$factors = [];
|
|
|
|
$limit = (pyjslib_int(pow($n, 0.5)) + 1);
|
|
|
|
foreach ($smallprimes as $checker) {
|
|
|
|
if (($checker > $limit)) {
|
|
|
|
break;
|
|
|
|
}
|
|
|
|
while ((($n % $checker) == 0)) {
|
|
|
|
$factors[] = $checker;
|
2016-06-25 00:01:38 +02:00
|
|
|
$n = floor($n / $checker);
|
2016-06-23 23:51:08 +02:00
|
|
|
$limit = (pyjslib_int(pow($n, 0.5)) + 1);
|
|
|
|
if (($checker > $limit)) {
|
|
|
|
break;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
if (($n < 2)) {
|
|
|
|
return $factors;
|
|
|
|
}
|
|
|
|
while (($n > 1)) {
|
|
|
|
if (isprime($n)) {
|
|
|
|
$factors[] = $n;
|
|
|
|
break;
|
|
|
|
}
|
|
|
|
$factor = pollard_brent($n);
|
|
|
|
$factors->extend(primefactors($factor));
|
2016-06-25 00:01:38 +02:00
|
|
|
$n = floor($n / $factor);
|
2016-06-23 23:51:08 +02:00
|
|
|
}
|
|
|
|
if ($sort) {
|
|
|
|
$factors->sort();
|
|
|
|
}
|
2016-07-14 15:15:50 +02:00
|
|
|
|
2016-06-23 23:51:08 +02:00
|
|
|
return $factors;
|
|
|
|
}
|
2016-07-14 15:15:50 +02:00
|
|
|
function factorization($n)
|
|
|
|
{
|
2016-06-23 23:51:08 +02:00
|
|
|
$factors = [];
|
|
|
|
foreach (primefactors($n) as $p1) {
|
|
|
|
try {
|
2016-07-14 15:15:50 +02:00
|
|
|
$factors[$p1] += 1;
|
|
|
|
} catch (KeyError $e) {
|
2016-06-23 23:51:08 +02:00
|
|
|
$factors[$p1] = 1;
|
|
|
|
}
|
|
|
|
}
|
2016-07-14 15:15:50 +02:00
|
|
|
|
2016-06-23 23:51:08 +02:00
|
|
|
return $factors;
|
|
|
|
}
|
|
|
|
$totients = [];
|
2016-07-14 15:15:50 +02:00
|
|
|
function totient($n)
|
|
|
|
{
|
2016-06-23 23:51:08 +02:00
|
|
|
if (($n == 0)) {
|
|
|
|
return 1;
|
|
|
|
}
|
|
|
|
try {
|
|
|
|
return $totients[$n];
|
2016-07-14 15:15:50 +02:00
|
|
|
} catch (KeyError $e) {
|
2016-06-23 23:51:08 +02:00
|
|
|
}
|
|
|
|
$tot = 1;
|
|
|
|
foreach (factorization($n)->items() as list($p, $exp)) {
|
2016-07-14 15:15:50 +02:00
|
|
|
$tot *= (($p - 1) * pow($p, ($exp - 1)));
|
2016-06-23 23:51:08 +02:00
|
|
|
}
|
|
|
|
$totients[$n] = $tot;
|
2016-07-14 15:15:50 +02:00
|
|
|
|
2016-06-23 23:51:08 +02:00
|
|
|
return $tot;
|
|
|
|
}
|
2016-07-14 15:15:50 +02:00
|
|
|
function gcd($a, $b)
|
|
|
|
{
|
2016-06-23 23:51:08 +02:00
|
|
|
if (($a == $b)) {
|
|
|
|
return $a;
|
|
|
|
}
|
|
|
|
while (($b > 0)) {
|
2016-07-14 15:15:50 +02:00
|
|
|
list($a, $b) = [$b, ($a % $b)];
|
2016-06-23 23:51:08 +02:00
|
|
|
}
|
2016-07-14 15:15:50 +02:00
|
|
|
|
2016-06-23 23:51:08 +02:00
|
|
|
return $a;
|
|
|
|
}
|
2016-07-14 15:15:50 +02:00
|
|
|
function lcm($a, $b)
|
|
|
|
{
|
2016-06-26 12:59:15 +02:00
|
|
|
return floor(abs(($a * $b)) / gcd($a, $b));
|
2016-06-23 23:51:08 +02:00
|
|
|
}
|