Name for this Algorithm

algorithmsdivisibilityeuclidean-algorithm

I've managed to prove a bunch of properties about this algorithm that I came up with. I'm now curious to know it's name to see what other people have done.

Given a number in base b $$N_0 = b N_X + N_L, \qquad 0 \leq N_L \leq b-1 $$
And a divisor also in base b $$D = b D_X + D_L, \qquad 0 \leq D_L \leq b-1 $$

($N_L$ and $D_L$ are the last digits of N and D)

Then compute $$N_{i+1} = \left| \begin{matrix} N_X & N_L \\ D_X & D_L \end{matrix} \right|$$

I've shown this algorithm to a few lecturers at my uni and they don't know a name for it but suspect it is some sort of variation on the Euclidean Algorithm.

It can be used for divisibility tests (especially good in binary, if $D_L = 1$, or if $D$ is prime )

E.g. Does $9734$ divide by $31$

$N_X = 973$, $N_L = 4$

$D_X=3$, $D_L = 1$

$$N_0 = 9734$$

$$N_1 = \left| \begin{matrix} 973 & 4 \\ 3 & 1 \end{matrix} \right| = 961$$
$$N_2 = \left| \begin{matrix} 96 & 1 \\ 3 & 1 \end{matrix} \right| = 93$$
$$N_3 = \left| \begin{matrix} 9 & 3 \\ 3 & 1 \end{matrix} \right| = 0$$

If the algorithm hits $0$, then $D|N$ (Under certain conditions, counter example N=48, D=36)

Also, if $D_L = 1$ the quotient appears

$$N_1 = \left| \begin{matrix} 973 & \color{red}{4} \\ 3 & 1 \end{matrix} \right| = 961$$
$$N_2 = \left| \begin{matrix} 96 & \color{red}{1} \\ 3 & 1 \end{matrix} \right| = 93$$
$$N_3 = \left| \begin{matrix} 9 & \color{red}{3} \\ 3 & 1 \end{matrix} \right| = 0$$

$9734 = \color{red}{314} \times 31$

Any help is greatly appreciated! Thanks, Ben

Best Answer

It's the special case $m = 10c\!+\!d,\ n = 10a\!+\!b\ $ of the following

$\quad \bbox[5px,border:1px solid red]{(a,n) = 1\,\Rightarrow\, (m,n) = (am\!-\!cn,n)}\ \ $ by $\,(m,n)=(am,n) = (am\!-\!cn,n)$

So $\ (a,b)=1\,\Rightarrow\, (10c\!+\!d,10a\!+\!b) = (ad\!-\!bc,10a\!+\!b),\ $ exactly as you observed.

Remark $ $ Various well-known divisibity test rules are essentially special cases of this, e.g. the rule for divisibility by $7$ follows from the rule for $21$ given below (we negated $\,ad\! -\! bc = 2d\!-\!c\,$ below)

$\qquad\qquad \qquad\ \ \ \, (10c\!+\!d,10(2)\!+\!1) = (c\!-\!2d, 10(2)\!+\!1)$

${\rm similar\ to\ your}\ \ (10c\!+\!d,10(3)\!+\!1) = (c\!-\!3d, 10(3)\!+\!1)\ $ rule for $\,31$

Related Solutions

[Math] Computing Hermite Normal Form using Extended Euclidean Algorithm

I think there are some errors, either book errata, or transcription errors, in the stated algorithm. Firstly, $k$ is the pivot column, so it needs to be decremented (whenever nonzero) at the end of the first $i$ loop. Secondly, for each iteration in this loop, $j$ needs to be looped over from $k-1$ down to $1$, to clear out the lower triangular part, so it must be initialized inside the $i$ loop.

