experimental-projects/pseudoprimes/programs/lehmer_totient_weak_numbers_cached.pl

#!/usr/bin/perl

# Composite numbers n such that phi(n) divides p*(n - 1) for some prime factor p of n-1.
# https://oeis.org/A338998

# Known terms:
#   1729, 12801, 5079361, 34479361, 3069196417

use 5.020;
use strict;
use warnings;
use experimental qw(signatures);

use Storable;
use Math::GMPz;
use ntheory qw(:all);
use Math::Prime::Util::GMP;
use experimental qw(signatures);
use Math::Sidef  qw(trial_factor);
use List::Util   qw(uniq);
use POSIX        qw(ULONG_MAX);

my $carmichael_file = "cache/factors-carmichael.storable";
my $carmichael      = retrieve($carmichael_file);

#~ sub my_euler_phi ($factors) {    # assumes n is squarefree
#~ Math::Prime::Util::GMP::vecprod(map { ($_ < ~0) ? ($_ - 1) : Math::Prime::Util::GMP::subint($_, 1) } @$factors);
#~ }

sub my_euler_phi ($factors) {    # assumes n is squarefree

    state $t = Math::GMPz::Rmpz_init();
    state $u = Math::GMPz::Rmpz_init();

    Math::GMPz::Rmpz_set_ui($t, 1);

    foreach my $p (@$factors) {
        if ($p < ULONG_MAX) {
            Math::GMPz::Rmpz_mul_ui($t, $t, $p - 1);
        }
        else {
            Math::GMPz::Rmpz_set_str($u, $p, 10);
            Math::GMPz::Rmpz_sub_ui($u, $u, 1);
            Math::GMPz::Rmpz_mul($t, $t, $u);
        }
    }

    return $t;
}

sub is_smooth_over_prod ($n, $k) {

    state $g = Math::GMPz::Rmpz_init_nobless();
    state $t = Math::GMPz::Rmpz_init_nobless();

    Math::GMPz::Rmpz_set($t, $n);
    Math::GMPz::Rmpz_gcd($g, $t, $k);

    while (Math::GMPz::Rmpz_cmp_ui($g, 1) > 0) {
        Math::GMPz::Rmpz_remove($t, $t, $g);
        return 1 if Math::GMPz::Rmpz_cmp_ui($t, 1) == 0;
        Math::GMPz::Rmpz_gcd($g, $t, $g);
    }

    return 0;
}

my $t   = Math::GMPz::Rmpz_init();
my $nm1 = Math::GMPz::Rmpz_init();

while (my ($key, $value) = each %$carmichael) {

    # length($key) <= 30 or next;

    Math::GMPz::Rmpz_set_str($nm1, $key, 10);
    Math::GMPz::Rmpz_sub_ui($nm1, $nm1, 1);

    my @factors = split(' ', $value);
    my $phi     = my_euler_phi(\@factors);

    Math::GMPz::Rmpz_gcd($t, $phi, $nm1);
    Math::GMPz::Rmpz_divexact($t, $phi, $t);

    #~ if (Math::GMPz::Rmpz_divisible_p($nm1, $t)) {
    if (Math::GMPz::Rmpz_divisible_p($nm1, $t)) {
        if (is_prime(Math::GMPz::Rmpz_get_str($t, 10))) {
            say "\nFound term: $key\n";
        }
        else {
            say "$t -> ", $key;
        }
    }
}

__END__

# No terms > 2^64 are known.

# If we relax the problem such that k is allowed to be composite, then we have (k -> n):

256 -> 161053591977104597648385
4320 -> 1442607709754537310156481
276480 -> 121311396642375203328001
36664320 -> 18307378995828443136001
39243204 -> 522701579776994601985
164647296 -> 710498768604458891905
239368932 -> 166767515214109647553
239855616 -> 69331699150300274689
344198160 -> 26220447360756492481
470271360 -> 22577695684076482561
494534656 -> 9370403366634042163201
690069744 -> 26412091255782290298241
847888206 -> 508856401419300817312321
1008611352 -> 225751988643126523201
1022867120 -> 151320171400394261761
1270746144 -> 1973514672419436576001
1355878368 -> 25043558587894400823361
1713960000 -> 405475097425913520001
2595175200 -> 79919351614524024001
3936306240 -> 22695373452873769921
4414596480 -> 2819631734736201434649601
4766062840 -> 305022961592329915047961
5256049760 -> 1454162635820073015201
6723894240 -> 23152815997129790193601
9337407520 -> 588513459619254160321
12358835370 -> 271556534518852228430401
14172416640 -> 57199945443705577825921
14385833280 -> 8449460360571792219841
20042190000 -> 3976684712865158490001
22016352600 -> 43988320199517413227201
23021417472 -> 2484569421264140859210670081
35963136000 -> 3221275513469123100672001
42830346240 -> 38223250373318642652303361
47342232000 -> 2479639741736856792912001
100711683072 -> 172991335211656109051905
225735802880 -> 24454143142535092272994713601
299747952000 -> 3373971932766269202480001
521132748000 -> 24305672075484052148160001
638306611200 -> 244560890510134115626598401
648369446400 -> 21566311675749133571835648001
768266286228 -> 19674209294344913441970817
938064067200 -> 33803763411350231745638401
1050025132800 -> 4775814557195289074611201
1127001600000 -> 77755530778217753011200001
1231491998240 -> 1585372177726779981269921
1299266121600 -> 566296873170734780865792001
1712586456000 -> 18856002135543697675152001
3473510584320 -> 256551983459660482787450881
6546891384000 -> 60060884070685995775533024001
7062492272640 -> 169302315858010163156179353601
10390842000000 -> 224159025932793824796000001
11788391160000 -> 1836659845422904698410760001
14246648367828 -> 7429200724594400259443892735745
20751621203808 -> 6753985548147606751087157673985
28322297280000 -> 22584599278031754549573120001
43712316942336 -> 5423553691041141223181881345
53989739274240 -> 26468027537347458252574433281
144645490376640 -> 22031534549990632963014227521
194461205316000 -> 16267018441058472048178342860001
231593852160000 -> 280329366963716628966650880001
1004376125491200 -> 851213068769837939865599760537601