'********************************************** '* Solving a diophantian equation ax+by = c * '* (a,b,c,x,y are integer numbers) * '* ------------------------------------------ * '* Ref.: "Mathematiques en Turbo-Pascal By * '* M. Ducamp and A. Reverchon (2), * '* Eyrolles, Paris, 1988" [BIBLI 05]. * '* ------------------------------------------ * '* Sample run: * '* * '* SOLVING IN Z EQUATION AX + BY = C * '* * '* A = 3 * '* B = -2 * '* C = 7 * '* * '* Solutions are: * '* * '* X = 1 + 2*K * '* Y = -2 + 3*K * '* * '* BASIC version by J-P Moreau. * '* (www.jpmoreau.fr) * '********************************************** TRUE=1 FALSE=0 CLS print print " SOLVING IN Z EQUATION AX + BY = C" print input " A = ",a input " B = ",b input " C = ",c print ' iresult = Diophantian(a,b,c,x0,y0,p,q) gosub 2000 if iresult=TRUE then print print " Solutions are:" print print using " X=###### +######"; x0, abs(q);:print " *K" if p*q>0 then print using " Y=###### -######"; y0, abs(p);:print " *K" else print using " Y=###### +######"; y0, abs(p);:print " *K" end if else print " No solutions." end if print print END 'return greatest common divisor of two integer numbers 'FUNCTION GCD(aa,bb) 1000 aa = int(abs(aa)) bb = int(abs(bb)) if aa>1e10 or bb>1e10 then GCD=1: return end if if aa=0 or bb=0 then GCD=1: return end if if aa < bb then temp=aa: aa=bb: bb=temp end if 1010 r=aa-bb*int(aa/bb) aa=bb: bb=r if abs(r)>1e-10 then goto 1010 GCD = aa return '*********************************************************** '* Solving equation ax+by=c, a,b,c,x,y are integer numbers * '* ------------------------------------------------------- * '* INPUT: a,b,c coefficients of equation * '* OUTPUT: solutions are x0+kp and y0-kq, with k=0,1,2... * '* or k=-1,-2,-3... * '* The function returns TRUE if solutions exist (that is, * '* if the GCD of a,b is also a divisor of c). * '*********************************************************** 'Function Diophantian(a,b,c,x0,y0,p,q) 2000 iresult=FALSE if a=0 then return : if b=0 then return aa=a: bb=b 'send copies of a and b to subroutine 1000 gosub 1000 : pg=GCD a=a/pg: b=b/pg: c=c/pg if c<>INT(c) then return 'pg must be also a divisor of c x1=0: y2=0 ifound=FALSE 2010 y1=(c-a*x1)/b if y1=INT(y1) then x0=x1 : y0=y1 ifound=TRUE else x1=-x1: if x1>=0 then x1=x1+1 x2=(c-b*y2)/a if x2=INT(x2) then x0=x2: y0=y2: ifound=TRUE else y2=-y2: if y2>=0 then y2=y2+1 end if end if if ifound=FALSE then goto 2010 p=a: q=b iresult=TRUE return 'end of file diophan.bas