1+ # Extended Euclidean Algorithm
2+ #
3+ # The Extended Euclidean Algorithm not only finds the Greatest Common Divisor (GCD)
4+ # of two integers a and b, but also finds integers x and y such that:
5+ # ax + by = gcd(a, b) (Bézout's identity)
6+ #
7+ # This is particularly useful in modular arithmetic, RSA cryptography, and
8+ # finding modular multiplicative inverses.
9+ #
10+ # Time Complexity: O(log(min(a, b)))
11+ # Space Complexity: O(log(min(a, b))) due to recursion, O(1) for iterative version
12+ #
13+ # Input: Two integers a and b
14+ # Output: A list containing gcd, and coefficients x, y such that ax + by = gcd(a, b)
15+
16+ # Recursive implementation of Extended Euclidean Algorithm
17+ extended_gcd_recursive <- function (a , b ) {
18+ # Base case
19+ if (b == 0 ) {
20+ return (list (
21+ gcd = abs(a ),
22+ x = sign(a ), # Coefficient for a
23+ y = 0 # Coefficient for b
24+ ))
25+ }
26+
27+ # Recursive call
28+ result <- extended_gcd_recursive(b , a %% b )
29+
30+ # Update coefficients using the relation:
31+ # gcd(a, b) = gcd(b, a mod b)
32+ # If gcd(b, a mod b) = x1*b + y1*(a mod b)
33+ # Then gcd(a, b) = y1*a + (x1 - floor(a/b)*y1)*b
34+ x <- result $ y
35+ y <- result $ x - (a %/% b ) * result $ y
36+
37+ return (list (
38+ gcd = result $ gcd ,
39+ x = x ,
40+ y = y
41+ ))
42+ }
43+
44+ # Iterative implementation of Extended Euclidean Algorithm
45+ extended_gcd_iterative <- function (a , b ) {
46+ # Store original values for final adjustment
47+ orig_a <- a
48+ orig_b <- b
49+
50+ # Initialize coefficients
51+ old_x <- 1 ; x <- 0
52+ old_y <- 0 ; y <- 1
53+
54+ while (b != 0 ) {
55+ quotient <- a %/% b
56+
57+ # Update a and b
58+ temp <- b
59+ b <- a %% b
60+ a <- temp
61+
62+ # Update x coefficients
63+ temp <- x
64+ x <- old_x - quotient * x
65+ old_x <- temp
66+
67+ # Update y coefficients
68+ temp <- y
69+ y <- old_y - quotient * y
70+ old_y <- temp
71+ }
72+
73+ # Adjust signs based on original inputs
74+ if (orig_a < 0 ) {
75+ a <- - a
76+ old_x <- - old_x
77+ }
78+ if (orig_b < 0 ) {
79+ old_y <- - old_y
80+ }
81+
82+ return (list (
83+ gcd = abs(a ),
84+ x = old_x ,
85+ y = old_y
86+ ))
87+ }
88+
89+ # Function to find modular multiplicative inverse
90+ modular_inverse <- function (a , m ) {
91+ # ' Find modular multiplicative inverse of a modulo m
92+ # ' Returns x such that (a * x) ≡ 1 (mod m)
93+ # ' Only exists if gcd(a, m) = 1
94+
95+ result <- extended_gcd_iterative(a , m )
96+
97+ if (result $ gcd != 1 ) {
98+ return (NULL ) # Inverse doesn't exist
99+ }
100+
101+ # Make sure the result is positive
102+ inverse <- result $ x %% m
103+ if (inverse < 0 ) {
104+ inverse <- inverse + m
105+ }
106+
107+ return (inverse )
108+ }
109+
110+ # Function to solve linear Diophantine equation ax + by = c
111+ solve_diophantine <- function (a , b , c ) {
112+ # ' Solve the linear Diophantine equation ax + by = c
113+ # ' Returns NULL if no integer solutions exist
114+ # ' Returns one particular solution and the general solution pattern
115+
116+ result <- extended_gcd_iterative(a , b )
117+ gcd_ab <- result $ gcd
118+
119+ # Check if solution exists
120+ if (c %% gcd_ab != 0 ) {
121+ return (NULL ) # No integer solutions exist
122+ }
123+
124+ # Scale the coefficients
125+ scale <- c / gcd_ab
126+ x0 <- result $ x * scale
127+ y0 <- result $ y * scale
128+
129+ return (list (
130+ particular_solution = list (x = x0 , y = y0 ),
131+ general_solution = list (
132+ x_formula = paste0(x0 , " + " , b / gcd_ab , " *t" ),
133+ y_formula = paste0(y0 , " - " , a / gcd_ab , " *t" ),
134+ description = " where t is any integer"
135+ ),
136+ verification = a * x0 + b * y0 == c
137+ ))
138+ }
139+
140+ # Function to find all solutions in a given range
141+ find_diophantine_solutions_in_range <- function (a , b , c , x_min , x_max , y_min , y_max ) {
142+ # ' Find all integer solutions to ax + by = c in the given ranges
143+
144+ dioph_result <- solve_diophantine(a , b , c )
145+ if (is.null(dioph_result )) {
146+ return (NULL )
147+ }
148+
149+ x0 <- dioph_result $ particular_solution $ x
150+ y0 <- dioph_result $ particular_solution $ y
151+ gcd_ab <- extended_gcd_iterative(a , b )$ gcd
152+
153+ b_coeff <- b / gcd_ab
154+ a_coeff <- a / gcd_ab
155+
156+ # Find range of t values
157+ t_min_x <- ceiling((x_min - x0 ) / b_coeff )
158+ t_max_x <- floor((x_max - x0 ) / b_coeff )
159+ t_min_y <- ceiling((y0 - y_max ) / a_coeff )
160+ t_max_y <- floor((y0 - y_min ) / a_coeff )
161+
162+ t_min <- max(t_min_x , t_min_y )
163+ t_max <- min(t_max_x , t_max_y )
164+
165+ if (t_min > t_max ) {
166+ return (data.frame (x = integer(0 ), y = integer(0 ), t = integer(0 )))
167+ }
168+
169+ # Generate solutions
170+ t_values <- t_min : t_max
171+ x_values <- x0 + b_coeff * t_values
172+ y_values <- y0 - a_coeff * t_values
173+
174+ return (data.frame (x = x_values , y = y_values , t = t_values ))
175+ }
176+
177+ # Example usage and testing
178+ cat(" === Extended Euclidean Algorithm ===\n " )
179+
180+ # Test basic extended GCD
181+ cat(" Testing Extended GCD with a=240, b=46:\n " )
182+ result1 <- extended_gcd_iterative(240 , 46 )
183+ cat(" GCD:" , result1 $ gcd , " \n " )
184+ cat(" Coefficients: x =" , result1 $ x , " , y =" , result1 $ y , " \n " )
185+ cat(" Verification:" , 240 * result1 $ x + 46 * result1 $ y , " = GCD:" , result1 $ gcd , " \n " )
186+ cat(" Check:" , 240 * result1 $ x + 46 * result1 $ y == result1 $ gcd , " \n\n " )
187+
188+ # Compare recursive and iterative methods
189+ cat(" Comparing recursive vs iterative methods:\n " )
190+ result_rec <- extended_gcd_recursive(240 , 46 )
191+ result_iter <- extended_gcd_iterative(240 , 46 )
192+ cat(" Recursive - GCD:" , result_rec $ gcd , " , x:" , result_rec $ x , " , y:" , result_rec $ y , " \n " )
193+ cat(" Iterative - GCD:" , result_iter $ gcd , " , x:" , result_iter $ x , " , y:" , result_iter $ y , " \n " )
194+ cat(" Results match:" ,
195+ result_rec $ gcd == result_iter $ gcd &&
196+ result_rec $ x == result_iter $ x &&
197+ result_rec $ y == result_iter $ y , " \n\n " )
198+
199+ # Test modular multiplicative inverse
200+ cat(" === Modular Multiplicative Inverse ===\n " )
201+ cat(" Finding inverse of 7 modulo 26:\n " )
202+ inv <- modular_inverse(7 , 26 )
203+ if (! is.null(inv )) {
204+ cat(" 7^(-1) ≡" , inv , " (mod 26)\n " )
205+ cat(" Verification: 7 *" , inv , " ≡" , (7 * inv ) %% 26 , " (mod 26)\n " )
206+ } else {
207+ cat(" Inverse does not exist\n " )
208+ }
209+
210+ # Test case where inverse doesn't exist
211+ cat(" \n Testing case where inverse doesn't exist (a=6, m=9):\n " )
212+ inv2 <- modular_inverse(6 , 9 )
213+ if (is.null(inv2 )) {
214+ cat(" Inverse does not exist (as expected, since gcd(6,9) = 3 ≠ 1)\n " )
215+ }
216+
217+ # Test solving Diophantine equations
218+ cat(" \n === Linear Diophantine Equations ===\n " )
219+ cat(" Solving 25x + 9y = 7:\n " )
220+ dioph1 <- solve_diophantine(25 , 9 , 7 )
221+ if (! is.null(dioph1 )) {
222+ cat(" Particular solution: x =" , dioph1 $ particular_solution $ x ,
223+ " , y =" , dioph1 $ particular_solution $ y , " \n " )
224+ cat(" General solution:\n " )
225+ cat(" x =" , dioph1 $ general_solution $ x_formula , " \n " )
226+ cat(" y =" , dioph1 $ general_solution $ y_formula , " \n " )
227+ cat(" " , dioph1 $ general_solution $ description , " \n " )
228+ cat(" Verification:" , dioph1 $ verification , " \n " )
229+ } else {
230+ cat(" No integer solutions exist\n " )
231+ }
232+
233+ # Test equation with no solutions
234+ cat(" \n Solving 6x + 9y = 10 (should have no integer solutions):\n " )
235+ dioph2 <- solve_diophantine(6 , 9 , 10 )
236+ if (is.null(dioph2 )) {
237+ cat(" No integer solutions exist (as expected, since gcd(6,9)=3 does not divide 10)\n " )
238+ }
239+
240+ # Find solutions in a specific range
241+ cat(" \n === Solutions in Range ===\n " )
242+ cat(" Finding solutions to 25x + 9y = 7 where -10 ≤ x ≤ 10 and -10 ≤ y ≤ 10:\n " )
243+ range_solutions <- find_diophantine_solutions_in_range(25 , 9 , 7 , - 10 , 10 , - 10 , 10 )
244+ if (! is.null(range_solutions ) && nrow(range_solutions ) > 0 ) {
245+ print(range_solutions )
246+
247+ # Verify solutions
248+ cat(" Verification of solutions:\n " )
249+ for (i in 1 : nrow(range_solutions )) {
250+ x <- range_solutions $ x [i ]
251+ y <- range_solutions $ y [i ]
252+ result_check <- 25 * x + 9 * y
253+ cat(" 25 *" , x , " + 9 *" , y , " =" , result_check ,
254+ " (" , if (result_check == 7 ) " ✓" else " ✗" , " )\n " )
255+ }
256+ } else {
257+ cat(" No solutions found in the specified range\n " )
258+ }
259+
260+ # Applications example
261+ cat(" \n === Practical Applications ===\n " )
262+ cat(" Example: Making change with 3-cent and 5-cent coins\n " )
263+ cat(" Problem: Can we make exactly 14 cents? Find all ways.\n " )
264+ change_problem <- solve_diophantine(3 , 5 , 14 )
265+ if (! is.null(change_problem )) {
266+ cat(" Yes! Particular solution:" ,
267+ change_problem $ particular_solution $ x , " three-cent coins and" ,
268+ change_problem $ particular_solution $ y , " five-cent coins\n " )
269+
270+ # Find all non-negative solutions
271+ solutions_range <- find_diophantine_solutions_in_range(3 , 5 , 14 , 0 , 10 , 0 , 10 )
272+ if (nrow(solutions_range ) > 0 ) {
273+ cat(" All valid ways to make 14 cents:\n " )
274+ for (i in 1 : nrow(solutions_range )) {
275+ x <- solutions_range $ x [i ]
276+ y <- solutions_range $ y [i ]
277+ if (x > = 0 && y > = 0 ) {
278+ cat(x , " × 3¢ +" , y , " × 5¢ = 14¢\n " )
279+ }
280+ }
281+ }
282+ }
0 commit comments