Zsigmondy numbers

Zsigmondy numbers n to a, b, are the greatest divisor of aⁿ - bⁿ that is coprime to a^m - b^m for all positive integers m < n.

E.G.

Suppose we set a = 2 and b = 1. (Zs(n,2,1))

For each n, find the divisors of aⁿ - bⁿ and return the largest that is coprime to all of a^m - b^m, where m is each of the positive integers 1 to n - 1.

When n = 4, 2⁴ - 1⁴ = 15. The divisors of 15 are 1, 3, 5, and 15.

For m = 1, 2, 3 we get 2-1, 2²-1², 2³-1³, or 1, 3, 7.

The divisors of 15 that are coprime to each are 5 and 1, (1 is always included).

The largest coprime divisor is 5, so Zs(4,2,1) = 5.

When n = 6, 2⁶ - 1⁶ = 63; its divisors are 1, 3, 7, 9, 21, 63. The largest divisor coprime to all of 1, 3, 7, 15, 31 is 1, (1 is always included), so Zs(6,2,1) = 1.

If a particular aⁿ - bⁿ is prime, then Zs(n,a,b) will be equal to that prime. 2⁵ - 1⁵ = 31 so Zs(5,2,1) = 31.

Task

Write a general routine (function, procedure, whatever) to find the Zsigmondy number sequence given a set of radices.
Use that routine to generate the first several elements, (at least 10), for the following radix sets.

(2,1)
(3,1)
(4,1)
(5,1)
(6,1)
(7,1)
(3,2)
(5,3)
(7,3)
(7,5)

See also

OEIS:A064078 - Zsigmondy numbers for a = 2, b = 1 OEIS:A064077 - Zsigmondy numbers for a = 3, b = 1 OEIS:A064080 - Zsigmondy numbers for a = 4, b = 1 OEIS:A064081 - Zsigmondy numbers for a = 5, b = 1 OEIS:A064082 - Zsigmondy numbers for a = 6, b = 1 OEIS:A064083 - Zsigmondy numbers for a = 7, b = 1 OEIS:A109325 - Zsigmondy numbers for a = 3, b = 2 OEIS:A109347 - Zsigmondy numbers for a = 5, b = 3 OEIS:A109348 - Zsigmondy numbers for a = 7, b = 3 OEIS:A109349 - Zsigmondy numbers for a = 7, b = 5

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <vector>

std::vector<uint64_t> divisors(uint64_t n) {
    std::vector<uint64_t> result{1};
    uint64_t power = 2;
    for (; (n & 1) == 0; power <<= 1, n >>= 1)
        result.push_back(power);
    for (uint64_t p = 3; p * p <= n; p += 2) {
        size_t size = result.size();
        for (power = p; n % p == 0; power *= p, n /= p)
            for (size_t i = 0; i != size; ++i)
                result.push_back(power * result[i]);
    }
    if (n > 1) {
        size_t size = result.size();
        for (size_t i = 0; i != size; ++i)
            result.push_back(n * result[i]);
    }
    sort(result.begin(), result.end());
    return result;
}

uint64_t ipow(uint64_t base, uint64_t exp) {
    if (exp == 0)
        return 1;
    if ((exp & 1) == 0)
        return ipow(base * base, exp >> 1);
    return base * ipow(base * base, (exp - 1) >> 1);
}

uint64_t zsigmondy(uint64_t n, uint64_t a, uint64_t b) {
    auto p = ipow(a, n) - ipow(b, n);
    auto d = divisors(p);
    if (d.size() == 2)
        return p;
    std::vector<uint64_t> dms(n - 1);
    for (uint64_t m = 1; m < n; ++m)
        dms[m - 1] = ipow(a, m) - ipow(b, m);
    for (auto i = d.rbegin(); i != d.rend(); ++i) {
        uint64_t z = *i;
        if (all_of(dms.begin(), dms.end(),
                   [z](uint64_t x) { return std::gcd(x, z) == 1; }))
            return z;
    }
    return 1;
}

void test(uint64_t a, uint64_t b) {
    std::cout << "Zsigmondy(n, " << a << ", " << b << "):\n";
    for (uint64_t n = 1; n <= 20; ++n) {
        std::cout << zsigmondy(n, a, b) << ' ';
    }
    std::cout << "\n\n";
}

int main() {
    test(2, 1);
    test(3, 1);
    test(4, 1);
    test(5, 1);
    test(6, 1);
    test(7, 1);
    test(3, 2);
    test(5, 3);
    test(7, 3);
    test(7, 5);
}

Output:

Zsigmondy(n, 2, 1):
1 3 7 5 31 1 127 17 73 11 2047 13 8191 43 151 257 131071 19 524287 41 

Zsigmondy(n, 3, 1):
2 1 13 5 121 7 1093 41 757 61 88573 73 797161 547 4561 3281 64570081 703 581130733 1181 

Zsigmondy(n, 4, 1):
3 5 7 17 341 13 5461 257 1387 41 1398101 241 22369621 3277 49981 65537 5726623061 4033 91625968981 61681 

Zsigmondy(n, 5, 1):
4 3 31 13 781 7 19531 313 15751 521 12207031 601 305175781 13021 315121 195313 190734863281 5167 4768371582031 375601 

