|
%I #39 May 17 2020 06:18:17
%S 1,1,1,1,5,1,1,1,1,5,1,1,13,1,5,1,17,1,1,5,1,1,1,1,25,13,1,1,29,5,1,1,
%T 1,17,5,1,37,1,13,5,41,1,1,1,5,1,1,1,1,25,17,13,53,1,5,1,1,29,1,5,61,
%U 1,1,1,65,1,1,17,1,5,1,1,73,37,25,1,1,13,1,5,1,41,1,1,85,1,29,1
%N a(n) is the product of primes (with multiplicity) of form 4*k+1 that divide n.
%C Completely multiplicative with a(p) = p if p = 4k+1 and a(p) = 1 otherwise. - _Tom Edgar_, Mar 05 2015
%H Alois P. Heinz, <a href="/A170818/b170818.txt">Table of n, a(n) for n = 1..10000</a>
%H A. Tripathi, <a href="http://www.fq.math.ca/Papers1/46_47-4/Tripathi.pdf">On Pythagorean triples containing a fixed integer</a>, Fib. Q., 46/47 (2008/2009), 331-340.
%H <a href="/index/Di#divseq">Index to divisibility sequences</a>
%F a(n) = n/A072438(n). - _Michel Marcus_, Mar 05 2015
%p a:= n-> mul(`if`(irem(i[1], 4)=1, i[1]^i[2], 1), i=ifactors(n)[2]):
%p seq(a(n), n=1..100); # _Alois P. Heinz_, Jun 09 2014
%t a[n_] := Product[{p, e} = pe; If[Mod[p, 4] == 1, p^e, 1], {pe, FactorInteger[n]}];
%t Array[a, 100] (* _Jean-François Alcover_, May 29 2019 *)
%o (PARI) a(n)=my(f=factor(n)); prod(i=1,#f~,if(f[i,1]%4>1,1,f[i,1])^f[i,2]) \\ _Charles R Greathouse IV_, Jun 28 2015
%o (Python)
%o from sympy import factorint, prod
%o def a072438(n):
%o f = factorint(n)
%o return 1 if n == 1 else prod(i**f[i] for i in f if i % 4 != 1)
%o def a(n): return n//a072438(n) # _Indranil Ghosh_, May 08 2017
%Y Cf. A170817-A170819, A097706, A083025, A072438, A286361.
%K nonn,mult
%O 1,5
%A _N. J. A. Sloane_, Dec 22 2009
|
