Fermat pseudoprimes

From Rosetta Code
Task
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


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++

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#

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

Java

import java.util.ArrayList;
import java.util.List;
import java.util.stream.Collectors;

public final class FermatPseudoprimes {

	public static void main(String[] args) {
		final int limit = 50_000;

		for ( int base = 1; base <= 20; base++ ) {
			int count = 0;
			List<Integer> firstTwenty = new ArrayList<Integer>();
			for ( int number = 4; number <= limit; number++ ) {
				if ( isFermatPseudoprime(base, number) ) {
					count += 1;
					if ( count <= 20 ) {
						firstTwenty.addLast(number);
					}
				}
			}
			
			System.out.println("Base " + base + ":");
			System.out.println("    The number of Fermat pseudoprimes up to " + limit + ": " + count);
			System.out.println("    The first 20: " +
				firstTwenty.stream().map(String::valueOf).collect(Collectors.joining(" ")));
		}
	}	

	private static boolean isFermatPseudoprime(int base, int number) {
	    return ! isPrime(number) && moduloPower(base, number - 1, number) == 1;
	}
	
	private static long moduloPower(long base, long exponent, long modulus) {
	    if ( modulus == 1 ) {
	        return 0;
	    }
	    
	    base %= modulus;
	    long result = 1;	    
	    while ( exponent > 0 ) {
	        if ( ( exponent & 1 ) == 1 ) {
	            result = ( result * base ) % modulus;
	        }
	        base = ( base * base ) % modulus;
	        exponent >>= 1;
	    }
	    
	    return result;
	}

	private static boolean isPrime(int number) {
	    if ( number < 2 ) {
	        return false;
	    }
	    if ( number % 2 == 0 ) {
	        return number == 2;
	    }
	    if ( number % 3 == 0 ) {
	        return number == 3;
	    }
	    
	    for ( int p = 5; p * p <= number; p += 4 ) {
	        if ( number % p == 0 ) {
	            return false;
	        }
	        p += 2;
	        if ( number % p == 0 ) {
	            return false;
	        }
	    }

	    return true;
	}
	
}
Output:
Base 1:
    The number of Fermat pseudoprimes up to 50000: 44866
    The first 20: 4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32
Base 2:
    The number of Fermat pseudoprimes up to 50000: 55
    The first 20: 341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 4681 5461 6601 7957 8321
Base 3:
    The number of Fermat pseudoprimes up to 50000: 53
    The first 20: 91 121 286 671 703 949 1105 1541 1729 1891 2465 2665 2701 2821 3281 3367 3751 4961 5551 6601
Base 4:
    The number of Fermat pseudoprimes up to 50000: 111
    The first 20: 15 85 91 341 435 451 561 645 703 1105 1247 1271 1387 1581 1695 1729 1891 1905 2047 2071
Base 5:
    The number of Fermat pseudoprimes up to 50000: 54
    The first 20: 4 124 217 561 781 1541 1729 1891 2821 4123 5461 5611 5662 5731 6601 7449 7813 8029 8911 9881
Base 6:
    The number of Fermat pseudoprimes up to 50000: 74
    The first 20: 35 185 217 301 481 1105 1111 1261 1333 1729 2465 2701 2821 3421 3565 3589 3913 4123 4495 5713
Base 7:
    The number of Fermat pseudoprimes up to 50000: 49
    The first 20: 6 25 325 561 703 817 1105 1825 2101 2353 2465 3277 4525 4825 6697 8321 10225 10585 10621 11041
Base 8:
    The number of Fermat pseudoprimes up to 50000: 150
    The first 20: 9 21 45 63 65 105 117 133 153 231 273 341 481 511 561 585 645 651 861 949
Base 9:
    The number of Fermat pseudoprimes up to 50000: 113
    The first 20: 4 8 28 52 91 121 205 286 364 511 532 616 671 697 703 946 949 1036 1105 1288
Base 10:
    The number of Fermat pseudoprimes up to 50000: 65
    The first 20: 9 33 91 99 259 451 481 561 657 703 909 1233 1729 2409 2821 2981 3333 3367 4141 4187
Base 11:
    The number of Fermat pseudoprimes up to 50000: 61
    The first 20: 10 15 70 133 190 259 305 481 645 703 793 1105 1330 1729 2047 2257 2465 2821 4577 4921
Base 12:
    The number of Fermat pseudoprimes up to 50000: 91
    The first 20: 65 91 133 143 145 247 377 385 703 1045 1099 1105 1649 1729 1885 1891 2041 2233 2465 2701
Base 13:
    The number of Fermat pseudoprimes up to 50000: 68
    The first 20: 4 6 12 21 85 105 231 244 276 357 427 561 1099 1785 1891 2465 2806 3605 5028 5149
