/* Diophantine equation solver for AX + BY = C - May 1995 */ /* Euclidean algorithm and basic solution from Byte Magazine, March 1986 by Bob Kurosaka */ /* REXX version by Claude VIE, with final calculations added */ If Arg()>0 Then Parse Upper Arg arg1 '(' opt . Else arg1='' det = 0 /* Display details if 1 */ If arg1 = '?' | arg1 = '/?' Then Do Call help Exit End If arg1 == '' Then Do Say "Enter equation: ax + by = c (or ? or /?)" Say "for example, enter 154X + 69Y = 5000 or 154X - 69Y = 5000" Pull rep If rep = '?' | rep = '/?' Then Do Call help Exit End equation = Strip(rep) End Else equation = Strip(arg1) If equation = '' Then Exit Else Say "Equation is" equation /* parms verification */ neq=Countstr('=',equation); nx=Countstr('X',equation); ny=Countstr('Y',equation) If neq \= 1 | nx \= 1 | ny \= 1 | Countstr('(',equation) \= 0 | Countstr(')',equation) \= 0 Then Call err "Equation non valid, non unique number of X ("nx") or Y ("ny") or = ("neq") or unwanted parentheses" Parse var equation a_x 'X' b_y 'Y' fil '=' c_xy If a_x='' Then a_x=1; b_y=Strip(b_y); If b_y='+' | b_y='-' Then b_y=b_y'1' If \(Datatype(a_x,'W') & Datatype(b_y,'W') & Datatype(c_xy,'W')) | fil \= '' Then Call err 'Integer values not found or unknown filler found: a >'a_x'<' 'b >'b_y 'c >'c_xy'<' '>'fil'<' If a_x = 0 | b_y = 0 Then Call err "a or b is zero" a_x=a_x+0; b_y=b_y+0; c_xy=c_xy+0; /* format numbers */ /* End of verif */ /* Euclidean algorithm for finding GCD */ swap_xy = 0 /* initialize the terms for the algorithm */ If Abs(a_x) >= Abs(b_y) Then Do; dividend = a_x; divisor = b_y; End Else Do; dividend = b_y; divisor = a_x; swap_xy = 1;End /* 1st iteration */ quotient = dividend % divisor /* be sure to trunc, not round negatives */ If det Then Say "==> truncation (not rounding) of negatives:" quotient Trunc(dividend/divisor) remainder = dividend // divisor If det Then Do; r=dividend - divisor * quotient Say "==> remainder verification:" remainder r; End /* x1 is the ongoing count of x, y1 is the ongoing count of y. You can keep track of all ongoing counts by using only the previous two values for x or y, so we need only x1, x2, x3 and y1, y2, y3. */ x1 = 1; y1 = -quotient /* If either a or b is even multiple of the other, then either x or y will equal 1 while the other equals 0. */ If remainder = 0 Then Do; x2 = 0; y2 = 1; End /* End of 1st iteration */ Else Do /*$*/ /* 2nd iteration */ dividend = divisor; divisor = remainder quotient = dividend % divisor remainder = dividend // divisor x2 = - quotient * x1; y2 = 1 - quotient * y1 /* If a GCD is found on the second itereration of the euclidean algorithm, then x = x1, y = y1. In all subsequent cases, x = x2, y = y2 */ If remainder = 0 Then Do; x2 = x1; y2 = y1; End /* end of 2nd iteration */ /* The first two iterations are the only ones that do not follow the pattern: x(n) = x(n-2) - quotient * x(n-1) y(n) = y(n-2) - quotient * y(n-1) */ Else Do While remainder \= 0 /* 3rd iteration and more */ dividend = divisor; divisor = remainder quotient = dividend % divisor remainder = dividend // divisor If remainder = 0 Then leave x3 = x1 - quotient * x2; y3 = y1 - quotient * y2 x1 = x2; x2 = x3; y1 = y2; y2 = y3 End /* 3rd iteration and more */ End /*$*/ /* Calculate the basic solution for ax + by = c from GCD results, which have provided ax + by = d basic solution */ d = divisor; e = c_xy / d If c_xy/d \= c_xy%d Then Say "No integer solutions" Else Do If swap_xy Then Do; fil = x2; x2 = y2; y2 = fil; End If b_y>=0 Then equ = a_x'x' '+' b_y'y' '=' c_xy Else equ = a_x'x' '-' Abs(b_y)'y' '=' c_xy Say Say "The basic solution to the Diophantine equation," equ "is:" Say "x =" x2 * e " y =" y2 * e " (the GCD of" a_x "and" b_y "is" Abs(d)")." Say fil = "for all integer values 'n'." Say "The parametric equations for all integer answers are:" If b_y/d>0 Then fil2 = "x =" x2*e "+" b_y/d"n and" Else fil2 = "x =" x2*e b_y/d"n and" If a_x/d<0 Then Say fil2 "y =" y2*e "+" Abs(a_x/d)"n" fil Else Say fil2 "y =" y2*e "-" Abs(a_x/d)"n" fil Say /* Calculate the positive values for x and y */ nx0 = (-x2*e)/(b_y/d) If a_x/d<0 Then ny0 = (-y2*e)/(Abs(a_x/d)) Else ny0 = (y2*e)/(Abs(a_x/d)) If nx0 > ny0 Then Do; If det Then Say "==>" ny0 '<= N <=' nx0; i1 = ny0; i2 = nx0; End Else Do; If det Then Say "==>" nx0 '<= n <=' ny0; i1 = nx0; i2 = ny0; End /* Calculate now the integer values */ If i1>0 Then i1 = Trunc(i1-1); Else i1 = Trunc(i1-1) If i2>0 Then i2 = Trunc(i2+2); Else i2 = Trunc(i2+1) If det Then Say "==>" i1 '<= n <=' i2 If det Then Say Say "The positive answers are:" Say Right("n",8) Right("x",10) Right("y",10) Right("equation",40) Do i=i1 To i2 x = (x2 * e) + (b_y/d*i) y = (y2 * e) - (a_x/d*i) If b_y>=0 Then equ = '(' a_x 'x' x ') + (' b_y 'x' y ') =' c_xy Else equ = '(' a_x 'x' x ') - (' Abs(b_y) 'x' y ') =' c_xy If x>0 & y>0 Then Say Right(i,8) Right(x,10) Right(y,10) Right(equ,40) Else If det Then Say Right(i,8) Right(x,10) Right(y,10) "<==" End Say End Exit err: Parse arg msg Say msg Call help Exit Return help: Say "Syntax: DioEqSolver [equation]" Say Say "This program solves equations of the form AX + BY = C," Say "where A, B, C, X and Y are all integer values." Say "Enter your equation as shown in the general form." Say "For example, enter 154X + 69Y = 5000 or 154X - 69Y = 5000." Say "Do no enter negative coefficients with parentheses;" Say "that is, do NOT enter 154X + (-69Y) = 5000." Say Say "The program will not work properly for the degenerate case" Say "where either A or B is 0." Say Return