iframe-proxy | Sunbelt Computer Software

History

131 lines (91 loc) · 2.97 KB

Raw

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

# 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.

# GCD ( Greatest Common Divisor ) or HCF ( Highest Common Factor )

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 % greatest_common_divisor(a, b) == 0

) # greatest_common_divisor(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) # Initial value

d = greatest_common_divisor(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 greatest_common_divisor(a, b):

"""

>>> greatest_common_divisor(7,5)

Note : In number theory, two integers a and b are said to be relatively prime,

mutually prime, or co-prime if the only positive integer (factor) that

divides both of them is 1 i.e., gcd(a,b) = 1.

>>> greatest_common_divisor(121, 11)

"""

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):

"""

>>> extended_gcd(10, 6)

(2, -1, 2)

>>> extended_gcd(7, 5)

(1, -2, 3)

"""

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__":

from doctest import testmod

testmod(name="diophantine", verbose=True)

testmod(name="diophantine_all_soln", verbose=True)

testmod(name="extended_gcd", verbose=True)

testmod(name="greatest_common_divisor", verbose=True)

Sunbelt Computer Software

PL/B Language Development and Support

Navigation Menu

Sunbelt Computer Software

PL/B Language Development and Support

FilesExpand file tree

diophantine_equation.py

Latest commit

History

diophantine_equation.py

File metadata and controls