Base 14:
    The number of Fermat pseudoprimes up to 50000: 69
    The first 20: 15 39 65 195 481 561 781 793 841 985 1105 1111 1541 1891 2257 2465 2561 2665 2743 3277
Base 15:
    The number of Fermat pseudoprimes up to 50000: 42
    The first 20: 14 341 742 946 1477 1541 1687 1729 1891 1921 2821 3133 3277 4187 6541 6601 7471 8701 8911 9073
Base 16:
    The number of Fermat pseudoprimes up to 50000: 145
    The first 20: 15 51 85 91 255 341 435 451 561 595 645 703 1105 1247 1261 1271 1285 1387 1581 1687
Base 17:
    The number of Fermat pseudoprimes up to 50000: 63
    The first 20: 4 8 9 16 45 91 145 261 781 1111 1228 1305 1729 1885 2149 2821 3991 4005 4033 4187
Base 18:
    The number of Fermat pseudoprimes up to 50000: 98
    The first 20: 25 49 65 85 133 221 323 325 343 425 451 637 931 1105 1225 1369 1387 1649 1729 1921
Base 19:
    The number of Fermat pseudoprimes up to 50000: 93
    The first 20: 6 9 15 18 45 49 153 169 343 561 637 889 905 906 1035 1105 1629 1661 1849 1891
Base 20:
    The number of Fermat pseudoprimes up to 50000: 66
    The first 20: 21 57 133 231 399 561 671 861 889 1281 1653 1729 1891 2059 2413 2501 2761 2821 2947 3059

jq

Adapted from Wren

Works with gojq, the Go implementation of jq

gojq supports infinite-precision integer arithmetic and is thus suitable for the given tasks, including the stretch task: the illustration below includes thresholds up to 100,000.

### Generic functions

def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .;

def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in); 

def is_prime:
  . as $n
  | if ($n < 2)         then false
    elif ($n % 2 == 0)  then $n == 2
    elif ($n % 3 == 0)  then $n == 3
    elif ($n % 5 == 0)  then $n == 5
    elif ($n % 7 == 0)  then $n == 7
    elif ($n % 11 == 0) then $n == 11
    elif ($n % 13 == 0) then $n == 13
    elif ($n % 17 == 0) then $n == 17
    elif ($n % 19 == 0) then $n == 19
    else sqrt as $s
    | 23
    | until( . > $s or ($n % . == 0); . + 2)
    | . > $s
    end;

# The input to the power $exp modulo $mod, which is assumed to be positive
def modPow($exp; $mod):
  { base: (. % $mod),
    exp: $exp,
    r: 1
  }
  | until( .exp == 0;
      if .base == 0 then .r = 0 | .exp = 0
      else if .exp % 2 == 1 then .r = (.r * .base) % $mod end
      | .exp |= ((. / 2) | trunc)
      | .base |= (. * .) % $mod
      end )
  | .r;

### Fermat Pseudoprimes

def isFermatPseudoprime($a):
  if is_prime then false
  else . as $x
  | ($a | modPow($x-1; $x)) == 1
  end;


def display($n):
  "First \($n) Fermat pseudoprimes:",
  (range(1; $n+1) as $a
   | { x: 2,  count: 0 }
   | until (.count >= 20;
       if .x | isFermatPseudoprime($a)
       then .first20 += [.x]
       | .count += 1
       end
       | .x += 1 )
   | "Base \($a|lpad(2)): \(.first20|map(lpad(5))|join(" "))" ) ;

def task($limits):
  def heading: "\nCount <= \($limits | map(lpad(6)) | join(" "))";

  heading as $heading
  | $heading,
  "-" * ($heading|length - 1),
  (range(1; 21) as $a
   | "Base \($a|lpad(2)): " +
     ([foreach $limits[] as $limit ( {x: 2,  count: 0 };
         until (.x > $limit;
           if .x | isFermatPseudoprime($a) then .count += 1  end
           | .x += 1 ) )
       | .count
       | lpad(6) ]
     | join(" ") ) ) ;

display(20),
task( [12000, 25000, 50000, 100000] )
Output:

As for Wren.

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

PascalABC.NET

Translation of: Nim
function IsPrime(n: int64): boolean;
begin
  if (n = 2) or (n = 3) then Result := true
  else if (n <= 1) or ((n mod 2) = 0) or ((n mod 3) = 0) then Result := false
  else
  begin
    var i := 5;
    Result := False;
    while i <= trunc(sqrt(n)) do
    begin
      if ((n mod i) = 0) or ((n mod (i + 2)) = 0) then exit;
      i += 6;
    end;
    Result := True;
  end;
end;

