This is closer. Since we only insist on equivalence modulo integers this example is in a sense too trivial-- one can solve it directly. I will propose the following modest alteration on the right hand sides.

sol = {-(1/4), -(7/8), -(5/16)};
mats = {{{-1, -1, 0}, {-31, -32, -32}, {30, 31, 
     30}}, {{-121, -109, -120}, {-101, -91, -100}, {210, 189, 208}}};
vecs = {{17/8, 143/4, -44}, {65/8, 1289/8, (137/8)}};
results = # . sol & /@ mats
results - vecs

(* Out[552]= {{9/8, 183/4, -44}, {1305/8, 1089/8, -(2263/8)}}

Out[553]= {{-1, 10, 0}, {155, -25, -300}} *)

So sol does what is required mod Z.

Yes, it is easy to get the same results as you obtained with GAP. They are (in WL notation) {{3/4, 1/8, 3/16}, {3/4, 1/8, 11/16 ]}}. Recall that I showed that first one. Also I mentioned the null vectors are the same as in a prior example. They were/are {1,0,0}, {0,1,0}, and {0,0,1/2}.

As with linear systems in general, any solution can be cast as a specific solution plus a null vector. You are interested in solutions modulo Z and, since the matrices had all integer entries, any solution can trivially be modified by an integer. This means it makes sense to restrict to solutions lying in the unit cube. One such, as noted, is {3/4, 1/8, 3/16}. All others in this restricted range will of necessity be found by adding multiples of the third null vector, {0,0,1/2}. Quite obviously the only other restricted solution will be {3/4, 1/8, 3/16} + {0,0,1/2}, or {3/4, 1/8, 11/16}. So that's how we get the second GAP solution.

One can of course put this reasoning into code, using the specific solution and null vectors. I'll show for the example at hand.

