login
Number of coincidence site modules of index 10n+1 with 10-fold symmetry in plane.
1

%I #21 Apr 05 2020 09:33:56

%S 1,4,0,4,4,0,4,4,0,0,4,0,8,4,0,4,0,0,4,4,0,4,0,0,4,4,0,4,4,0,0,4,0,4,

%T 16,0,0,0,0,0,4,0,4,4,0,16,4,0,0,4,0,0,4,0,4,0,0,4,0,0,4,0,0,4,4,0,4,

%U 16,0,4,4,0,0,0,0,4,4,0,16,0,0,4,4,0,0,0,0,0,4,0,0,4,0,0,4,0,8,4,0,4

%N Number of coincidence site modules of index 10n+1 with 10-fold symmetry in plane.

%C The Dirichlet g.f. is Sum_{n>=0} a(n+1)/(1+10n)^s. - _R. J. Mathar_, Jul 16 2010

%D M. Baake, "Solution of coincidence problem...", in R. V. Moody, ed., Math. of Long-Range Aperiodic Order, Kluwer 1997, pp. 9-44.

%D Baake, M. and P. A. B. Pleasants. "The coincidence problem for crystals and quasicrystals." Aperiodic, vol. 94, pp. 25-29. 1995.

%F Dirichlet series: Product ((1+p^(-s))/(1-p^(-s)))^2 (p=1 mod 5).

%p read("transforms") : maxOrd := 1000 :

%p ZetaNum := proc(p,nmax,f) local n ; L := [1,seq(0,n=2..p-1),f,seq(0,n=p+1..nmax)] ; end proc:

%p Zeta := proc(p,nmax,f) local L,e; L := [1,seq(0,n=2..nmax)] ; for e from 1 do if p^e > nmax then break; else L := subsop(p^e=f^e,L) ; end if; end do: L ; end proc:

%p Zetap := [1,seq(0,n=2..maxOrd)] : for e from 1 to maxOrd by 5 do if isprime(e) then ZetaNum(e,maxOrd,1) ; Zetap := DIRICHLET(Zetap,%) ; Zeta(e,maxOrd,1) ; Zetap := DIRICHLET(Zetap,%) ; end if; end do:

%p Zetap := DIRICHLET(Zetap,Zetap) ; seq( Zetap[1+10*e],e=0..(nops(Zetap)-1)/10) ;

%p # _R. J. Mathar_, Jul 16 2010

%t did[m_, n_] := If[Mod[m, n] == 0, 1, 0];

%t DIRICHLET[a_List, b_List] := Module[{c = {}, i, s, d}, For[i = 1, i <= Min[Length[a], Length[b]], i++, s = 0; For[d = 1, d <= i, d++, If[did[i, d] == 1, s = s + a[[d]] b[[i/d]]]]; c = Append[c, s]]; c];

%t maxOrd = 1000;

%t zetaNum [p_, nmax_, f_] := Module[{n}, L = Join[{1}, Table[0, {n, 2, p - 1}], {f}, Table[0, {n, p + 1, nmax}]]];

%t zeta[p_, nmax_, f_] := Module[{L, e}, L = Join[{1}, Table[0, {n, 2, nmax}]]; For[e = 1, True, e++, If [p^e > nmax, Break[], L = ReplacePart[L, p^e -> f^e]]]; L];

%t zetap = Join[{1}, Table[0, {n, 2, maxOrd}]]; For[e = 1, e <= maxOrd, e += 5, If[PrimeQ[e], ze = zetaNum[e, maxOrd, 1];

%t zetap = DIRICHLET[zetap, ze]; ze = zeta[e, maxOrd, 1];

%t zetap = DIRICHLET[zetap, ze]]];

%t zetap = DIRICHLET[zetap, zetap];

%t Table[zetap[[1 + 10 e]], {e, 0, (Length[zetap] - 1)/10}] (* _Jean-François Alcover_, Apr 05 2020, after _R. J. Mathar_ *)

%o (PARI)

%o M=1200

%o t4=direuler(p=2,M,(1+(p%5<2)*X)) \\ p == 0 or 1 mod 5

%o t5=direuler(p=2,M,1/(1+(p%5<1)*X)) \\ p == 0 mod 5

%o t6=dirmul(t4,t5) \\ p == 1 mod 5

%o t7=direuler(p=2,M,1/(1-(p%5<2)*X))

%o t8=direuler(p=2,M,(1-(p%5<1)*X))

%o t9=dirmul(t7,t8)

%o t10=dirmul(t6,t9)

%o t10b=dirmul(t10,t10)

%o t11=vector(40,n,t10b[10*n+1]) \\ (and then prepend 1)

%o \\ _N. J. A. Sloane_, Nov 15 2019

%K nonn,easy

%O 1,2

%A _N. J. A. Sloane_

%E More terms from _R. J. Mathar_, Jul 16 2010