Zsigmondy(n, 6, 1):
5 7 43 37 311 31 55987 1297 46873 1111 72559411 1261 2612138803 5713 1406371 1679617 3385331888947 46441 121871948002099 1634221 

Zsigmondy(n, 7, 1):
6 1 19 25 2801 43 137257 1201 39331 2101 329554457 2353 16148168401 102943 4956001 2882401 38771752331201 117307 1899815864228857 1129901 

Zsigmondy(n, 3, 2):
1 5 19 13 211 7 2059 97 1009 11 175099 61 1586131 463 3571 6817 129009091 577 1161737179 4621 

Zsigmondy(n, 5, 3):
2 1 49 17 1441 19 37969 353 19729 421 24325489 481 609554401 10039 216001 198593 381405156481 12979 9536162033329 288961 

Zsigmondy(n, 7, 3):
4 5 79 29 4141 37 205339 1241 127639 341 494287399 2041 24221854021 82573 3628081 2885681 58157596211761 109117 2849723505777919 4871281 

Zsigmondy(n, 7, 5):
2 3 109 37 6841 13 372709 1513 176149 1661 964249309 1801 47834153641 75139 3162961 3077713 115933787267041 30133 5689910849522509 3949201

J

Implementation:

dp=: -/@:(^/) NB. (a,b) dp n  is (a^n)-(b^n)
Zsigmondy=: dp (-.,)&.:q: (dp 1 }. i.)

In other words, 2 1 dp 1 2 3 4 is 1 3 7 15. And coprime is sequence difference (not set difference, since prime factors may repeat) under factorization (generate the sequence of prime factors for each number, remove any common factors and form the product of any that remain).

Task examples:

   Zs=: Zsigmondy"1 0&(1x+i.20) NB. first 20 in sequence
   Zs 2 1
1 3 7 5 31 1 127 17 73 11 2047 13 8191 43 151 257 131071 19 524287 41
   Zs 3 1
2 1 13 5 121 7 1093 41 757 61 88573 73 797161 547 4561 3281 64570081 703 581130733 1181
   Zs 4 1
3 5 7 17 341 13 5461 257 1387 41 1398101 241 22369621 3277 49981 65537 5726623061 4033 91625968981 61681
   Zs 5 1
4 3 31 13 781 7 19531 313 15751 521 12207031 601 305175781 13021 315121 195313 190734863281 5167 4768371582031 375601
   Zs 6 1
5 7 43 37 311 31 55987 1297 46873 1111 72559411 1261 2612138803 5713 1406371 1679617 3385331888947 46441 121871948002099 1634221
   Zs 7 1
6 1 19 25 2801 43 137257 1201 39331 2101 329554457 2353 16148168401 102943 4956001 2882401 38771752331201 117307 1899815864228857 1129901
   Zs 3 2
1 5 19 13 211 7 2059 97 1009 11 175099 61 1586131 463 3571 6817 129009091 577 1161737179 4621
   Zs 5 3
2 1 49 17 1441 19 37969 353 19729 421 24325489 481 609554401 10039 216001 198593 381405156481 12979 9536162033329 288961
   Zs 7 3
4 5 79 29 4141 37 205339 1241 127639 341 494287399 2041 24221854021 82573 3628081 2885681 58157596211761 109117 2849723505777919 4871281
   Zs 7 5
2 3 109 37 6841 13 372709 1513 176149 1661 964249309 1801 47834153641 75139 3162961 3077713 115933787267041 30133 5689910849522509 3949201

jq

Adapted from Wren

Works with: jq

Works with gojq, the Go implementation of jq, and with fq.

The integer arithmetic precision of the C implementation of jq is limited, but is sufficient for the specific set of tasks defined in this entry.

Generic Functions

# Input: an integer
def isPrime:
  . as $n
  | if   ($n < 2)       then false
    elif ($n % 2 == 0)  then $n == 2
    elif ($n % 3 == 0)  then $n == 3
    else 5
    | until( . <= 0;
        if .*. > $n then -1
	elif ($n % . == 0) then 0
        else . + 2
        |  if ($n % . == 0) then 0
           else . + 4
           end
        end)
    | . == -1
    end;

# divisors as an unsorted stream
def divisors:
  if . == 1 then 1
  else . as $n
  | label $out
  | range(1; $n) as $i
  | ($i * $i) as $i2
  | if $i2 > $n then break $out
    else if $i2 == $n
         then $i
         elif ($n % $i) == 0
         then $i, ($n/$i)
         else empty
	 end
    end
  end;

def gcd(a; b):
  # subfunction expects [a,b] as input
  # i.e. a ~ .[0] and b ~ .[1]
  def rgcd: if .[1] == 0 then .[0]
         else [.[1], .[0] % .[1]] | rgcd
         end;
  [a,b] | rgcd;

# To take advantage of gojq's arbitrary-precision integer arithmetic:
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in);

The Zsigmondy Function

