-
-
Notifications
You must be signed in to change notification settings - Fork 47k
Chinese Remainder Theorem | Diophantine Equation | Modular Division #1248
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Changes from 6 commits
f810adf
2b7a31e
41a9414
72e9db4
33de1b1
df00bf0
2b5c654
2492b52
8037edb
7241693
67f4fee
b241701
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,70 @@ | ||
# Chinese Remainder Theorem: | ||
|
||
# If GCD(a,b) = 1, then for any remainder ra modulo a and any remainder rb modulo b there exists integer n, | ||
# such that n = ra (mod a) and n = ra(mod b). If n1 and n2 are two such integers, then n1=n2(mod ab) | ||
|
||
# Algorithm : | ||
|
||
# 1. Use extended euclid algorithm to find x,y such that a*x + b*y = 1 | ||
# 2. Take n = ra*by + rb*ax | ||
|
||
# import testmod for testing our function | ||
from doctest import testmod | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Let's move this line down just below line 68... |
||
|
||
|
||
# Extended Euclid | ||
def extended_euclid(a, b): | ||
if b == 0: | ||
return (1, 0) | ||
(x, y) = extended_euclid(b, a % b) | ||
k = a // b | ||
return (y, x - k * y) | ||
|
||
|
||
# Uses ExtendedEuclid to find inverses | ||
def chinese_remainder_theorem(n1, r1, n2, r2): | ||
""" | ||
>>> chinese_remainder_theorem(5,1,7,3) | ||
31 | ||
|
||
Explanation : 31 is the smallest number such that | ||
(i) When we divide it by 5, we get remainder 1 | ||
(ii) When we divide it by 7, we get remainder 3 | ||
|
||
>>> chinese_remainder_theorem(6,1,4,3) | ||
14 | ||
""" | ||
(x, y) = extended_euclid(n1, n2) | ||
m = n1 * n2 | ||
n = r2 * x * n1 + r1 * y * n2 | ||
return ((n % m + m) % m) | ||
|
||
|
||
# ----------SAME SOLUTION USING InvertModulo instead ExtendedEuclid---------------- | ||
|
||
# This function find the inverses of a i.e., a^(-1) | ||
def invert_modulo(a, n): | ||
(b, x) = extended_euclid(a, n) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. No doctest? |
||
if b < 0: | ||
b = (b % n + n) % n | ||
return b | ||
|
||
|
||
# Same a above using InvertingModulo | ||
def chinese_remainder_theorem2(n1, r1, n2, r2): | ||
""" | ||
>>> chinese_remainder_theorem2(5,1,7,3) | ||
31 | ||
|
||
>>> chinese_remainder_theorem2(6,1,4,3) | ||
14 | ||
""" | ||
x, y = invert_modulo(n1, n2), invert_modulo(n2, n1) | ||
m = n1 * n2 | ||
n = r2 * x * n1 + r1 * y * n2 | ||
return (n % m + m) % m | ||
|
||
|
||
if __name__ == '__main__': | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. from doctest import testmod |
||
testmod(name='chinese_remainder_theorem', verbose=True) | ||
testmod(name='chinese_remainder_theorem2', verbose=True) |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,102 @@ | ||
# Diophantine Equation : Given integers a,b,c ( at least one of a and b != 0), the diophantine equation | ||
# a*x + b*y = c has a solution (where x and y are integers) iff gcd(a,b) divides c. | ||
|
||
# import testmod for testing our function | ||
from doctest import testmod | ||
|
||
|
||
def diophantine(a, b, c): | ||
""" | ||
>>> diophantine(10,6,14) | ||
(-7.0, 14.0) | ||
|
||
>>> diophantine(391,299,-69) | ||
(9.0, -12.0) | ||
|
||
But above equation has one more solution i.e., x = -4, y = 5. | ||
That's why we need diophantine all solution function. | ||
|
||
""" | ||
|
||
assert c % gcd(a, b) == 0 # gcd(a,b) function implemented below | ||
|
||
(d, x, y) = extended_gcd(a, b) # extended_gcd(a,b) function implemented below | ||
r = c / d | ||
return (r * x, r * y) | ||
|
||
|
||
# Lemma : if n|ab and gcd(a,n) = 1, then n|b. | ||
|
||
# Finding All solutions of Diophantine Equations: | ||
|
||
# Theorem : Let gcd(a,b) = d, a = d*p, b = d*q. If (x0,y0) is a solution of Diophantine Equation a*x + b*y = c. | ||
# a*x0 + b*y0 = c, then all the solutions have the form a(x0 + t*q) + b(y0 - t*p) = c, where t is an arbitrary integer. | ||
|
||
# n is the number of solution you want, n = 2 by default | ||
|
||
def diophantine_all_soln(a, b, c, n=2): | ||
""" | ||
>>> diophantine_all_soln(10, 6, 14) | ||
-7.0 14.0 | ||
-4.0 9.0 | ||
|
||
>>> diophantine_all_soln(10, 6, 14, 4) | ||
-7.0 14.0 | ||
-4.0 9.0 | ||
-1.0 4.0 | ||
2.0 -1.0 | ||
|
||
>>> diophantine_all_soln(391, 299, -69, n = 4) | ||
9.0 -12.0 | ||
22.0 -29.0 | ||
35.0 -46.0 | ||
48.0 -63.0 | ||
|
||
""" | ||
(x0, y0) = diophantine(a, b, c) | ||
d = gcd(a, b) | ||
p = a // d | ||
q = b // d | ||
|
||
for i in range(n): | ||
x = x0 + i * q | ||
y = y0 - i * p | ||
print(x, y) | ||
|
||
|
||
# Euclid's Lemma : d divides a and b, if and only if d divides a-b and b | ||
|
||
# Euclid's Algorithm | ||
|
||
def gcd(a, b): | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Please expand the acronym gcd |
||
if a < b: | ||
a, b = b, a | ||
|
||
while a % b != 0: | ||
a, b = b, a % b | ||
|
||
return b | ||
|
||
|
||
# Extended Euclid's Algorithm : If d divides a and b and d = a*x + b*y for integers x and y, then d = gcd(a,b) | ||
|
||
|
||
def extended_gcd(a, b): | ||
assert a >= 0 and b >= 0 | ||
|
||
if b == 0: | ||
d, x, y = a, 1, 0 | ||
else: | ||
(d, p, q) = extended_gcd(b, a % b) | ||
x = q | ||
y = p - q * (a // b) | ||
|
||
assert a % d == 0 and b % d == 0 | ||
assert d == a * x + b * y | ||
|
||
return (d, x, y) | ||
|
||
|
||
if __name__ == '__main__': | ||
testmod(name='diophantine', verbose=True) | ||
testmod(name='diophantine_all_soln', verbose=True) |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,106 @@ | ||
# Modular Division | ||
# An efficient algorithm for dividing b by a modulo n. | ||
|
||
# Given three integers a, b, and n, such that gcd(a,n)=1 and n>1, the algorithm should return an integer x such that | ||
# 0≤x≤n−1, and b/a=x(modn) (that is, b=ax(modn)). | ||
|
||
# Theorem: | ||
# a has a multiplicative inverse modulo n iff gcd(a,n) = 1 | ||
|
||
|
||
# This find x = b*a^(-1) mod n | ||
# Uses ExtendedEuclid to find the inverse of a | ||
|
||
# Import testmod for testing our function | ||
from doctest import testmod | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Move down to main() |
||
|
||
|
||
def modular_division(a, b, n): | ||
""" | ||
>>> modular_division(4,8,5) | ||
2 | ||
|
||
>>> modular_division(3,8,5) | ||
1 | ||
|
||
>>> modular_division(4, 11, 5) | ||
4 | ||
|
||
""" | ||
assert n > 1 and a > 0 and gcd(a, n) == 1 | ||
(d, t, s) = extended_gcd(n, a) # Implemented below | ||
x = (b * s) % n | ||
return x | ||
|
||
|
||
# This function find the inverses of a i.e., a^(-1) | ||
def invert_modulo(a, n): | ||
(b, x) = extended_euclid(a, n) # Implemented below | ||
if b < 0: | ||
b = (b % n + n) % n | ||
return b | ||
|
||
|
||
# ------------------ Finding Modular division using invert_modulo ------------------- | ||
|
||
# This function used the above inversion of a to find x = (b*a^(-1))mod n | ||
def modular_division2(a, b, n): | ||
""" | ||
>>> modular_division2(4,8,5) | ||
2 | ||
|
||
>>> modular_division2(3,8,5) | ||
1 | ||
|
||
>>> modular_division2(4, 11, 5) | ||
4 | ||
|
||
""" | ||
s = invert_modulo(a, n) | ||
x = (b * s) % n | ||
return x | ||
|
||
|
||
# Extended Euclid's Algorithm : If d divides a and b and d = a*x + b*y for integers x and y, then d = gcd(a,b) | ||
|
||
|
||
def extended_gcd(a, b): | ||
assert a >= 0 and b >= 0 | ||
|
||
if b == 0: | ||
d, x, y = a, 1, 0 | ||
else: | ||
(d, p, q) = extended_gcd(b, a % b) | ||
x = q | ||
y = p - q * (a // b) | ||
|
||
assert a % d == 0 and b % d == 0 | ||
assert d == a * x + b * y | ||
|
||
return (d, x, y) | ||
|
||
|
||
# Extended Euclid | ||
def extended_euclid(a, b): | ||
if b == 0: | ||
return (1, 0) | ||
(x, y) = extended_euclid(b, a % b) | ||
k = a // b | ||
return (y, x - k * y) | ||
|
||
# Euclid's Lemma : d divides a and b, if and only if d divides a-b and b | ||
# Euclid's Algorithm | ||
|
||
def gcd(a, b): | ||
if a < b: | ||
a, b = b, a | ||
|
||
while a % b != 0: | ||
a, b = b, a % b | ||
|
||
return b | ||
|
||
|
||
if __name__ == '__main__': | ||
testmod(name='modular_division', verbose=True) | ||
testmod(name='modular_division2', verbose=True) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Please expand the acronym GCD and try to keep the lines to 88 characters or less.