|
1 | 1 | from __future__ import annotations |
2 | 2 |
|
3 | | -from maths.greatest_common_divisor import greatest_common_divisor |
| 3 | +def gcd(a: int, b: int) -> int: |
| 4 | + while b: |
| 5 | + a, b = b, a % b |
| 6 | + return a |
4 | 7 |
|
| 8 | +def extended_gcd(a: int, b: int) -> tuple[int, int, int]: |
| 9 | + if b == 0: |
| 10 | + return a, 1, 0 |
| 11 | + d, x, y = extended_gcd(b, a % b) |
| 12 | + return d, y, x - (a // b) * y |
5 | 13 |
|
6 | 14 | def diophantine(a: int, b: int, c: int) -> tuple[float, float]: |
7 | | - """ |
8 | | - Diophantine Equation : Given integers a,b,c ( at least one of a and b != 0), the |
9 | | - diophantine equation a*x + b*y = c has a solution (where x and y are integers) |
10 | | - iff greatest_common_divisor(a,b) divides c. |
11 | | -
|
12 | | - GCD ( Greatest Common Divisor ) or HCF ( Highest Common Factor ) |
13 | | -
|
14 | | - >>> diophantine(10,6,14) |
15 | | - (-7.0, 14.0) |
16 | | -
|
17 | | - >>> diophantine(391,299,-69) |
18 | | - (9.0, -12.0) |
19 | | -
|
20 | | - But above equation has one more solution i.e., x = -4, y = 5. |
21 | | - That's why we need diophantine all solution function. |
22 | | -
|
23 | | - """ |
24 | | - |
25 | | - assert ( |
26 | | - c % greatest_common_divisor(a, b) == 0 |
27 | | - ) # greatest_common_divisor(a,b) is in maths directory |
28 | | - (d, x, y) = extended_gcd(a, b) # extended_gcd(a,b) function implemented below |
| 15 | + d = gcd(a, b) |
| 16 | + assert c % d == 0 |
| 17 | + x, y = extended_gcd(a, b)[1:] |
29 | 18 | r = c / d |
30 | | - return (r * x, r * y) |
31 | | - |
| 19 | + return r * x, r * y |
32 | 20 |
|
33 | 21 | def diophantine_all_soln(a: int, b: int, c: int, n: int = 2) -> None: |
34 | | - """ |
35 | | - Lemma : if n|ab and gcd(a,n) = 1, then n|b. |
36 | | -
|
37 | | - Finding All solutions of Diophantine Equations: |
38 | | -
|
39 | | - Theorem : Let gcd(a,b) = d, a = d*p, b = d*q. If (x0,y0) is a solution of |
40 | | - Diophantine Equation a*x + b*y = c. a*x0 + b*y0 = c, then all the |
41 | | - solutions have the form a(x0 + t*q) + b(y0 - t*p) = c, |
42 | | - where t is an arbitrary integer. |
43 | | -
|
44 | | - n is the number of solution you want, n = 2 by default |
45 | | -
|
46 | | - >>> diophantine_all_soln(10, 6, 14) |
47 | | - -7.0 14.0 |
48 | | - -4.0 9.0 |
49 | | -
|
50 | | - >>> diophantine_all_soln(10, 6, 14, 4) |
51 | | - -7.0 14.0 |
52 | | - -4.0 9.0 |
53 | | - -1.0 4.0 |
54 | | - 2.0 -1.0 |
55 | | -
|
56 | | - >>> diophantine_all_soln(391, 299, -69, n = 4) |
57 | | - 9.0 -12.0 |
58 | | - 22.0 -29.0 |
59 | | - 35.0 -46.0 |
60 | | - 48.0 -63.0 |
61 | | -
|
62 | | - """ |
63 | | - (x0, y0) = diophantine(a, b, c) # Initial value |
64 | | - d = greatest_common_divisor(a, b) |
65 | | - p = a // d |
66 | | - q = b // d |
67 | | - |
| 22 | + x0, y0 = diophantine(a, b, c) |
| 23 | + p, q = a // gcd(a, b), b // gcd(a, b) |
68 | 24 | for i in range(n): |
69 | | - x = x0 + i * q |
70 | | - y = y0 - i * p |
71 | | - print(x, y) |
72 | | - |
73 | | - |
74 | | -def extended_gcd(a: int, b: int) -> tuple[int, int, int]: |
75 | | - """ |
76 | | - Extended Euclid's Algorithm : If d divides a and b and d = a*x + b*y for integers |
77 | | - x and y, then d = gcd(a,b) |
78 | | -
|
79 | | - >>> extended_gcd(10, 6) |
80 | | - (2, -1, 2) |
81 | | -
|
82 | | - >>> extended_gcd(7, 5) |
83 | | - (1, -2, 3) |
84 | | -
|
85 | | - """ |
86 | | - assert a >= 0 |
87 | | - assert b >= 0 |
88 | | - |
89 | | - if b == 0: |
90 | | - d, x, y = a, 1, 0 |
91 | | - else: |
92 | | - (d, p, q) = extended_gcd(b, a % b) |
93 | | - x = q |
94 | | - y = p - q * (a // b) |
95 | | - |
96 | | - assert a % d == 0 |
97 | | - assert b % d == 0 |
98 | | - assert d == a * x + b * y |
99 | | - |
100 | | - return (d, x, y) |
101 | | - |
| 25 | + print(x0 + i * q, y0 - i * p) |
102 | 26 |
|
103 | 27 | if __name__ == "__main__": |
104 | 28 | from doctest import testmod |
105 | | - |
106 | | - testmod(name="diophantine", verbose=True) |
107 | | - testmod(name="diophantine_all_soln", verbose=True) |
108 | | - testmod(name="extended_gcd", verbose=True) |
109 | | - testmod(name="greatest_common_divisor", verbose=True) |
| 29 | + testmod() |
0 commit comments