function powMod(a, n, m: int64): int64;
begin
  a := a mod m;
  if a > 0 then
  begin
    result := 1;
    while n > 0 do
    begin
      if (n and 1) <> 0 then
        result := (result * a) mod m;
      n := n shr 1;
      a := (a * a) mod m;
    end;
  end;
end;

function isFermatPseudoprime(x, a: integer): boolean;
begin
  if isPrime(x) then result := false
  else result := powMod(a, x - 1, x) = 1
end;

begin
  for var a := 1 to 20 do
  begin
    var count := 0;
    var first20 := new List<integer>;
    for var x := 1 to 50_000 do
      if isFermatPseudoprime(x, a) then
      begin
        count += 1;
        if count <= 20 then first20.add(x)
      end;
    println('Base', a, ':');
    println('  Number of Fermat pseudoprimes up to 50_000:', count);
    println('  First 20:', first20);
  end;
end.
Output:
See output of Nim

Perl

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

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)

SETL

program fermat_pseudoprimes;
    init count := 50000;
    init prime := sieve(count);

    loop for a in [1..20] do
        p := fermat_pseudoprimes(a, count);
        print("a=" + lpad(str a,2) + ": #" + lpad(str #p, 5)
                 + "; first: " + str p(1..20));
    end loop;

    proc fermat_pseudoprimes(a, maxval);
        return [x in [2..maxval] | is_fermat_pseudoprime(x, a)];
    end proc;
    
    proc is_fermat_pseudoprime(x, a);
        return not prime(x) and (a**(x-1)-1) mod x = 0;
    end proc;

    proc sieve(n);
        pr := [True] * n;
        pr(1) := False;
        loop for p in [2..floor sqrt n] do
            loop for c in [p*p, p*p+p..n] do
                pr(c) := False;
            end loop;
        end loop;
        return pr;
    end proc;
end program;
Output:
a= 1: #44866; first: [4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32]
a= 2: #   55; first: [341 561 645 1105 1387 1729 1905 2047 2465 2701 2821 3277 4033 4369 4371 4681 5461 6601 7957 8321]
a= 3: #   53; first: [91 121 286 671 703 949 1105 1541 1729 1891 2465 2665 2701 2821 3281 3367 3751 4961 5551 6601]
a= 4: #  111; first: [15 85 91 341 435 451 561 645 703 1105 1247 1271 1387 1581 1695 1729 1891 1905 2047 2071]
a= 5: #   54; first: [4 124 217 561 781 1541 1729 1891 2821 4123 5461 5611 5662 5731 6601 7449 7813 8029 8911 9881]
a= 6: #   74; first: [35 185 217 301 481 1105 1111 1261 1333 1729 2465 2701 2821 3421 3565 3589 3913 4123 4495 5713]
a= 7: #   49; first: [6 25 325 561 703 817 1105 1825 2101 2353 2465 3277 4525 4825 6697 8321 10225 10585 10621 11041]
a= 8: #  150; first: [9 21 45 63 65 105 117 133 153 231 273 341 481 511 561 585 645 651 861 949]
a= 9: #  113; first: [4 8 28 52 91 121 205 286 364 511 532 616 671 697 703 946 949 1036 1105 1288]
a=10: #   65; first: [9 33 91 99 259 451 481 561 657 703 909 1233 1729 2409 2821 2981 3333 3367 4141 4187]
a=11: #   61; first: [10 15 70 133 190 259 305 481 645 703 793 1105 1330 1729 2047 2257 2465 2821 4577 4921]
a=12: #   91; first: [65 91 133 143 145 247 377 385 703 1045 1099 1105 1649 1729 1885 1891 2041 2233 2465 2701]
a=13: #   68; first: [4 6 12 21 85 105 231 244 276 357 427 561 1099 1785 1891 2465 2806 3605 5028 5149]
a=14: #   69; first: [15 39 65 195 481 561 781 793 841 985 1105 1111 1541 1891 2257 2465 2561 2665 2743 3277]
a=15: #   42; first: [14 341 742 946 1477 1541 1687 1729 1891 1921 2821 3133 3277 4187 6541 6601 7471 8701 8911 9073]
a=16: #  145; first: [15 51 85 91 255 341 435 451 561 595 645 703 1105 1247 1261 1271 1285 1387 1581 1687]
a=17: #   63; first: [4 8 9 16 45 91 145 261 781 1111 1228 1305 1729 1885 2149 2821 3991 4005 4033 4187]
a=18: #   98; first: [25 49 65 85 133 221 323 325 343 425 451 637 931 1105 1225 1369 1387 1649 1729 1921]
a=19: #   93; first: [6 9 15 18 45 49 153 169 343 561 637 889 905 906 1035 1105 1629 1661 1849 1891]
a=20: #   66; first: [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

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