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

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

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

Phix

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,"Zsigmony(1..%d,%d,%d): %s\n",{lim,a,b,res})
end for

Output:

As per the Wren entry, 7¹⁹ exceeds 2⁵³ and hence 32-bit is limited to 18 entries when a is 7.

Zsigmony(1..20,2,1): 1 3 7 5 31 1 127 17 73 11 89 13 8191 43 151 257 131071 19 524287 41
Zsigmony(1..20,3,1): 2 1 13 5 121 7 1093 41 757 61 88573 73 797161 547 4561 3281 64570081 703 581130733 1181
Zsigmony(1..20,4,1): 3 5 7 17 341 13 5461 257 1387 41 1398101 241 22369621 3277 49981 65537 5726623061 4033 91625968981 61681
Zsigmony(1..20,5,1): 1 3 31 13 781 7 19531 313 15751 521 12207031 601 305175781 13021 315121 195313 190734863281 5167 4768371582031 375601
Zsigmony(1..20,6,1): 5 7 43 37 311 31 55987 1297 46873 1111 72559411 1261 2612138803 5713 1406371 1679617 3385331888947 46441 121871948002099 1634221
Zsigmony(1..20,7,1): 1 1 19 25 2801 43 137257 1201 39331 2101 329554457 2353 16148168401 102943 4956001 2882401 38771752331201 117307 1899815864228857 1129901
Zsigmony(1..20,3,2): 1 5 19 13 211 7 71 97 1009 11 7613 61 29927 463 3571 6817 129009091 577 745181 4621
Zsigmony(1..20,5,3): 2 1 49 17 1441 19 37969 353 19729 421 24325489 481 609554401 10039 216001 198593 381405156481 12979 9536162033329 288961
Zsigmony(1..20,7,3): 1 5 79 29 4141 37 205339 1241 127639 341 494287399 2041 24221854021 82573 3628081 2885681 58157596211761 109117 2849723505777919 4871281
Zsigmony(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, * × $a … *;
    my @bexp = 1, $b, * × $b … *;
    (1..∞).map: -> $n {
        my @divisors = divisors(@aexp[$n] - @bexp[$n]).sort: -*;
        @divisors.first: -> $d { all (1..^$n).map: -> $m { (@aexp[$m] - @bexp[$m]) gcd $d == 1 } }
    }
}

for '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)
  -> $name, $seq { say "\n$name:\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

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 - probably due to intermediate calculations exceeding Wren's safe integer limit of 2^53, though I'll need to investigate further to find the exact cause.

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