%I #40 Dec 18 2022 09:11:45
%S 1,10,21,20,11,2730,1,680,1197,550,23,5460,1,290,7161,1360,1,5757570,
%T 1,45100,6321,230,47,185640,11,530,3591,580,59,283933650,1,2720,32361,
%U 10,781,840605220,1,10,1659,1533400,83,23830170,1,40940,408177,470,1,36014160,1,277750,2163,1060,107,1882725390
%N Let (n)_p denote the exponent of prime p in the prime power factorization of n. Then a(n) is defined by the formulas a(1)=1; for n >= 2, (a(n))_2 = (n)_2, (a(n))_3 = (n)_3 and, for p >= 5, (a(n))_p = 1 + ((2n)/(p-1))_p if p-1|2*n, and (a(n))_p = 0 otherwise.
%C a(n)=1 iff n has form 6n+-1 and, if d >= 5 is a divisor of n, then 2*d+1 is not prime. The places of 1's form sequence A045979.
%C If p is an odd prime and p^n is the side length of the odd leg of a primitive Pythagorean triangle (PPT) it constrains the other leg and hypotenuse to be (p^(2n)-1)/2 and (p^(2n)+1)/2 and the area to be (p^n-1)p^n(p^n+1)/4. Now consider the term (p^n-1)p^n(p^n+1): it must at least be divisible by 24 for all odd primes p because the area of a PPT is divisible by 6 (see A127922 for n=1). a(n) equals the common divisor of the term (p^n-1)p^n(p^n+1)/24 for all odd primes p. - _Frank M Jackson_, Dec 09 2017
%H Frank M Jackson, <a href="/A202318/b202318.txt">Table of n, a(n) for n = 1..10000</a>
%F a(n) = (1/24)*b(2n+1)/b(2n-1), where b(n) = A053657(n).
%F a(p) = A002445(p)/6, for prime p >= 5.
%F a(n) = numerator of e^(real(lim_{s -> 1} (zeta(s)*(zeta(-1)^(s-1) - zeta(-(2*n-1))^(s-1))))). - _Mats Granvik_, Feb 05 2016
%F a(n) = A036283(n)/6. - _Hugo Pfoertner_, Dec 18 2022
%e Let n=6. Since 2*6+1=13 is prime, the max p that should be considered is 13. We have
%e (a(6))_2 = (a(6))_3 = 1,
%e (a(6))_5 = (12/4)_5 + 1 = 1,
%e (a(6))_7 = (12/6)_7 + 1 = 1,
%e (a(6))_13 = (12/12)_13 + 1 = 1.
%e Thus a(6) = 2*3*5*7*13 = 2730.
%t Table[Numerator[Exp[Re[Limit[Zeta[s] (Zeta[-1]^(s - 1) - Zeta[-(2*n - 1)]^(s - 1)), s -> 1]]]], {n, 1, 54}] (* _Mats Granvik_, Feb 05 2016 *)
%t Table[(lst=Table[p=Prime[m+1]; (p^n-1)p^n(p^n+1), {m, 1, 10}]; GCD@@lst/24), {n, 1, 100}] (* _Frank M Jackson_, Dec 09 2017 *)
%t a[n_] := Product[p^Sum[Floor[(n-1)/((p-1) p^k)], {k, 0, n}], {p, Prime[Range[n]]}]; Array[a[2#+1]/(24 a[2#-1]) &, 100] (* using _Jean-François Alcover_'s program A053657 *)(* _Frank M Jackson_, Dec 16 2017 *)
%o (PARI) a(n) = {my(r = 1); forprime(p=2, 2*n+1, if (p<=3, r *= p^valuation(n, p), if (! (2*n % (p-1)), r *= p^(1+valuation((2*n)/(p-1), p))););); r;} \\ _Michel Marcus_, Feb 06 2016
%Y Cf. A053657, A045979, A002445, A127922.
%K nonn
%O 1,2
%A _Vladimir Shevelev_ and _Peter J. C. Moses_, Dec 16 2011