Explanation for the Euclid Extended Algorithm

algorithmseuclidean-algorithmnumber theoryself-learning

Please this is how the code for the extended-euclid algorithm was implemented in the book Introduction to Algorithm (Chapter 31 page 937

EXTENDED-EUCLID(a,b)

if b == 0

return (a,1,0)

else (d_1,x_1,y_1) = EXTENDED-EUCLID(b,a mod b)

$(d,x,y) = (d_1,y_1,x_1 – \lfloor \frac{a}{b} \rfloor)$

return(d,x,y)

Then I tried to implement it in python

def extended_euclid(a,b):
    if b == 0:
        return (a,1,0)
    else:
        (d_1,x_1,y_1) = extended_euclid(b, a % b)
        (d,x,y) = (d_1,y_1,x_1 - (a//b) * y_1)
        return (d,x,y)

But when I run the function with input (99,78) it returns (3, -11, 14)

But the explanation the book gave about how to to find the EXTENDED-EUCLID(99,78) using the above algorithm was

Step 1: The algorithm returns (3,-11,14)

Step 2: The algorithm returns (3,3,-11)

Step 3: The algorithm returns (3,-2,3)

Step 4: The algorithm returns (3,1,-2)

Step 5: The algorithm returns (3,0,1)

Step 6: The algorithm returns (3,1,1)

So I was expecting the final output from the code I wrote to be (3,1,1). Please can someone explain why I am not getting similar answer ?

Best Answer

The steps are given in reverse order. The $(3,1,1)$ should be $(3,1,0)$. You can see what is going on by using print() function before returning.

def extended_euclid(a,b):
    if b == 0:
        print([a,b],[a,1,0]);
        return (a,1,0)
    else:
        (d_1,x_1,y_1) = extended_euclid(b, a % b)
        (d,x,y) = (d_1,y_1,x_1 - (a//b) * y_1)
        print([a,b],[d,x,y]);
        return (d,x,y)

If you execute extended_euclid(99,78) then the resulting output is

[3, 0] [3, 1, 0]
[6, 3] [3, 0, 1]
[15, 6] [3, 1, -2]
[21, 15] [3, -2, 3]
[78, 21] [3, 3, -11]
[99, 78] [3, -11, 14]

while the function returns [3,-11,14]. This means that $$ 3 = (-11)\cdot 99 + (14)\cdot 78 = -1089 + 1092 $$ is the gcd $\,3\,$ expressed as a linear combination of the two numbers $\,99,78.$

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!)?

[Math] Knuth’s mastermind algorithm

Let's say $S_i$ is an element of $S$, the possible values of the winning code.

What you're asking of each $S_i$ at the step you're puzzling about is this:

If the answer were $S_i$, would I have gotten the response I got with my current guess?

If the answer is no, then $S_i$ cannot be the code, because you would have gotten a different response. If the answer is yes, then you keep that around for another turn, because it's consistent with the response you got. In other words, you can't eliminate it as impossible.

Looking at the example in the other thread:

You guess $1122$ and get one black peg (one peg color correct and in the right position) and one white peg (one peg color correct but in the wrong position).

The step of reducing the list involves asking, for every $S_i$, "Would I have gotten one black peg and one white peg with my guess of $1122$ if the code were $S_i$?"

OK, so let's do a few:

If the code were $1111$ I would have gotten two black pegs ($BB$) with my guess of $1122$, which is not the same as one black peg and one white peg ($BW$). (Represent this by $F(1122, 1111) = BB$.) So, I remove $1111$ from the list of potential solutions.
$F(1122, 1112) = BBB \neq BW \to \text{Remove 1112 from S}$
$F(1122, 1113) = BB \neq BW \to \text{Remove 1113 from S}$
$F(1122, 1114) = BB \neq BW \to \text{Remove 1114 from S}$
$F(1122, 1115) = BB \neq BW \to \text{Remove 1115 from S}$

And so forth.

However,

$F(1122, 1314) = BW = BW \to \text{Keep 1314 in S}$

Hope this helps.

Best Answer

Related Solutions

[Math] Computing Hermite Normal Form using Extended Euclidean Algorithm

[Math] Knuth’s mastermind algorithm

Related Question