NYCJUG/linearDiophantineEquations

Solving Linear Diophantine equations

A Linear Diophantine equation has the form $ax+by=c$ where the solutions (a, b, c) are restricted to positive integers. The linear form of the Diophantine equation potentially has infinitely many solutions.

The solver following has been copied directly from an old issue of APL Quote-Quad (APL Quote-Quad: the Early Years, APL Press, Palo Alto, 1982, pp. 448-449). It was submitted by Raymond O. Folse, an assistant professor of Computer Science at Nicholls State University, Thibodaux, Louisiana, USA. Images of the original submission follow.

Here is a J translation File:DiophantineSolver.ijs :

NB.* diophantineSolver: solve Diophantine equations.
NB. Algorithm from APL Quote-Quad
NB. Raymond O. Folse, Asst. Prof. of Computer Science, Nicholls State
NB. University, Thibodaux, Louisiana 70301 U.S.A.

SHOWERRMSG=: 1
NB.* diophantineSolve: determine X and Y satisfying Diophantine equation
NB. (A1*X)+(A2Y)=A3 where A1, A2, and A3 are positive integers.  The values
NB. of X, Y are expressed as functions of T, where T is an integer.
diophantineSolve=: 3 : 0
   if. 3~:#y do.
       if. SHOWERRMSG do.
           smoutput 'Must enter 3 coefficients A,B,C for AX + BY = C.'
       end.
       return.
   end.
   k=. egcd 2{.y           NB. GCD of A & B, A & B%gcd
   if. 0~:r=. (0{k)|2{y do.
       if. SHOWERRMSG do.
           smoutput 'No solution for ','.',~":y
       end.
       return.
   end.
   xx=. (q*1{k),(2{k)*q=. <.(2{y)%0{k
   xx,.1 _1*(1 0{y)%0{k
)

NB.* egcd: compute Greatest Common Divisor Z1 of A, B and integers Z2, Z3
NB. such that Z1 = (A*Z2)+(B*Z3).
egcd=: 3 : 0
   'a b'=. y
   m=. a,1 0
   n=. b,0 1
   while. 0~:0{n do. 'm n'=. n;m-n*<.(0{m)%0{n end.
   m
)

Note elimination of awkward temporary by use of multiple assignment in "while" loop.

Even better J versions of this Euclidean GCD solver can be found here. And, of course, J has a primitive for gcd: +..

diophantineExplain=: 0 : 0
   ]soln=. diophantineSolve 7 3 100     NB. Equation 7X + 3Y = 100 has solution
 100  3                                 NB.  X =  100 + 3T
_200 _7                                 NB.  Y = _200 - 7T
   7 3 +/ . * soln +/ . * 1 1           NB. EG for T=1:
100
   7 3 +/ . * soln +/ . * 1 2           NB. EG for T=2:
100
   7 3 +/ . * soln +/ . * 1 3           NB. EG for T=3:
100
   7 3 +/ . * soln +/ . * 1 99          NB. EG for T=99:
100

NB. More examples:

   ]soln=. diophantineSolve 2 3 5
_5  3
 5 _2
   (<2 3)+/ . *&>(<soln)+/ . *&.>1,&.>>:i.10      NB. EG for 1-10
5 5 5 5 5 5 5 5 5 5

   ]soln=. diophantineSolve 4 6 7                 NB. Unsolvable
No solution for 4 6 7.
   $soln
0 0

   ]soln=. diophantineSolve 18 15 12
 4  5
_4 _6
   (<18 15)+/ . *&>(<soln)+/ . *&.>1,&.>i:10      NB. EG for _10 to 10
12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
)

findSolvableDEs=: 3 : 0
   'npi ni mf'=. y    NB. #/iteration, # iterations, max factor
   goodeq=. 0 2$a:
   SHOWERRMSG=: 0
   while. 0<:ni=. <:ni do.
       nn=. <"1 ?mf$~npi,3
       ss=. diophantineSolve&.>nn
       issolved=. (0~:#&>ss)*.2=#&>$&.>ss
       goodeq=. goodeq,issolved#ss,.nn
   end.
   goodeq
)

You can make an interesting picture with the complete solution set over >:i.100: