%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=input('input reference divide vector, [R1 R2 ...] \n');
k=length(R);
fprintf('%g loops.  R = ',k)
for i = 1:k
    fprintf('%g   ',R(i))
end
%particular solution to Eq. (18), which forms basis for others.
%   Start procedure in Appendix B
x=zeros(1,k);y=zeros(1,k);z=zeros(1,k);
for i = 2:k
    [g,z(i),w(i)]=gcd(R(i),prod(R(1:i-1)));
    if g ~= 1
        fprintf('\nTwo divider ratios are divisible by same integer.  Re-enter\n')
        break
    end %if
end %for
end %while
v(1)=prod(z(2:k));
for r=2:k-1
    v(r)=prod(z(r+1:k))*w(r);
end
v(k)=w(k);
%    End procedure in Appendix B
for j=1:k
    fprintf('\nz(%g)=%g, ',j,v(j));
end
sum = 0;
for j=1:k
    sum = v(j)/R(j)+ sum;
end
prodp = prod(R);
uno = sum*prodp;
fprintf('\nsum of (z/R)''s * product of R''s is %g.  It should be 1. \n',uno)
%  Solution of x(1)/R(1) + x(2)/R(2) + ... = a/[R(1)*R(2)*...] for any a, algorithm
%      under Eq. (18).
a=0;
while a<= prodp && a>=-prodp
fprintf('\ninput a between %10.0f and %10.0f (outside range to quit): ',-prodp,prodp);
a=input('');
if a>prodp || a<-prodp
    break
elseif a == -prodp
    x = zeros(2:k); x(1) = -R(1);
    fprintf('x(1)=%g, others = 0\n',x(1));
elseif a == prodp
    x = zeros(2:k); x(1) = R(1);
    fprintf('x(1)=%g, others = 0\n',x(1));
else
    q=0;
    for i=1:k
    y(i) = mod(a*v(i),R(i));
    q = y(i)/R(i) + q ;
    end
    q = q - a/prodp;
    for j = 1:k
        if j<=q+.1
            x(j) = y(j)-R(j);
        else
            x(j) = y(j);
        end
        fprintf('x(%g) = %g\n',j,x(j))
    end
end
sum = 0;
for j=1:k
    sum = x(j)/R(j)+ sum;
end
uno = sum*prodp;
fprintf('sum of (x/R)''s * product of R''s is %10.0f.\n',uno)
end