|
%I #59 Sep 27 2023 13:41:31
%S 0,1,9,36,100,229,473,910,1648,2795,4469,6818,10032,14315,19907,27190,
%T 36502,48233,62803,80736,102550,128847,160271,197516,241314,292737,
%U 352591,421764,501204,592257,696281,814450,948112,1098607,1267367
%N Number of cubes that can be formed from the points of a cubical grid of n X n X n points.
%C Skew cubes are allowed.
%H Baitian Li, <a href="/A098928/b098928.txt">Table of n, a(n) for n = 1..10000</a> (terms 1..101 from E. J. Ionascu and R. A. Obando)
%H E. J. Ionascu and R. A. Obando, <a href="http://arxiv.org/abs/1003.4569">Counting all cubes in {0,1,...,n}^3</a>, arXiv:1003.4569 [math.NT], 2010.
%H Eugen J. Ionascu and Andrei Markov, <a href="http://dx.doi.org/10.1016/j.jnt.2010.07.008">Platonic solids in Z^3</a>, Journal of Number Theory, Volume 131, Issue 1, January 2011, Pages 138-145.
%H Eugen J. Ionascu and R. A. Obando, <a href="http://www.emis.de/journals/INTEGERS/papers/a9self/a9self.Abstract.html">Cubes in {0,1,...,N}^3</a>, INTEGERS, 12A (2012), #A9. - From _N. J. A. Sloane_, Feb 05 2013
%H I. Larrosa, <a href="http://faculty.missouristate.edu/l/lesreid/POW08_03.html">SMSU Problem Corner</a>.
%H Baitian Li, <a href="/A098928/a098928.txt">C++ program for A098928</a>
%e For n = 3 there are 8 cubes of volume 1 and 1 cube of volume 8; thus a(3)=9. - _José María Grau Ribas_, Mar 15 2014
%e a(6)=229 because we can place 15^2 cubes in a 6 X 6 X 6 cubical grid with their edges parallel to the faces of the grid, plus 4 cubes of edge 3 with a vertex in each face of the lattice and the other two vertices on a diagonal.
%t Needs["Quaternions`"];
%t (* Initialize variables *)
%t R = 20;
%t NN = 1010;
%t (* Quaternion operations *)
%t test[q_Quaternion] :=
%t Module[{unit, res, a, b, c, u, v, w, p},
%t If[Round[Norm[q]] > R, Return[]];
%t If[q == Quaternion[0, 0, 0, 0], Return[]];
%t unit = Quaternion[0, 1, 0, 0];
%t res = q ** unit ** Conjugate[q];
%t a = Abs[res[[2]]] + Abs[res[[3]]] + Abs[res[[4]]];
%t unit = Quaternion[0, 0, 1, 0];
%t res = q ** unit ** Conjugate[q];
%t b = Abs[res[[2]]] + Abs[res[[3]]] + Abs[res[[4]]];
%t unit = Quaternion[0, 0, 0, 1];
%t res = q ** unit ** Conjugate[q];
%t c = Abs[res[[2]]] + Abs[res[[3]]] + Abs[res[[4]]];
%t For[i = 1, i <= (R - 1)/Max[a, b, c], i++,
%t If[SquareFreeQ[i], {u = a*i;
%t v = b*i;
%t w = c*i;
%t p = Max[u, v, w] + 1;
%t coe[[p + 1, 4]] += (1);
%t coe[[p + 1, 3]] -= (u + v + w);
%t coe[[p + 1, 2]] += (u*v + v*w + w*u);
%t coe[[p + 1, 1]] -= (u*v*w)}]]];
%t (* Set up coefficient matrix *)
%t coe = ConstantArray[0, {NN, 4}];
%t (* Loop through quaternions *)
%t rt = Ceiling[Sqrt[R]] + 1;
%t For[s = -rt, s <= rt, s++,
%t For[x = -rt, x <= rt, x++,
%t For[y = -rt, y <= rt, y++,
%t For[z = -rt, z <= rt, z++, test[Quaternion[s, x, y, z]];
%t test[Quaternion[s + 0.5, x + 0.5, y + 0.5, z + 0.5]];]]]];
%t newCoe = coe;
%t newCoe[[2 ;; ;; 2]] = coe[[2 ;; ;; 2]]/2;
%t (* Calculate and output results *)
%t For[i = 2, i <= R + 1, i++, ans = 0;
%t For[j = 4, j >= 1, j--, newCoe[[i, j]] += newCoe[[i - 1, j]];
%t ans = ans*(i - 1) + newCoe[[i, j]];
%t ];
%t Print[i - 1, " ", ans/24];];
%t (* _Haomin Yang_, Aug 29 2023 *)
%o (C++) // see link above
%Y Cf. A103158.
%Y Cf. A000537 (without skew cubes), A002415 (number of squares with corners on an n X n grid), A108279, A102698.
%K nonn
%O 1,3
%A _Ignacio Larrosa Cañestro_, Oct 19 2004, Sep 29 2009
%E Edited by _Ray Chandler_, Apr 05 2010
%E Further edited by _N. J. A. Sloane_, Mar 31 2016
| 