Path: news.mathworks.com!not-for-mail
From: "Bruno Luong" <b.luong@fogale.fr>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Largest Common Factor
Date: Thu, 1 May 2008 09:10:06 +0000 (UTC)
Organization: FOGALE nanotech
Lines: 53
Message-ID: <fvc1de$s6q$1@fred.mathworks.com>
References: <fv8j8a$2j8$1@fred.mathworks.com> <19640430.1209534092123.JavaMail.jakarta@nitrogen.mathforum.org>  <48d14e27-7130-48a7-9c9c-ddf9058363dc@f63g2000hsf.googlegroups.com>
Reply-To: "Bruno Luong" <b.luong@fogale.fr>
NNTP-Posting-Host: webapp-02-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1209633006 28890 172.30.248.37 (1 May 2008 09:10:06 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Thu, 1 May 2008 09:10:06 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 390839
Xref: news.mathworks.com comp.soft-sys.matlab:466084


Greg Heath <heath@alumni.brown.edu> wrote in message
<48d14e27-7130-48a7-9c9c-ddf9058363dc@f63g2000hsf.googlegroups.com>...
> 
> There seems to be something in gcd that is missing in
> Bruno and Rogers code ... or does one have to choose
> a more appropriate value for tol? I tried a few, but half of
> the time there were errors of O(1e-16):
> 

Hi Greg,

No I d'ont believe the accuracy issue lies on tol, but
rather on roundoff of successive division (in mod and rem).

One quick fix it adjusting the common factor "g" at the end
of the loop as following:

k=round(a/g);
g=mean(a./k);

%%%%%%%%
%% Here is the full test for my algo

clear
tol =1e-10;
format long e

z0 =randperm(128)';
h0 =sort(z0(1:32));

g00 = 0.4567890123456789;
g0 = tol*round(g00/tol);
dg = g0-g00;
a = g0*h0;

% The engine!
aworking=a;
while 1
    [g imin]=min(aworking);
    irest=[1:imin-1 imin+1:length(aworking)];
    arest=mod(aworking(irest),g);
    arest(abs(arest)<tol)=nan; % lock
    if any(abs(arest))
        aworking(irest)=arest;
    else
        break
    end
end
k=round(a/g);
g=mean(a./k);
disp(g)

% Bruno