Fermat pseudoprimes
You are encouraged to solve this task according to the task description, using any language you may know.
A Fermat pseudoprime is a positive composite integer that passes the Fermat primality test.
Fermat's little theorem states that if p is prime and a is coprime to p, then ap−1 − 1 is divisible by p.
For an integer a > 1, if a composite integer x evenly divides ax−1 − 1, then x is called a Fermat pseudoprime to base a.
Fermat pseudoprimes to base 2 are sometimes called Sarrus numbers or Poulet numbers.
Fermat pseudoprimes can be found to any positive integer base. When using a base integer a = 1, this method returns all composite numbers.
- Task
For base integers a of 1 through 20:
- Find the count of pseudoprimes up to and including 12,000.
- Show the first 20 pseudoprimes.
- Stretch
- Extend the count threshold out to 25,000, 50,000 or higher.
- See also
- Wikipedia: Fermat pseudoprime
- OEIS:A002808 - Composite numbers
- OEIS:A001567 - Fermat pseudoprimes to base 2
- OEIS:A005935 - Fermat pseudoprimes to base 3
- OEIS:A020136 - Fermat pseudoprimes to base 4
- OEIS:A005936 - Fermat pseudoprimes to base 5
- OEIS:A005937 - Fermat pseudoprimes to base 6
- OEIS:A005938 - Fermat pseudoprimes to base 7
- OEIS:A020137 - Fermat pseudoprimes to base 8
- OEIS:A020138 - Fermat pseudoprimes to base 9
- OEIS:A005939 - Fermat pseudoprimes to base 10
- OEIS:A020139 - Fermat pseudoprimes to base 11
- OEIS:A020140 - Fermat pseudoprimes to base 12
- OEIS:A020141 - Fermat pseudoprimes to base 13
- OEIS:A020142 - Fermat pseudoprimes to base 14
- OEIS:A020143 - Fermat pseudoprimes to base 15
- OEIS:A020144 - Fermat pseudoprimes to base 16
- OEIS:A020145 - Fermat pseudoprimes to base 17
- OEIS:A020146 - Fermat pseudoprimes to base 18
- OEIS:A020147 - Fermat pseudoprimes to base 19
- OEIS:A020148 - Fermat pseudoprimes to base 20
- Related
ALGOL 68
BEGIN # find some Fermat pseudo primes: x is a Fermat pseudoprime in base a #
# if a^(x-1) - 1 is divisible by x and x is not prime #
# iterative Greatest Common Divisor routine, returns the gcd of m and n #
PROC gcd = ( INT m, n )INT:
BEGIN
INT a := ABS m, b := ABS n;
WHILE b /= 0 DO
INT new a = b;
b := a MOD b;
a := new a
OD;
a
END # gcd # ;
[ 0 : 50 000 ]BOOL prime;
prime[ 0 ] := prime[ 1 ] := FALSE;
prime[ 2 ] := TRUE;
FOR i FROM 3 BY 2 TO UPB prime DO prime[ i ] := TRUE OD;
FOR i FROM 4 BY 2 TO UPB prime DO prime[ i ] := FALSE OD;
FOR i FROM 3 BY 2 TO ENTIER sqrt( UPB prime ) DO
IF prime[ i ] THEN
FOR s FROM i * i BY i + i TO UPB prime DO prime[ s ] := FALSE OD
FI
OD;
FOR base FROM 1 TO 20 DO
[ 1 : 20 ]INT fp;
INT count := 0;
# x from 3 as 1 is neither composite nor prime and 2 is prime #
IF base = 1 THEN
# 1^(x-1) is 1 for all x, so all composites are #
# fermat pseudo primes in base 1 #
FOR x FROM 3 TO UPB prime DO
IF NOT prime[ x ] THEN
IF ( count +:= 1 ) <= UPB fp THEN
fp[ count ] := x
FI
FI
OD
ELSE
# must test base^(x-1) mod x #
FOR x FROM 3 BY IF ODD base THEN 1 ELSE 2 FI TO UPB prime DO
IF NOT prime[ x ] THEN
IF gcd( x, base ) = 1 THEN
# have a composite x co-prime to the base #
INT base to x minus 1 mod x := 1;
FOR p TO x - 1
WHILE
base to x minus 1 mod x *:= base MODAB x;
IF base to x minus 1 mod x = 1
THEN # the power is now 1 mod x, #
# if the power divides x - 1 then base^(x-1) #
# will be 1 mod x otherwise it won't be #
IF ( x - 1 ) MOD p /= 0 THEN
base to x minus 1 mod x := 0
FI;
FALSE
ELSE TRUE # keep checking the powers #
FI
DO SKIP OD;
IF base to x minus 1 mod x = 1 THEN
# have a composite x that divides base^(x-1)-1 #
IF ( count +:= 1 ) <= UPB fp THEN
fp[ count ] := x
FI
FI
FI
FI
OD
FI;
print( ( "base ", whole( base, -2 ), " to ", whole( UPB prime, 0 )
, " total: ", whole( count, -6 ), ", first: "
)
);
FOR f FROM LWB fp TO IF count < UPB fp THEN count ELSE UPB fp FI DO
print( ( " ", whole( fp[ f ], 0 ) ) )
OD;
print( ( newline ) )
OD
END- Output:
base 1 to 50000 total: 44866, first: 4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32 base 2 to 50000 total: 55, first: 341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 4681 5461 6601 7957 8321 base 3 to 50000 total: 53, first: 91 121 286 671 703 949 1105 1541 1729 1891 2465 2665 2701 2821 3281 3367 3751 4961 5551 6601 base 4 to 50000 total: 111, first: 15 85 91 341 435 451 561 645 703 1105 1247 1271 1387 1581 1695 1729 1891 1905 2047 2071 base 5 to 50000 total: 54, first: 4 124 217 561 781 1541 1729 1891 2821 4123 5461 5611 5662 5731 6601 7449 7813 8029 8911 9881 base 6 to 50000 total: 74, first: 35 185 217 301 481 1105 1111 1261 1333 1729 2465 2701 2821 3421 3565 3589 3913 4123 4495 5713 base 7 to 50000 total: 49, first: 6 25 325 561 703 817 1105 1825 2101 2353 2465 3277 4525 4825 6697 8321 10225 10585 10621 11041 base 8 to 50000 total: 150, first: 9 21 45 63 65 105 117 133 153 231 273 341 481 511 561 585 645 651 861 949 base 9 to 50000 total: 113, first: 4 8 28 52 91 121 205 286 364 511 532 616 671 697 703 946 949 1036 1105 1288 base 10 to 50000 total: 65, first: 9 33 91 99 259 451 481 561 657 703 909 1233 1729 2409 2821 2981 3333 3367 4141 4187 base 11 to 50000 total: 61, first: 10 15 70 133 190 259 305 481 645 703 793 1105 1330 1729 2047 2257 2465 2821 4577 4921 base 12 to 50000 total: 91, first: 65 91 133 143 145 247 377 385 703 1045 1099 1105 1649 1729 1885 1891 2041 2233 2465 2701 base 13 to 50000 total: 68, first: 4 6 12 21 85 105 231 244 276 357 427 561 1099 1785 1891 2465 2806 3605 5028 5149 base 14 to 50000 total: 69, first: 15 39 65 195 481 561 781 793 841 985 1105 1111 1541 1891 2257 2465 2561 2665 2743 3277 base 15 to 50000 total: 42, first: 14 341 742 946 1477 1541 1687 1729 1891 1921 2821 3133 3277 4187 6541 6601 7471 8701 8911 9073 base 16 to 50000 total: 145, first: 15 51 85 91 255 341 435 451 561 595 645 703 1105 1247 1261 1271 1285 1387 1581 1687 base 17 to 50000 total: 63, first: 4 8 9 16 45 91 145 261 781 1111 1228 1305 1729 1885 2149 2821 3991 4005 4033 4187 base 18 to 50000 total: 98, first: 25 49 65 85 133 221 323 325 343 425 451 637 931 1105 1225 1369 1387 1649 1729 1921 base 19 to 50000 total: 93, first: 6 9 15 18 45 49 153 169 343 561 637 889 905 906 1035 1105 1629 1661 1849 1891 base 20 to 50000 total: 66, first: 21 57 133 231 399 561 671 861 889 1281 1653 1729 1891 2059 2413 2501 2761 2821 2947 3059
C++
#include <cstdint>
#include <iomanip>
#include <iostream>
uint64_t modpow(uint64_t base, uint64_t exp, uint64_t mod) {
if (mod == 1)
return 0;
uint64_t result = 1;
base %= mod;
for (; exp > 0; exp >>= 1) {
if ((exp & 1) == 1)
result = (result * base) % mod;
base = (base * base) % mod;
}
return result;
}
bool is_prime(uint64_t n) {
if (n < 2)
return false;
if (n % 2 == 0)
return n == 2;
if (n % 3 == 0)
return n == 3;
for (uint64_t p = 5; p * p <= n; p += 4) {
if (n % p == 0)
return false;
p += 2;
if (n % p == 0)
return false;
}
return true;
}
bool is_fermat_pseudoprime(uint64_t a, uint64_t x) {
return !is_prime(x) && modpow(a, x - 1, x) == 1;
}
int main() {
std::cout << "First 20 Fermat pseudoprimes:\n";
for (uint64_t a = 1; a <= 20; ++a) {
std::cout << "Base " << std::setw(2) << a << ": ";
int count = 0;
for (uint64_t x = 4; count < 20; ++x) {
if (is_fermat_pseudoprime(a, x)) {
++count;
std::cout << std::setw(5) << x << ' ';
}
}
std::cout << '\n';
}
const uint64_t limits[] = {12000, 25000, 50000, 100000};
std::cout << "\nCount <= ";
for (uint64_t limit : limits) {
std::cout << std::setw(6) << limit << ' ';
}
std::cout << "\n------------------------------------\n";
for (uint64_t a = 1; a <= 20; ++a) {
std::cout << "Base " << std::setw(2) << a << ": ";
int count = 0;
uint64_t x = 4;
for (uint64_t limit : limits) {
for (; x <= limit; ++x) {
if (is_fermat_pseudoprime(a, x))
++count;
}
std::cout << std::setw(6) << count << ' ';
}
std::cout << '\n';
}
}
- Output:
First 20 Fermat pseudoprimes: Base 1: 4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32 Base 2: 341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 4681 5461 6601 7957 8321 Base 3: 91 121 286 671 703 949 1105 1541 1729 1891 2465 2665 2701 2821 3281 3367 3751 4961 5551 6601 Base 4: 15 85 91 341 435 451 561 645 703 1105 1247 1271 1387 1581 1695 1729 1891 1905 2047 2071 Base 5: 4 124 217 561 781 1541 1729 1891 2821 4123 5461 5611 5662 5731 6601 7449 7813 8029 8911 9881 Base 6: 35 185 217 301 481 1105 1111 1261 1333 1729 2465 2701 2821 3421 3565 3589 3913 4123 4495 5713 Base 7: 6 25 325 561 703 817 1105 1825 2101 2353 2465 3277 4525 4825 6697 8321 10225 10585 10621 11041 Base 8: 9 21 45 63 65 105 117 133 153 231 273 341 481 511 561 585 645 651 861 949 Base 9: 4 8 28 52 91 121 205 286 364 511 532 616 671 697 703 946 949 1036 1105 1288 Base 10: 9 33 91 99 259 451 481 561 657 703 909 1233 1729 2409 2821 2981 3333 3367 4141 4187 Base 11: 10 15 70 133 190 259 305 481 645 703 793 1105 1330 1729 2047 2257 2465 2821 4577 4921 Base 12: 65 91 133 143 145 247 377 385 703 1045 1099 1105 1649 1729 1885 1891 2041 2233 2465 2701 Base 13: 4 6 12 21 85 105 231 244 276 357 427 561 1099 1785 1891 2465 2806 3605 5028 5149 Base 14: 15 39 65 195 481 561 781 793 841 985 1105 1111 1541 1891 2257 2465 2561 2665 2743 3277 Base 15: 14 341 742 946 1477 1541 1687 1729 1891 1921 2821 3133 3277 4187 6541 6601 7471 8701 8911 9073 Base 16: 15 51 85 91 255 341 435 451 561 595 645 703 1105 1247 1261 1271 1285 1387 1581 1687 Base 17: 4 8 9 16 45 91 145 261 781 1111 1228 1305 1729 1885 2149 2821 3991 4005 4033 4187 Base 18: 25 49 65 85 133 221 323 325 343 425 451 637 931 1105 1225 1369 1387 1649 1729 1921 Base 19: 6 9 15 18 45 49 153 169 343 561 637 889 905 906 1035 1105 1629 1661 1849 1891 Base 20: 21 57 133 231 399 561 671 861 889 1281 1653 1729 1891 2059 2413 2501 2761 2821 2947 3059 Count <= 12000 25000 50000 100000 ------------------------------------ Base 1: 10561 22237 44866 90407 Base 2: 25 38 55 78 Base 3: 25 38 53 78 Base 4: 50 75 111 153 Base 5: 22 35 54 73 Base 6: 31 46 74 104 Base 7: 21 32 49 73 Base 8: 76 110 150 218 Base 9: 55 83 113 164 Base 10: 35 53 65 90 Base 11: 30 41 61 89 Base 12: 35 60 91 127 Base 13: 31 51 68 91 Base 14: 33 51 69 96 Base 15: 22 31 42 70 Base 16: 69 102 145 200 Base 17: 31 44 63 83 Base 18: 46 69 98 134 Base 19: 48 70 93 121 Base 20: 35 49 66 101
F#
This task uses Extensible Prime Generator (F#)
// Fermat pseudoprimes. Nigel Galloway: August 17th., 2022
let fp(a:int)=let a=bigint a in primesI()|>Seq.pairwise|>Seq.collect(fun(n,g)->seq{for n in n+1I..g-1I do if bigint.ModPow(a,n-1I,n)=1I then yield n})
{1..20}|>Seq.iter(fun n->printf $"Base %2d{n} - Up to 50000: %5d{fp n|>Seq.takeWhile((>=)50000I)|>Seq.length} First 20: ("; fp n|>Seq.take 20|>Seq.iter(printf "%A "); printfn ")")
- Output:
Base 1 - Up to 50000: 44866 First 20: (4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32 ) Base 2 - Up to 50000: 55 First 20: (341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 4681 5461 6601 7957 8321 ) Base 3 - Up to 50000: 53 First 20: (91 121 286 671 703 949 1105 1541 1729 1891 2465 2665 2701 2821 3281 3367 3751 4961 5551 6601 ) Base 4 - Up to 50000: 111 First 20: (15 85 91 341 435 451 561 645 703 1105 1247 1271 1387 1581 1695 1729 1891 1905 2047 2071 ) Base 5 - Up to 50000: 54 First 20: (4 124 217 561 781 1541 1729 1891 2821 4123 5461 5611 5662 5731 6601 7449 7813 8029 8911 9881 ) Base 6 - Up to 50000: 74 First 20: (35 185 217 301 481 1105 1111 1261 1333 1729 2465 2701 2821 3421 3565 3589 3913 4123 4495 5713 ) Base 7 - Up to 50000: 49 First 20: (6 25 325 561 703 817 1105 1825 2101 2353 2465 3277 4525 4825 6697 8321 10225 10585 10621 11041 ) Base 8 - Up to 50000: 150 First 20: (9 21 45 63 65 105 117 133 153 231 273 341 481 511 561 585 645 651 861 949 ) Base 9 - Up to 50000: 113 First 20: (4 8 28 52 91 121 205 286 364 511 532 616 671 697 703 946 949 1036 1105 1288 ) Base 10 - Up to 50000: 65 First 20: (9 33 91 99 259 451 481 561 657 703 909 1233 1729 2409 2821 2981 3333 3367 4141 4187 ) Base 11 - Up to 50000: 61 First 20: (10 15 70 133 190 259 305 481 645 703 793 1105 1330 1729 2047 2257 2465 2821 4577 4921 ) Base 12 - Up to 50000: 91 First 20: (65 91 133 143 145 247 377 385 703 1045 1099 1105 1649 1729 1885 1891 2041 2233 2465 2701 ) Base 13 - Up to 50000: 68 First 20: (4 6 12 21 85 105 231 244 276 357 427 561 1099 1785 1891 2465 2806 3605 5028 5149 ) Base 14 - Up to 50000: 69 First 20: (15 39 65 195 481 561 781 793 841 985 1105 1111 1541 1891 2257 2465 2561 2665 2743 3277 ) Base 15 - Up to 50000: 42 First 20: (14 341 742 946 1477 1541 1687 1729 1891 1921 2821 3133 3277 4187 6541 6601 7471 8701 8911 9073 ) Base 16 - Up to 50000: 145 First 20: (15 51 85 91 255 341 435 451 561 595 645 703 1105 1247 1261 1271 1285 1387 1581 1687 ) Base 17 - Up to 50000: 63 First 20: (4 8 9 16 45 91 145 261 781 1111 1228 1305 1729 1885 2149 2821 3991 4005 4033 4187 ) Base 18 - Up to 50000: 98 First 20: (25 49 65 85 133 221 323 325 343 425 451 637 931 1105 1225 1369 1387 1649 1729 1921 ) Base 19 - Up to 50000: 93 First 20: (6 9 15 18 45 49 153 169 343 561 637 889 905 906 1035 1105 1629 1661 1849 1891 ) Base 20 - Up to 50000: 66 First 20: (21 57 133 231 399 561 671 861 889 1281 1653 1729 1891 2059 2413 2501 2761 2821 2947 3059 ) Real: 00:00:00.632
FreeBASIC
'#include "isprime.bas"
Function powMod (a As Uinteger, n As Uinteger, m As Uinteger) As Uinteger
' Return "a^n mod m".
Dim As Uinteger a1 = a Mod m
Dim As Uinteger n1 = n
Dim As Uinteger result
If a1 > 0 Then
result = 1
While n1 > 0
If (n1 And 1) <> 0 Then result = (result * a1) Mod m
n1 Shr= 1
a1 = (a1 * a1) Mod m
Wend
End If
Return result
End Function
Function isFermatPseudoprime(x As Uinteger, a As Uinteger) As Boolean
Return Iif(isPrime(x), False, powMod(a, x - 1, x) = 1)
End Function
Dim As Uinteger a, b, limite = 50000
Print "First 20 Fermat pseudoprimes:"
For a = 1 To 20
Dim As Uinteger total = 0
Dim As Uinteger x = 2
Dim As Uinteger primer20(1 To 20)
For b = 1 To limite
While x <= b
If isFermatPseudoprime(x, a) Then
total += 1
If total <= 20 Then primer20(total) = x
End If
x += 1
Wend
Next
Print Using "Base ## to ##### total: #####_, first: "; a; limite; total;
For b = 1 To 20
Print Using "######"; primer20(b);
Next b
Print
Next a
Sleep- Output:
First 20 Fermat pseudoprimes: Base 1 to 50000 total: 44866, first: 4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32 Base 2 to 50000 total: 55, first: 341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 4681 5461 6601 7957 8321 Base 3 to 50000 total: 53, first: 91 121 286 671 703 949 1105 1541 1729 1891 2465 2665 2701 2821 3281 3367 3751 4961 5551 6601 Base 4 to 50000 total: 111, first: 15 85 91 341 435 451 561 645 703 1105 1247 1271 1387 1581 1695 1729 1891 1905 2047 2071 Base 5 to 50000 total: 54, first: 4 124 217 561 781 1541 1729 1891 2821 4123 5461 5611 5662 5731 6601 7449 7813 8029 8911 9881 Base 6 to 50000 total: 74, first: 35 185 217 301 481 1105 1111 1261 1333 1729 2465 2701 2821 3421 3565 3589 3913 4123 4495 5713 Base 7 to 50000 total: 49, first: 6 25 325 561 703 817 1105 1825 2101 2353 2465 3277 4525 4825 6697 8321 10225 10585 10621 11041 Base 8 to 50000 total: 150, first: 9 21 45 63 65 105 117 133 153 231 273 341 481 511 561 585 645 651 861 949 Base 9 to 50000 total: 113, first: 4 8 28 52 91 121 205 286 364 511 532 616 671 697 703 946 949 1036 1105 1288 Base 10 to 50000 total: 65, first: 9 33 91 99 259 451 481 561 657 703 909 1233 1729 2409 2821 2981 3333 3367 4141 4187 Base 11 to 50000 total: 61, first: 10 15 70 133 190 259 305 481 645 703 793 1105 1330 1729 2047 2257 2465 2821 4577 4921 Base 12 to 50000 total: 91, first: 65 91 133 143 145 247 377 385 703 1045 1099 1105 1649 1729 1885 1891 2041 2233 2465 2701 Base 13 to 50000 total: 68, first: 4 6 12 21 85 105 231 244 276 357 427 561 1099 1785 1891 2465 2806 3605 5028 5149 Base 14 to 50000 total: 69, first: 15 39 65 195 481 561 781 793 841 985 1105 1111 1541 1891 2257 2465 2561 2665 2743 3277 Base 15 to 50000 total: 42, first: 14 341 742 946 1477 1541 1687 1729 1891 1921 2821 3133 3277 4187 6541 6601 7471 8701 8911 9073 Base 16 to 50000 total: 145, first: 15 51 85 91 255 341 435 451 561 595 645 703 1105 1247 1261 1271 1285 1387 1581 1687 Base 17 to 50000 total: 63, first: 4 8 9 16 45 91 145 261 781 1111 1228 1305 1729 1885 2149 2821 3991 4005 4033 4187 Base 18 to 50000 total: 98, first: 25 49 65 85 133 221 323 325 343 425 451 637 931 1105 1225 1369 1387 1649 1729 1921 Base 19 to 50000 total: 93, first: 6 9 15 18 45 49 153 169 343 561 637 889 905 906 1035 1105 1629 1661 1849 1891 Base 20 to 50000 total: 66, first: 21 57 133 231 399 561 671 861 889 1281 1653 1729 1891 2059 2413 2501 2761 2821 2947 3059
J
fermat=. {{1 = x (y&|@^) <: y}}"0
(>: i. 20) ([ ,. fermat/ (# , 20&{.)@# ]) (#~ 0&p:) >: i. 50000
- Output:
1 44866 4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32 2 55 341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 4681 5461 6601 7957 8321 3 53 91 121 286 671 703 949 1105 1541 1729 1891 2465 2665 2701 2821 3281 3367 3751 4961 5551 6601 4 111 15 85 91 341 435 451 561 645 703 1105 1247 1271 1387 1581 1695 1729 1891 1905 2047 2071 5 54 4 124 217 561 781 1541 1729 1891 2821 4123 5461 5611 5662 5731 6601 7449 7813 8029 8911 9881 6 74 35 185 217 301 481 1105 1111 1261 1333 1729 2465 2701 2821 3421 3565 3589 3913 4123 4495 5713 7 49 6 25 325 561 703 817 1105 1825 2101 2353 2465 3277 4525 4825 6697 8321 10225 10585 10621 11041 8 150 9 21 45 63 65 105 117 133 153 231 273 341 481 511 561 585 645 651 861 949 9 113 4 8 28 52 91 121 205 286 364 511 532 616 671 697 703 946 949 1036 1105 1288 10 65 9 33 91 99 259 451 481 561 657 703 909 1233 1729 2409 2821 2981 3333 3367 4141 4187 11 61 10 15 70 133 190 259 305 481 645 703 793 1105 1330 1729 2047 2257 2465 2821 4577 4921 12 91 65 91 133 143 145 247 377 385 703 1045 1099 1105 1649 1729 1885 1891 2041 2233 2465 2701 13 68 4 6 12 21 85 105 231 244 276 357 427 561 1099 1785 1891 2465 2806 3605 5028 5149 14 69 15 39 65 195 481 561 781 793 841 985 1105 1111 1541 1891 2257 2465 2561 2665 2743 3277 15 42 14 341 742 946 1477 1541 1687 1729 1891 1921 2821 3133 3277 4187 6541 6601 7471 8701 8911 9073 16 145 15 51 85 91 255 341 435 451 561 595 645 703 1105 1247 1261 1271 1285 1387 1581 1687 17 63 4 8 9 16 45 91 145 261 781 1111 1228 1305 1729 1885 2149 2821 3991 4005 4033 4187 18 98 25 49 65 85 133 221 323 325 343 425 451 637 931 1105 1225 1369 1387 1649 1729 1921 19 93 6 9 15 18 45 49 153 169 343 561 637 889 905 906 1035 1105 1629 1661 1849 1891 20 66 21 57 133 231 399 561 671 861 889 1281 1653 1729 1891 2059 2413 2501 2761 2821 2947 3059
Julia
using Primes
ispseudo(n, base) = !isprime(n) && BigInt(base)^(n - 1) % n == 1
for b in 1:20
pseudos = filter(n -> ispseudo(n, b), 1:50000)
println("Base ", lpad(b, 2), " up to 50000: ", lpad(length(pseudos), 5), " First 20: ", pseudos[1:20])
end
- Output:
Base 1 up to 50000: 44866 First 20: [4, 6, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32] Base 2 up to 50000: 55 First 20: [341, 561, 645, 1105, 1387, 1729, 1905, 2047, 2465, 2701, 2821, 3277, 4033, 4369, 4371, 4681, 5461, 6601, 7957, 8321] Base 3 up to 50000: 53 First 20: [91, 121, 286, 671, 703, 949, 1105, 1541, 1729, 1891, 2465, 2665, 2701, 2821, 3281, 3367, 3751, 4961, 5551, 6601] Base 4 up to 50000: 111 First 20: [15, 85, 91, 341, 435, 451, 561, 645, 703, 1105, 1247, 1271, 1387, 1581, 1695, 1729, 1891, 1905, 2047, 2071] Base 5 up to 50000: 54 First 20: [4, 124, 217, 561, 781, 1541, 1729, 1891, 2821, 4123, 5461, 5611, 5662, 5731, 6601, 7449, 7813, 8029, 8911, 9881] Base 6 up to 50000: 74 First 20: [35, 185, 217, 301, 481, 1105, 1111, 1261, 1333, 1729, 2465, 2701, 2821, 3421, 3565, 3589, 3913, 4123, 4495, 5713] Base 7 up to 50000: 49 First 20: [6, 25, 325, 561, 703, 817, 1105, 1825, 2101, 2353, 2465, 3277, 4525, 4825, 6697, 8321, 10225, 10585, 10621, 11041] Base 8 up to 50000: 150 First 20: [9, 21, 45, 63, 65, 105, 117, 133, 153, 231, 273, 341, 481, 511, 561, 585, 645, 651, 861, 949] Base 9 up to 50000: 113 First 20: [4, 8, 28, 52, 91, 121, 205, 286, 364, 511, 532, 616, 671, 697, 703, 946, 949, 1036, 1105, 1288] Base 10 up to 50000: 65 First 20: [9, 33, 91, 99, 259, 451, 481, 561, 657, 703, 909, 1233, 1729, 2409, 2821, 2981, 3333, 3367, 4141, 4187] Base 11 up to 50000: 61 First 20: [10, 15, 70, 133, 190, 259, 305, 481, 645, 703, 793, 1105, 1330, 1729, 2047, 2257, 2465, 2821, 4577, 4921] Base 12 up to 50000: 91 First 20: [65, 91, 133, 143, 145, 247, 377, 385, 703, 1045, 1099, 1105, 1649, 1729, 1885, 1891, 2041, 2233, 2465, 2701] Base 13 up to 50000: 68 First 20: [4, 6, 12, 21, 85, 105, 231, 244, 276, 357, 427, 561, 1099, 1785, 1891, 2465, 2806, 3605, 5028, 5149] Base 14 up to 50000: 69 First 20: [15, 39, 65, 195, 481, 561, 781, 793, 841, 985, 1105, 1111, 1541, 1891, 2257, 2465, 2561, 2665, 2743, 3277] Base 15 up to 50000: 42 First 20: [14, 341, 742, 946, 1477, 1541, 1687, 1729, 1891, 1921, 2821, 3133, 3277, 4187, 6541, 6601, 7471, 8701, 8911, 9073] Base 16 up to 50000: 145 First 20: [15, 51, 85, 91, 255, 341, 435, 451, 561, 595, 645, 703, 1105, 1247, 1261, 1271, 1285, 1387, 1581, 1687] Base 17 up to 50000: 63 First 20: [4, 8, 9, 16, 45, 91, 145, 261, 781, 1111, 1228, 1305, 1729, 1885, 2149, 2821, 3991, 4005, 4033, 4187] Base 18 up to 50000: 98 First 20: [25, 49, 65, 85, 133, 221, 323, 325, 343, 425, 451, 637, 931, 1105, 1225, 1369, 1387, 1649, 1729, 1921] Base 19 up to 50000: 93 First 20: [6, 9, 15, 18, 45, 49, 153, 169, 343, 561, 637, 889, 905, 906, 1035, 1105, 1629, 1661, 1849, 1891] Base 20 up to 50000: 66 First 20: [21, 57, 133, 231, 399, 561, 671, 861, 889, 1281, 1653, 1729, 1891, 2059, 2413, 2501, 2761, 2821, 2947, 3059]
Mathematica/Wolfram Language
ispseudo[n_, base_] := (!PrimeQ[n]) && (PowerMod[base, n - 1, n] == 1)
Do[
pseudos = Select[Range[1, 50000], ispseudo[#, b] &];
Print["Base ", b, " up to 50000: ", Length[pseudos], " First 20: ", Take[pseudos, 20]],
{b, 1, 20}
]
- Output:
Base 1 up to 50000: 44866 First 20: {4, 6, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32}
Base 2 up to 50000: 55 First 20: {341, 561, 645, 1105, 1387, 1729, 1905, 2047, 2465, 2701, 2821, 3277, 4033, 4369, 4371, 4681, 5461, 6601, 7957, 8321}
Base 3 up to 50000: 53 First 20: {91, 121, 286, 671, 703, 949, 1105, 1541, 1729, 1891, 2465, 2665, 2701, 2821, 3281, 3367, 3751, 4961, 5551, 6601}
Base 4 up to 50000: 111 First 20: {15, 85, 91, 341, 435, 451, 561, 645, 703, 1105, 1247, 1271, 1387, 1581, 1695, 1729, 1891, 1905, 2047, 2071}
Base 5 up to 50000: 54 First 20: {4, 124, 217, 561, 781, 1541, 1729, 1891, 2821, 4123, 5461, 5611, 5662, 5731, 6601, 7449, 7813, 8029, 8911, 9881}
Base 6 up to 50000: 74 First 20: {35, 185, 217, 301, 481, 1105, 1111, 1261, 1333, 1729, 2465, 2701, 2821, 3421, 3565, 3589, 3913, 4123, 4495, 5713}
Base 7 up to 50000: 49 First 20: {6, 25, 325, 561, 703, 817, 1105, 1825, 2101, 2353, 2465, 3277, 4525, 4825, 6697, 8321, 10225, 10585, 10621, 11041}
Base 8 up to 50000: 150 First 20: {9, 21, 45, 63, 65, 105, 117, 133, 153, 231, 273, 341, 481, 511, 561, 585, 645, 651, 861, 949}
Base 9 up to 50000: 113 First 20: {4, 8, 28, 52, 91, 121, 205, 286, 364, 511, 532, 616, 671, 697, 703, 946, 949, 1036, 1105, 1288}
Base 10 up to 50000: 65 First 20: {9, 33, 91, 99, 259, 451, 481, 561, 657, 703, 909, 1233, 1729, 2409, 2821, 2981, 3333, 3367, 4141, 4187}
Base 11 up to 50000: 61 First 20: {10, 15, 70, 133, 190, 259, 305, 481, 645, 703, 793, 1105, 1330, 1729, 2047, 2257, 2465, 2821, 4577, 4921}
Base 12 up to 50000: 91 First 20: {65, 91, 133, 143, 145, 247, 377, 385, 703, 1045, 1099, 1105, 1649, 1729, 1885, 1891, 2041, 2233, 2465, 2701}
Base 13 up to 50000: 68 First 20: {4, 6, 12, 21, 85, 105, 231, 244, 276, 357, 427, 561, 1099, 1785, 1891, 2465, 2806, 3605, 5028, 5149}
Base 14 up to 50000: 69 First 20: {15, 39, 65, 195, 481, 561, 781, 793, 841, 985, 1105, 1111, 1541, 1891, 2257, 2465, 2561, 2665, 2743, 3277}
Base 15 up to 50000: 42 First 20: {14, 341, 742, 946, 1477, 1541, 1687, 1729, 1891, 1921, 2821, 3133, 3277, 4187, 6541, 6601, 7471, 8701, 8911, 9073}
Base 16 up to 50000: 145 First 20: {15, 51, 85, 91, 255, 341, 435, 451, 561, 595, 645, 703, 1105, 1247, 1261, 1271, 1285, 1387, 1581, 1687}
Base 17 up to 50000: 63 First 20: {4, 8, 9, 16, 45, 91, 145, 261, 781, 1111, 1228, 1305, 1729, 1885, 2149, 2821, 3991, 4005, 4033, 4187}
Base 18 up to 50000: 98 First 20: {25, 49, 65, 85, 133, 221, 323, 325, 343, 425, 451, 637, 931, 1105, 1225, 1369, 1387, 1649, 1729, 1921}
Base 19 up to 50000: 93 First 20: {6, 9, 15, 18, 45, 49, 153, 169, 343, 561, 637, 889, 905, 906, 1035, 1105, 1629, 1661, 1849, 1891}
Base 20 up to 50000: 66 First 20: {21, 57, 133, 231, 399, 561, 671, 861, 889, 1281, 1653, 1729, 1891, 2059, 2413, 2501, 2761, 2821, 2947, 3059}
Nim
import std/[strformat, strutils]
proc powMod*(a, n, m: int): int =
## Return "a^n mod m".
var a = a mod m
var n = n
if a > 0:
result = 1
while n > 0:
if (n and 1) != 0:
result = (result * a) mod m
n = n shr 1
a = (a * a) mod m
func isPrime(n: Natural): bool =
## Return true if "n" is prime.
if n < 2: return false
if (n and 1) == 0: return n == 2
if n mod 3 == 0: return n == 3
var k = 5
var delta = 2
while k * k <= n:
if n mod k == 0: return false
inc k, delta
delta = 6 - delta
result = true
func isFermatPseudoprime(x, a: int): bool =
## Return true is "x" is a Fermat pseudoprime to base "a".
if x.isPrime: return false
result = powMod(a, x - 1, x) == 1
const Lim = 50_000
for a in 1..20:
var count = 0
var first20: seq[int]
for x in 1..Lim:
if x.isFermatPseudoprime(a):
inc count
if count <= 20:
first20.add x
echo &"Base {a}:"
echo &" Number of Fermat pseudoprimes up to {insertSep($Lim)}: {count}"
echo &" First 20: {first20.join(\" \")}"
- Output:
Base 1: Number of Fermat pseudoprimes up to 50_000: 44866 First 20: 4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32 Base 2: Number of Fermat pseudoprimes up to 50_000: 55 First 20: 341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 4681 5461 6601 7957 8321 Base 3: Number of Fermat pseudoprimes up to 50_000: 53 First 20: 91 121 286 671 703 949 1105 1541 1729 1891 2465 2665 2701 2821 3281 3367 3751 4961 5551 6601 Base 4: Number of Fermat pseudoprimes up to 50_000: 111 First 20: 15 85 91 341 435 451 561 645 703 1105 1247 1271 1387 1581 1695 1729 1891 1905 2047 2071 Base 5: Number of Fermat pseudoprimes up to 50_000: 54 First 20: 4 124 217 561 781 1541 1729 1891 2821 4123 5461 5611 5662 5731 6601 7449 7813 8029 8911 9881 Base 6: Number of Fermat pseudoprimes up to 50_000: 74 First 20: 35 185 217 301 481 1105 1111 1261 1333 1729 2465 2701 2821 3421 3565 3589 3913 4123 4495 5713 Base 7: Number of Fermat pseudoprimes up to 50_000: 49 First 20: 6 25 325 561 703 817 1105 1825 2101 2353 2465 3277 4525 4825 6697 8321 10225 10585 10621 11041 Base 8: Number of Fermat pseudoprimes up to 50_000: 150 First 20: 9 21 45 63 65 105 117 133 153 231 273 341 481 511 561 585 645 651 861 949 Base 9: Number of Fermat pseudoprimes up to 50_000: 113 First 20: 4 8 28 52 91 121 205 286 364 511 532 616 671 697 703 946 949 1036 1105 1288 Base 10: Number of Fermat pseudoprimes up to 50_000: 65 First 20: 9 33 91 99 259 451 481 561 657 703 909 1233 1729 2409 2821 2981 3333 3367 4141 4187 Base 11: Number of Fermat pseudoprimes up to 50_000: 61 First 20: 10 15 70 133 190 259 305 481 645 703 793 1105 1330 1729 2047 2257 2465 2821 4577 4921 Base 12: Number of Fermat pseudoprimes up to 50_000: 91 First 20: 65 91 133 143 145 247 377 385 703 1045 1099 1105 1649 1729 1885 1891 2041 2233 2465 2701 Base 13: Number of Fermat pseudoprimes up to 50_000: 68 First 20: 4 6 12 21 85 105 231 244 276 357 427 561 1099 1785 1891 2465 2806 3605 5028 5149 Base 14: Number of Fermat pseudoprimes up to 50_000: 69 First 20: 15 39 65 195 481 561 781 793 841 985 1105 1111 1541 1891 2257 2465 2561 2665 2743 3277 Base 15: Number of Fermat pseudoprimes up to 50_000: 42 First 20: 14 341 742 946 1477 1541 1687 1729 1891 1921 2821 3133 3277 4187 6541 6601 7471 8701 8911 9073 Base 16: Number of Fermat pseudoprimes up to 50_000: 145 First 20: 15 51 85 91 255 341 435 451 561 595 645 703 1105 1247 1261 1271 1285 1387 1581 1687 Base 17: Number of Fermat pseudoprimes up to 50_000: 63 First 20: 4 8 9 16 45 91 145 261 781 1111 1228 1305 1729 1885 2149 2821 3991 4005 4033 4187 Base 18: Number of Fermat pseudoprimes up to 50_000: 98 First 20: 25 49 65 85 133 221 323 325 343 425 451 637 931 1105 1225 1369 1387 1649 1729 1921 Base 19: Number of Fermat pseudoprimes up to 50_000: 93 First 20: 6 9 15 18 45 49 153 169 343 561 637 889 905 906 1035 1105 1629 1661 1849 1891 Base 20: Number of Fermat pseudoprimes up to 50_000: 66 First 20: 21 57 133 231 399 561 671 861 889 1281 1653 1729 1891 2059 2413 2501 2761 2821 2947 3059
Perl
use v5.36;
use ntheory 'is_prime';
sub expmod ($a, $b, $n) {
my $c = 1;
do {
($c *= $a) %= $n if $b % 2;
($a *= $a) %= $n;
} while ($b = int $b/2);
$c
}
my $threshold = 50000;
say "For each base: # of terms up to $threshold, first 20 displayed";
for my $b (1..20) {
my @pseudo = grep { !is_prime($_) && (1 == expmod $b, $_ - 1, $_) } 1..$threshold;
printf "base %2d: %5d (%s)\n", $b, $#pseudo, join ' ', @pseudo[1..20];
}
- Output:
base 1: 44866 (4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32) base 2: 55 (341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 4681 5461 6601 7957 8321) base 3: 53 (91 121 286 671 703 949 1105 1541 1729 1891 2465 2665 2701 2821 3281 3367 3751 4961 5551 6601) base 4: 111 (15 85 91 341 435 451 561 645 703 1105 1247 1271 1387 1581 1695 1729 1891 1905 2047 2071) base 5: 54 (4 124 217 561 781 1541 1729 1891 2821 4123 5461 5611 5662 5731 6601 7449 7813 8029 8911 9881) base 6: 74 (35 185 217 301 481 1105 1111 1261 1333 1729 2465 2701 2821 3421 3565 3589 3913 4123 4495 5713) base 7: 49 (6 25 325 561 703 817 1105 1825 2101 2353 2465 3277 4525 4825 6697 8321 10225 10585 10621 11041) base 8: 150 (9 21 45 63 65 105 117 133 153 231 273 341 481 511 561 585 645 651 861 949) base 9: 113 (4 8 28 52 91 121 205 286 364 511 532 616 671 697 703 946 949 1036 1105 1288) base 10: 65 (9 33 91 99 259 451 481 561 657 703 909 1233 1729 2409 2821 2981 3333 3367 4141 4187) base 11: 61 (10 15 70 133 190 259 305 481 645 703 793 1105 1330 1729 2047 2257 2465 2821 4577 4921) base 12: 91 (65 91 133 143 145 247 377 385 703 1045 1099 1105 1649 1729 1885 1891 2041 2233 2465 2701) base 13: 68 (4 6 12 21 85 105 231 244 276 357 427 561 1099 1785 1891 2465 2806 3605 5028 5149) base 14: 69 (15 39 65 195 481 561 781 793 841 985 1105 1111 1541 1891 2257 2465 2561 2665 2743 3277) base 15: 42 (14 341 742 946 1477 1541 1687 1729 1891 1921 2821 3133 3277 4187 6541 6601 7471 8701 8911 9073) base 16: 145 (15 51 85 91 255 341 435 451 561 595 645 703 1105 1247 1261 1271 1285 1387 1581 1687) base 17: 63 (4 8 9 16 45 91 145 261 781 1111 1228 1305 1729 1885 2149 2821 3991 4005 4033 4187) base 18: 98 (25 49 65 85 133 221 323 325 343 425 451 637 931 1105 1225 1369 1387 1649 1729 1921) base 19: 93 (6 9 15 18 45 49 153 169 343 561 637 889 905 906 1035 1105 1629 1661 1849 1891) base 20: 66 (21 57 133 231 399 561 671 861 889 1281 1653 1729 1891 2059 2413 2501 2761 2821 2947 3059)
Phix
with javascript_semantics include mpfr.e function fermat_pseudoprime(integer x, base) if is_prime(x) then return false end if mpz z = mpz_init(x), a = mpz_init(base) mpz_powm_ui(z, a, x-1, z) return mpz_cmp_si(z,1) == 0 end function constant limits = {12000, 25000, 50000, 100000}, sp = repeat('-',53) printf(1,"Base <%s first 20 %s> <=%s\n",{sp,sp,join(limits,fmt:="%6d")}) for base=1 to 20 do integer count = 0, nlx = 1, nl = limits[1] sequence first20 = {}, counts = repeat(0,length(limits)) for x=2 to limits[$] do if fermat_pseudoprime(x, base) then if count<20 then first20 &= x end if count += 1 end if if x=nl then counts[nlx] = count if nlx<length(limits) then nlx += 1 nl = limits[nlx] end if end if end for string s = join(first20,fmt:="%5d"), c = join(counts,fmt:="%6d") printf(1,"%2d: %s %s\n", {base, s, c}) end for
- Output:
Base <----------------------------------------------------- first 20 -----------------------------------------------------> <= 12000 25000 50000 100000 1: 4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32 10561 22237 44866 90407 2: 341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 4681 5461 6601 7957 8321 25 38 55 78 3: 91 121 286 671 703 949 1105 1541 1729 1891 2465 2665 2701 2821 3281 3367 3751 4961 5551 6601 25 38 53 78 4: 15 85 91 341 435 451 561 645 703 1105 1247 1271 1387 1581 1695 1729 1891 1905 2047 2071 50 75 111 153 5: 4 124 217 561 781 1541 1729 1891 2821 4123 5461 5611 5662 5731 6601 7449 7813 8029 8911 9881 22 35 54 73 6: 35 185 217 301 481 1105 1111 1261 1333 1729 2465 2701 2821 3421 3565 3589 3913 4123 4495 5713 31 46 74 104 7: 6 25 325 561 703 817 1105 1825 2101 2353 2465 3277 4525 4825 6697 8321 10225 10585 10621 11041 21 32 49 73 8: 9 21 45 63 65 105 117 133 153 231 273 341 481 511 561 585 645 651 861 949 76 110 150 218 9: 4 8 28 52 91 121 205 286 364 511 532 616 671 697 703 946 949 1036 1105 1288 55 83 113 164 10: 9 33 91 99 259 451 481 561 657 703 909 1233 1729 2409 2821 2981 3333 3367 4141 4187 35 53 65 90 11: 10 15 70 133 190 259 305 481 645 703 793 1105 1330 1729 2047 2257 2465 2821 4577 4921 30 41 61 89 12: 65 91 133 143 145 247 377 385 703 1045 1099 1105 1649 1729 1885 1891 2041 2233 2465 2701 35 60 91 127 13: 4 6 12 21 85 105 231 244 276 357 427 561 1099 1785 1891 2465 2806 3605 5028 5149 31 51 68 91 14: 15 39 65 195 481 561 781 793 841 985 1105 1111 1541 1891 2257 2465 2561 2665 2743 3277 33 51 69 96 15: 14 341 742 946 1477 1541 1687 1729 1891 1921 2821 3133 3277 4187 6541 6601 7471 8701 8911 9073 22 31 42 70 16: 15 51 85 91 255 341 435 451 561 595 645 703 1105 1247 1261 1271 1285 1387 1581 1687 69 102 145 200 17: 4 8 9 16 45 91 145 261 781 1111 1228 1305 1729 1885 2149 2821 3991 4005 4033 4187 31 44 63 83 18: 25 49 65 85 133 221 323 325 343 425 451 637 931 1105 1225 1369 1387 1649 1729 1921 46 69 98 134 19: 6 9 15 18 45 49 153 169 343 561 637 889 905 906 1035 1105 1629 1661 1849 1891 48 70 93 121 20: 21 57 133 231 399 561 671 861 889 1281 1653 1729 1891 2059 2413 2501 2761 2821 2947 3059 35 49 66 101
Python
Translated from the Julia version.
from sympy import isprime
def ispseudo(n, base):
return not isprime(n) and pow(base, n - 1, n) == 1
for b in range(1, 21):
pseudos = [n for n in range(1, 50001) if ispseudo(n, b)]
print(f"Base {str(b).rjust(2)} up to 50000: {str(len(pseudos)).rjust(5)} First 20: {pseudos[:20]}")
Output:
Base 1 up to 50000: 44866 First 20: [4, 6, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32] Base 2 up to 50000: 55 First 20: [341, 561, 645, 1105, 1387, 1729, 1905, 2047, 2465, 2701, 2821, 3277, 4033, 4369, 4371, 4681, 5461, 6601, 7957, 8321] Base 3 up to 50000: 53 First 20: [91, 121, 286, 671, 703, 949, 1105, 1541, 1729, 1891, 2465, 2665, 2701, 2821, 3281, 3367, 3751, 4961, 5551, 6601] Base 4 up to 50000: 111 First 20: [15, 85, 91, 341, 435, 451, 561, 645, 703, 1105, 1247, 1271, 1387, 1581, 1695, 1729, 1891, 1905, 2047, 2071] Base 5 up to 50000: 54 First 20: [4, 124, 217, 561, 781, 1541, 1729, 1891, 2821, 4123, 5461, 5611, 5662, 5731, 6601, 7449, 7813, 8029, 8911, 9881] Base 6 up to 50000: 74 First 20: [35, 185, 217, 301, 481, 1105, 1111, 1261, 1333, 1729, 2465, 2701, 2821, 3421, 3565, 3589, 3913, 4123, 4495, 5713] Base 7 up to 50000: 49 First 20: [6, 25, 325, 561, 703, 817, 1105, 1825, 2101, 2353, 2465, 3277, 4525, 4825, 6697, 8321, 10225, 10585, 10621, 11041] Base 8 up to 50000: 150 First 20: [9, 21, 45, 63, 65, 105, 117, 133, 153, 231, 273, 341, 481, 511, 561, 585, 645, 651, 861, 949] Base 9 up to 50000: 113 First 20: [4, 8, 28, 52, 91, 121, 205, 286, 364, 511, 532, 616, 671, 697, 703, 946, 949, 1036, 1105, 1288] Base 10 up to 50000: 65 First 20: [9, 33, 91, 99, 259, 451, 481, 561, 657, 703, 909, 1233, 1729, 2409, 2821, 2981, 3333, 3367, 4141, 4187] Base 11 up to 50000: 61 First 20: [10, 15, 70, 133, 190, 259, 305, 481, 645, 703, 793, 1105, 1330, 1729, 2047, 2257, 2465, 2821, 4577, 4921] Base 12 up to 50000: 91 First 20: [65, 91, 133, 143, 145, 247, 377, 385, 703, 1045, 1099, 1105, 1649, 1729, 1885, 1891, 2041, 2233, 2465, 2701] Base 13 up to 50000: 68 First 20: [4, 6, 12, 21, 85, 105, 231, 244, 276, 357, 427, 561, 1099, 1785, 1891, 2465, 2806, 3605, 5028, 5149] Base 14 up to 50000: 69 First 20: [15, 39, 65, 195, 481, 561, 781, 793, 841, 985, 1105, 1111, 1541, 1891, 2257, 2465, 2561, 2665, 2743, 3277] Base 15 up to 50000: 42 First 20: [14, 341, 742, 946, 1477, 1541, 1687, 1729, 1891, 1921, 2821, 3133, 3277, 4187, 6541, 6601, 7471, 8701, 8911, 9073] Base 16 up to 50000: 145 First 20: [15, 51, 85, 91, 255, 341, 435, 451, 561, 595, 645, 703, 1105, 1247, 1261, 1271, 1285, 1387, 1581, 1687] Base 17 up to 50000: 63 First 20: [4, 8, 9, 16, 45, 91, 145, 261, 781, 1111, 1228, 1305, 1729, 1885, 2149, 2821, 3991, 4005, 4033, 4187] Base 18 up to 50000: 98 First 20: [25, 49, 65, 85, 133, 221, 323, 325, 343, 425, 451, 637, 931, 1105, 1225, 1369, 1387, 1649, 1729, 1921] Base 19 up to 50000: 93 First 20: [6, 9, 15, 18, 45, 49, 153, 169, 343, 561, 637, 889, 905, 906, 1035, 1105, 1629, 1661, 1849, 1891] Base 20 up to 50000: 66 First 20: [21, 57, 133, 231, 399, 561, 671, 861, 889, 1281, 1653, 1729, 1891, 2059, 2413, 2501, 2761, 2821, 2947, 3059]
Quackery
eratosthenes and isprime are defined at Sieve of Eratosthenes#Quackery.
**mod is defined at Modular exponentiation#Quackery.
Count threshold extended to 50000.
50000 eratosthenes
[ over 2 < iff
[ 2drop false ]
done
over isprime iff
[ 2drop false ]
done
over 1 - rot **mod
1 = ] is fpp ( n n --> b )
20 times
[ say "Base: "i^ 1+ dup echo cr
temp put
[]
50000 times
[ i^ 1+ temp share
fpp if [ i^ 1+ join ] ]
temp release
dup size echo
say " Fermat Pseudoprimes up to 50000"
cr
20 split drop
witheach
[ echo
i^ 19 < if sp ]
say "..."
cr cr ]- Output:
Base: 1 44866 Fermat Pseudoprimes up to 50000 4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32... Base: 2 55 Fermat Pseudoprimes up to 50000 341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 4681 5461 6601 7957 8321... Base: 3 53 Fermat Pseudoprimes up to 50000 91 121 286 671 703 949 1105 1541 1729 1891 2465 2665 2701 2821 3281 3367 3751 4961 5551 6601... Base: 4 111 Fermat Pseudoprimes up to 50000 15 85 91 341 435 451 561 645 703 1105 1247 1271 1387 1581 1695 1729 1891 1905 2047 2071... Base: 5 54 Fermat Pseudoprimes up to 50000 4 124 217 561 781 1541 1729 1891 2821 4123 5461 5611 5662 5731 6601 7449 7813 8029 8911 9881... Base: 6 74 Fermat Pseudoprimes up to 50000 35 185 217 301 481 1105 1111 1261 1333 1729 2465 2701 2821 3421 3565 3589 3913 4123 4495 5713... Base: 7 49 Fermat Pseudoprimes up to 50000 6 25 325 561 703 817 1105 1825 2101 2353 2465 3277 4525 4825 6697 8321 10225 10585 10621 11041... Base: 8 150 Fermat Pseudoprimes up to 50000 9 21 45 63 65 105 117 133 153 231 273 341 481 511 561 585 645 651 861 949... Base: 9 113 Fermat Pseudoprimes up to 50000 4 8 28 52 91 121 205 286 364 511 532 616 671 697 703 946 949 1036 1105 1288... Base: 10 65 Fermat Pseudoprimes up to 50000 9 33 91 99 259 451 481 561 657 703 909 1233 1729 2409 2821 2981 3333 3367 4141 4187... Base: 11 61 Fermat Pseudoprimes up to 50000 10 15 70 133 190 259 305 481 645 703 793 1105 1330 1729 2047 2257 2465 2821 4577 4921... Base: 12 91 Fermat Pseudoprimes up to 50000 65 91 133 143 145 247 377 385 703 1045 1099 1105 1649 1729 1885 1891 2041 2233 2465 2701... Base: 13 68 Fermat Pseudoprimes up to 50000 4 6 12 21 85 105 231 244 276 357 427 561 1099 1785 1891 2465 2806 3605 5028 5149... Base: 14 69 Fermat Pseudoprimes up to 50000 15 39 65 195 481 561 781 793 841 985 1105 1111 1541 1891 2257 2465 2561 2665 2743 3277... Base: 15 42 Fermat Pseudoprimes up to 50000 14 341 742 946 1477 1541 1687 1729 1891 1921 2821 3133 3277 4187 6541 6601 7471 8701 8911 9073... Base: 16 145 Fermat Pseudoprimes up to 50000 15 51 85 91 255 341 435 451 561 595 645 703 1105 1247 1261 1271 1285 1387 1581 1687... Base: 17 63 Fermat Pseudoprimes up to 50000 4 8 9 16 45 91 145 261 781 1111 1228 1305 1729 1885 2149 2821 3991 4005 4033 4187... Base: 18 98 Fermat Pseudoprimes up to 50000 25 49 65 85 133 221 323 325 343 425 451 637 931 1105 1225 1369 1387 1649 1729 1921... Base: 19 93 Fermat Pseudoprimes up to 50000 6 9 15 18 45 49 153 169 343 561 637 889 905 906 1035 1105 1629 1661 1849 1891... Base: 20 66 Fermat Pseudoprimes up to 50000 21 57 133 231 399 561 671 861 889 1281 1653 1729 1891 2059 2413 2501 2761 2821 2947 3059...
Raku
use List::Divvy;
for 1..20 -> $base {
my @pseudo = lazy (2..*).hyper.grep: { !.is-prime && (1 == expmod $base, $_ - 1, $_) }
my $threshold = 100_000;
say $base.fmt("Base %2d - Up to $threshold: ") ~ (+@pseudo.&upto: $threshold).fmt('%5d')
~ " First 20: " ~ @pseudo[^20].gist
}
Base 1 - Up to 100000: 90407 First 20: (4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32) Base 2 - Up to 100000: 78 First 20: (341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 4681 5461 6601 7957 8321) Base 3 - Up to 100000: 78 First 20: (91 121 286 671 703 949 1105 1541 1729 1891 2465 2665 2701 2821 3281 3367 3751 4961 5551 6601) Base 4 - Up to 100000: 153 First 20: (15 85 91 341 435 451 561 645 703 1105 1247 1271 1387 1581 1695 1729 1891 1905 2047 2071) Base 5 - Up to 100000: 73 First 20: (4 124 217 561 781 1541 1729 1891 2821 4123 5461 5611 5662 5731 6601 7449 7813 8029 8911 9881) Base 6 - Up to 100000: 104 First 20: (35 185 217 301 481 1105 1111 1261 1333 1729 2465 2701 2821 3421 3565 3589 3913 4123 4495 5713) Base 7 - Up to 100000: 73 First 20: (6 25 325 561 703 817 1105 1825 2101 2353 2465 3277 4525 4825 6697 8321 10225 10585 10621 11041) Base 8 - Up to 100000: 218 First 20: (9 21 45 63 65 105 117 133 153 231 273 341 481 511 561 585 645 651 861 949) Base 9 - Up to 100000: 164 First 20: (4 8 28 52 91 121 205 286 364 511 532 616 671 697 703 946 949 1036 1105 1288) Base 10 - Up to 100000: 90 First 20: (9 33 91 99 259 451 481 561 657 703 909 1233 1729 2409 2821 2981 3333 3367 4141 4187) Base 11 - Up to 100000: 89 First 20: (10 15 70 133 190 259 305 481 645 703 793 1105 1330 1729 2047 2257 2465 2821 4577 4921) Base 12 - Up to 100000: 127 First 20: (65 91 133 143 145 247 377 385 703 1045 1099 1105 1649 1729 1885 1891 2041 2233 2465 2701) Base 13 - Up to 100000: 91 First 20: (4 6 12 21 85 105 231 244 276 357 427 561 1099 1785 1891 2465 2806 3605 5028 5149) Base 14 - Up to 100000: 96 First 20: (15 39 65 195 481 561 781 793 841 985 1105 1111 1541 1891 2257 2465 2561 2665 2743 3277) Base 15 - Up to 100000: 70 First 20: (14 341 742 946 1477 1541 1687 1729 1891 1921 2821 3133 3277 4187 6541 6601 7471 8701 8911 9073) Base 16 - Up to 100000: 200 First 20: (15 51 85 91 255 341 435 451 561 595 645 703 1105 1247 1261 1271 1285 1387 1581 1687) Base 17 - Up to 100000: 83 First 20: (4 8 9 16 45 91 145 261 781 1111 1228 1305 1729 1885 2149 2821 3991 4005 4033 4187) Base 18 - Up to 100000: 134 First 20: (25 49 65 85 133 221 323 325 343 425 451 637 931 1105 1225 1369 1387 1649 1729 1921) Base 19 - Up to 100000: 121 First 20: (6 9 15 18 45 49 153 169 343 561 637 889 905 906 1035 1105 1629 1661 1849 1891) Base 20 - Up to 100000: 101 First 20: (21 57 133 231 399 561 671 861 889 1281 1653 1729 1891 2059 2413 2501 2761 2821 2947 3059)
Sidef
func generate_fermat_psp(base, upto, n) {
if (base == 1) {
return(n.by { .is_composite }, upto.composite_count)
}
var psp = []
for k in (1..Inf) {
break if (k.pn_primorial > upto)
psp << k.fermat_psp(base, 1, upto)...
}
return(psp.sort.first(n), psp.len)
}
var upto = 1e7
for base in (1..20) {
var (psp, count) = generate_fermat_psp(base, upto, 20)
printf("Base %2d - up to #{upto}: %7d First 20: %s\n", base, count, psp)
}
- Output:
Base 1 - up to 10000000: 9335420 First 20: [4, 6, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32] Base 2 - up to 10000000: 750 First 20: [341, 561, 645, 1105, 1387, 1729, 1905, 2047, 2465, 2701, 2821, 3277, 4033, 4369, 4371, 4681, 5461, 6601, 7957, 8321] Base 3 - up to 10000000: 760 First 20: [91, 121, 286, 671, 703, 949, 1105, 1541, 1729, 1891, 2465, 2665, 2701, 2821, 3281, 3367, 3751, 4961, 5551, 6601] Base 4 - up to 10000000: 1347 First 20: [15, 85, 91, 341, 435, 451, 561, 645, 703, 1105, 1247, 1271, 1387, 1581, 1695, 1729, 1891, 1905, 2047, 2071] Base 5 - up to 10000000: 745 First 20: [4, 124, 217, 561, 781, 1541, 1729, 1891, 2821, 4123, 5461, 5611, 5662, 5731, 6601, 7449, 7813, 8029, 8911, 9881] Base 6 - up to 10000000: 895 First 20: [35, 185, 217, 301, 481, 1105, 1111, 1261, 1333, 1729, 2465, 2701, 2821, 3421, 3565, 3589, 3913, 4123, 4495, 5713] Base 7 - up to 10000000: 659 First 20: [6, 25, 325, 561, 703, 817, 1105, 1825, 2101, 2353, 2465, 3277, 4525, 4825, 6697, 8321, 10225, 10585, 10621, 11041] Base 8 - up to 10000000: 1993 First 20: [9, 21, 45, 63, 65, 105, 117, 133, 153, 231, 273, 341, 481, 511, 561, 585, 645, 651, 861, 949] Base 9 - up to 10000000: 1418 First 20: [4, 8, 28, 52, 91, 121, 205, 286, 364, 511, 532, 616, 671, 697, 703, 946, 949, 1036, 1105, 1288] Base 10 - up to 10000000: 766 First 20: [9, 33, 91, 99, 259, 451, 481, 561, 657, 703, 909, 1233, 1729, 2409, 2821, 2981, 3333, 3367, 4141, 4187] Base 11 - up to 10000000: 695 First 20: [10, 15, 70, 133, 190, 259, 305, 481, 645, 703, 793, 1105, 1330, 1729, 2047, 2257, 2465, 2821, 4577, 4921] Base 12 - up to 10000000: 1091 First 20: [65, 91, 133, 143, 145, 247, 377, 385, 703, 1045, 1099, 1105, 1649, 1729, 1885, 1891, 2041, 2233, 2465, 2701] Base 13 - up to 10000000: 750 First 20: [4, 6, 12, 21, 85, 105, 231, 244, 276, 357, 427, 561, 1099, 1785, 1891, 2465, 2806, 3605, 5028, 5149] Base 14 - up to 10000000: 817 First 20: [15, 39, 65, 195, 481, 561, 781, 793, 841, 985, 1105, 1111, 1541, 1891, 2257, 2465, 2561, 2665, 2743, 3277] Base 15 - up to 10000000: 628 First 20: [14, 341, 742, 946, 1477, 1541, 1687, 1729, 1891, 1921, 2821, 3133, 3277, 4187, 6541, 6601, 7471, 8701, 8911, 9073] Base 16 - up to 10000000: 1749 First 20: [15, 51, 85, 91, 255, 341, 435, 451, 561, 595, 645, 703, 1105, 1247, 1261, 1271, 1285, 1387, 1581, 1687] Base 17 - up to 10000000: 763 First 20: [4, 8, 9, 16, 45, 91, 145, 261, 781, 1111, 1228, 1305, 1729, 1885, 2149, 2821, 3991, 4005, 4033, 4187] Base 18 - up to 10000000: 1161 First 20: [25, 49, 65, 85, 133, 221, 323, 325, 343, 425, 451, 637, 931, 1105, 1225, 1369, 1387, 1649, 1729, 1921] Base 19 - up to 10000000: 932 First 20: [6, 9, 15, 18, 45, 49, 153, 169, 343, 561, 637, 889, 905, 906, 1035, 1105, 1629, 1661, 1849, 1891] Base 20 - up to 10000000: 850 First 20: [21, 57, 133, 231, 399, 561, 671, 861, 889, 1281, 1653, 1729, 1891, 2059, 2413, 2501, 2761, 2821, 2947, 3059]
Wren
import "./math" for Int
import "./gmp" for Mpz
import "./fmt" for Fmt
var one = Mpz.one
var isFermatPseudoprime = Fn.new { |x, a|
if (Int.isPrime(x)) return false
var bx = Mpz.from(x)
a = Mpz.from(a)
return a.modPow(x-1, bx) == one
}
System.print("First 20 Fermat pseudoprimes:")
for (a in 1..20) {
var count = 0
var x = 2
var first20 = List.filled(20, 0)
while (count < 20) {
if (isFermatPseudoprime.call(x, a)) {
first20[count] = x
count = count + 1
}
x = x + 1
}
Fmt.print("Base $2d: $5d", a, first20)
}
var limits = [12000, 25000, 50000, 100000]
var heading = Fmt.swrite("\nCount <= $6d", limits)
System.print(heading)
System.print("-" * (heading.count - 1))
for (a in 1..20) {
Fmt.write("Base $2d: ", a)
var x = 2
var count = 0
for (limit in limits) {
while (x <= limit) {
if (isFermatPseudoprime.call(x, a)) count = count + 1
x = x + 1
}
Fmt.write("$6d ", count)
}
System.print()
}
- Output:
First 20 Fermat pseudoprimes: Base 1: 4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32 Base 2: 341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 4681 5461 6601 7957 8321 Base 3: 91 121 286 671 703 949 1105 1541 1729 1891 2465 2665 2701 2821 3281 3367 3751 4961 5551 6601 Base 4: 15 85 91 341 435 451 561 645 703 1105 1247 1271 1387 1581 1695 1729 1891 1905 2047 2071 Base 5: 4 124 217 561 781 1541 1729 1891 2821 4123 5461 5611 5662 5731 6601 7449 7813 8029 8911 9881 Base 6: 35 185 217 301 481 1105 1111 1261 1333 1729 2465 2701 2821 3421 3565 3589 3913 4123 4495 5713 Base 7: 6 25 325 561 703 817 1105 1825 2101 2353 2465 3277 4525 4825 6697 8321 10225 10585 10621 11041 Base 8: 9 21 45 63 65 105 117 133 153 231 273 341 481 511 561 585 645 651 861 949 Base 9: 4 8 28 52 91 121 205 286 364 511 532 616 671 697 703 946 949 1036 1105 1288 Base 10: 9 33 91 99 259 451 481 561 657 703 909 1233 1729 2409 2821 2981 3333 3367 4141 4187 Base 11: 10 15 70 133 190 259 305 481 645 703 793 1105 1330 1729 2047 2257 2465 2821 4577 4921 Base 12: 65 91 133 143 145 247 377 385 703 1045 1099 1105 1649 1729 1885 1891 2041 2233 2465 2701 Base 13: 4 6 12 21 85 105 231 244 276 357 427 561 1099 1785 1891 2465 2806 3605 5028 5149 Base 14: 15 39 65 195 481 561 781 793 841 985 1105 1111 1541 1891 2257 2465 2561 2665 2743 3277 Base 15: 14 341 742 946 1477 1541 1687 1729 1891 1921 2821 3133 3277 4187 6541 6601 7471 8701 8911 9073 Base 16: 15 51 85 91 255 341 435 451 561 595 645 703 1105 1247 1261 1271 1285 1387 1581 1687 Base 17: 4 8 9 16 45 91 145 261 781 1111 1228 1305 1729 1885 2149 2821 3991 4005 4033 4187 Base 18: 25 49 65 85 133 221 323 325 343 425 451 637 931 1105 1225 1369 1387 1649 1729 1921 Base 19: 6 9 15 18 45 49 153 169 343 561 637 889 905 906 1035 1105 1629 1661 1849 1891 Base 20: 21 57 133 231 399 561 671 861 889 1281 1653 1729 1891 2059 2413 2501 2761 2821 2947 3059 Count <= 12000 25000 50000 100000 ------------------------------------ Base 1: 10561 22237 44866 90407 Base 2: 25 38 55 78 Base 3: 25 38 53 78 Base 4: 50 75 111 153 Base 5: 22 35 54 73 Base 6: 31 46 74 104 Base 7: 21 32 49 73 Base 8: 76 110 150 218 Base 9: 55 83 113 164 Base 10: 35 53 65 90 Base 11: 30 41 61 89 Base 12: 35 60 91 127 Base 13: 31 51 68 91 Base 14: 33 51 69 96 Base 15: 22 31 42 70 Base 16: 69 102 145 200 Base 17: 31 44 63 83 Base 18: 46 69 98 134 Base 19: 48 70 93 121 Base 20: 35 49 66 101