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

The two paragraphs plus code that follow show how to proceed when no solution is found on the first try.

Hi Daniel,

I notice that in the last step of your approach there are some difficulties and problems as follows:

• You said the following:

The catch is that we might not obtain such a vector in the HNF.

This means that the method may miss solutions.

• In this step, there are a lot of empirical judgment and trial-and-error, which make this approach difficult to implement programmatically.

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.

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.

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}} *)

Hi Daniel,

First, I appreciate your wonderful analysis and detailed solving process. I need some time to fully understand your method and specific details.

However, based on the general characteristics of the problem, I still have the following questions:

1. For a problem like this, if there are solutions, then there should be infinitely many solutions, how do I find and represent all of them in a reasonable form?
2. What are the characteristics and properties of the set of all solutions?
3. How to quickly determine if there is no solution?

$a \times b \mod d = a \times (b + c) \mod d \implies a \times c \mod d = 0 \implies c = k \times d$, where $k$ is an integer. This property entirely explains your identity that:

Mod[mat2 . #&/@sol,1]==Mod[mat2 . (#+rnd)&/@sol,1]
Out[] = True

Because all elements in rnd are integers of the form $k \times 1$. Is vec2 meant to be distraction or did you want that in the identity somewhere? Possible dimension mismatch?

You say that :

mat2 . sol = vec2

but the dimensions of mat2, sol, and vec2 are inconsistent:

$$(3\times3)\cdot(3\times 3) = (3 \times 3 )\neq1 \times 3.$$

I thought you may have made a mistake, meaning by sol in this equation, one row from the sol example you gave, but none of the three dot products equal to vec2.

Hi Zhao,

I think I understand the problem. What has been elusive thus far is how to solve it. I'll try some ,ore and post a response if I get anywhere.

---edit ---

I posted a notebook that describes how to proceed. I did not elaborate too much but here is the outline.

First I show a not-quite-standard way to set up a linear system to solve using the reduced echelon form.

Then I show how it readily extends to find, when one exists, a common solution to two (or more) such systems.

The next step is to extend the idea to solving over integers modulo some given integer. Here we do not assume the given integer is prime.

The important last step is to adjust this technology to finding a rational solution that is correct modulo the integers. This solves the particular example I had posed, which in turn was a minor change to one you had provided.

--- end edit ---

I used the following option for Solve in the previous post:

Modulus->Integers 

However, it should be changed to the following option, which has been corrected in the previous post:

Modulus->1

But, even so, WL still does not support related processing.

You are right. This is the exact mathematical description of the problem, as shown below:

results - vecs == 0, (mod Z^n), where n is the dimension of these vectors.

But how to solve this problem in WL? I tried the following method, but WL does not support it at all:

In[51]:= m=Array[v,3]
Solve[ConstantArray[ConstantArray[0,Length[vecs[[1]]]], Length[mats]]==(# . m & /@ mats), m, Modulus->1]

Out[51]= {v[1], v[2], v[3]}

During evaluation of In[51]:= Solve::rmod: Modulus 1 must be a non-negative integer other than 1.

Out[52]= Solve[{{0, 0, 0}, {0, 0, 
    0}} == {{-v[1] - v[2], -31 v[1] - 32 v[2] - 32 v[3], 
    30 v[1] + 31 v[2] + 30 v[3]}, {-121 v[1] - 109 v[2] - 
     120 v[3], -101 v[1] - 91 v[2] - 100 v[3], 
    210 v[1] + 189 v[2] + 208 v[3]}}, {v[1], v[2], v[3]}, 
 Modulus -> 1]

See the following example:

In[86]:= 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={{9/8, 183/4, -44}, {1305/8, 1089/8, -(2263/8)}};
# . sol&/@mats
%==vecs

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

Out[90]= True

I really don't know how to solve this problem in Wolfram. Here is a GAP-based solution to the first matrix and its corresponding vector given in the above example:

gap> mat:=[ [ -2, 0, 0 ], [ 0, -2, 0 ], [ 1, 1, 0 ] ];
[ [ -2, 0, 0 ], [ 0, -2, 0 ], [ 1, 1, 0 ] ]
gap> vec1:=[ -23/8, 17/8, -9/8 ];
[ -23/8, 17/8, -9/8 ]
gap> vec2:=List(vec1, x -> FractionModOne(x));
[ 1/8, 1/8, 7/8 ]
gap> # false means acting on left, a.k.a,
gap> # compute the set of solutions of the equation
gap> # M * x = b  (mod Z).
gap> sol1:=SolveInhomEquationsModZ(mat, vec1, false);
[ [ [ 15/16, 15/16, 0 ], [ 7/16, 7/16, 0 ] ], [ [ 0, 0, 1 ] ] ]
gap> sol2:=SolveInhomEquationsModZ(mat, vec2, false);
[ [ [ 15/16, 15/16, 0 ], [ 7/16, 7/16, 0 ] ], [ [ 0, 0, 1 ] ] ]
gap> sol1=sol2;
true
gap> # These soloutions are all on the integer lattice spanned by vec2:
gap> List(Union(sol1), x -> VectorModL(mat * x, [vec2] ));
[ [ 0, 0, 0 ], [ 0, 0, 7 ], [ 0, 0, 15 ] ]