Algebra Solutions: Program to compute the greatest common divisor of two integers

Write a computer program to determine the greatest common divisor $d$ of two integers $a$ and $b$ and to express $d$ as $ax + by$ for some integers $x$ and $y$ .

I will be writing my programs in Haskell.
First the gcd program. I have to admit that I wrote my own gcd function, but was disappointed by the performance when compared to the library gcd that is built in to Haskell. My own code took about four times as much time and twice as much space. After tweaking for a while I got the performance to twice the time. I finally peeked at the Prelude definition to see what it did differently, to find out that the only major difference was the use of “rem” (for integer division remainder) instead of “mod”. FYI: rem is considerably more efficient than mod. I do disagree with the library function on one point, however; I feel very strongly that $\mathsf{gcd}(0,0)$ should be 0 and not undefined, as the Haskell specification says. First of all, it satisfies the definition rather trivially. We have $0|0$ since $0 \cdot k = 0$ for all integers $k$ , and if $k|0$ , then $k|0$ . Second, this has the side effect of making gcd idempotent on $\mathbb{Z}$ . Perhaps that’s why the Haskellers don’t like it.
Anyway, here is the code for gcd.

{- Compute the GCD of a and b. -}

gcd :: Integer -> Integer -> Integer

gcd a 0 = abs a

gcd a b = gcd b (a `rem` b)

This function uses about 15 reductions and 20 cells per step. (“Reduction” and “cell” are the fancy words Haskell uses to measure time and space complexity, respectively.) By the way, I found this number by computing the gcd of consecutive Fibonacci numbers; I forget the reason, but Fibonacci numbers give the worst case time complexity (for their size) for this implementation of the gcd function.
Now for the second program. What this program needs to find are two integers which satisfy Bezout’s identity for $a$ and $b$ . For this reason I’ll call the function “bezout”. Here is my code.

swap (a,b) = (b,a)

{- Find integers (q,r) so that a = qb + r and 0 <= r < |b|. -}

divAlg :: Integer -> Integer -> (Integer, Integer)

divAlg _ 0 = error "division by zero"

divAlg a b = (quot a b, rem a b)

{- Find integers (x,y) so that ax + by = gcd(a,b). -}

bezout :: Integer -> Integer -> (Integer, Integer)

bezout a b = ((signum a)*x,(signum b)*y)

 where (x,y) = bezout' (abs a) (abs b)

        where bezout' _ 0 = (1,0)

              bezout' _ 1 = (0,1)

              bezout' a b = (y, x - q*y)

               where (x,y) = bezout' b r

                     (q,r) = divAlg a b

The bezout function went through several iterations while I picked out bugs. To give us some confidence that this function does what we think it does, we can run a simple test.

{- Verify the correctness of bezout on some input a b. -}

bezoutTest :: Integer -> Integer -> Bool

bezoutTest a b = (a*x) + (b*y) == d

 where d = gcd a b; (x,y) = bezout a b

Without going into too much detail, Haskell has a library called QuickCheck that takes a function that returns either True or False, generates a bunch of random input, and makes sure that the function returns True for all the input. To run the tests, simply run one of the following commands in your favorite Haskell interpreter.

import Test.QuickCheck

{- Simple test on 100 random inputs. -}

quickCheck bezoutTest

{- Display each test input before checking it. -}

verboseCheck bezoutTest

{- With this command you can supply the number of tests to run. -}

check (defaultConfig { configMaxTest = 10000}) bezoutTest

Algebra Solutions

Program to compute the greatest common divisor of two integers

No comments:

Post a Comment