inverse_znorder.sf

#!/usr/bin/ruby

# Find all the integers d such that znorder(b, d) = n, for given integers (n,b).

# See also:
#   https://oeis.org/A226477
#   https://oeis.org/A345319

func inverse_znorder(n, b=10) {
    divisors(b**n - 1).grep {|d| znorder(b, d) == n }
}

# Example:
#   find Fermat pseudoprimes to base 2

gather {
    var base = 2
    for n in (1..50) {
        inverse_znorder(n, base).each {|k|
            take(k) if (k.is_psp(base) && k.is_composite)
        }
    }
}.sort.uniq.say

__END__
[341, 561, 645, 1105, 1387, 1729, 1905, 2047, 2701, 3277, 4033, 4369, 4681, 5461, 7957, 10261, 13981, 14491, 15709, 15841, 16705, 18705, 23001, 23377, 35333, 42799, 49981, 55245, 60787, 63973, 65281, 72885, 107185, 126217, 150851, 158369, 162193, 188057, 188461, 215265, 233017, 256999, 266305, 294409, 423793, 452051, 486737, 617093, 647089, 688213, 743665, 838861, 1016801, 1052929, 1082401, 1398101, 1549411, 1620385, 1826203, 1840357, 2113665, 2304167, 2434651, 2528921, 2940337, 2953711, 4025905, 4188889, 4670029, 6236473, 6242685, 6973057, 7295851, 8388607, 8727391, 9588151, 9995671, 10610063, 11242465, 13757653, 14709241, 15732721, 16843009, 17895697, 22369621, 25831585, 29593159, 33193117, 34603041, 72135505, 74658629, 78526729, 102134113, 109118791, 131421541, 179223265, 231927781, 271682651, 285212689, 369269551, 390512785, 509033161, 536870911, 555046097, 611097401, 708621217, 905040953, 994133479, 1090519105, 1227133513, 1418724681, 1831258601, 2709434065, 2954187413, 4384806451, 5726623061, 9972894583, 13177667449, 17384656705, 20408568497, 31184907679, 59862819377, 68736258049, 68853957121, 78536544841, 91625968981, 137438953471, 233009086681, 262815104305, 499069107643, 1100586419201, 1172812402961, 2199023255551, 8796093022207, 23456248059221, 140737488355327]