OFFSET
1,3
COMMENTS
a(n) = 2*A060594(n) for n = 5, 10, 13, 15, 16, 17, 20, 25, 26, 29, .... This subsequence, which contains all the primes of form 4k+1, seems to be asymptotic to 2n.
Shadow transform of A123865. - Michel Marcus, Jun 06 2013
LINKS
T. D. Noe, Table of n, a(n) for n = 1..1000
Steven R. Finch, Quartic and Octic Characters Modulo n, arXiv:0907.4894 [math.NT], 2009.
Steven Finch, Greg Martin and Pascal Sebah, Roots of unity and nullity modulo n, Proc. Amer. Math. Soc., Vol. 138, No. 8 (2010), pp. 2729-2743.
Lorenz Halbeisen and Norbert Hungerbühler, Number theoretic aspects of a combinatorial function, Notes on Number Theory and Discrete Mathematics 5(4) (1999), 138-150. (ps, pdf); see Definition 7 for the shadow transform.
Lorenz Halbeisen and Norbert Hungerbühler, Number theoretic aspects of a combinatorial function, Notes on Number Theory and Discrete Mathematics 5(4) (1999), 138-150; see Definition 7 for the shadow transform.
N. J. A. Sloane, Transforms.
FORMULA
Sum_{k=1..n} a(k) seems to be asymptotic to C*n*Log(n) with C>1.4...(when Sum_{k=1..n} A060594(k) is asymptotic to C/2*n*Log(n)).
Multiplicative with a(p^e) = p^min(e-1, 3) if p = 2, 4 if p == 1 (mod 4), 2 if p == 3 (mod 4). - David W. Wilson, Jun 09 2005
In fact, Sum_{k=1..n} a(k) is asymptotic to c*n*log(n)^2 where 2*c=0.190876.... - Steven Finch, Aug 12 2009 [c = 7/(2*Pi^3) * Product_{p prime == 1 (mod 4)} (1 - 4/(p+1)^2) = 0.0954383605... (Finch et al., 2010). - Amiram Eldar, Mar 26 2021]
MAPLE
a:= n-> add(`if`(irem(j^4-1, n)=0, 1, 0), j=0..n-1):
seq(a(n), n=1..120); # Alois P. Heinz, Jun 06 2013
# alternative
A073103 := proc(n)
local a, pf, p, r;
a := 1 ;
for pf in ifactors(n)[2] do
p := op(1, pf);
r := op(2, pf);
if p = 2 then
a := a*p^min(r-1, 3) ;
else
if modp(p, 4) = 1 then
a := 4*a ;
else
a := 2*a ;
end if;
end if;
end do:
a ;
end proc: # R. J. Mathar, Mar 02 2015
MATHEMATICA
a[n_] := Sum[If[Mod[j^4-1, n] == 0, 1, 0], {j, 0, n-1}]; Table[a[n], {n, 1, 120}] (* Jean-François Alcover, Jun 12 2015, after Alois P. Heinz *)
f[2, e_] := 2^Min[e-1, 3]; f[p_, e_] := If[Mod[p, 4] == 1, 4, 2]; a[1] = 1; a[n_] := Times @@ f @@@ FactorInteger[n]; Array[a, 100] (* Amiram Eldar, Sep 17 2020 *)
PROG
(PARI) a(n)=sum(i=1, n, if((i^4-1)%n, 0, 1))
(PARI) a(n)=my(f=factor(n)); prod(i=1, #f~, if(f[i, 1]==2, 2^min(f[i, 2]-1, 3), if(f[i, 1]%4==1, 4, 2))) \\ Charles R Greathouse IV, Mar 02 2015
(Python)
from math import prod
from sympy import primefactors
def A073103(n): return (1<<min(s-1, 3) if (s:=(~n & n-1).bit_length()) else 1)*prod(1<<(5-(p&3)>>1) for p in primefactors(n>>s)) # Chai Wah Wu, Oct 26 2022
CROSSREFS
KEYWORD
easy,nonn,mult
AUTHOR
Benoit Cloitre, Aug 19 2002
STATUS
approved