%Reference R. Sotiriadis, "Diophantine Frequency Synthesis", IEEE Trans.
%    Ultrsonics, Ferro. & Freq. Control, Nov. 2006, pp. 1988-1998
%    Synthesize a/[R(1)* R(2)] = x(1)/R(1) + x(2)/R(2).  Get x's given a.
clear;
g = 0;
while g ~= 1
    R(1) =input('\ninput first reference divider ratio: ');
    R(2) =input('input second reference divider ratio: ');
    [g,z2,z1]=gcd(R(1),R(2)); %Greatest common divisor (g)should be 1.  R(1)*z2 + R(2)*z1 = g
    if g ~= 1
        fprintf('Divider ratios both divisible by %g.  Re-enter\n',g)
    end %if
end %while
%  Particular solution, Eq. (18)<- in reference
fprintf('R(1)=%g, R(2)=%g, g=%g, z(1)=%g, z(2)=%g\n',R(1),R(2),g,z1,z2);
prod = R(1)*R(2);
a=0;
%  Solution of x(1)/R(1) + x(2)/R(2) = a/[R(1)*R(2)] for any a, method under Eq. (18) in reference.
while a<= prod && a>=-prod
fprintf('input a between %g and %g: ',-prod,prod);
a=input('');
if a>prod || a<-prod
    break
elseif a == -prod
    fprintf('a = %g, x(1)=%g, x(2)=0\n',-prod,-R(1));
elseif a == prod
    fprintf('a = %g, x(1)=%g, x(2)=0\n',prod,R(1));
else
    y(1) = mod(a*z1,R(1));
    y(2) = mod(a*z2,R(2));
    q = y(1)/R(1) + y(2)/R(2) - a/prod;
    for k = 1:2
        if k<=q+.1
            x(k) = y(k)-R(k);
        else
            x(k) = y(k);
        end
    end
    %For Test -> fprintf('a = %g, y(1) = %g, y(2) = %g, q = %g, x(1) = %g, x(2} = %g\n',a,y(1),y(2),q, x(1),x(2))
    fprintf('a = %g, x(1) = %g, x(2} = %g\n',a,x(1),x(2)) 
    fprintf('a = %g to verify \n',x(1)*R(2)+ x(2)*R(1))
end
end