# Input: $n
def zs($a; $b):
  . as $n
  | (($a|power($n)) - ($b|power($n))) as $dn
  | if ($dn|isPrime) then $dn
    else ([$dn|divisors]|sort) as $divs
    | [range(1; $n) as $m | ($a|power($m)) - ($b|power($m))] as $dms
    | first( range( $divs|length-1; 0 ; -1) as $i
        | if (all($dms[]; gcd(.; $divs[$i]) == 1 )) then $divs[$i] else empty end)
       // 1
  end;

The Task

def abs:
  [ [2, 1], [3, 1], [4, 1], [5, 1], [6, 1], [7, 1], [3, 2], [5, 3], [7, 3], [7, 5] ];

abs[] as [$a, $b]
| (if ($a == 7 and $b != 3) then 18 else 20 end) as $lim
  | "Zsigmony(\($a), \($b)) - first \($lim) terms:",
     [range(1; $lim+1) | zs($a; $b)]

Invocation: gojq -nrcf zsigmondy-numbers.jq

Output:

Zsigmony(2, 1) - first 20 terms:
[1,3,7,5,31,1,127,17,73,11,2047,13,8191,43,151,257,131071,19,524287,41]
Zsigmony(3, 1) - first 20 terms:
[2,1,13,5,121,7,1093,41,757,61,88573,73,797161,547,4561,3281,64570081,703,581130733,1181]
Zsigmony(4, 1) - first 20 terms:
[3,5,7,17,341,13,5461,257,1387,41,1398101,241,22369621,3277,49981,65537,5726623061,4033,91625968981,61681]
Zsigmony(5, 1) - first 20 terms:
[4,3,31,13,781,7,19531,313,15751,521,12207031,601,305175781,13021,315121,195313,190734863281,5167,4768371582031,375601]
Zsigmony(6, 1) - first 20 terms:
[5,7,43,37,311,31,55987,1297,46873,1111,72559411,1261,2612138803,5713,1406371,1679617,3385331888947,46441,121871948002099,1634221]
Zsigmony(7, 1) - first 18 terms:
[6,1,19,25,2801,43,137257,1201,39331,2101,329554457,2353,16148168401,102943,4956001,2882401,38771752331201,117307]
Zsigmony(3, 2) - first 20 terms:
[1,5,19,13,211,7,2059,97,1009,11,175099,61,1586131,463,3571,6817,129009091,577,1161737179,4621]
Zsigmony(5, 3) - first 20 terms:
[2,1,49,17,1441,19,37969,353,19729,421,24325489,481,609554401,10039,216001,198593,381405156481,12979,9536162033329,288961]
Zsigmony(7, 3) - first 20 terms:
[4,5,79,29,4141,37,205339,1241,127639,341,494287399,2041,24221854021,82573,3628081,2885681,58157596211761,109117,2849723505777919,4871281]
Zsigmony(7, 5) - first 18 terms:
[2,3,109,37,6841,13,372709,1513,176149,1661,964249309,1801,47834153641,75139,3162961,3077713,115933787267041,30133]

Julia

""" Rosetta code task: rosettacode.org/wiki/Zsigmondy_numbers """

using Primes

function divisors(n)
    f = [one(n)]
    for (p,e) in factor(n)
        f = reduce(vcat, [f*p^j for j in 1:e], init=f)
    end
    return length(f) == 1 ? [one(n), n] : sort!(f)
end

function Zs(n, a, b)
    @assert a > b
    dexpms = map(i -> a^i - b^i, 1:n-1)
    dexpn = a^n - b^n
    return maximum(reduce(vcat, filter(d -> all(k -> gcd(k, d) == 1, dexpms), divisors(dexpn))))
end

tests = [(2, 1), (3, 1), (4, 1), (5, 1), (6, 1), (7, 1), (3, 2), (5, 3), (7, 3), (7, 5)]
for (a, b) in tests
    println("\nZsigmondy(n, $a, $b): ", join([Zs(n, a, b) for n in 1:20], ", "))
end

Output:

Zsigmondy(n, 2, 1): 1, 3, 7, 5, 31, 1, 127, 17, 73, 11, 2047, 13, 8191, 43, 151, 257, 131071, 19, 524287, 41

Zsigmondy(n, 3, 1): 2, 1, 13, 5, 121, 7, 1093, 41, 757, 61, 88573, 73, 797161, 547, 4561, 3281, 64570081, 703, 581130733, 1181

Zsigmondy(n, 4, 1): 3, 5, 7, 17, 341, 13, 5461, 257, 1387, 41, 1398101, 241, 22369621, 3277, 49981, 65537, 5726623061, 4033, 91625968981, 61681

Zsigmondy(n, 5, 1): 4, 3, 31, 13, 781, 7, 19531, 313, 15751, 521, 12207031, 601, 305175781, 13021, 315121, 195313, 190734863281, 5167, 4768371582031, 375601

Zsigmondy(n, 6, 1): 5, 7, 43, 37, 311, 31, 55987, 1297, 46873, 1111, 72559411, 1261, 2612138803, 5713, 1406371, 1679617, 3385331888947, 46441, 121871948002099, 1634221

Zsigmondy(n, 7, 1): 6, 1, 19, 25, 2801, 43, 137257, 1201, 39331, 2101, 329554457, 2353, 16148168401, 102943, 4956001, 2882401, 38771752331201, 117307, 1899815864228857, 1129901

