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;
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 );
INT step = IF ODD an minus bn THEN 2 ELSE 1 FI;
WHILE ( d +:= step ) <= 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}
Arturo
zsigmondy: function [a b num][
map num 'n ->
first select.first reverse factors (a^n)-(b^n) 'x ->
every? n-1 'm ->
one? gcd @[x (a^m)-(b^m)]
]
n: 20
loop [[2 1][3 1][4 1][5 1][6 1][7 1][3 2][5 3][7 3][7 5]] 'pair [
[a,b]: pair
print ["Zsigmondy(" n "," a "," b "):" join.with:", " zsigmondy a b n "\n"]
]
- 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
C
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
uint64_t ipow(uint64_t n, uint64_t p) {
if (n == 1 || p == 0) return 1;
uint64_t x = ipow(n, p>>1);
return p&1 ? x*x*n : x*x;
}
uint64_t gcd(uint64_t a, uint64_t b) {
return b ? gcd(b, a%b) : a;
}
bool all_coprime(uint64_t a, uint64_t b, uint64_t d, uint64_t n) {
for (uint64_t m = 1; m < n; m++) {
uint64_t dm = ipow(a,m)-ipow(b,m);
if (gcd(dm, d) != 1) return false;
}
return true;
}
uint64_t zsigmondy(uint64_t n, uint64_t a, uint64_t b) {
uint64_t dn = ipow(a,n) - ipow(b,n);
uint64_t maxdiv = 0;
for (uint64_t d = 1; d*d <= dn; d++) {
if (dn % d != 0) continue;
if (all_coprime(a, b, d, n))
maxdiv = d > maxdiv ? d : maxdiv;
uint64_t dnd = dn/d;
if (all_coprime(a, b, dnd, n))
maxdiv = dnd > maxdiv ? dnd : maxdiv;
};
return maxdiv;
}
void zsig_row(uint64_t a, uint64_t b) {
printf("zsigmondy(n, %lu, %lu):\n", a, b);
for (uint64_t n = 1; n <= 18; n++) {
printf("%lu ", zsigmondy(n, a, b));
}
printf("\n");
}
int main(void) {
uint64_t pairs[][2] = {
{2, 1}, {3, 1}, {4, 1}, {5, 1}, {6, 1}, {7, 1},
{3, 2}, {5, 3}, {7, 3}, {7, 5}
};
for (size_t pair=0; pair<sizeof(pairs)/sizeof(*pairs); pair++) {
zsig_row(pairs[pair][0], pairs[pair][1]);
}
return 0;
}
- 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
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
C#
internal class Program
{
private static void Main(string[] args)
{
//record Pair(int a, int b) { }
List<KeyValuePair<int, int>> pairs =
[
KeyValuePair.Create(2, 1),
KeyValuePair.Create(3, 1),
KeyValuePair.Create(4, 1),
KeyValuePair.Create(5, 1),
KeyValuePair.Create(6, 1),
KeyValuePair.Create(7, 1),
KeyValuePair.Create(3, 2),
KeyValuePair.Create(5, 3),
KeyValuePair.Create(7, 3),
KeyValuePair.Create(7, 5),
];
foreach (KeyValuePair<int, int> pair in pairs)
{
Console.WriteLine("Zsigmondy(n, " + pair.Key + ", " + pair.Value + ")");
for (int n = 1; n <= 18; n++)
{
Console.Write(ZsigmondyNumber(n, pair.Key, pair.Value) + " ");
}
Console.WriteLine();
}
}
private static long ZsigmondyNumber(int n, int a, int b)
{
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.RemoveAll(d => GCD(dm, d) > 1);
}
return divisors[divisors.Count() - 1];
}
private static List<long> Divisors(long number)
{
List<long> result = new List<long>();
for (long d = 1; d * d <= number; d++)
{
if (number % d == 0)
{
result.Add(d);
if (d * d < number)
{
result.Add(number / d);
}
}
}
result.Sort();
return result;
}
private static bool 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
CLU
isqrt = proc (n: int) returns (int)
x0: int := n/2
if x0=0 then return(n) end
x1: int := (x0 + n/x0)/2
while x1 < x0 do
x0, x1 := x1, (x1 + n/x1)/2
end
return(x0)
end isqrt
divisors = iter (n: int) yields (int)
for d: int in int$from_to(1, isqrt(n)) do
if n//d ~= 0 then continue end
yield(d)
if n/d ~= d then yield(n/d) end
end
end divisors
gcd = proc (a, b: int) returns (int)
if b=0 then return(a) else return(gcd(b, a//b)) end
end gcd
zsigmondy = proc (n, a, b: int) returns (int)
dn: int := a**n - b**n
dms: array[int] := array[int]$predict(1, n-1)
for m: int in int$from_to(1, n-1) do
array[int]$addh(dms, a**m - b**m)
end
maxdiv: int := 0
for d: int in divisors(dn) do
begin
for m: int in int$from_to(1, n-1) do
if gcd(dms[m], d) ~= 1 then exit not_coprime end
end
if d > maxdiv then maxdiv := d end
end except when not_coprime: end
end
return(maxdiv)
end zsigmondy
zsig_row = proc (po: stream, a, b: int)
stream$putl(po, "zsigmondy(n, " || int$unparse(a) || ", " || int$unparse(b) || "):")
for n: int in int$from_to(1, 18) do
stream$puts(po, int$unparse(zsigmondy(n, a, b)) || " ")
end
stream$putl(po, "")
end zsig_row
start_up = proc ()
pair = struct[a, b: int]
po: stream := stream$primary_output()
pairs: sequence[pair] := sequence[pair]$[
pair${a:2, b:1}, pair${a:3, b:1}, pair${a:4, b:1}, pair${a:5, b:1},
pair${a:6, b:1}, pair${a:7, b:1}, pair${a:3, b:2}, pair${a:5, b:3},
pair${a:7, b:3}, pair${a:7, b:5}
]
for p: pair in sequence[pair]$elements(pairs) do
zsig_row(po, p.a, p.b)
end
end start_up- 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
Cowgol
The biggest integer type that Cowgol supports is 32 bits, so this program prints the values for n from 1 upto 12.
include "cowgol.coh";
sub ipow(base: uint32, exp: uint32): (result: uint32) is
result := 1;
loop
if (exp & 1) == 1 then
result := result * base;
end if;
exp := exp >> 1;
if exp == 0 then
break;
end if;
base := base * base;
end loop;
end sub;
sub gcd(a: uint32, b: uint32): (result: uint32) is
while b>0 loop
result := a;
a := b;
b := result % b;
end loop;
result := a;
end sub;
sub max(a: uint32, b: uint32): (result: uint32) is
if a>b then
result := a;
else
result := b;
end if;
end sub;
sub all_coprime(a: uint32, b: uint32, d: uint32, n: uint32): (result: uint8) is
result := 1;
var m: uint32 := 1;
while m<n loop
var dm: uint32 := ipow(a,m) - ipow(b,m);
if (gcd(dm, d) != 1) then
result := 0;
break;
end if;
m := m+1;
end loop;
end sub;
sub zsigmondy(n: uint32, a: uint32, b: uint32): (maxdiv: uint32) is
var dn: uint32 := ipow(a,n) - ipow(b,n);
var d: uint32 := 1;
maxdiv := 0;
while d*d <= dn loop
if dn % d == 0 then
if all_coprime(a, b, d, n) != 0 then
maxdiv := max(maxdiv, d);
end if;
var dnd := dn/d;
if all_coprime(a, b, dnd, n) != 0 then
maxdiv := max(maxdiv, dnd);
end if;
end if;
d := d + 1;
end loop;
end sub;
sub zsig_row(a: uint32, b: uint32) is
print("zsigmondy(n, ");
print_i32(a);
print(", ");
print_i32(b);
print("):");
print_nl();
var n: uint32 := 1;
while n <= 12 loop
print_i32(zsigmondy(n, a, b));
print_char(' ');
n := n + 1;
end loop;
print_nl();
end sub;
record Pair is
a: uint32;
b: uint32;
end record;
var pairs: Pair[] := {
{2, 1}, {3, 1}, {4, 1}, {5, 1}, {6, 1}, {7, 1},
{3, 2}, {5, 3}, {7, 3}, {7, 5}
};
var pair: @indexof pairs := 0;
while pair < @sizeof pairs loop
zsig_row(pairs[pair].a, pairs[pair].b);
pair := pair + 1;
end loop;- Output:
zsigmondy(n, 2, 1): 1 3 7 5 31 1 127 17 73 11 2047 13 zsigmondy(n, 3, 1): 2 1 13 5 121 7 1093 41 757 61 88573 73 zsigmondy(n, 4, 1): 3 5 7 17 341 13 5461 257 1387 41 1398101 241 zsigmondy(n, 5, 1): 4 3 31 13 781 7 19531 313 15751 521 12207031 601 zsigmondy(n, 6, 1): 5 7 43 37 311 31 55987 1297 46873 1111 72559411 1261 zsigmondy(n, 7, 1): 6 1 19 25 2801 43 137257 1201 39331 2101 329554457 9962347 zsigmondy(n, 3, 2): 1 5 19 13 211 7 2059 97 1009 11 175099 61 zsigmondy(n, 5, 3): 2 1 49 17 1441 19 37969 353 19729 421 24325489 481 zsigmondy(n, 7, 3): 4 5 79 29 4141 37 205339 1241 127639 341 494287399 59740867 zsigmondy(n, 7, 5): 2 3 109 37 6841 13 372709 1513 176149 1661 964249309 14838431
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
FreeBASIC
Function GCD(n As Ulongint, d As Ulongint) As Ulongint
Return Iif(d = 0, n, GCD(d, n Mod d))
End Function
Sub SortInt(array() As Ulongint)
Dim As Integer i, j
Dim As Integer lb = Lbound(array), ub = Ubound(array)
For i = lb To ub - 1
For j = i + 1 To ub
If array(i) > array(j) Then Swap array(i), array(j)
Next j
Next i
End Sub
Sub divisors(n As Ulongint, result() As Ulongint)
Dim As Integer i, size
Redim result(0)
result(0) = 1
If n > 1 Then
Redim Preserve result(1)
result(1) = n
End If
For i As Ulongint = 2 To Sqr(n)
If n Mod i = 0 Then
Redim Preserve result(Ubound(result) + 1)
result(Ubound(result)) = i
If i * i <> n Then
Redim Preserve result(Ubound(result) + 1)
result(Ubound(result)) = n \ i
End If
End If
Next
SortInt(result())
End Sub
Function ipow(base_ As Ulongint, exp_ As Ulongint) As Ulongint
Dim As Ulongint result = 1
While exp_ > 0
If (exp_ And 1) Then result *= base_
base_ *= base_
exp_ Shr= 1
Wend
Return result
End Function
Function zsigmondy(n As Ulongint, a As Ulongint, b As Ulongint) As Ulongint
Dim As Ulongint p = ipow(a, n) - ipow(b, n)
Dim As Ulongint d()
divisors(p, d())
For i As Integer = Ubound(d) To 0 Step -1
Dim As Ulongint z = d(i)
Dim As Integer all_coprime = 1
For m As Ulongint = 1 To n - 1
If GCD(ipow(a, m) - ipow(b, m), z) <> 1 Then
all_coprime = 0
Exit For
End If
Next
If all_coprime Then Return z
Next
Return 1
End Function
Dim As Uinteger avs(1 To 10, 1 To 2) = {{2,1}, {3,1}, {4,1}, {5,1}, {6,1}, {7,1}, {3,2}, {5,3}, {7,3}, {7,5}}
For ab As Integer = 1 To Ubound(avs)
Dim As Uinteger a = avs(ab,1), b = avs(ab,2)
Dim As Uinteger lim = 20
Print Using "Zsigmondy(n, #_, #) - first ## terms:"; a; b; lim
Print " "; Zsigmondy(lim, a, b);
Print Chr(10)
Next ab
Sleep
- Output:
Zsigmondy(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 Zsigmondy(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 Zsigmondy(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 Zsigmondy(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 Zsigmondy(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 Zsigmondy(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 Zsigmondy(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 Zsigmondy(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 Zsigmondy(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 Zsigmondy(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
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}}
Miranda
main = [Stdout (zsig_row a b) | (a,b)<-pairs]
where pairs = [(2,1),(3,1),(4,1),(5,1),(6,1),(7,1),
(3,2),(5,3),(7,3),(7,5)]
zsig_row :: num->num->[char]
zsig_row a b = lay [title, show zsigs, ""]
where title = "zsigmondy n " ++ show a ++ " " ++ show b ++ ":"
zsigs = [zsigmondy n a b | n<-[1..18]]
zsigmondy :: num->num->num->num
zsigmondy n a b = max [d | d<-divisors dn; and [dm $gcd d = 1 | dm<-dms]]
where dn = a^n - b^n
dms = [a^m - b^m | m <- [1..n-1]]
divisors :: num->[num]
divisors n = concat [unique [d, n div d] | d<-[1..entier (sqrt n)]; n mod d=0]
gcd :: num->num->num
gcd a 0 = a
gcd a b = gcd b (a mod b)
unique :: [*]->[*]
unique = unique' []
where unique' acc [] = acc
unique' acc (a:as) = unique' acc as, if a $in acc
= unique' (a:acc) as, otherwise
in :: *->[*]->bool
in a [] = False
in a (a:as) = True
in a (b:as) = a $in as- 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]
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
Refal
$ENTRY Go {
= <Each (ZsigRow) <Pairs>>;
};
Pairs {
= (2 1) (3 1) (4 1) (5 1) (6 1) (7 1) (3 2) (5 3) (7 3) (7 5);
};
ZsigRow {
s.A s.B =
<Prout 'zsigmondy(n, ' <Symb s.A> ', ' <Symb s.B> '):'>
<Prout <Each (Zsigmondy s.A s.B) <Iota 1 18>>>
<Prout>;
};
Zsigmondy {
s.A s.B s.N,
<ZsD s.A s.B s.N>: s.DN,
<Divs s.DN>: e.Divs,
<Each (ZsD s.A s.B) <Iota 1 <Sub s.N 1>>>: e.DMs,
<FilterAllCoprime (e.Divs) e.DMs>: e.1 s.MaxDiv
= s.MaxDiv;
};
ZsD {
s.A s.B s.N = <Sub <Pow s.A s.N> <Pow s.B s.N>>;
};
FilterAllCoprime {
() e.DMs = ;
(s.D e.Divs) e.DMs, <AllCoprime s.D e.DMs>: {
True = s.D <FilterAllCoprime (e.Divs) e.DMs>;
False = <FilterAllCoprime (e.Divs) e.DMs>;
};
};
AllCoprime {
s.D = True;
s.D s.DM e.DMs, <Gcd s.D s.DM>: {
1 = <AllCoprime s.D e.DMs>;
s.N = False;
};
};
Isqrt {
s.N s.X0,
<Div <Add s.X0 <Div s.N s.X0>> 2>: s.X1,
<Compare s.X1 s.X0>: {
'-' = <Isqrt s.N s.X1>;
s.C = s.X0;
};
s.N, <Div s.N 2>: {
0 = s.N;
s.X0 = <Isqrt s.N s.X0>;
};
};
Iota {
s.E s.E = s.E;
s.S s.E, <Compare s.S s.E>: '+' = ;
s.S s.E = s.S <Iota <Add 1 s.S> s.E>;
};
Each {
(e.F) = ;
(e.F) s.X e.Xs = <Mu e.F s.X> <Each (e.F) e.Xs>;
(e.F) (e.X) e.Xs = <Mu e.F e.X> <Each (e.F) e.Xs>;
};
Gcd {
s.A 0 = s.A;
s.A s.B = <Gcd s.B <Mod s.A s.B>>;
};
Divs {
s.N s.D s.Max, <Compare s.D s.Max>: '+' = ;
s.N s.D s.Max,
<Mod s.N s.D>: {
0, <Div s.N s.D>: {
s.D = s.D <Divs s.N <Add 1 s.D> s.Max>;
s.D2 = s.D <Divs s.N <Add 1 s.D> s.Max> s.D2;
};
s.M = <Divs s.N <Add 1 s.D> s.Max>;
};
s.N = <Divs s.N 1 <Isqrt s.N>>;
};
Pow {
s.N 0 = 1;
s.N s.P, <Divmod s.P 2>: {
(s.P2) 0, <Pow s.N s.P2>: s.X = <Mul s.X s.X>;
(s.P2) 1, <Pow s.N s.P2>: s.X = <Mul s.N <Mul s.X s.X>>;
};
};- 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
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