Zsigmondy numbers
You are encouraged to solve this task according to the task description, using any language you may know.
Zsigmondy numbers n to a, b, are the greatest divisor of an - bn that is coprime to am - bm 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 an - bn and return the largest that is coprime to all of am - bm, where m is each of the positive integers 1 to n - 1.
When n = 4, 24 - 14 = 15. The divisors of 15 are 1, 3, 5, and 15.
For m = 1, 2, 3 we get 2-1, 22-12, 23-13, 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, 26 - 16 = 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 an - bn is prime, then Zs(n,a,b) will be equal to that prime. 25 - 15 = 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
ALGOL 68
64-bit integers are needed to show some of the "larger" a, b values, adjust the INTEGER mode and isqrt procedure accordingly, if necessary for your compiler/interpreter.
BEGIN # find some Zsigmondy numbers #
# integral mode that has precision to suit the task - adjust to suit #
MODE INTEGER = LONG INT;
# returns the square root of n, use the sqrt routine appropriate for #
# the INTEGER mode #
PROC isqrt = ( INTEGER n )INTEGER: ENTIER long sqrt( n );
# returns the gcd of m and n #
PRIO GCD = 1;
OP GCD = ( INTEGER m, n )INTEGER:
BEGIN
INTEGER a := ABS m, b := ABS n;
WHILE b /= 0 DO
INTEGER new a = b;
b := a MOD b;
a := new a
OD;
a
END # GCD # ;
# returns TRUE if all m are coprime to n, FALSE otherwise #
PRIO ALLCOPRIME = 1;
OP ALLCOPRIME = ( []INTEGER m, INTEGER n )BOOL:
BEGIN
BOOL coprime := TRUE;
INTEGER abs n = ABS n;
FOR i FROM LWB m TO UPB m WHILE coprime := ( m[ i ] GCD n ) = 1 DO SKIP OD;
coprime
END # ALLCOPRIME # ;
# returns the sequence of Zsigmondy numbers of n to a, b #
PROC zsigmondy = ( INT n, a, b )[]INTEGER:
BEGIN
[ 1 : n - 1 ]INTEGER am minus bm;
INTEGER am := 1, bm := 1;
FOR m TO n - 1 DO
am minus bm[ m ] := ( am *:= a ) - ( bm *:= b )
OD;
INTEGER an := 1, bn := 1;
[ 1 : n ]INTEGER z;
FOR i TO n DO
INTEGER an minus bn = ( an *:= a ) - ( bn *:= b );
INTEGER largest divisor := 1;
IF am minus bm[ 1 : i - 1 ] ALLCOPRIME an minus bn THEN
largest divisor := an minus bn
ELSE
INTEGER d := 1;
INTEGER max d := isqrt( an minus bn );
WHILE ( d +:= 1 ) <= max d AND largest divisor <= max d DO
IF an minus bn MOD d = 0 THEN
# d is a divisor of a^n - b^n #
IF d > largest divisor THEN
IF am minus bm[ 1 : i - 1 ] ALLCOPRIME d THEN largest divisor := d FI
FI;
INTEGER other d = an minus bn OVER d;
IF other d > d AND other d > largest divisor THEN
IF am minus bm[ 1 : i - 1 ] ALLCOPRIME other d THEN largest divisor := other d FI
FI
FI
OD
FI;
z[ i ] := largest divisor
OD;
z
END # zsigmondy # ;
[,]INT tests = ( ( 2, 1 ), ( 3, 1 ), ( 4, 1 ), ( 5, 1 ), ( 6, 1 )
, ( 7, 1 ), ( 3, 2 ), ( 5, 3 ), ( 7, 3 ), ( 7, 5 )
);
FOR i FROM LWB tests TO UPB tests DO
INT a = tests[ i, 1 ], b = tests[ i, 2 ];
[]INTEGER z = zsigmondy( 20, a, b );
print( ( "zsigmondy( 20, ", whole( a, 0 ), ", ", whole( b, 0 ), " ):", newline, " " ) );
FOR z pos FROM LWB z TO UPB z DO print( ( " ", whole( z[ z pos ], 0 ) ) ) OD;
print( ( newline, newline ) )
OD
END- Output:
zsigmondy( 20, 2, 1 ):
1 3 7 5 31 1 127 17 73 11 2047 13 8191 43 151 257 131071 19 524287 41
zsigmondy( 20, 3, 1 ):
2 1 13 5 121 7 1093 41 757 61 88573 73 797161 547 4561 3281 64570081 703 581130733 1181
zsigmondy( 20, 4, 1 ):
3 5 7 17 341 13 5461 257 1387 41 1398101 241 22369621 3277 49981 65537 5726623061 4033 91625968981 61681
zsigmondy( 20, 5, 1 ):
4 3 31 13 781 7 19531 313 15751 521 12207031 601 305175781 13021 315121 195313 190734863281 5167 4768371582031 375601
zsigmondy( 20, 6, 1 ):
5 7 43 37 311 31 55987 1297 46873 1111 72559411 1261 2612138803 5713 1406371 1679617 3385331888947 46441 121871948002099 1634221
zsigmondy( 20, 7, 1 ):
6 1 19 25 2801 43 137257 1201 39331 2101 329554457 2353 16148168401 102943 4956001 2882401 38771752331201 117307 1899815864228857 1129901
zsigmondy( 20, 3, 2 ):
1 5 19 13 211 7 2059 97 1009 11 175099 61 1586131 463 3571 6817 129009091 577 1161737179 4621
zsigmondy( 20, 5, 3 ):
2 1 49 17 1441 19 37969 353 19729 421 24325489 481 609554401 10039 216001 198593 381405156481 12979 9536162033329 288961
zsigmondy( 20, 7, 3 ):
4 5 79 29 4141 37 205339 1241 127639 341 494287399 2041 24221854021 82573 3628081 2885681 58157596211761 109117 2849723505777919 4871281
zsigmondy( 20, 7, 5 ):
2 3 109 37 6841 13 372709 1513 176149 1661 964249309 1801 47834153641 75139 3162961 3077713 115933787267041 30133 5689910849522509 3949201
Amazing Hopper
#include <basico.h>
#proto zsigmondy(_X_,_Y_,_Z_)
#define TOPE 11
algoritmo
decimales '0'
matrices(a,b)
'2,3,4,5,6,7,3,5,7,7' enlistar en 'a'
'1,1,1,1,1,1,2,3,3,5' enlistar en 'b'
imprimir( "Zsigmondy Numbers\n\n")
iterar para ( k=1, #(k<=10), ++k)
imprimir ("[11,",#(a[k]),",",#(b[k]),"] {")
iterar para( i=1, #(i<=TOPE), ++i )
imprimir( _zsigmondy(i, #(a[k]), #(b[k])),\
solo si ( #(i<TOPE), ","))
siguiente
"}\n", imprimir
siguiente
terminar
subrutinas
zsigmondy(i,a,b)
dn=0
si ' #( dn :=( a^i - b^i ) ), es primo? '
tomar 'dn'
sino
divisores=0
guardar 'divisores de (dn)' en 'divisores'
iterar para( m=1, #(m<i), ++m )
remover según( #(gcd((a^m-b^m),divisores )>1),\
divisores)
siguiente
tomar 'último elemento de (divisores)'
fin si
retornar
- Output:
Zsigmondy Numbers
[11,2,1] {1,3,7,5,31,1,127,17,73,11,2047}
[11,3,1] {2,1,13,5,121,7,1093,41,757,61,88573}
[11,4,1] {3,5,7,17,341,13,5461,257,1387,41,1398101}
[11,5,1] {4,3,31,13,781,7,19531,313,15751,521,12207031}
[11,6,1] {5,7,43,37,311,31,55987,1297,46873,1111,72559411}
[11,7,1] {6,1,19,25,2801,43,137257,1201,39331,2101,329554457}
[11,3,2] {1,5,19,13,211,7,2059,97,1009,11,175099}
[11,5,3] {2,1,49,17,1441,19,37969,353,19729,421,24325489}
[11,7,3] {4,5,79,29,4141,37,205339,1241,127639,341,494287399}
[11,7,5] {2,3,109,37,6841,13,372709,1513,176149,1661,964249309}
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
EasyLang
func isprim num .
if num < 2
return 0
.
i = 2
while i <= sqrt num
if num mod i = 0
return 0
.
i += 1
.
return 1
.
proc sort . d[] .
for i = 1 to len d[] - 1
for j = i + 1 to len d[]
if d[j] < d[i]
swap d[j] d[i]
.
.
.
.
proc divisors num . res[] .
res[] = [ ]
d = 1
while d <= sqrt num
if num mod d = 0
res[] &= d
if d <> sqrt num
res[] &= num / d
.
.
d += 1
.
sort res[]
.
func gcd a b .
while b <> 0
h = b
b = a mod b
a = h
.
return a
.
func zsigmo n a b .
dn = pow a n - pow b n
if isprim dn = 1
return dn
.
divisors dn divs[]
for m = 1 to n - 1
dms[] &= pow a m - pow b m
.
for i = len divs[] downto 1
d = divs[i]
for m = 1 to n - 1
if gcd dms[m] d <> 1
break 1
.
.
if m = n
return d
.
.
return 1
.
proc test . .
test[][] = [ [ 2 1 ] [ 3 1 ] [ 4 1 ] [ 5 1 ] [ 6 1 ] [ 7 1 ] [ 3 2 ] [ 5 3 ] [ 7 3 ] [ 7 5 ] ]
for i to len test[][]
a = test[i][1]
b = test[i][2]
write "Zsigmondy(n, " & a & ", " & b & "):"
for n = 1 to 10
write " " & zsigmo n a b
.
print ""
.
.
test
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
Java
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
public final class ZsigmondyNumbers {
public static void main(String[] args) {
record Pair(int a, int b) {}
List<Pair> pairs = List.of( new Pair(2, 1), new Pair(3, 1), new Pair(4, 1), new Pair(5, 1),
new Pair(6, 1), new Pair(7, 1), new Pair(3, 2), new Pair(5, 3), new Pair(7, 3), new Pair(7, 5) );
for ( Pair pair : pairs ) {
System.out.println("Zsigmondy(n, " + pair.a + ", " + pair.b + ")");
for ( int n = 1; n <= 18; n++ ) {
System.out.print(zsigmondyNumber(n, pair.a, pair.b) + " ");
}
System.out.println(System.lineSeparator());
}
}
private static long zsigmondyNumber(int n, int a, int b) {
final long dn = (long) ( Math.pow(a, n) - Math.pow(b, n) );
if ( isPrime(dn) ) {
return dn;
}
List<Long> divisors = divisors(dn);
for ( int m = 1; m < n; m++ ) {
long dm = (long) ( Math.pow(a, m) - Math.pow(b, m) );
divisors.removeIf( d -> gcd(dm, d) > 1 );
}
return divisors.get(divisors.size() - 1);
}
private static List<Long> divisors(long number) {
List<Long> result = new ArrayList<Long>();
for ( long d = 1; d * d <= number; d++ ) {
if ( number % d == 0 ) {
result.add(d);
if ( d * d < number ) {
result.add(number / d);
}
}
}
Collections.sort(result);
return result;
}
private static boolean isPrime(long number) {
if ( number < 2 ) {
return false;
}
if ( number % 2 == 0 ) {
return number == 2;
}
if ( number % 3 == 0 ) {
return number == 3;
}
int delta = 2;
long k = 5;
while ( k * k <= number ) {
if ( number % k == 0 ) {
return false;
}
k += delta;
delta = 6 - delta;
}
return true;
}
private static long gcd(long a, long b) {
while ( b != 0 ) {
long temp = a; a = b; b = temp % b;
}
return a;
}
}
- Output:
Zsigmondy(n, 2, 1) 1 3 7 5 31 1 127 17 73 11 2047 13 8191 43 151 257 131071 19 Zsigmondy(n, 3, 1) 2 1 13 5 121 7 1093 41 757 61 88573 73 797161 547 4561 3281 64570081 703 Zsigmondy(n, 4, 1) 3 5 7 17 341 13 5461 257 1387 41 1398101 241 22369621 3277 49981 65537 5726623061 4033 Zsigmondy(n, 5, 1) 4 3 31 13 781 7 19531 313 15751 521 12207031 601 305175781 13021 315121 195313 190734863281 5167 Zsigmondy(n, 6, 1) 5 7 43 37 311 31 55987 1297 46873 1111 72559411 1261 2612138803 5713 1406371 1679617 3385331888947 46441 Zsigmondy(n, 7, 1) 6 1 19 25 2801 43 137257 1201 39331 2101 329554457 2353 16148168401 102943 4956001 2882401 38771752331201 117307 Zsigmondy(n, 3, 2) 1 5 19 13 211 7 2059 97 1009 11 175099 61 1586131 463 3571 6817 129009091 577 Zsigmondy(n, 5, 3) 2 1 49 17 1441 19 37969 353 19729 421 24325489 481 609554401 10039 216001 198593 381405156481 12979 Zsigmondy(n, 7, 3) 4 5 79 29 4141 37 205339 1241 127639 341 494287399 2041 24221854021 82573 3628081 2885681 58157596211761 109117 Zsigmondy(n, 7, 5) 2 3 109 37 6841 13 372709 1513 176149 1661 964249309 1801 47834153641 75139 3162961 3077713 115933787267041 30133
jq
Adapted from Wren
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(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:
ClearAll[zsigmondy]
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
data2 = MapThread[
zsigmondy, {{2, 3, 4, 5, 6, 7, 3, 5, 7, 7}, {1, 1, 1, 1, 1, 1, 2,
3, 3, 5}, {Range[10], Range[10], Range[10], Range[10], Range[10],
Range[10], Range[10], Range[10], Range[10], Range[10]}}];
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"];
Grid[Transpose@{data1, data2}~
Prepend~{"pair of numbers", "Zsigmondy numbers"},
Dividers -> {All, {1 -> True, 2 -> True}}]
output:
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}}
Nim
import std/[algorithm, math, strformat]
func isPrime(n: Natural): bool =
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 divisors(n: Positive): seq[int] =
for d in 1..sqrt(n.toFloat).int:
if n mod d == 0:
result.add d
if n div d != d:
result.add n div d
result.sort()
func zs(n, a, b: Positive): int =
let dn = a^n - b^n
if dn.isPrime: return dn
var divs = dn.divisors
for m in 1..<n:
let dm = a^m - b^m
for i in countdown(divs.high, 1):
if gcd(dm, divs[i]) != 1:
divs.delete i
result = divs[^1]
const N = 15
for (a, b) in [(2, 1), (3, 1), (4, 1), (5, 1), (6, 1), (7, 1), (3, 2), (5, 3), (7, 3), (7, 5)]:
echo &"Zsigmondy(n, {a}, {b}) – First {N} terms:"
for n in 1..N:
stdout.write zs(n, a, b), ' '
echo '\n'
- Output:
Zsigmondy(n, 2, 1) – First 15 terms: 1 3 7 5 31 1 127 17 73 11 2047 13 8191 43 151 Zsigmondy(n, 3, 1) – First 15 terms: 2 1 13 5 121 7 1093 41 757 61 88573 73 797161 547 4561 Zsigmondy(n, 4, 1) – First 15 terms: 3 5 7 17 341 13 5461 257 1387 41 1398101 241 22369621 3277 49981 Zsigmondy(n, 5, 1) – First 15 terms: 4 3 31 13 781 7 19531 313 15751 521 12207031 601 305175781 13021 315121 Zsigmondy(n, 6, 1) – First 15 terms: 5 7 43 37 311 31 55987 1297 46873 1111 72559411 1261 2612138803 5713 1406371 Zsigmondy(n, 7, 1) – First 15 terms: 6 1 19 25 2801 43 137257 1201 39331 2101 329554457 2353 16148168401 102943 4956001 Zsigmondy(n, 3, 2) – First 15 terms: 1 5 19 13 211 7 2059 97 1009 11 175099 61 1586131 463 3571 Zsigmondy(n, 5, 3) – First 15 terms: 2 1 49 17 1441 19 37969 353 19729 421 24325489 481 609554401 10039 216001 Zsigmondy(n, 7, 3) – First 15 terms: 4 5 79 29 4141 37 205339 1241 127639 341 494287399 2041 24221854021 82573 3628081 Zsigmondy(n, 7, 5) – First 15 terms: 2 3 109 37 6841 13 372709 1513 176149 1661 964249309 1801 47834153641 75139 3162961
Perl
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
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, 719 exceeds 253 (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
Rust
use itertools::{max, all};
use gcd::Gcd;
use divisors::get_divisors;
fn zsigmondy(a: u64, b: u64, n: u64) -> u64 {
assert!(a > b);
let dexpms: Vec<u64> = (1..(n as u32)).map(|i| a.pow(i) - b.pow(i)).collect();
let dexpn: u64 = a.pow(n as u32) - b.pow(n as u32);
let mut divisors = get_divisors(dexpn).to_vec(); // get_divisors(n) does not include 1 and n
divisors.append(&mut [1, dexpn].to_vec()); // so add those
let z = divisors.iter().filter(|d| all(dexpms.clone(), |k| Gcd::gcd(k, **d) == 1));
return *max(z).unwrap();
}
fn main() {
for (name, a, b) in [
("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,),] {
println!("\n{name}:");
for n in 1..21 {
print!("{} ", zsigmondy(a, b, n));
}
}
}
- 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
SETL
program zsigmondy;
pairs := [[2, 1], [3, 1], [4, 1], [5, 1], [6, 1], [7, 1],
[3, 2], [5, 3], [7, 3], [7, 5]];
loop for [a, b] in pairs do
print("Zsigmondy(n, " + str a + ", " + str b + ")");
loop for n in [1..18] do
putchar(str zs(n, a, b) + " ");
end loop;
print;
end loop;
proc zs(n,a,b);
dn := a**n - b**n;
divs := divisors(dn);
if #divs = 1 then return dn; end if;
loop for m in [1..n-1] do
dm := a**m - b**m;
divs -:= {d : d in divs | gcd(dm, d) > 1};
end loop;
return max/divs;
end proc;
proc divisors(n);
ds := {d : d in {1..floor(sqrt(n))} | n mod d=0};
return ds + {n div d : d in ds};
end proc;
proc gcd(a,b);
loop while b/=0 do
[a, b] := [b, a mod b];
end loop;
return a;
end proc;
end program;- Output:
Zsigmondy(n, 2, 1) 1 3 7 5 31 1 127 17 73 11 2047 13 8191 43 151 257 131071 19 Zsigmondy(n, 3, 1) 2 1 13 5 121 7 1093 41 757 61 88573 73 797161 547 4561 3281 64570081 703 Zsigmondy(n, 4, 1) 3 5 7 17 341 13 5461 257 1387 41 1398101 241 22369621 3277 49981 65537 5726623061 4033 Zsigmondy(n, 5, 1) 4 3 31 13 781 7 19531 313 15751 521 12207031 601 305175781 13021 315121 195313 190734863281 5167 Zsigmondy(n, 6, 1) 5 7 43 37 311 31 55987 1297 46873 1111 72559411 1261 2612138803 5713 1406371 1679617 3385331888947 46441 Zsigmondy(n, 7, 1) 6 1 19 25 2801 43 137257 1201 39331 2101 329554457 2353 16148168401 102943 4956001 2882401 38771752331201 117307 Zsigmondy(n, 3, 2) 1 5 19 13 211 7 2059 97 1009 11 175099 61 1586131 463 3571 6817 129009091 577 Zsigmondy(n, 5, 3) 2 1 49 17 1441 19 37969 353 19729 421 24325489 481 609554401 10039 216001 198593 381405156481 12979 Zsigmondy(n, 7, 3) 4 5 79 29 4141 37 205339 1241 127639 341 494287399 2041 24221854021 82573 3628081 2885681 58157596211761 109117 Zsigmondy(n, 7, 5) 2 3 109 37 6841 13 372709 1513 176149 1661 964249309 1801 47834153641 75139 3162961 3077713 115933787267041 30133
Sidef
func zsigmondy (a, b, c) {
var aexp = (0..c -> map {|k| a**k })
var bexp = (0..c -> map {|k| b**k })
var z = []
for n in (1 .. c) {
divisors(aexp[n] - bexp[n]).flip.each {|d|
1 ..^ n -> all {|k|
aexp[k] - bexp[k] -> is_coprime(d)
} && (z << d; break)
}
}
return z
}
for name,a,b in ([
'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,
].slices(3)) {
say ("\n#{name}:\n", 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]
Wren
Normal arithmetic
Shows the first 20 terms for each radix set apart from [7, 1], [7, 3] 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.
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 18 terms: 4 5 79 29 4141 37 205339 1241 127639 341 494287399 2041 24221854021 82573 3628081 2885681 58157596211761 109117 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
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