Zsigmondy(n, 3, 2): 1, 5, 19, 13, 211, 7, 2059, 97, 1009, 11, 175099, 61, 1586131, 463, 3571, 6817, 129009091, 577, 1161737179, 4621

Zsigmondy(n, 5, 3): 2, 1, 49, 17, 1441, 19, 37969, 353, 19729, 421, 24325489, 481, 609554401, 10039, 216001, 198593, 381405156481, 12979, 9536162033329, 288961

Zsigmondy(n, 7, 3): 4, 5, 79, 29, 4141, 37, 205339, 1241, 127639, 341, 494287399, 2041, 24221854021, 82573, 3628081, 2885681, 58157596211761, 109117, 2849723505777919, 4871281

Zsigmondy(n, 7, 5): 2, 3, 109, 37, 6841, 13, 372709, 1513, 176149, 1661, 964249309, 1801, 47834153641, 75139, 3162961, 3077713, 115933787267041, 30133, 5689910849522509, 3949201

Mathematica/Wolfram Language

The Function: Return the Zsigmondy-number to a,b,n integer: Program:

Attributes[zsigmondy] = Listable;
zsigmondy[a_Integer, b_Integer, 1] := a - b /; a >= b;
zsigmondy[a_Integer, b_Integer, n_Integer] := (
    hatvanyok = a^Range[n] - b^Range[n];
   kishatvany = a^Range[n - 1] - b^Range[n - 1];
   biggestelement = Part[hatvanyok, n];
   divisors = Divisors[biggestelement];
   divisorpairs = Join @@ (Thread[{#, kishatvany}] & /@ divisors);
   coprimepairs = Select[divisorpairs, CoprimeQ[#[[1]], #[[2]]] &];
   firstelement = coprimepairs[[All, 1]];
   revesedelements = ReverseSort[firstelement];
   Commonest[revesedelements, 1]) /; a >= b

Attributes[zsigmondy] = Listable;
zsigmondy[a_Integer, b_Integer, 1] := a - b /; a >= b;
zsigmondy[a_Integer, b_Integer, n_Integer] := (
    hatvanyok = a^Range[n] - b^Range[n];
   kishatvany = a^Range[n - 1] - b^Range[n - 1];
   biggestelement = Part[hatvanyok, n];
   divisors = Divisors[biggestelement];
   divisorpairs = Join @@ (Thread[{#, kishatvany}] & /@ divisors);
   coprimepairs = Select[divisorpairs, CoprimeQ[#[[1]], #[[2]]] &];
   firstelement = coprimepairs[[All, 1]];
   revesedelements = ReverseSort[firstelement];
   Commonest[revesedelements, 1]) /; a >= b

l1 = Table[zsigmondy[2, 1, k], {k, Range[10]}];
l2 = Table[zsigmondy[3, 1, k], {k, Range[10]}];
l3 = Table[zsigmondy[4, 1, k], {k, Range[10]}];
l4 = Table[zsigmondy[5, 1, k], {k, Range[10]}];
l5 = Table[zsigmondy[6, 1, k], {k, Range[10]}];
l6 = Table[zsigmondy[7, 1, k], {k, Range[10]}];
l7 = Table[zsigmondy[3, 2, k], {k, Range[10]}];
l8 = Table[zsigmondy[5, 3, k], {k, Range[10]}];
l9 = Table[zsigmondy[7, 3, k], {k, Range[10]}];
l10 = Table[zsigmondy[7, 5, k], {k, Range[10]}];
matrix = {{l1}, {l2}, {l3}, {l4}, {l5}, {l6}, {l7}, {l8}, {l9}, {l10}};
data1 = List["Zsigmondy:2,1,n", "Zsigmondy:3,1,n", "Zsigmondy:4,1,n", 
   "Zsigmondy:5,1,n", "Zsigmondy:6,1,n", "Zsigmondy:7,1,n", 
   "Zsigmondy:3,2,n", "Zsigmondy:5,3,n", "Zsigmondy:7,3,n", 
   "Zsigmondy:7,5,n"];
data2 = matrix;

Grid[Transpose@{data1, data2}~
  Prepend~{"pair of numbers", "Zsigmondy numbers"}, 
 Dividers -> {All, {1 -> True, 2 -> True}}]

output:

Grid[{{"pair of numbers", 
   "Zsigmondy numbers"}, {"Zsigmondy:2,1,n", {{{1}, {3}, {7}, {5}, \
{31}, {1}, {127}, {17}, {73}, {11}}}}, 
     {"Zsigmondy:3,1,n", {{{2}, {1}, {13}, {5}, {121}, {7}, {1093}, \
{41}, {757}, {61}}}}, 
     {"Zsigmondy:4,1,n", {{{3}, {5}, {7}, {17}, {341}, {13}, {5461}, \
{257}, {1387}, {41}}}}, 
     {"Zsigmondy:5,1,n", {{{4}, {3}, {31}, {13}, {781}, {7}, {19531}, \
{313}, {15751}, {521}}}}, 
     {"Zsigmondy:6,1,n", {{{5}, {7}, {43}, {37}, {311}, {31}, \
{55987}, {1297}, {46873}, {1111}}}}, 
     {"Zsigmondy:7,1,n", {{{6}, {1}, {19}, {25}, {2801}, {43}, \
{137257}, {1201}, {39331}, {2101}}}}, 
     {"Zsigmondy:3,2,n", {{{1}, {5}, {19}, {13}, {211}, {7}, {2059}, \
{97}, {1009}, {11}}}}, 
     {"Zsigmondy:5,3,n", {{{2}, {1}, {49}, {17}, {1441}, {19}, \
{37969}, {353}, {19729}, {421}}}}, 
     {"Zsigmondy:7,3,n", {{{4}, {5}, {79}, {29}, {4141}, {37}, \
{205339}, {1241}, {127639}, {341}}}}, 
     {"Zsigmondy:7,5,n", {{{2}, {3}, {109}, {37}, {6841}, {13}, \
{372709}, {1513}, {176149}, {1661}}}}}, 
 Dividers -> {All, {1 -> True, 2 -> True}}, 
   ItemSize -> {Automatic, Automatic}]

Perl

Translation of: Raku

Library: ntheory

use v5.36;
use experimental 'for_list';
use ntheory <gcd divisors vecall>;

sub Zsigmondy ($A, $B, $c) {
    my @aexp = map { $A ** $_ } 0..$c;
    my @bexp = map { $B ** $_ } 0..$c;
    my @z;
    for my $n (1..$c) {
        for my $d (sort { $b <=> $a } divisors $aexp[$n] - $bexp[$n]) {
            push @z, $d and last if vecall { gcd($aexp[$_] - $bexp[$_], $d) == 1 } 1..$n-1
        }
        return @z if 20 == @z
    }
}

for my($name,$a,$b) (
    'A064078: Zsigmondy(n,2,1)', 2,1,
    'A064079: Zsigmondy(n,3,1)', 3,1,
    'A064080: Zsigmondy(n,4,1)', 4,1,
    'A064081: Zsigmondy(n,5,1)', 5,1,
    'A064082: Zsigmondy(n,6,1)', 6,1,
    'A064083: Zsigmondy(n,7,1)', 7,1,
    'A109325: Zsigmondy(n,3,2)', 3,2,
    'A109347: Zsigmondy(n,5,3)', 5,3,
    'A109348: Zsigmondy(n,7,3)', 7,3,
    'A109349: Zsigmondy(n,7,5)', 7,5,
) {
    say "\n$name:\n" . join ' ', Zsigmondy($a, $b, 20)
}

Output:

A064078: Zsigmondy(n,2,1):
1 3 7 5 31 1 127 17 73 11 2047 13 8191 43 151 257 131071 19 524287 41

A064079: Zsigmondy(n,3,1):
2 1 13 5 121 7 1093 41 757 61 88573 73 797161 547 4561 3281 64570081 703 581130733 1181

A064080: Zsigmondy(n,4,1):
3 5 7 17 341 13 5461 257 1387 41 1398101 241 22369621 3277 49981 65537 5726623061 4033 91625968981 61681

A064081: Zsigmondy(n,5,1):
4 3 31 13 781 7 19531 313 15751 521 12207031 601 305175781 13021 315121 195313 190734863281 5167 4768371582031 375601

A064082: Zsigmondy(n,6,1):
5 7 43 37 311 31 55987 1297 46873 1111 72559411 1261 2612138803 5713 1406371 1679617 3385331888947 46441 121871948002099 1634221

A064083: Zsigmondy(n,7,1):
6 1 19 25 2801 43 137257 1201 39331 2101 329554457 2353 16148168401 102943 4956001 2882401 38771752331201 117307 1899815864228857 1129901

A109325: Zsigmondy(n,3,2):
1 5 19 13 211 7 2059 97 1009 11 175099 61 1586131 463 3571 6817 129009091 577 1161737179 4621

A109347: Zsigmondy(n,5,3):
2 1 49 17 1441 19 37969 353 19729 421 24325489 481 609554401 10039 216001 198593 381405156481 12979 9536162033329 288961

A109348: Zsigmondy(n,7,3):
4 5 79 29 4141 37 205339 1241 127639 341 494287399 2041 24221854021 82573 3628081 2885681 58157596211761 109117 2849723505777919 4871281

A109349: Zsigmondy(n,7,5):
2 3 109 37 6841 13 372709 1513 176149 1661 964249309 1801 47834153641 75139 3162961 3077713 115933787267041 30133 5689910849522509 3949201

Phix

Translation of: Wren

with javascript_semantics
function Zsigmondy(integer n, a, b)
    atom dn = power(a,n) - power(b,n)
    if is_prime(dn) then return dn end if
    sequence dms = repeat(0,n-1)
    for m=1 to n-1 do
        dms[m] = power(a,m)-power(b,m)
    end for
    for d in reverse(factors(dn)) do
        for m=1 to n-1 do
            if gcd(dms[m],d)!=1 then exit end if
            if m=n-1 then return d end if
        end for
    end for
    return 1
end function

for ab in {{2,1},{3,1},{4,1},{5,1},{6,1},{7,1},{3,2},{5,3},{7,3},{7,5}} do
    integer {a,b} = ab, lim = iff(machine_bits()=32 and a=7?18:20)
    string res = join(apply(true,Zsigmondy,{tagset(lim),a,b}),fmt:="%d")
    printf(1,"Zsigmondy(1..%d,%d,%d): %s\n",{lim,a,b,res})
end for

Output:

As per Wren, 7¹⁹ exceeds 2⁵³ (plus is_prime/factors actually crash, as they should when precision is exceeded) and hence 32-bit limits itself to 18 entries when a=7.

Zsigmondy(1..20,2,1): 1 3 7 5 31 1 127 17 73 11 89 13 8191 43 151 257 131071 19 524287 41
Zsigmondy(1..20,3,1): 2 1 13 5 121 7 1093 41 757 61 88573 73 797161 547 4561 3281 64570081 703 581130733 1181
Zsigmondy(1..20,4,1): 3 5 7 17 341 13 5461 257 1387 41 1398101 241 22369621 3277 49981 65537 5726623061 4033 91625968981 61681
Zsigmondy(1..20,5,1): 1 3 31 13 781 7 19531 313 15751 521 12207031 601 305175781 13021 315121 195313 190734863281 5167 4768371582031 375601
Zsigmondy(1..20,6,1): 5 7 43 37 311 31 55987 1297 46873 1111 72559411 1261 2612138803 5713 1406371 1679617 3385331888947 46441 121871948002099 1634221
Zsigmondy(1..20,7,1): 1 1 19 25 2801 43 137257 1201 39331 2101 329554457 2353 16148168401 102943 4956001 2882401 38771752331201 117307 1899815864228857 1129901
Zsigmondy(1..20,3,2): 1 5 19 13 211 7 71 97 1009 11 7613 61 29927 463 3571 6817 129009091 577 745181 4621
Zsigmondy(1..20,5,3): 2 1 49 17 1441 19 37969 353 19729 421 24325489 481 609554401 10039 216001 198593 381405156481 12979 9536162033329 288961
Zsigmondy(1..20,7,3): 1 5 79 29 4141 37 205339 1241 127639 341 494287399 2041 24221854021 82573 3628081 2885681 58157596211761 109117 2849723505777919 4871281
Zsigmondy(1..20,7,5): 2 3 109 37 6841 13 372709 1513 176149 1661 964249309 1801 47834153641 75139 3162961 3077713 115933787267041 30133 5689910849522509 3949201

Python

''' Rosetta code task: rosettacode.org/wiki/Zsigmondy_numbers '''

from math import gcd
from sympy import divisors


def zsig(num, aint, bint):
    ''' Get Zs(n, a, b) in task. '''
    assert aint > bint
    dexpms = [aint**i - bint**i for i in range(1, num)]
    dexpn = aint**num - bint**num
    return max([d for d in divisors(dexpn) if all(gcd(k, d) == 1 for k in dexpms)])


tests = [(2, 1), (3, 1), (4, 1), (5, 1), (6, 1),
         (7, 1), (3, 2), (5, 3), (7, 3), (7, 5)]
for (a, b) in tests:
    print(f'\nZsigmondy(n, {a}, {b}):', ', '.join(
        [str(zsig(n, a, b)) for n in range(1, 21)]))

Output:

Zsigmondy(n, 2, 1): 1, 3, 7, 5, 31, 1, 127, 17, 73, 11, 2047, 13, 8191, 43, 151, 257, 131071, 19, 524287, 41

Zsigmondy(n, 3, 1): 2, 1, 13, 5, 121, 7, 1093, 41, 757, 61, 88573, 73, 797161, 547, 4561, 3281, 64570081, 703, 581130733, 1181

Zsigmondy(n, 4, 1): 3, 5, 7, 17, 341, 13, 5461, 257, 1387, 41, 1398101, 241, 22369621, 3277, 49981, 65537, 5726623061, 4033, 91625968981, 61681

Zsigmondy(n, 5, 1): 4, 3, 31, 13, 781, 7, 19531, 313, 15751, 521, 12207031, 601, 305175781, 13021, 315121, 195313, 190734863281, 5167, 4768371582031, 375601

Zsigmondy(n, 6, 1): 5, 7, 43, 37, 311, 31, 55987, 1297, 46873, 1111, 72559411, 1261, 2612138803, 5713, 1406371, 1679617, 3385331888947, 46441, 121871948002099, 1634221

Zsigmondy(n, 7, 1): 6, 1, 19, 25, 2801, 43, 137257, 1201, 39331, 2101, 329554457, 2353, 16148168401, 102943, 4956001, 2882401, 38771752331201, 117307, 1899815864228857, 1129901

Zsigmondy(n, 3, 2): 1, 5, 19, 13, 211, 7, 2059, 97, 1009, 11, 175099, 61, 1586131, 463, 3571, 6817, 129009091, 577, 1161737179, 4621

Zsigmondy(n, 5, 3): 2, 1, 49, 17, 1441, 19, 37969, 353, 19729, 421, 24325489, 481, 609554401, 10039, 216001, 198593, 381405156481, 12979, 9536162033329, 288961

Zsigmondy(n, 7, 3): 4, 5, 79, 29, 4141, 37, 205339, 1241, 127639, 341, 494287399, 2041, 24221854021, 82573, 3628081, 2885681, 58157596211761, 109117, 2849723505777919, 4871281

Zsigmondy(n, 7, 5): 2, 3, 109, 37, 6841, 13, 372709, 1513, 176149, 1661, 964249309, 1801, 47834153641, 75139, 3162961, 3077713, 115933787267041, 30133, 5689910849522509, 3949201

Raku

First twenty elements of each.

use Prime::Factor;

sub Zsigmondy ($a, $b) {
    my @aexp = 1, * × $a … *;
    my @bexp = 1, * × $b … *;
    (1..∞).map: -> $n {
        (@aexp[$n] - @bexp[$n]).&divisors.sort(-*).first: -> $d {
            all (1..^$n).map: { (@aexp[$_] - @bexp[$_]) gcd $d == 1 }
        }
    }
}

for '064078', (2,1), '064079', (3,1), '064080', (4,1), '064081', (5,1), '064082', (6,1),
    '064083', (7,1), '109325', (3,2), '109347', (5,3), '109348', (7,3), '109349', (7,5)
  -> $oeis, $seq {
    say "\nA$oeis: Zsigmondy(n,$seq[0],$seq[1]):\n" ~ Zsigmondy(|$seq)[^20]
}

Output:

A064078: Zsigmondy(n,2,1):
1 3 7 5 31 1 127 17 73 11 2047 13 8191 43 151 257 131071 19 524287 41

A064079: Zsigmondy(n,3,1):
2 1 13 5 121 7 1093 41 757 61 88573 73 797161 547 4561 3281 64570081 703 581130733 1181

A064080: Zsigmondy(n,4,1):
3 5 7 17 341 13 5461 257 1387 41 1398101 241 22369621 3277 49981 65537 5726623061 4033 91625968981 61681

A064081: Zsigmondy(n,5,1):
4 3 31 13 781 7 19531 313 15751 521 12207031 601 305175781 13021 315121 195313 190734863281 5167 4768371582031 375601

A064082: Zsigmondy(n,6,1):
5 7 43 37 311 31 55987 1297 46873 1111 72559411 1261 2612138803 5713 1406371 1679617 3385331888947 46441 121871948002099 1634221

A064083: Zsigmondy(n,7,1):
6 1 19 25 2801 43 137257 1201 39331 2101 329554457 2353 16148168401 102943 4956001 2882401 38771752331201 117307 1899815864228857 1129901

A109325: Zsigmondy(n,3,2):
1 5 19 13 211 7 2059 97 1009 11 175099 61 1586131 463 3571 6817 129009091 577 1161737179 4621

A109347: Zsigmondy(n,5,3):
2 1 49 17 1441 19 37969 353 19729 421 24325489 481 609554401 10039 216001 198593 381405156481 12979 9536162033329 288961

A109348: Zsigmondy(n,7,3):
4 5 79 29 4141 37 205339 1241 127639 341 494287399 2041 24221854021 82573 3628081 2885681 58157596211761 109117 2849723505777919 4871281

A109349: Zsigmondy(n,7,5):
2 3 109 37 6841 13 372709 1513 176149 1661 964249309 1801 47834153641 75139 3162961 3077713 115933787267041 30133 5689910849522509 3949201

Wren

Normal arithmetic

Library: Wren-math

Library: Wren-fmt

Shows the first 20 terms for each radix set apart from [7, 1] and [7, 5] where I've had to limit the number of terms to 18 for now. This is because the 19th term is being calculated incorrectly due to 7^19 exceeding Wren's safe integer limit of 2^53. It's only by good fortune that [7, 3] works correctly.

import "./math" for Int
import "./fmt" for Fmt

var zs = Fn.new { |n, a, b|
    var dn = a.pow(n) - b.pow(n)
    if (Int.isPrime(dn)) return dn
    var divs = Int.divisors(dn)
    var dms = (1...n).map { |m| a.pow(m) - b.pow(m) }.toList
    for (i in divs.count-1..1) {
        if (dms.all { |dm| Int.gcd(dm, divs[i]) == 1 }) return divs[i]
    }
    return 1
}

var abs = [ [2, 1], [3, 1], [4, 1], [5, 1], [6, 1], [7, 1], [3, 2], [5, 3], [7, 3], [7, 5] ]
for (ab in abs) {
    var a = ab[0]
    var b = ab[1]
    var lim = 20
    if (a == 7 && b != 3) lim = 18
    System.print("Zsigmony(n, %(a), %(b)) - first %(lim) terms:")
    Fmt.print("$d", (1..lim).map { |n| zs.call(n, a, b) }.toList)
    System.print()
}

Output:

Zsigmony(n, 2, 1) - first 20 terms:
1 3 7 5 31 1 127 17 73 11 2047 13 8191 43 151 257 131071 19 524287 41

Zsigmony(n, 3, 1) - first 20 terms:
2 1 13 5 121 7 1093 41 757 61 88573 73 797161 547 4561 3281 64570081 703 581130733 1181

Zsigmony(n, 4, 1) - first 20 terms:
3 5 7 17 341 13 5461 257 1387 41 1398101 241 22369621 3277 49981 65537 5726623061 4033 91625968981 61681

Zsigmony(n, 5, 1) - first 20 terms:
4 3 31 13 781 7 19531 313 15751 521 12207031 601 305175781 13021 315121 195313 190734863281 5167 4768371582031 375601

Zsigmony(n, 6, 1) - first 20 terms:
5 7 43 37 311 31 55987 1297 46873 1111 72559411 1261 2612138803 5713 1406371 1679617 3385331888947 46441 121871948002099 1634221

Zsigmony(n, 7, 1) - first 18 terms:
6 1 19 25 2801 43 137257 1201 39331 2101 329554457 2353 16148168401 102943 4956001 2882401 38771752331201 117307

Zsigmony(n, 3, 2) - first 20 terms:
1 5 19 13 211 7 2059 97 1009 11 175099 61 1586131 463 3571 6817 129009091 577 1161737179 4621

Zsigmony(n, 5, 3) - first 20 terms:
2 1 49 17 1441 19 37969 353 19729 421 24325489 481 609554401 10039 216001 198593 381405156481 12979 9536162033329 288961

Zsigmony(n, 7, 3) - first 20 terms:
4 5 79 29 4141 37 205339 1241 127639 341 494287399 2041 24221854021 82573 3628081 2885681 58157596211761 109117 2849723505777919 4871281

Zsigmony(n, 7, 5) - first 18 terms:
2 3 109 37 6841 13 372709 1513 176149 1661 964249309 1801 47834153641 75139 3162961 3077713 115933787267041 30133

BigInt

Library: Wren-big

Library: Wren-seq

However, we can deal with integers of any size by switching to BigInt.

import "./big" for BigInt
import "./seq" for Lst
import "./fmt" for Fmt

var divisors = Fn.new { |n|
    var factors = BigInt.primeFactors(n)
    var divs = [BigInt.one]
    for (p in factors) {
        for (i in 0...divs.count) divs.add(divs[i]*p)
    }
    return Lst.prune(divs.sort { |i, j| i >= j })
}

var zs = Fn.new { |n, a, b|
    a = BigInt.new(a)
    b = BigInt.new(b)
    var dn = a.pow(n) - b.pow(n)
    if (dn.isPrime) return dn
    var divs = divisors.call(dn)
    var dms = (1...n).map { |m| a.pow(m) - b.pow(m) }.toList
    for (div in divs) {
        if (dms.all { |dm| BigInt.gcd(dm, div) == 1 }) return div
    }
    return BigInt.one
}

var abs = [ [2, 1], [3, 1], [4, 1], [5, 1], [6, 1], [7, 1], [3, 2], [5, 3], [7, 3], [7, 5] ]
var lim = 20
for (ab in abs) {
    var a = ab[0]
    var b = ab[1]
    System.print("Zsigmony(n, %(a), %(b)) - first %(lim) terms:")
    Fmt.print("$i", (1..lim).map { |n| zs.call(n, a, b) }.toList)
    System.print()
}

Output:

Zsigmony(n, 2, 1) - first 20 terms:
1 3 7 5 31 1 127 17 73 11 2047 13 8191 43 151 257 131071 19 524287 41

Zsigmony(n, 3, 1) - first 20 terms:
2 1 13 5 121 7 1093 41 757 61 88573 73 797161 547 4561 3281 64570081 703 581130733 1181

Zsigmony(n, 4, 1) - first 20 terms:
3 5 7 17 341 13 5461 257 1387 41 1398101 241 22369621 3277 49981 65537 5726623061 4033 91625968981 61681

Zsigmony(n, 5, 1) - first 20 terms:
4 3 31 13 781 7 19531 313 15751 521 12207031 601 305175781 13021 315121 195313 190734863281 5167 4768371582031 375601

Zsigmony(n, 6, 1) - first 20 terms:
5 7 43 37 311 31 55987 1297 46873 1111 72559411 1261 2612138803 5713 1406371 1679617 3385331888947 46441 121871948002099 1634221

Zsigmony(n, 7, 1) - first 20 terms:
6 1 19 25 2801 43 137257 1201 39331 2101 329554457 2353 16148168401 102943 4956001 2882401 38771752331201 117307 1899815864228857 1129901

Zsigmony(n, 3, 2) - first 20 terms:
1 5 19 13 211 7 2059 97 1009 11 175099 61 1586131 463 3571 6817 129009091 577 1161737179 4621

Zsigmony(n, 5, 3) - first 20 terms:
2 1 49 17 1441 19 37969 353 19729 421 24325489 481 609554401 10039 216001 198593 381405156481 12979 9536162033329 288961

Zsigmony(n, 7, 3) - first 20 terms:
4 5 79 29 4141 37 205339 1241 127639 341 494287399 2041 24221854021 82573 3628081 2885681 58157596211761 109117 2849723505777919 4871281

Zsigmony(n, 7, 5) - first 20 terms:
2 3 109 37 6841 13 372709 1513 176149 1661 964249309 1801 47834153641 75139 3162961 3077713 115933787267041 30133 5689910849522509 3949201