I am also skeptical of whether special handling of xgcd is necessary or even helpful, since it always returns a "positive definite" $d$ (an admitted abuse of terminology), in the sense that it is (or should be) always positive given $(a,b)\ne(0,0)$. Because of this, we also do not need the $b < 0$ test (and correction) in the final (partial upper triangular) reductions. In fact, it seems rather artificial to require $v \le 0$ and $u > 0$ of xgcd. What about $(a,b)=(3,1)$ and $(-1,3)$? Also, your pseudocode for the xgcd wrapper only alters its output if $u$ or $v$ is zero, and I think the loop could be avoided by solving some linear equations (taking quotients and a minimum).

However, I don't have Cohen's book, and don't know how or why he requires $|u|$ and $|v|$ to be minimal. Does he mean that $\max\{|u|,|v|\}$ must be minimal? I admit, I haven't thought deeply about it. But I think it is a fact that xgcd already does a good job of producing a minimal solution.

I have corrected the algorithm, as I think it should be, and implemented it in sage. After testing described below, I also made a correction.

m = Matrix(ZZ,4,4,\
    [[10,-2,-2, 4]\
    ,[ 0,-2,-1,-6]\
    ,[ 1, 2, 3,-2]\
    ,[ 0, 1,-1, 0]]).transpose(); m

def hnf(m,show=False):
    nr = m.nrows()
    nc = m.ncols()
    if show:
        print 'hnf: '+str(nr)+'x'+str(nc)
    assert(nr == nc)
    c = m.columns()
    M = Matrix(ZZ,nc,nr,c).transpose()
    if show:
        print M
    def r(n):
        s = range(n)
        s.reverse()
        return s
    k = nc - 1
    pivots = []
    for i in r(nr):
        for j in r(k):
            ck = c[k]; a = ck[i]
            cj = c[j]; b = cj[i]
            if b != 0:
                d,u,v = xgcd(a,b)
                if show:
                    print 'row '+str(i+1)+', cols '+str(j+1)+','+str(k+1)\
                        , 'xgcd:',u,'(',a,') + ',v,'(',b, ') = ',d
                ad = a/d
                bd = b/d
                c[k] = ck*u  + cj*v  # if a != d:
                c[j] = cj*ad - ck*bd
                M = Matrix(ZZ,nc,nr,c).transpose()
                if show:
                    print M
        if c[k][i] != 0:
            pivots.append((i,k))
            if c[k][i] < 0:  # this is necessary for the last pivot column,
                c[k] = -c[k] # since it cannot utilize xgcd to "normalize"
            k = k - 1
    for (i,k) in pivots:
        b = c[k][i]; assert(b > 0)
        if show:
            print 'pivot('+str(i+1)+','+str(k+1)+')='+str(b)
        for j in [k+1..nc-1]:
            q = floor(c[j][i] / b)
            c[j] -= c[k]*q
        M = Matrix(ZZ,nc,nr,c).transpose()
        if show:
            print M
    pivots.reverse()
    #if show:
    #    print 'pivots:', pivots
    return M, pivots

h,p = hnf(m,True)

Here is the output (for your test case):

