1# Copyright (C) 2020 Atsushi Togo
2# All rights reserved.
3#
4# This file is part of phonopy.
5#
6# Redistribution and use in source and binary forms, with or without
7# modification, are permitted provided that the following conditions
8# are met:
9#
10# * Redistributions of source code must retain the above copyright
11#   notice, this list of conditions and the following disclaimer.
12#
13# * Redistributions in binary form must reproduce the above copyright
14#   notice, this list of conditions and the following disclaimer in
15#   the documentation and/or other materials provided with the
16#   distribution.
17#
18# * Neither the name of the phonopy project nor the names of its
19#   contributors may be used to endorse or promote products derived
20#   from this software without specific prior written permission.
21#
22# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
23# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
24# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
25# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
26# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
27# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
28# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
29# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
30# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
32# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
33# POSSIBILITY OF SUCH DAMAGE.
34
35import numpy as np
36
37
38#
39# Smith normal form for 3x3 integer matrix
40# This code is maintained at https://github.com/atztogo/snf3x3.
41#
42class SNF3x3(object):
43    """Smith normal form for 3x3 matrix
44
45    Input 3x3 interger matrix A is transformed to 3x3 diagonal interger
46    matrix (D) by two 3x3 interger matrices P and Q, which is written as
47
48        D = PAQ
49
50    with abs(det(A)) = det(D), det(P) = 1, det(Q) = sgn(det(A)).
51
52    The algorithm is implemented following
53    https://en.wikipedia.org/wiki/Smith_normal_form,
54    but the diagonal elements don't follow the rule of it.
55
56    Usage
57    -----
58    snf = SNF3x3(int_mat)
59    snf.run()
60
61    Attributes
62    ----------
63    D, P, Q : ndarray
64        These 3x3 interger matrices explained above.
65        shape=(3, 3), dtype='intc'
66
67    """
68
69    def __init__(self, A):
70        """
71
72        Parameters
73        ----------
74        A : array_like
75            3x3 interger matrix.
76            shape=(3, 3)
77
78        """
79
80        self._A_orig = np.array(A, dtype='int_', order='C')
81        self._A = np.array(A, dtype='int_', order='C')
82        self._Ps = []
83        self._Qs = []
84        self._L = []
85        self._P = None
86        self._Q = None
87        self._D = None
88        self._attempt = 0
89
90    def run(self):
91        for i in self:
92            pass
93
94        A = np.dot(self._P, np.dot(self._A_orig, self._Q))
95        assert (A == self._A).all()
96
97    def __iter__(self):
98        return self
99
100    def __next__(self):
101        self._attempt += 1
102        if self._first():
103            if self._second():
104                self._finalize()
105                raise StopIteration
106        return self._attempt
107
108    def next(self):
109        self.__next__()
110
111    @property
112    def A(self):
113        return self._A
114
115    @property
116    def D(self):
117        return self._D
118
119    @property
120    def P(self):
121        return self._P
122
123    @property
124    def Q(self):
125        return self._Q
126
127    def _first(self):
128        self._first_one_loop()
129        A = self._A
130        if A[1, 0] == 0 and A[2, 0] == 0:
131            return True
132        elif A[1, 0] % A[0, 0] == 0 and A[2, 0] % A[0, 0] == 0:
133            self._first_finalize()
134            self._Ps += self._L
135            self._L = []
136            return True
137        else:
138            return False
139
140    def _first_one_loop(self):
141        self._first_column()
142        self._Ps += self._L
143        self._L = []
144        self._A[:] = self._A.T
145        self._first_column()
146        self._Qs += self._L
147        self._L = []
148        self._A[:] = self._A.T
149
150    def _first_column(self):
151        i = self._search_first_pivot()
152        if i > 0:
153            self._swap_rows(0, i)
154        if i < 0:
155            raise RuntimeError("Determinant is 0.")
156
157        if self._A[1, 0] != 0:
158            self._zero_first_column(1)
159        if self._A[2, 0] != 0:
160            self._zero_first_column(2)
161
162    def _zero_first_column(self, j):
163        A = self._A
164        r, s, t = xgcd([A[0, 0], A[j, 0]])
165        self._set_zero(0, j, A[0, 0], A[j, 0], r, s, t)
166
167    def _search_first_pivot(self):
168        for i in range(3):  # column index
169            if self._A[i, 0] != 0:
170                return i
171        return -1
172
173    def _first_finalize(self):
174        """Set zeros along the first colomn except for A[0, 0]
175
176        This is possible only when A[1,0] and A[2,0] are dividable by A[0,0].
177
178        """
179
180        A = self._A
181        L = np.eye(3, dtype='int_')
182        L[1, 0] = -A[1, 0] // A[0, 0]
183        L[2, 0] = -A[2, 0] // A[0, 0]
184        self._L.append(L.copy())
185        self._A[:] = np.dot(L, self._A)
186
187    def _second(self):
188        """Find Smith normal form for Right-low 2x2 matrix"""
189
190        self._second_one_loop()
191        A = self._A
192        if A[2, 1] == 0:
193            return True
194        elif A[2, 1] % A[1, 1] == 0:
195            self._second_finalize()
196            self._Ps += self._L
197            self._L = []
198            return True
199        else:
200            return False
201
202    def _second_one_loop(self):
203        self._second_column()
204        self._Ps += self._L
205        self._L = []
206        self._A[:] = self._A.T
207        self._second_column()
208        self._Qs += self._L
209        self._L = []
210        self._A[:] = self._A.T
211
212    def _second_column(self):
213        """Right-low 2x2 matrix
214
215        Assume elements in first row and column are all zero except for A[0,0].
216
217        """
218
219        if self._A[1, 1] == 0 and self._A[2, 1] != 0:
220            self._swap_rows(1, 2)
221
222        if self._A[2, 1] != 0:
223            self._zero_second_column()
224
225    def _zero_second_column(self):
226        A = self._A
227        r, s, t = xgcd([A[1, 1], A[2, 1]])
228        self._set_zero(1, 2, A[1, 1], A[2, 1], r, s, t)
229
230    def _second_finalize(self):
231        """Set zero at A[2, 1]
232
233        This is possible only when A[2,1] is dividable by A[1,1].
234
235        """
236
237        A = self._A
238        L = np.eye(3, dtype='int_')
239        L[2, 1] = -A[2, 1] // A[1, 1]
240        self._L.append(L.copy())
241        self._A[:] = np.dot(L, self._A)
242
243    def _finalize(self):
244        for i in range(3):
245            if self._A[i, i] < 0:
246                self._flip_sign_row(i)
247                self._Ps += self._L
248                self._L = []
249        self._finalize_sort()
250        self._finalize_disturb(0, 1)
251        self._first()
252        self._finalize_sort()
253        self._finalize_disturb(1, 2)
254        self._second()
255        self._set_PQ()
256
257    def _finalize_sort(self):
258        if self._A[0, 0] > self._A[1, 1]:
259            self._swap_diag_elems(0, 1)
260        if self._A[1, 1] > self._A[2, 2]:
261            self._swap_diag_elems(1, 2)
262        if self._A[0, 0] > self._A[1, 1]:
263            self._swap_diag_elems(0, 1)
264
265    def _finalize_disturb(self, i, j):
266        if self._A[j, j] % self._A[i, i] != 0:
267            self._A[:] = self._A.T
268            self._disturb_rows(i, j)
269            self._Qs += self._L
270            self._L = []
271            self._A[:] = self._A.T
272
273    def _disturb_rows(self, i, j):
274        L = np.eye(3, dtype='int_')
275        L[i, i] = 1
276        L[i, j] = 1
277        L[j, i] = 0
278        L[j, j] = 1
279        self._L.append(L.copy())
280        self._A[:] = np.dot(L, self._A)
281
282    def _swap_diag_elems(self, i, j):
283        self._swap_rows(i, j)
284        self._Ps += self._L
285        self._L = []
286        self._A[:] = self._A.T
287        self._swap_rows(i, j)
288        self._Qs += self._L
289        self._L = []
290        self._A[:] = self._A.T
291
292    def _swap_rows(self, i, j):
293        """Swap i and j rows
294
295        As the side effect, determinant flips.
296
297        """
298
299        L = np.eye(3, dtype='int_')
300        L[i, i] = 0
301        L[j, j] = 0
302        L[i, j] = 1
303        L[j, i] = 1
304        self._L.append(L.copy())
305        self._A[:] = np.dot(L, self._A)
306
307    def _flip_sign_row(self, i):
308        """Multiply -1 for all elements in row"""
309
310        L = np.eye(3, dtype='int_')
311        L[i, i] = -1
312        self._L.append(L.copy())
313        self._A[:] = np.dot(L, self._A)
314
315    def _set_zero(self, i, j, a, b, r, s, t):
316        """Let A[i, j] be zero based on Bezout's identity
317
318           [ii ij]
319           [ji jj] is a (k,k) minor of original 3x3 matrix.
320
321        """
322
323        L = np.eye(3, dtype='int_')
324        L[i, i] = s
325        L[i, j] = t
326        L[j, i] = -b // r
327        L[j, j] = a // r
328        self._L.append(L.copy())
329        self._A[:] = np.dot(L, self._A)
330
331    def _set_PQ(self):
332        P = np.eye(3, dtype='int_')
333        for _P in self._Ps:
334            P = np.dot(_P, P)
335        Q = np.eye(3, dtype='int_')
336        for _Q in self._Qs:
337            Q = np.dot(Q, _Q.T)
338
339        if self._det(P) < 0:
340            P *= -1
341            Q *= -1
342
343        self._P = P
344        self._Q = Q
345        self._D = self._A.copy()
346
347    def _det(self, m):
348        return (m[0, 0] * (m[1, 1] * m[2, 2] - m[1, 2] * m[2, 1])
349                + m[0, 1] * (m[1, 2] * m[2, 0] - m[1, 0] * m[2, 2])
350                + m[0, 2] * (m[1, 0] * m[2, 1] - m[1, 1] * m[2, 0]))
351
352
353def xgcd(vals):
354    _xgcd = Xgcd(vals)
355    return _xgcd.run()
356
357
358class Xgcd(object):
359    def __init__(self, vals):
360        self._vals = np.array(vals, dtype='intc')
361
362    def run(self):
363        r0, r1 = self._vals
364        s0 = 1
365        s1 = 0
366        t0 = 0
367        t1 = 1
368        for i in range(1000):
369            r0, r1, s0, s1, t0, t1 = self._step(r0, r1, s0, s1, t0, t1)
370            if r1 == 0:
371                break
372
373        assert r0 == self._vals[0] * s0 + self._vals[1] * t0
374
375        self._rst = np.array([r0, s0, t0], dtype='intc')
376
377        return self._rst
378
379    def _step(self, r0, r1, s0, s1, t0, t1):
380        q, m = divmod(r0, r1)
381        if m < 0:
382            if r1 > 0:
383                m += r1
384                q -= 1
385            if r1 < 0:
386                m -= r1
387                q += 1
388        r2 = m
389        s2 = s0 - q * s1
390        t2 = t0 - q * t1
391        return r1, r2, s1, s2, t1, t2
392
393    def __str__(self):
394        v = self._vals
395        r, s, t = self._rst
396        return "%d = %d * (%d) + %d * (%d)" % (r, v[0], s, v[1], t)
397