EGCD development

Functions in this file are NOT designed to handle negative integers.  They do not work correctly if their arguments
are negative.  Presumably this would be fairly easy to fix.

At each step in the extended Euclidean algorithm, we need to keep track of the following items:

a:  the greater of the two numbers of which we take the GCD
b:  the smaller of the two numbers of which we take the GCD
d:  the larger of the two previous remainders computed
e:  the smaller of the two previous remainders computed
x1, y1 :  numbers such that a*x1+b*y1 = d
x2,y2:   numbers such that a*x2 + b*y2 = e.

At each step, we divide d by e, obtaining quotient q = floor(d/e) and remainder r = d mod e.

d will be replaced by e at the next step, and e by r.

x1, y1 will be replaced by x2 and y2 (because d is replaced by e).

r = d - qe = (a*x1+b*y1)-q(a*x2+b*y2) = (x1-q*y1)*a + (x2-q*y2)b

so the new values for x2 and y2 need to be x1-q*y1, x2-q*y2.

The in-class version kept a,b as part of the data structure, but this turns out not to be necessary (we use a,b
to state the loop invariant, but they play no direct role in the calculation).

All of this motivates the following definition:
> onestep := (d,e,x1,y1,x2,y2) -> (e,d mod e,x2,y2,x1-floor(d/e)*x2,y1-floor(d/e)*y2);

  onestep := (d, e, x1, y1, x2, y2) ->

        (e, d mod e, x2, y2, x1 - floor(d/e) x2, y1 - floor(d/e) y2)

The following is a test of onestep
> onestep(120,32,1,0,0,1);

                         32, 24, 0, 1, 1, -3

> onestep(%);

                         24, 8, 1, -3, -1, 4

The EGCD algorithm starts with (a,b,1,0,0,1) and continues until e goes evenly into d; at that point e is the gcd.
> loop := (d,e,x1,y1,x2,y2)->if d mod e = 0 then (x2,y2) #extra data discarded on termination
> else loop(onestep(d,e,x1,y1,x2,y2)) end if;
> 

  loop := proc(d, e, x1, y1, x2, y2)
option operator, arrow;
    if d mod e = 0 then x2, y2
    else loop(onestep(d, e, x1, y1, x2, y2))
    end if
end proc

A test follows.
> loop(120,32,1,0,0,1);

                                -1, 4

Now for the real EGCD algorithm.  As noted above, we start with (a,b,1,0,0,1).  If b>a we compute EGCD(b,a) and reverse the coefficients in the result.
> reverse:=(a,b)->(b,a);

                     reverse := (a, b) -> (b, a)

> EGCD:=(a,b) -> if b>a then reverse(EGCD(b,a)) else loop(a,b,1,0,0,1) end if;

  EGCD := proc(a, b)
option operator, arrow;
    if a < b then reverse(EGCD(b, a))
    else loop(a, b, 1, 0, 0, 1)
    end if
end proc

A test follows.
> EGCD(120,32);

                                -1, 4

> EGCD(32,120);

                                4, -1

Now we define a function which computes modular inverses.
> first_item:=(a,b)->a;

                      first_item := (a, b) -> a

> modinverse:=(a,n)->if a>=n then modinverse(a mod n,n)
> else if gcd(a,n)<>1 then "Error"
> else first_item(EGCD(a,n)) mod n end if end if;

  modinverse := proc(a, n)
option operator, arrow;
    if n <= a then modinverse(a mod n, n)
    else
        if gcd(a, n) <> 1 then "Error"
        else first_item(EGCD(a, n)) mod n
        end if
    end if
end proc

Now the inevitable tests:
> modinverse(7,10);

                                  3

> modinverse(17,10);

                                  3

> modinverse(5,10);

                               "Error"

> modinverse(1154646356462543627,1947386476427150937158216782157165716);

                 713181439101404684559360074474141803

> (1/1154646356462543627) mod 1947386476427150937158216782157165716;

                 713181439101404684559360074474141803

> 

                                -1, 1

> 

                                -5, 2

> 

                                 4, 1

>