extended algorithm of Euclides and the modular inverse of an element

bz_euclides.go

// available at https://play.golang.org/p/COP6ErBkD47
package main

import (
	"fmt"
)

func extended_algorithm_of_euclides(a, b int64) (d, s, t int64) {
	var (
		rx = a
		ry = b
		sx = int64(1)
		sy = int64(0)
		tx = int64(0)
		ty = int64(1)
		q int64
	)
	for ry != 0 {
		// rx = ry*q + rz
		// => q = rx/ry   &&   rz = rx - ry*q
		// 1. compute q
		q = rx / ry
		// 2. compute rx, ry for this iteration
		rx, ry = ry, rx - ry*q
		// 3. compute sx, sy for this iteration
		sx, sy = sy, sx - sy*q
		// 3. compute tx, ty for this iteration
		tx, ty = ty, tx - ty*q
	}
	return rx, sx, tx
}

func test_extended_algorithm_of_euclides(a, b int64) {
	d, s, t := extended_algorithm_of_euclides(a, b)
	fmt.Printf("found that for a=%d and b=%d, that ggd(a,b)=%d, and that (s, t) = (%d, %d)\n", a, b, d, s, t)
	r := a*s + b*t
	if a%d != 0 {
		panic(fmt.Sprintf("expected d=%d to be a divisor of a=%d", d, a))
	}
	if b%d != 0 {
		panic(fmt.Sprintf("expected d=%d to be a divisor of b=%d", d, b))
	}
	if d != r {
		panic(fmt.Sprintf("d := %d != %d = %d*%d + %d*%d = a*s + b*t", d, r, a, s, b, t))
	}
}

func main() {
	test_extended_algorithm_of_euclides(126, 35)
	test_extended_algorithm_of_euclides(35, 126)

	// the algorithm still works if b is a divisor of a,
	// so not sure why the theorem in my text book suggests that it doesn't,
	// from my own analysis this is what I suspected and my Go implementation seems
	// to confirm it. Right now I am pretty puzzled...  
	test_extended_algorithm_of_euclides(8, 4)
	test_extended_algorithm_of_euclides(4, 8)
	
	// keeps working even with very big numbers
	test_extended_algorithm_of_euclides(252, 369)
	test_extended_algorithm_of_euclides(4390394392932135, 30459304949303034)
	
}

bze_mod_inverse.go

// available at https://play.golang.org/p/8G0Q4K_Q7iu
package main

import (
	"fmt"
)

func extended_algorithm_of_euclides(a, b int64) (d, s, t int64) {
	var (
		rx = a
		ry = b
		sx = int64(1)
		sy = int64(0)
		tx = int64(0)
		ty = int64(1)
		q int64
	)
	for ry != 0 {
		// rx = ry*q + rz
		// => q = rx/ry   &&   rz = rx - ry*q
		// 1. compute q
		q = rx / ry
		// 2. compute rx, ry for this iteration
		rx, ry = ry, rx - ry*q
		// 3. compute sx, sy for this iteration
		sx, sy = sy, sx - sy*q
		// 3. compute tx, ty for this iteration
		tx, ty = ty, tx - ty*q
	}
	return rx, sx, tx
}

func test_extended_algorithm_of_euclides(a, b int64) {
	d, s, t := extended_algorithm_of_euclides(a, b)
	fmt.Printf("found that for a=%d and b=%d, that ggd(a,b)=%d, and that (s, t) = (%d, %d)\n", a, b, d, s, t)
	r := a*s + b*t
	if a%d != 0 {
		panic(fmt.Sprintf("expected d=%d to be a divisor of a=%d", d, a))
	}
	if b%d != 0 {
		panic(fmt.Sprintf("expected d=%d to be a divisor of b=%d", d, b))
	}
	if d != r {
		panic(fmt.Sprintf("d := %d != %d = %d*%d + %d*%d = a*s + b*t", d, r, a, s, b, t))
	}
}

func modular_inverse(a, m int64) (int64, bool) {
	d, x, _ := extended_algorithm_of_euclides(a, m)
	if d != 1 {
		return 0, false
	}
	return x % m, true
}

func norm_modular_inverse(a, m int64) (int64, bool) {
	i, ok := modular_inverse(a, m)
	if !ok {
		return 0, false
	}
	if i > 0 {
		return i%m, true
	}
	return m + (i%m), true
}

func test_modular_inverse(a, m, expected int64) {
	i, ok := modular_inverse(a, m)
	if expected != 0 {
		if !ok {
			panic(fmt.Sprintf("expected an inverse for a=%d modulo m=%d but found none", a, m))
		}
		if i != expected && (m+i)%m != expected {
			panic(fmt.Sprintf("expected %d as inverse for a=%d modulo m=%d, but found %d as the inverse for a modulo m", expected, a, m, i))
		}
		fmt.Printf("confirmed that %d is the inverse of %d modulo %d\n", expected, a, m)
	} else if ok {
		panic(fmt.Sprintf("expected no inverse for a=%d modulo m=%d but found %d as inverse of a modulo m", a, m, i))
	} else {
		fmt.Printf("confirmed that %d modulo %d has no inverse\n", a, m)
	}
}

func main() {
	test_extended_algorithm_of_euclides(126, 35)
	test_extended_algorithm_of_euclides(35, 126)
	test_extended_algorithm_of_euclides(8, 4)
	test_extended_algorithm_of_euclides(4, 8)
	test_extended_algorithm_of_euclides(252, 369)
	test_extended_algorithm_of_euclides(4390394392932135, 30459304949303034)
	
	test_modular_inverse(1, 2, 1)
	test_modular_inverse(2, 2, 0)
	test_modular_inverse(1, 4, 1)
	test_modular_inverse(2, 4, 0)
	test_modular_inverse(3, 4, -1)
	test_modular_inverse(3, 4, 3)
	
	seen := map[int64]struct{}{}
	for i := int64(1); i < 37; i++ {
		ai, ok := norm_modular_inverse(i, 37)
		if !ok {
			panic(fmt.Sprintf("found no inverse for %d while one was expected modulo 37", i))
		}
		if ai <= 0 || ai >= 37 {
			panic(fmt.Sprintf("found an inverse %d for %d modulo 37, outside the allowed exclusive range of ]0, 37[", ai, i))
		}
		if aii, ok := seen[ai]; ok {
			panic(fmt.Sprintf("already found an inverse %d for %d modulo 37, %d is a duplicate ?!", aii, i, ai))
		}
		seen[ai] = struct{}{}
		fmt.Printf("%d is the inverse of %d module 37\n", ai, i)
	}
	if len(seen) != 36 {
		panic(fmt.Sprintf("found %d inverses modulo 37, while 36 was expected", len(seen)))
	}
	fmt.Println("Found 36 inverses for 36 elements modulo 37, as expected")
}