hnf: 4x4
[10  0  1  0]
[-2 -2  2  1]
[-2 -1  3 -1]
[ 4 -6 -2  0]
row 4, cols 3,4 xgcd: 0 ( 0 ) +  -1 ( -2 ) =  2
[10  0  0 -1]
[-2 -2  1 -2]
[-2 -1 -1 -3]
[ 4 -6  0  2]
row 4, cols 2,4 xgcd: 1 ( 2 ) +  0 ( -6 ) =  2
[ 10  -3   0  -1]
[ -2  -8   1  -2]
[ -2 -10  -1  -3]
[  4   0   0   2]
row 4, cols 1,4 xgcd: 1 ( 2 ) +  0 ( 4 ) =  2
[ 12  -3   0  -1]
[  2  -8   1  -2]
[  4 -10  -1  -3]
[  0   0   0   2]
row 3, cols 2,3 xgcd: -1 ( -1 ) +  0 ( -10 ) =  1
[12  3  0 -1]
[ 2 18 -1 -2]
[ 4  0  1 -3]
[ 0  0  0  2]
row 3, cols 1,3 xgcd: 1 ( 1 ) +  0 ( 4 ) =  1
[12  3  0 -1]
[ 6 18 -1 -2]
[ 0  0  1 -3]
[ 0  0  0  2]
row 2, cols 1,2 xgcd: 0 ( 18 ) +  1 ( 6 ) =  6
[33 12  0 -1]
[ 0  6 -1 -2]
[ 0  0  1 -3]
[ 0  0  0  2]
pivot(4,4)=2
[33 12  0 -1]
[ 0  6 -1 -2]
[ 0  0  1 -3]
[ 0  0  0  2]
pivot(3,3)=1
[33 12  0 -1]
[ 0  6 -1 -5]
[ 0  0  1  0]
[ 0  0  0  2]
pivot(2,2)=6
[33 12 12 11]
[ 0  6  5  1]
[ 0  0  1  0]
[ 0  0  0  2]
pivot(1,1)=33
[33 12 12 11]
[ 0  6  5  1]
[ 0  0  1  0]
[ 0  0  0  2]

Sage also has a hermite_form function for matrices,

var('s')
for s in ['pari', 'pari0', 'pari4', 'padic', 'default', 'ntl']:
    print str(s)
    print m.hermite_form(algorithm=s)

