experimental-projects/oeis-research/Thomas Ordowski/12 June/prog.pl

#!/usr/bin/perl


#~ Consider the following sequences:

#~ Let a(n) be the smallest prime p such that n^p == n (mod q), for n = 0,1,2,...
#~ (*) where q is the smallest prime > p such that p-1 | q-1.
#~ (**) where q is the next prime after p.

#~ (*) 2, 2, 11, 2, 2, 3, 2, 2, 5, 2, ... Is this sequence periodic?
#~ (**) 2, 2, 113, 2, 2, 3, 2, 2, 5, 2, ... Is this sequence bounded?

#~ Have a nice day!

#~ Corrected set: A = {9, 15, 21, 25, 33, 39, 49, 51, 57, 65, 69, 85, 87, 91, 93, 111, 121, 123, 129, 133, 141, 145, 159, 169, 177, 183, 185, 201, 205, 213, 217, 219, 237, 249, 259, 265, 267, 289, 291, 301, 303, 305, 309, 321, 327, 339, 341, 361, 365, 381, 393, 411, 417, 427, 445, 447, 451, 453, 469, 471, 481, 485, 489, 501, 505, 511, 519, 529, 537, 543, 545, 553, 561, 565, 573, 579, 591, 597, 633, 645, 669, 671, 679, 681, 685, 687, 699, 703, 717, 721, 723, 745, 753, 763, 771, 781, 785, 789, 793, 807, 813, 831, 841, 843, 849, 865, 879, 889, 905, 921, 933, 939, 949, 951, 961, 965, 973, 985, 993, 1011, 1041, 1047, 1057, 1059, 1065, 1077, 1099, 1101, 1111, 1119, 1137, 1141, 1145, 1149, 1165, 1167, 1191, 1203, 1205, 1227, 1247, 1257, 1261, 1263, 1267, 1285, 1293, 1299, 1317, 1329, 1345, 1347, 1351, 1369, 1371, 1383, 1385, 1387, 1389, 1393, 1401, 1405, 1417, 1437, 1441, 1461, 1465, 1473, 1477, 1497, 1509, 1527, 1541, 1561, 1563, 1565, 1569, 1585, 1603, 1623, 1641, 1649, 1661, 1671, 1681, 1685, 1687, 1689, 1707, 1713, 1729}

#~ LCM(A) = 19974189529267233805048246587648609684002517045072696325440852333740211305925469232124401559514495489578549049794060699217424449904662437604539103637539339468499990506897971814716522960950895471907013503223284820116877757184319629345000591535688675

use 5.014;
use ntheory qw(:all);

sub find_q {
    my ($p) = @_;

    foreach my $k(2..1e10) {
        if (is_prime(($p-1)*$k + 1)) {
            return (($p-1)*$k + 1);
        }
    }
}

sub a {     # p-1 | q-1
    my ($n) = @_;

    my $iter = prime_iterator();

    for (my $p = $iter->(); ; $p = $iter->()) {
        my $q = find_q($p);

        if (powmod($n, $p, $q) == ($n % $q)) {
            return $p; #($p, $q);
        }
    }
}

sub b {     # nextprime (**)
    my ($n) = @_;

    my $iter = prime_iterator();

    for (my $p = $iter->(); ; $p = $iter->()) {
        my $q = next_prime($p);

        if (powmod($n, $p, $q) == ($n % $q)) {
            return $p;
        }
    }
}

sub c {
    my ($n) = @_;

    for(my $k = 3; $k <= 1729; $k += 2) {

        next if is_prime($k);

        my $t = powmod($n, ($k+1)>>1, $k);

        if($t == ($n%$k) or $t == (-$n)%$k) {
            return $k;
        }
    }
}


use Math::GMPz;
my $period = Math::GMPz->new("19974189529267233805048246587648609684002517045072696325440852333740211305925469232124401559514495489578549049794060699217424449904662437604539103637539339468499990506897971814716522960950895471907013503223284820116877757184319629345000591535688675");


foreach my $k(1..1e6) {
    my $x = c($k);
    my $y = c($k + $period);

    if ($x != $y) {
        die "non period: $x != $y for k = $k";
    }
}

