# Swinging Primes
# A Maple code companion to
# http://www.luschny.de/math/primes/SwingingPrimes.html
# ================================
# Peter Luschny, 2009/07/01
# Licensed under a Creative
# Commons Attribution 3.0 License
#================================

swing := proc(n) option remember;
if n = 0 then 1 elif irem(n, 2) = 1
then swing(n-1)*n else 4*swing(n-1)/n fi end:

A74 := proc(f,n)
select(isprime, map(x -> f(x)+1,[$1..n]));
select(isprime, map(x -> f(x)-1,[$1..n]));
sort(convert(convert(%%,set) union convert(%,set),list)) end:

A75 := proc(f,n)
select(isprime, map(x -> f(x)+1,[$1..n])) end:

A76 := proc(f,n)
select(isprime, map(x ->f(x)-1,[$0..n]));
sort(%) end:

A77 := proc(f,n)
select(x -> isprime(f(x)+1),[$0..n]) end:

A78 := proc(f,n)
select(x -> isprime(f(x)-1),[$0..n]) end:

A79 := proc(f,n)
select(k -> isprime(f(k)+1),[$0..n]);
select(isprime,%); end:

A80 := proc(f,n)
select(k -> isprime(f(k)-1),[$0..n]);
select(isprime,%); end:

A81 := proc(f,n)
select(isprime,[$2..n]);
select(isprime, map(x -> f(x)+1,%)) end:

A82 := proc(f,n)
select(isprime,[$2..n]);
select(isprime, map(x -> f(x)-1,%)) end:

A83 := proc(f,n)
select(s->isprime(s) and isprime(s-2),
map(k->f(k)+1,[$4..n])) end;

A84 := proc(n)
select(isprime,[$2..n]);
select(isprime, map(x -> x! - 1,%)) end:

A85 := proc(n) local i; mul(swing(i),i=0..n) end:
A86 := proc(n) local i; mul(A85(i),i=0..n) end:

A87 := proc(n) local i;
mul(i,i=map(swing,numtheory[divisors](n))) end:

A88 := n -> A87(n)/swing(n):
A89 := n -> n!/A87(n):
A90 := n -> n!/swing(n):

A91 := proc(n) local i;
add(i,i=map(swing,numtheory[divisors](n))) end:

##################################################

A74(factorial,26);
A74(swing,46);

# [2, 3, 5, 7, 23, 719, 5039, 39916801, 479001599,
# [2, 3, 5, 7, 19, 29, 31, 71, 139, 251, 631, 3433,

A75(factorial,28);
A75(swing,68);

# [2, 3, 7, 39916801, 10888869450418352160768000001]
# [2, 3, 7, 7, 31, 71, 631, 3433, 51481, 2704157,

A76(factorial,28);
A76(swing,58);

# [5, 23, 719, 5039, 479001599, 87178291199]
# [5, 5, 19, 29, 139, 251, 12011, 48619, 51479,

A77(factorial,120);
A77(swing,100);

# [0, 1, 2, 3, 11, 27, 37, 41, 73, 77, 116]
# [0, 1, 2, 3, 4, 5, 8, 9, 14, 15, 24, 27, 31,

A78(factorial,22);
A78(swing,110);

# [3, 4, 6, 7, 12, 14]
# [3, 4, 5, 6, 7, 10, 13, 15, 18, 30, 35, 39, 41,

A79(factorial,126);
A79(swing,910);

# [2, 3, 11, 37, 41, 73]
# [2, 3, 5, 31, 67, 139, 631]

A80(factorial,100);
A80(swing,240);

# [3, 7]
# [3, 5, 7, 13, 41, 47, 83, 137, 151, 229]

A81(factorial,40);
A81(swing,88);

# [3, 7, 39916801, 13763753091226345046315979581580902400000001]
# [3, 7, 31, 4808643121, 483701705079089804581]

A82(factorial,40);
A82(swing,68);

# [5, 5039]
# [5, 29, 139, 12011, 5651707681619, 386971244197199]

A83(factorial,40);
A83(swing,88);

# [5]
# [5, 7, 31, 51481, 1580132580471901]

A84(20);

# [5, 5039]

seq(A85(i),i=0..9);

# 1, 1, 2, 12, 72, 2160, 43200, 6048000, 423360000,

seq(A86(i),i=0..9);

# 1, 1, 2, 24, 1728, 3732480, 161243136000,

seq(A87(i),i=0..12);
seq(A88(i),i=0..12);
seq(A89(i),i=0..12);
seq(A90(i),i=0..12);

# 1, 1, 2, 6, 12, 30, 240, 140, 840, 3780, 15120,
# 1, 1, 1, 1, 2, 1, 12, 1, 12, 6, 60, 1, 240
# 1, 1, 1, 1, 2, 4, 3, 36, 48, 96, 240, 14400, 
# 1, 1, 1, 1, 4, 4, 36, 36, 576, 576, 14400, 14400, 

seq(A91(i),i=0..12);

# 0, 1, 3, 7, 9, 31, 29, 141, 79, 637, 285, 2773, 

##################################################

# Wilson type sequences

WQ := proc(f,r,n) map(p->(f(p-1)+r(p))/p,
select(isprime,[$1..n])) end:

WR := proc(f,r,n) map(p->(f(p-1)+r(p))/p mod p,
select(isprime,[$1..n])) end:

WS := proc(f,r,n)
select(p->(f(p-1)+r(p)) mod p = 0,[$2..n]);
select(q -> not isprime(q),%) end:

WP := proc(f,r,n)
select(p->(f(p-1)+r(p))/p mod p = 0, 
select(isprime,[$1..n])) end:

A163209 := n -> WS(swing,p->(-1)^iquo(p+2,2),n):
A163210 := n -> WQ(swing,p->(-1)^iquo(p+2,2),n):
A163211 := n -> select(isprime,A163210(n)):
A163212 := n -> select(isprime,WQ(factorial,p->1,n)):
A163213 := n -> WR(swing,p->(-1)^iquo(p+2,2),n): 
A002068 := n -> WR(factorial,p->1,n):
A007619 := n -> WQ(factorial,p->1,n):

A163209(6000);
A163210(30);
A163211(30);
A163212(30);
A163213(30);
A002068(30);
A007619(20);

WP(factorial,p->1, 600);
WP(swing,p->(-1)^iquo(p+2,2), 3600);