but each implementation produces the matrix below $$ \left(\matrix{ 2 & 0 & 2 & 27 \\ 0 & 1 & 1 & 64 \\ 0 & 0 & 3 & 45 \\ 0 & 0 & 0 & 66 }\right) $$ which, as you point out in a comment, is arrived at by iterating over $i$ increasing, i.e. by pivoting from the top (left) downwards, rather than decreasing from the bottom (right) upwards. We can verify this directly using the transformation $\phi:m\mapsto rm^tr=(rmr)^t$, for $r$ the (orthogonal) matrix $$ r= \left(\begin{array}{rrrr} 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \end{array}\right), $$ which transposes $m$ by both diagonals. Our (Cohen's) HNF (which pivots from the bottom right upward) and sage's (which forms pivots from the top left downwards) can be obtained from each other by pre- and post- transforming by $\phi:$

r = Matrix(ZZ,4,4,[\
    floor(1/(1+n%3))*\
 (1-floor(1/(1+n%15)))\
    for n in range(16)])
print r
phi = lambda x: (r*x*r).transpose()
print phi(m).hermite_form() == phi(h)
phi(phi(m).hermite_form())

By the way, if you don't mind a little excursion, the group of transformations of square matrices generated under composition by the (self-inverse) maps $\color{red}{t}:m\mapsto m^t$ (transpose) and $\color{blue}{r}:m\mapsto rm$ (row reversal) or $\color{green}{c}:m\mapsto mr$ (column reversal) is in fact isomorphic to the dihedral group $D_4$ of the square. We can visualize this group by its action on a fixed matrix, e.g. the $2\times2$ matrix $e=\left(\matrix{1&2\\3&4}\right)$. The dihedral group of a square as vertices of a cube

Each transformation is represented by, or identified with, (a path from $e$ to) a vertex. The order two elements consist of the surface edges and interior diagonals, while the rotations consist of the face diagonals, and their cosets partition the vertices into two intersecting tetrahedra. This is also a commutative diagram, and thus can be used to verify such formulas as $rc=cr$, $rt=tc$, and $ct=tr$ as well as $\phi=rtc=trc$ from above.

Anyway after some testing

ranks = [0,0,0,0,0]
passed = 0
failed = 0
for n in range(10^2):
    m = random_matrix(ZZ,4)
    h, p = hnf(phi(m))
    if phi(h) == m.hermite_form():
        ranks[len(p)] += 1
        passed += 1
    else:
        print 'Test FAIL!'
        print 'm:'; print m
        print 'm"=phi(m):'; print phi(m)
        print 'hnf(m"):'; print h
        print 'hnf(m")":'; print phi(h)
        print 'm.hnf():'; print m.hermite_form()
        failed += 1
        break
print n+1,'tests performed:',passed,'passed,',failed,'failed.'
print 'rank counts:',ranks

I found that I needed to sometimes correct the sign of the first pivot column, and incorporated this into the hnf function above (along with making the print output optional but off by default). In any case, in sage, you can just define $\phi$ and use phi(phi(m).hermite_form()); after all, why should we add a seventh implementation (except to learn!)?

Sorting numbers with a special set of permutations

Here is an algorithm that can sort in exactly $d$ steps ($d$ permutations). It uses only the so-called $2$nd subset, where all operations bring the last entry (the $d$) to the front, while allowing you to "shoot" the front entry (the $0$) anywhere.

My key idea is that, the power to move $0$ anywhere makes this similar to how one can sort a hand of cards. I.e. maintain a sorted sub-sequence, and at every step take a new card and insert it into the right place in the sorted sub-sequence, growing the sorted sub-sequence by one card. After inserting $d$ cards you'd be done. It took me a few tries to get all the details right, to map this onto the $2$nd subset.

Here is my first attempt:

Initialization: Color the $2$nd entry (i.e. position $1$ in your numbering scheme) in the input array red. (This step is conceptual - no permutation is actually done.)

Invariant: In between turns, the red entries will always start at the second position (position $1$), be contiguous, and be sorted.

Every turn: Pick the $1$st entry (position $0$), which will be not red. Color it red, and insert it into the existing red sub-array at the right place.

Here is an example: We will sort $FCBDGEA$.

$$ \begin{align} FCBDGEA & = F\color{red}{C}BDGEA \, \, \, \, \, \, \text{(initialization)}\\ & \rightarrow A\color{red}{CF}BDGE \\ & \rightarrow E\color{red}{ACF}BDG \\ & \rightarrow G\color{red}{ACEF}BD \\ & \rightarrow D\color{red}{ACEFG}B \\ & \rightarrow B\color{red}{ACDEFG} \\ & \rightarrow \color{red}{GABCDEF} \\ \end{align} $$

Clearly (e.g. by the invariant), after $d$ turns (i.e. $d$ permutations) we end up with the sorted array, but starting at the $2$nd position. I.e., the result is one left-shift from the desired answer. Left-shift is not allowed, but right-shift is in the $2$nd subset, so we can do $d$ more right-shifts and arrive at the desired answer $ABCDEFG$. Therefore this first-attempt runs in $2d$ steps.

My second attempt fixes the "bug" by declaring that the minimum element (here, $A$) be treated as the maximum. Except for this modified sort order, the algorithm is identical to the first attempt. The same example:

$$ \begin{align} FCBDGEA & = F\color{red}{C}BDGEA \, \, \, \, \, \, \text{(initialization)}\\ & \rightarrow A\color{red}{CF}BDGE \\ & \rightarrow E\color{red}{CFA}BDG \, \, \, \, \, \, \text{(treat $A$ as maximum)}\\ & \rightarrow G\color{red}{CEFA}BD \\ & \rightarrow D\color{red}{CEFGA}B \\ & \rightarrow B\color{red}{CDEFGA} \\ & \rightarrow \color{red}{ABCDEFG} \\ \end{align} $$

Just like in the first attempt, after $d$ steps (permutations) we end up with the "sorted" array starting at the $2$nd position - except now "sorted" means according to the "modified" sort order, i.e. $BCDEFGA$. Clearly, this means the array starting at the front will be sorted according to the unmodified order. I.e. this second attempt algorithm sorts the array in exactly $d$ steps (permutations), as claimed.

Best Answer

Related Solutions

[Math] Computing Hermite Normal Form using Extended Euclidean Algorithm

Sorting numbers with a special set of permutations

Related Question