__END__
#~ foreach my $k(1..1e8) {
    #~ if (c($k) == 1247) {
        #~ say "c($k) = 1247";
        #~ last;
    #~ }
#~ }
#~ __END__

#say c(2153846);
#~ use Math::GMPz;
#~ say c(Math::GMPz->new("560500033"));
#~ say c(Math::GMPz->new("13668555257082559489"));
#~ say c(Math::GMPz->new("9498733884118309169"));
#~ say c(Math::GMPz->new("9621022700125202321"));


#~ __END__

my %seen;
use Math::GMPz;

my %known;
@known{

9, 15, 21, 25, 33, 39, 49, 51, 57, 65, 69, 85, 87, 91, 93, 111, 121, 123, 129, 133, 141, 145, 159, 169, 177, 183, 185, 201, 205, 213, 217, 219, 237, 249, 259, 265, 267, 289, 291, 301, 303, 305, 309, 321, 327, 339, 341, 361, 365, 381, 393, 411, 417, 427, 445, 447, 451, 453, 469, 471, 481, 485, 489, 501, 505, 511, 519, 529, 537, 543, 545, 553, 561, 565, 573, 579, 591, 597, 633, 645, 669, 671, 679, 681, 685, 687, 699, 703, 717, 721, 723, 745, 753, 763, 771, 781, 785, 789, 793, 807, 831, 841, 849, 865, 879, 889, 905, 921, 933, 939, 949, 951, 961, 965, 973, 985, 1041, 1057, 1059, 1065, 1099, 1101, 1111, 1137, 1141, 1145, 1149, 1203, 1247, 1261, 1263, 1267, 1299, 1317, 1329, 1345, 1369, 1385, 1387, 1393, 1417, 1441, 1461, 1477, 1541, 1561, 1565, 1585, 1661, 1685, 1687, 1713, 1729,

1681,
1689,
1383,
1119,
1437,
1191,
1649,
1011,
1371,
1671,
1167,
1285,
1293,
843,
1077,
993,
 1473,
 813,
 1047,
 1509,
 1257,
 1389,
 1165,
 1351,
 1227,
 1465,
 1347,
 1405,
 1527,
 1569,
 1707,
 1497,
 1563,
 1401,
 1603,
 1641,
 1205,
 1623,



    } = ();


foreach my $k(1..1e8) {

    if (not exists $known{c($k)}) {
        say "Found: a($k) = ", c($k);
        $known{c($k)} = $k;
    }
}

__END__

use Math::AnyNum;
    say c(Math::AnyNum->new("13668555257082559489"));

    exit;

    #~ say join ', ', keys %known;
    #~ exit;

    #~ foreach my $n(1..1e8) {
    #~ if (not exists $known{c($n)}) {
        #~ say "Found: $n -- ", c($n);
        #~ $known{c($n)} = $n;
    #~ }
#~ }

#~ __END__

    #say scalar keys %known;
#    %known = ();

#~ foreach my $n(1e8..1e9) {
    #~ if (not exists $known{c($n)}) {
        #~ say "Found: $n -- ", c($n);
        #~ $known{c($n)} = $n;
    #~ }
#~ }

#~ __END__

@seen{keys %known} = ();
#foreach my $k(1..2e6) {
#    $seen{c($k)}   = $k;
#}

while (<>) {

    next if /^#/;
    next if !/\S/;

    my $n = (split(' '))[-1];
    $n || next;

    #if ($n > ~0) {
        $n = Math::GMPz->new("$n");
  #  }

    #say "Testing: $n";

    my $t = c($n);

    if (not exists $seen{$t}) {
        say "New: a($n) = $t";
        $seen{$t} = $n;
    }
}




#foreach my $k(sort {$a <=> $b }keys %seen) {
#    say "a($seen{$k}) = $k";
#}

say "There are ", scalar(keys %seen), " unique values";

foreach my $k(keys %seen) {
    if (not exists $known{$k}) {
        say "New: $k  (for $seen{$k})";
    }
}

say join ', ', sort {$a <=> $b} keys %seen;

__END__

New: 119  (for 13668555257082559489)
New: 1267  (for 560500033)
New: 35  (for 9498733884118309169)
New: 27  (for 9621022700125202321)