sol = {3/4, 1/8, 3/16};
nulls = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1/2}};
vars = Array[c, Length[nulls]];
vals = Array[d, Length[sol]];
vals /. Solve[
  Flatten[{Thread[sol + vars . nulls == vals], 
    Map[0 <= # <= 1 &, vals], Element[vars, Integers]}], 
  Join[vars, vals]]

(* Out[193]= {{3/4, 1/8, 3/16}, {3/4, 1/8, 11/16}} *)

In order:

(1) Just like working over a field, solutions are comprised of a specific solution plus arbitrary combinations over Z of null solutions. Those latter are also found in the HNF.

nullsZ = 
 Cases[hhnew, {Repeated[0, {Length[bigvec] + 1}], 
    rest__} :> {rest}[[1 ;; Length[mats[[1, 1]]]]]]

(* Out[576]= {{16, 0, 0}, {0, 16, 0}, {0, 0, 8}} *)

nullvecs = nullsZ/(mult*lcm)

(* Out[577]= {{1, 0, 0}, {0, 1, 0}, {0, 0, 1/2}} *)

The first two are "obvious" insofar as we can add any integer to components of the specific solution found and still have a solution modulo 1. The third null vector shows we can add half integers to the third component. For example:

In[578]:= mats[[1]] . (mod1soln + nullvecs[[3]]) - vecs[[1]]

(* Out[578]= {-3, -85, 91} *)

(2) See (1).

(3) I do not know of a quick way offhand. I think you'd have to do the solve process as shown, and see if there are no rows in the HNF that have annihilated that bigvec. In the terminology I used, that would mean keeprow came back empty or in the form Missing[...].

Are you sure there is such an x? Maybe for a start work out what it could be if you use only the first two matrices and corresponding vectors.

Now, I try to obtain the nullvecs computed by you, which can be generated in GAP as follows:

gap> mat1:= [ [-210, -210, -220], [ -221, -222, -232], [410, 411, 430 ] ];
[ [ -210, -210, -220 ], [ -221, -222, -232 ], [ 410, 411, 430 ] ]
gap> mat2:= [ [-132, -121, -132], [ -120, -110, -120], [240, 220, 240] ];
[ [ -132, -121, -132 ], [ -120, -110, -120 ], [ 240, 220, 240 ] ]
gap> mats:=[mat1, mat2];
[ [ [ -210, -210, -220 ], [ -221, -222, -232 ], [ 410, 411, 430 ] ], [ [ -132, -121, -132 ], [ -120, -110, -120 ], [ 240, 220, 240 ] ] ]
gap> LoadPackage("fr");
true
gap> # M * x = b
gap> d:=DimensionsMat(M)[1];
6
gap> nullvecs:=SolveInhomEquationsModZ(M, 0*[1..d], false )[1];
[ [ 0, 0, 0 ], [ 0, 0, 1/2 ] ]
gap> # TransposedMat(M)* x = b
gap> # M_mn * x_n1 = b_m1
gap> d:=DimensionsMat(TransposedMat(M))[1];
3
gap> nullvecs:=SolveHomEquationsModZ(TransposedMat(M))[1];
[ [ 0, 0, 0 ], [ 0, 0, 1/2 ] ]
gap> # or
gap> nullvecs:=SolutionMatMod1(TransposedMat(M), 0*[1..d]);
[ [ 0, 0, 0 ], [ 0, 0, 1/2 ] ]

Then I try to compute with your WL based method as follows, but failed:

In[88]:= 
mat1 = {{-210, -210, -220}, {-221, -222, -232}, {410, 411, 430}};
mat2 = {{-132, -121, -132}, {-120, -110, -120}, {240, 220, 240}};
mats = {mat1, mat2};
bigmat = Apply[Join[##, 2] &, Map[Transpose, mats]];

In[92]:= vec1 = {-27, -28, 105/2};
vec2 = {-47/8, -29/4, 25/2};
vecs = {vec1, vec2};
bigvec = vecs // Flatten;

In[102]:= lcm = LCM @@ Denominator[bigvec];
{uu, hh} = HermiteDecomposition[bigmat];
keeprow = 
  FirstCase[
   hh, {Repeated[0, {Length[bigvec]}], nz : Except[0], rest__} :> {nz,
      rest}];
mult = keeprow // First

nullsZ = 
 Cases[hh, {Repeated[0, {Length[bigvec] + 1}], 
    rest__} :> {rest}[[1 ;; Length[mats[[1, 1]]]]]]
nullvecs = nullsZ/(mult*lcm)


Out[105]= "NotFound"

Out[106]= {}

Out[107]= {}

Any further tips/comments/remarks will be appreciated.

Regards, Zhao

Dear Zhao,

I found in another thread (Pari/GP bulletin board it seems) that you also asked about a different matrix and vector pair.

mat1 = {{-210, -210, -220}, {-221, -222, -232}, {410, 411, 430}};
mat2 = {{-132, -121, -132}, {-120, -110, -120}, {240, 220, 240}};
vec1 = {-27, -28, 105/2};
vec2 = {-47/8, -29/4, 25/2};

Using the same code as I showed in this thread I get (again) a common solution modulo Z of {3/4, 1/8, 3/16}. Also the null vectors are the same as with the prior example. Perhaps the two examples are related in your work. I am really not familiar with crystallographic groups so I'm a bit lost there.

Anyway I thought this might be useful to you. From what I could tell, the result in that other thread excluded the second matrix and vector. Possibly readers failed to understand that you are seeking a common solution to both pairs. Or maybe I misunderstand what is wanted.

I meant to ask: Do you have references for this type of solving? Or applications? I tried to look up this type of solving and got nowhere.

There is no trial-and-error guessing. The part ' mult=First[keeprow]' determines whether we have a solution, and, if not, what new multiplier we require in the denominators. It's actually a complete algorithm.

Sorry, my statement is unclear. I mean, in all the context mentioned here, it corresponds to the following meaning:

mat2 . #&/@sol

Is vec2 meant to be distraction or did you want that in the identity somewhere? Possible dimension mismatch?

I really don't understand why you say that. Can you give more specific examples?