AX + BY = C

[bplus]

For math fans:

Code: QB64: [Select]

1. _TITLE "Ax + By = C"
2. SCREEN _NEWIMAGE(1000, 600, 32)
3.  
4. ' Thanks BSpinoza, orig ref: http://jean-pierre.moreau.pagesperso-orange.fr/Basic/diophan_bas.txt
5. ' I am curious about the technique so I study while I translate to QB64.
6. ' 2018-08-16 OK seems to work now. But you sort of need a routine to translate the results of
7. ' Diophantian  x0, y0, p, q into equations for x, y when a solution is found.
8.  
9. '**********************************************
10. '* Solving a diophantian equation ax+by = c   *
11. '*      (a,b,c,x,y are integer numbers)       *
12. '* ------------------------------------------ *
13. '* Ref.: "Mathematiques en Turbo-Pascal By    *
14. '*        M. Ducamp and A. Reverchon (2),     *
15. '*        Eyrolles, Paris, 1988" [BIBLI 05].  *
16. '* ------------------------------------------ *
17. '* Sample run:                                *
18. '*                                            *
19. '* SOLVING IN Z EQUATION AX + BY = C          *
20. '*                                            *
21. '*    A = 3                                   *
22. '*    B = -2                                  *
23. '*    C = 7                                   *
24. '*                                            *
25. '* Solutions are:                             *
26. '*                                            *
27. '*   X =    1 +     2*K                       *
28. '*   Y =   -2 +     3*K                       *
29. '*                                            *
30. '*              BASIC version by J-P Moreau.  *
31. '*                   (www.jpmoreau.fr)        *
32. '**********************************************
33.  
34.  
35. DEFINT A-Z
36.  
37.  
38. WHILE NOT _KEYDOWN(27)
39.     CLS
40.     PRINT
41.     PRINT " SOLVING EQUATION A*X + B*Y = C"
42.     PRINT
43.     INPUT "   A = ", a
44.     INPUT "   B = ", b
45.     INPUT "   C = ", c
46.     IF a = 0 OR b = 0 THEN END
47.     solnFound = FALSE: x0 = 0: y0 = 0: p = 0: q = 0
48.     Diophantian a, b, c, x0, y0, p, q, solnFound
49.  
50.     'now to translate these numbers!
51.     IF p * q > 0 THEN s$ = "-": m = -1 ELSE s$ = "+": m = 1
52.  
53.     IF solnFound = -1 THEN
54.         PRINT
55.         PRINT " Solutions are:"
56.         PRINT
57.         PRINT "   X = "; x0; " + "; ABS(q); " * K"
58.         PRINT "   Y = "; y0; " "; s$; " "; ABS(p); " * K"
59.         PRINT: PRINT " Check: Plug-in 10 integers for K:": PRINT
60.         FOR k = -5 TO 5
61.             x = x0 + ABS(q) * k
62.             y = y0 + m * ABS(p) * k
63.             PRINT " For K = "; k;
64.             LOCATE CSRLIN, 20: PRINT " X, Y is "; x; " , "; y;
65.             LOCATE CSRLIN, 50: PRINT a; " * "; x; " + "; b; " * "; y; " = "; a * x + b * y,
66.             LOCATE CSRLIN, 85: IF a * x + b * y = c THEN PRINT " Yes, C = "; c ELSE PRINT " No, C = "; c
67.         NEXT
68.     ELSE
69.         PRINT " No solutions."
70.     END IF
71.     PRINT: PRINT "Press spacebar to go again, escape to quit."
72.     H& = _KEYHIT
73.     WHILE H& <> -32
74.         H& = _KEYHIT
75.         IF _KEYDOWN(27) THEN END
76.         _LIMIT 100
77.     WEND
78. WEND
79.  
80. '*************************************************************************************************
81. '* Solving equation a*x + b*y = c, a, b, c, x, y are integer numbers
82. '* ------------------------------------------------------- ---------------------------------------
83. '* INPUT:   a,b,c  coefficients of equation   note: these numbers won't be changed by sub
84. '*          x0, y0, p, q, solutionFound all zero'd for output
85. '* OUTPUT:  solutions are x0+kp and y0-kq, with k=0,1,2...
86. '*          or k=-1,-2,-3...
87. '* The solutionFound returns -1 if solutions exist ie if the GCD of a,b is also a divisor of c.
88. '*****************************************************************************************************
89. SUB Diophantian (aa, bb, cc, x0, y0, p, q, solutionFound)
90.     solutionFound = 0
91.     IF aa = 0 THEN EXIT SUB: IF bb = 0 THEN EXIT SUB
92.     a = aa: b = bb: c = cc 'make copies as these will be changed
93.     pg = GCD(a, b)
94.     a = a / pg: b = b / pg: c! = c / pg
95.     IF c! <> INT(c!) THEN EXIT SUB ELSE c = INT(c!) 'pg must be also a divisor of c
96.     x1 = 0: y2 = 0
97.     ifound = 0
98.     WHILE ifound = 0
99.         y1! = (c - a * x1) / b
100.         IF y1! <> INT(y1!) THEN
101.             x1 = -1 * x1: IF x1 >= 0 THEN x1 = x1 + 1
102.             x2! = (c - b * y2) / a
103.             IF x2! <> INT(x2!) THEN
104.                 y2 = -1 * y2: IF y2 >= 0 THEN y2 = y2 + 1
105.             ELSE
106.                 x0 = INT(x2!): y0 = y2: ifound = -1
107.             END IF
108.         ELSE
109.             x0 = x1: y0 = INT(y1!)
110.             ifound = -1
111.         END IF
112.     WEND
113.     p = a: q = b
114.     solutionFound = -1
115. END SUB
116.  
117. FUNCTION GCD (a, b)
118.     c = ABS(a): d = ABS(b) 'make copies as these will be changed
119.     IF c = 0 OR d = 0 THEN GCD = -1: EXIT FUNCTION 'signal error?
120.     WHILE c <> 0 AND d <> 0
121.         IF c > d THEN c = c MOD d ELSE d = d MOD c
122.     WEND
123.     GCD = c + d
124. END FUNCTION
125.