#my $max = 0;

#~ say c(0);

#~ __END__

my %seen;
foreach my $k(1..2e6) {
    undef $seen{c($k)};
    #my $t = c($k);
    #if ($t > $max) {
    #    $max = $t;
    #    say "a($k) = $t";
    #}
}

say join ', ', sort {$a <=> $b} keys %seen;

say scalar keys %seen;

# A = {9, 15, 21, 25, 33, 39, 49, 51, 57, 65, 69, 85, 87, 91, 93, 111, 121, 123, 129, 133, 141, 145, 159, 169, 177, 183, 185, 201, 205, 213, 217, 219, 237, 249, 259, 265, 267, 289, 291, 301, 303, 305, 309, 321, 327, 339, 341, 361, 365, 381, 393, 411, 417, 427, 445, 447, 451, 453, 469, 471, 481, 485, 489, 501, 505, 511, 519, 529, 537, 543, 545, 553, 561, 565, 573, 579, 591, 597, 633, 645, 671, 679, 681, 685, 687, 699, 703, 717, 721, 723, 745, 753, 763, 771, 781, 785, 789, 793, 831, 841, 843, 849, 865, 879, 889, 905, 921, 933, 939, 949, 951, 961, 965, 973, 985, 993, 1011, 1041, 1057, 1059, 1065, 1077, 1099, 1101, 1111, 1119, 1137, 1141, 1145, 1167, 1191, 1247, 1261, 1263, 1285, 1293, 1329, 1369, 1371, 1383, 1387, 1417, 1437, 1441, 1461, 1477, 1541, 1561, 1585, 1649, 1661, 1671, 1681, 1689, 1729}

__END__

use Math::GMPz;

my $period = lcm(map{find_q($_)}(2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 139, 149, 163, 167));
$period = Math::GMPz->new("$period");

foreach my $p(2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 139, 149, 163, 167) {
    say "Q($p) = ", find_q($p);
}

__END__

$period/=499;
$period/=739;
$period/=743;

say $period;
#~ exit;
#~ use Math::AnyNum qw(irand);


foreach my $n(1..100000) {

    #my $n = irand(factorial($k));

    my $x = a($n);
    my $y = a($n+$period);

    say "Testing: $n  ($x == $y)";

    if ($x != $y) {
        die "Not a period: $x != $y for n = $n";
    }
}

__END__
use Math::GMPz;

my $max = 0;

while (<>) {

    next if /^#/;
    next if !/\S/;

    my $n = (split(' '))[-1];
    $n || next;

    if ($n > ~0) {
        $n = Math::GMPz->new("$n");
    }

    #my ($p, $q) = a($n);
    my $x = a($n);
    my $y = a($n+$period);

    say "$x $y";

    if ($x != $y) {
        die "Period fail: $x != $y for n = $n";
    }

    #if ($p > $max) {
    #    $max = $p;
    #    say "a($n) = $p (with q = $q)";
    #}
}

__END__

foreach my $n(0..1e4) {
    say a($n);
}

__END__

#~ my %seen;
#~ foreach my $k(1..1e7) {
    #~ undef $seen{a($k)};
#~ }
#~ say join ', ', sort {$a <=> $b} keys %seen;


#~ __END__

#say join ', ', map{a($_)} 0..10;
#say join ', ', map{b($_)} 0..10;

#OUTER: foreach my $k(1..1e6) {
use Math::GMPz;
my $k = Math::GMPz->new("59135093664847200");

    say "Testing: $k";
    foreach my $n(1..1e6) {
        my $x = a($n);
        my $y = a($n+$k);

        if ($x != $y) {
            die "Not a period $x != $y for n = $n";
        }
    }

    say "Period seems to be $k";
 #   last;
#}

__END__

my $max = 163;

foreach my $k(1e7..1e8) {
    my ($t, $q) = a($k);
    if ($t > $max) {
        $max = $t;
        say "a($k) = $t (with q = $q)";
    }
}

__END__

my %seen;
foreach my $k(0..1e4) {
    undef $seen{b($k)};
}

say join ', ', sort {$a <=> $b} keys %seen;