Solution of Diophantine equation ax + by = c, by means of Euclidean algorithm

function [d,p,q,r,s] = euklidalg(a,b,c);
% Solution of Diophantine equation ax + by = c
% by means of Euclidean algorithm
% The greatest common divisor is the last nonzero remainder
% Back substitution:
% coefficients of the linear combination p, q
% r = -b/d, s = a/d
%
% test of equation correctness:
if isempty(find(a)) | isempty(find(b))
    display('Diophantine equation is not given correctly!')
    p=0; q=0; d=0; s=0; r=0;
    return
end
alfa=a; beta=b; rem=a;
% search for the greatest common divisor:
i=0;
n=abs(length(a)-length(b))+1;
if n == 1
    n=n+1;
end
% test of zero remainder:
while norm(rem,inf) > 100*eps
    i=i+1;
    % elimination the zero leading coefficients:
    ind=find(beta);
    beta=beta(ind(1):length(beta));
    % quotient and remainder:
    [quot,rem]=deconv(alfa,beta);
    i0=1+n-length(quot);
    % storing of quotients:
    qq(i,i0:n)=quot;
    % shift of polynomials:
    alfa=beta; beta=rem;
end
% recurrent computation of the coefficients
d=alfa; p=0; q=1; m=i-1;
for i=m:-1:1
    r=p; p=q
    % formal rearrangement for polynomial sum executing:
    rr=zeros(1,length(qq(i,:))+length(p)-1);
    rr(length(rr)-length(r)+1:length(rr))=r;
    % computation of further element of the sequence:
    q=rr-conv(qq(i,:),p);
end
% normalization of polynomial:
ind=find(q); q=q(ind(1):length(q));
% general solution of the reduced equation:
r=-deconv(b,d); s=deconv(a,d);
return

I want solve this:so my insert this:

a = [1, -1];
b = [0, 1, -2];
c = 1;

the result:

p =  1
d =  1
p =   1
q =   -1
r =     0    -1     2

but the correct is:

Have you any idea? or some Euclidean algorithm to solve diophatine.

[Torsten]

I do not even understand your question ... Could you be a bit more precise in explaining the problem ?

Given an integer (?) value for z, you want to find all integer values for x and y such that

(1-z^(-1))*x + z^(-1)*(1-2*z^(-1))*y = 1

holds ?

[Marek Hutta]

I want solve diophatine equation by Euclidean algorithm. I will get values for controller with algebraic theory.

[Torsten]

Is my interpretation of what you want correct or not ?

Given an integer (?) value for z, you want to find all integer values for x and y such that

(1-z^(-1))*x + z^(-1)*(1-2*z^(-1))*y = 1

holds ?

[Marek Hutta]

Yes.

[Torsten]

And how should the Euclean algorithm work for general z-dependent factors (1-z^(-1)) and z^(-1)*(1-2*z^(-1)) ?

[Marek Hutta]

yes

[Torsten]

This was a question. How ? I only know the Euclidean algorithm to work with explicitly given natural numbers, not with parametrized expressions in a variable z. Can you link to such an application ?

Answer 1

[Hassaan]

@Marek Hutta With a = [1, -1] and b = [0, 1, -2], it's unclear how these relate as a system of equations because of the differing lengths.

-----------------------------------------------------------------------------------------------------------------------------------------------------

If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.

It's important to note that the advice and code are based on limited information and meant for educational purposes. Users should verify and adapt the code to their specific needs, ensuring compatibility and adherence to ethical standards.

Professional Interests

• Technical Services and Consulting
• Embedded Systems | Firmware Developement | Simulations
• Electrical and Electronics Engineering

Feel free to contact me.

[Marek Hutta]

And i use my values:

a = [1, -1];
b = [0, 1, -2];
c = 1;

.

a = [1,-1];
b = [0,1,-2];
c = 1;
[x, y] = solveDiophantine(a, b, c);
fprintf('Particular solution: x = %d, y = %d\n', x, y);

Arrays have incompatible sizes for this operation.

Error in solveDiophantine>gcdExtended (line 28)
        [g, x1, y1] = gcdExtended(b, mod(a, b));

Error in solveDiophantine (line 5)
    [g, p, q] = gcdExtended(a, b);

Related documentation

[Marek Hutta]

and how i define it? like matrix with same dimension?

Answer 2

[Andika]

Hi, Marek.

If you have Symbolic Math Toolbox installed, you can use the gcd function with 3 output arguments. The syntax [G,C,D] = gcd(A,B,X) finds the greatest common divisor of A and B, and also returns the Bézout coefficients, C and D, such that G = A*C + B*D. Here, X is the free parameter for the polynomials in A and B. Note that the polynomials must be linear combinations of positive integer powers.

For your specific example, you can try the code below.

syms z zp
a = 1 - zp;         % zp represents 1/z
b = (1-2*zp)*zp;
c = 1;

[g,u,v] = gcd(a,b,zp);

if mod(c,g) == 0    % check if the gcd g divides c, or c/gcd(a,b) must be integer
    cp = c/g;
    disp('Particular solution:')
    x1 = u*cp       % one of the particular solution
    y1 = v*cp       % one of the particular solution
    
    disp('General solution:')
    syms k integer
    x = x1 + b/g*k      % find the general solution
    y = y1 - a/g*k      % find the general solution
    
    % Confirm the general solution is correct
    tf = isAlways(simplify(a*x + b*y) == c)

    % Substitute zp to 1/z
    x = subs(x,zp,1/z)
    y = subs(y,zp,1/z)
else
    disp('This Diophantine equation does not have a solution.')
end

Particular solution:

x1 = 

y1 = 

General solution:

x = 

y = 

tf = logical
   1

x = 

y =