Skip to content

Instantly share code, notes, and snippets.

View Hermann-SW's full-sized avatar

Hermann Stamm-Wilbrandt Hermann-SW

89 followers · 0 following
View GitHub Profile
@Hermann-SW
Hermann-SW / 23.gp
Last active March 5, 2026 15:18
There are no further Carmichael numbers N=2^a*3^b+1 below 10^70 (than 1729=2^6*3^3+1 and 46656=2^6*3^6+1)
is_carmichael_minus_1(f)={
n=factorback(f)+1;
v=[d+1|d<-divisors(n-1),n%(d+1)==0&&isprime(d+1)];
vecprod(v)==n; \\ Korselt's criterion
}
m=10^70;
{
for(a=1,oo,
if(2^a<=m,
@Hermann-SW
Hermann-SW / Car_n-1_3_prime_factors.gp
Last active March 8, 2026 07:39
Prime factorization of N-1 having exactly 3 prime factors, for Carchmichael numbers N ≤ 10^24
{[ [2, 4; 5, 1; 7, 1],
[2, 4; 3, 1; 23, 1],
[2, 5; 7, 1; 11, 1],
[2, 3; 3, 3; 7, 2],
[2, 2; 3, 2; 1777, 1],
[2, 3; 3, 2; 1753, 1],
[2, 4; 3, 3; 1733, 1],
[2, 8; 3, 2; 433, 1],
[2, 3; 3, 5; 557, 1],
[2, 3; 3, 3; 23, 3],
@Hermann-SW
Hermann-SW / Car_3.343dd.gp
Last active March 2, 2026 08:19
3 prime factor Charmichael number examples for roughly every 3 decimal digits up to 343
assert(b)=if(!(b),error());
factmul(f1,f2)=matreduce(matconcat([f1,f2]~));
factval(F)=vecprod([v[1]^v[2]|v<-F~]);
{Redu=[0,25,0,25,110,291,51,146,131,511,111,95,1121,2685,820,12481,16175,1866,
4500,11525,8960,441,390,14796,1280,1651,1730,24140,21226,18555,43391,3716,2980,
46701,38580,15450,5560,19445,14376,83660,32560,7516,5060,23806,57806,44636,
28985,73445,60936,55146,91400,82190,54255,8016,25591,71945,259946,147035,11301,
3375,2371,18486,466191,436551,422806,6220,153406,493275,222755,1572896,453141,
5385,422511,663666,364225,84081,52590,916505,285466,827301,5671,137266,120160,
@Hermann-SW
Hermann-SW / 96522.gp
Created February 24, 2026 20:42
@neptune's largest known 96,522 decimal digits 3-Carmichael number, with single random base verification
\\ @Neptune's largest known 3-Carmichael number (96522 decimal digits):
\\ https://www.mersenneforum.org/node/22080/page3#post1066763
p = 3*(5752211*43#/2-1)^1069/2+1;
q = 3*(5752211*43#/2-1)^1069+1;
r = 3*((5752211*43#/2-1)^1069+(5752211*43#/2-1)^2138)/1050650772710+1;
n = p*q*r;
print(#digits(n));
n1 = (n-1)/gcd(n-1,p-1);
@Hermann-SW
Hermann-SW / 3710.gp
Last active February 25, 2026 04:35
3710 decimal digits Carmichael number from 1989 Dubner paper table 1
\\ 3710 decimal digits Carmichael number from 1989 Dubner paper table 1:
\\ https://www.ams.org/journals/mcom/1989-53-187/S0025-5718-1989-0969484-8/S0025-5718-1989-0969484-8.pdf
\\
T=47#/2;A=41;C=141847;M=(T*C-1)^A/4;P=6*M+1;Q=12*M+1;X=123165;R=1+(P*Q-1)/X;
N=P*Q*R;
\\ known partial N-1 factorization: https://www.mersenneforum.org/node/1106127
{F=[2,41;3,1;11,1;13,1;19,1;29,1;31,1;37,1;41,2;43,1;47,1;59,41;79,41;83,1;
1709,1;3527,1;3691,1;16943,1;469793,41;1799411527,1;3463701403,1;
731646295847,1;9957992526379,41;677868618879887,1;278798236535678281,1;
61534897980248555544581,1;9929897004627382451681972907710143,1];}
@Hermann-SW
Hermann-SW / par2.gp
Created February 16, 2026 14:16
Parallel determination of primitive root of n-th Euclid prime
check(En, n, phi, g, i)={
parfor(i=1, n, lift(Mod(g, En)^(phi / prime(i))) == 1, r, if(r, return(0)));
1
}
euclid_prime_find_root(n) = {
my(En = vecprod(primes(n)) + 1, phi = En - 1);
for(g=2, En-1,
if(lift(kronecker(g, En))==-1,
@Hermann-SW
Hermann-SW / seq.gp
Created February 16, 2026 14:13
Sequential determination of primitive root of n-th Euclid prime
euclid_prime_find_root(n) = {
my(En = vecprod(primes(n)) + 1, phi = En - 1);
for(g=2, En-1,
if(lift(kronecker(g, En))==-1,
my(is_root = 1);
for(i=1, n,
if(lift(Mod(g, En)^(phi / prime(i))) == 1,
is_root = 0;
break;
@Hermann-SW
Hermann-SW / a11.3399714628553118047.gp
Created January 16, 2026 16:54
Compute number of primes of form x^2 + x + A for x <= 10^11 (for Jacobson and Williams A = 3399714628553118047)
export(A=3399714628553118047); a(n)=parsum(k=0, 10^n, isprime(k^2+k+A));
default(threadsizemax,1000*10^6); \\ default 8MB was too small
#
print(a(11));
@Hermann-SW
Hermann-SW / a11.A331876.gp
Created January 15, 2026 07:21
Script to compute a(11) for oeis sequence A331876 (Euler polynomial)
a(n)=parsum(k=0, 10^n, isprime(k^2 + k + 41));
default(threadsizemax,1000*10^6); \\ default 8MB was too small
#
print(a(11));
@Hermann-SW
Hermann-SW / a11.gp
Last active January 10, 2026 14:50
Script to verify a(11) for oeis sequence A392244
a(n)=parsum(k=1, 10^n, isprime(k^2+(k+1)^2));
default(threadsizemax,100*10^6); \\ default 8MB was too small
#
print(a(11));