Fermat pseudoprimes

From Rosetta Code
Fermat pseudoprimes is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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


Related


ALGOL 68[edit]

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++[edit]

Translation of: Wren
#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#[edit]

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

Julia[edit]

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]

Perl[edit]

Library: ntheory
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[edit]

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

Raku[edit]

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)

Wren[edit]

Library: Wren-math
Library: Wren-gmp
Library: Wren-fmt
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