Coordinate-symmetric convex polytopes with equal Erhart (quasi)-polynomials Recall that given a nondegenerate polytope $P \subset \mathbb{R}^n$ which is the convex set of some vectors with integral coordinates, the Erhart polynomial $p_P(t)$ a polynomial such that $p_P(t)$ for natural numbers $t$ is the number of integral lattice points in the $t$-fold dilation $tP$ of $P$. 
Do there exist two polytopes $P, P'$ which have equal Erhart polynomials, are preserved under Euclidean reflections along all coordinate hyperplanes $\{x_i = 0\} \subset \mathbb{R}^n$, and are not $GL(n, \mathbb{Z})$ equivalent, i.e. there is no integral matrix $A$ with $det(A) = \pm 1$ such that $AP = P'$?
In fact I would be happy if there were two polytopes which are convex hulls of rational vectors and are symmetric under reflection across all coordinate hyperplanes, which are not $GL(n, \mathbb{Z})$ equivalent but have the same Erhart quasi-polynomials (these are certain functions characterized in the same way as the Erhart polynomials using the same definition of their values on the positive natural numbers; they just cease to be polynomial functions in the rational case.) 
 A: How about the following two polygons, defined by their convex hulls:

{-6, -4}, {-6, 4}, {0, -8}, {0, 8}, {6, -4}, {6, 4}

and 

{-6, -2}, {-6, 2}, {-2, -8}, {-2, 8}, {2, -8}, {2, 8}, {6,-2}, {6, 2}


Both of these have (if my calculations are correct),
the Ehrhart polynomial $1 + 6 t + 36 t^2$.
Since the number of vertices are different,
well, there is no matrix taking the first to the second.
There are more examples of this type, the two sets of points

{{-8,-2},{-8,2},{-4,-8},{-4,8},{4,-8},{4,8},{8,-2},{8,2}}}

and

{{-8,-2},{-8,2},{-6,-6},{-6,6},{-2,-8},{-2,8},{0,-8},{0,8},{2,-8},{2,8},{6,-6},{6,6},{8,-2},{8,2}}


both have Ehrhart polynomial $1 + 8 t + 52 t^2$.
I used some  Mathematica code to randomly generate some points in a symmetric fashion, then compute convex hull.
Finally, the Ehrhart polynomial is just a matter of counting lattice points in dilations, followed by an interpolation. Then one can just look at the data.
RandomPolytope[] := Module[{pts, sympts, grid, f, region},
   pts = 1 + RandomInteger[6, {5, 2}];
   grid = Join @@ Table[{x, y}, {x, -30, 30}, {y, -30, 30}];
   sympts = 
    Union[Join @@ 
      Join[{{#1, #2}, {-#1, #2}, {#1, -#2}, {-#1, -#2}} & @@@ pts]];
   poly = Table[
     region = BoundaryMeshRegion[ConvexHullMesh[k sympts]];
     Length@Select[grid, RegionMember[region, #] &]
     , {k, 1, 2}];
   (Round@MeshCoordinates@region) -> 
    Expand[InterpolatingPolynomial[Prepend[poly, 1], k] /. 
      k -> k + 1]
   ];
data = Table[RandomPolytope[], {150}];

