|
%I #39 Mar 31 2018 14:48:03
%S 0,1,1,4,4,4,9,5,5,9,16,12,12,12,16,25,17,13,13,17,25,36,28,20,24,20,
%T 28,36,49,37,33,25,25,33,37,49,64,52,40,36,40,36,40,52,64,81,65,53,45,
%U 41,41,45,53,65,81,100,84,72,64,52,60,52,64,72,84,100
%N Square array a(p,q) = p^2 + q^2 - 2*p - 2*q + 2*gcd(p,q), p >= 1, q >= 1, read by antidiagonals.
%C In the Cartesian plane, let r(p,q) denote the rotation with center origin and angle associated to slope p/q (p: number of units upwards, q: number of units towards the right).
%C Let R(p,q) be the square of area p^2 + q^2 = R^2, with vertices (0,0), (0,R), (R,R), (R,0).
%C The natural unit squares (i.e., the (a,a+1) X (b,b+1) Cartesian products, with a and b integers) are transformed by r(p,q) into rotated unit squares.
%C a(p,q) is the number of rotated unit squares that fully land inside R(p,q).
%H Luc Rousseau, <a href="/A288797/a288797_example2.png">Illustration for p=6, q=4</a>
%F a(p,q) = p^2 + q^2 - 2*p - 2*q + 2*gcd(p,q).
%F a(p,1) = a(1,p) = (p-1)^2 = A000290(p-1).
%F a(p,p) = 2*p*(p-1) = 4*A000217(p-1).
%F a(p,p+1) = 2*p*(p-1)+1 = A001844(p-1).
%F a(p,p+2) = 2*p^2+2*gcd(2,p) = 2*p^2+3+(-1)^(p) = 4*A099392(p-1) = 4*A080827(p).
%F a(p,p+3) = 2*p^2+2*p+5+4*A079978(p) = 1+4*(1+A143101(p)).
%F a(p,2*p) = p*(5*p-4) = A051624(p).
%F a(p,3*p) = 2*p*(5*p-3) = 4*A000566(p).
%e Table begins:
%e 0 1 4 9 16 25 ...
%e 1 4 5 12 17 28 ...
%e 4 5 12 13 20 33 ...
%e 9 12 13 24 25 36 ...
%e 16 17 20 25 40 41 ...
%e 25 28 33 36 41 60 ...
%e ... ... ... ... ... ... ...
%t A[p_, q_] := p^2 + q^2 - 2*p - 2*q + 2*GCD[p, q];
%t (* or, checking without the formula: *)
%t okQ[{a_, b_}, p_, q_] := Module[{r2 = p^2 + q^2}, 0 <= a*q - b*p <= r2 && 0 <= a*p + b*q <= r2 && 0 <= a*q - b*p + q <= r2 && 0 <= a*p + b*q + p <= r2 && 0 <= a*q - (b + 1)*p <= r2 && 0 <= a*p + b*q + q <= r2 && 0 <= (a + 1)*q - (b + 1)*p <= r2 && 0 <= a*p + b*q + p + q <= r2];
%t A[p_, q_] := Module[{r}, r = Reduce[okQ[{a, b}, p, q], {a, b}, Integers]; If[r[[0]] === And, 1, Length[r]]];
%t Flatten[Table[A[p - q + 1, q], {p, 1, 11}, {q, 1, p}]] (* _Jean-François Alcover_, Jun 17 2017 *)
%o (PARI)
%o a(p,q)=p^2+q^2-2*p-2*q+2*gcd(p,q)
%o for(n=2,12,for(p=1,n-1,{q=n-p;print(a(p,q))}))
%Y Cf. A000217, A000290, A000566, A001844, A003989, A051624, A080827, A099392, A143101.
%K tabl,nonn,easy
%O 1,4
%A _Luc Rousseau_, Jun 16 2017
|
