login
A170826
a(n) = gcd(n^2, n!).
5
1, 2, 3, 8, 5, 36, 7, 64, 81, 100, 11, 144, 13, 196, 225, 256, 17, 324, 19, 400, 441, 484, 23, 576, 625, 676, 729, 784, 29, 900, 31, 1024, 1089, 1156, 1225, 1296, 37, 1444, 1521, 1600, 41, 1764, 43, 1936, 2025, 2116, 47, 2304, 2401, 2500, 2601, 2704, 53, 2916
OFFSET
1,2
COMMENTS
From Marco Ripà, Apr 17 2026: (Start)
For each prime number p, by Wilson's theorem, there exists (at least) an integer n >= 1 such that gcd(n^3+1, n!+1) = p (and ditto for gcd(n^2-1, n!+1)), whereas the only two integers 1 < n < 2500000 such that gcd(n^3-1, n!-1) <> 1 are 2283 and 138645 (it is worth noting that gcd(2283^3-1, 2283!-1) = gcd(138645^3-1, 138645!-1) = 140929 = 2283 + 138645 + 1 since, modulo 140929, 2283 and 138645 are the only nontrivial cube roots of unity so that 2283^3 == 1 (mod 140929) and also 138645^3 == 1 (mod 140929)).
On the other hand, there does not exist any integer n > 1 such that gcd(n^2-1, n-1!) <> 1 (see MathOverflow answer #510274), while for all integers 1 < n < 10^8, it has been verified that gcd(a(n), A002522(n)) = gcd(n!+1, n^2+1) = 1 (as reported in comments to MathOverflow question #509651). Furthermore, if gcd(n!+1, n^2+1) > 1 for some integer n > 1, then it is a prime greater than n that divides n^2+1 (see the mentioned MathOverflow thread, in Links). (End)
FORMULA
If n is prime then a(n) = n; otherwise, if n <> 4 then a(n) = n^2. - Zak Seidov, Dec 28 2009
a(n) = n!/A092043(n). - Johannes W. Meijer, Jun 04 2016
a(n) = n^2 / n^c(n), where c = A010051 for n >= 5. - Wesley Ivan Hurt, Nov 10 2023
MAPLE
GCDWITHFACTORIAL:=proc(a) local b, i, k:
if whattype(a) <> list then RETURN([]); fi:
b:=[]:
for i to nops(a) do b:=[op(b), gcd(a[i], i!)]: od;
RETURN(b);
end:
A170826 := proc(n): gcd(n^2, n!) end: seq(A170826(n), n=1..54); # Johannes W. Meijer, Jun 04 2016
MATHEMATICA
Table[GCD[n^2, n!], {n, 54}] (* Michael De Vlieger, Jun 05 2016 *)
PROG
(PARI) a(n)=if(isprime(n), n, if(n==4, 8, n^2)) \\ Charles R Greathouse IV, Feb 01 2013
(Python)
from math import prod
from sympy import factorint
from sympy.ntheory.factor_ import digits
def A170826(n): return prod(p**min(e<<1, (n-sum(digits(n, p)[1:]))//(p-1)) for p, e, in factorint(n).items()) # Chai Wah Wu, May 05 2026
CROSSREFS
KEYWORD
nonn,easy
AUTHOR
N. J. A. Sloane, Dec 27 2009
STATUS
approved