Path: news.mathworks.com!newsfeed-00.mathworks.com!newscon02.news.prodigy.net!prodigy.net!news.glorb.com!postnews.google.com!f63g2000hsf.googlegroups.com!not-for-mail
From: Greg Heath <heath@alumni.brown.edu>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Largest Common Factor
Date: Wed, 30 Apr 2008 18:31:17 -0700 (PDT)
Organization: http://groups.google.com
Lines: 69
Message-ID: <48d14e27-7130-48a7-9c9c-ddf9058363dc@f63g2000hsf.googlegroups.com>
References: <fv8j8a$2j8$1@fred.mathworks.com> <19640430.1209534092123.JavaMail.jakarta@nitrogen.mathforum.org> 
NNTP-Posting-Host: 69.141.173.117
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: quoted-printable
X-Trace: posting.google.com 1209605478 15944 127.0.0.1 (1 May 2008 01:31:18 GMT)
X-Complaints-To: groups-abuse@google.com
NNTP-Posting-Date: Thu, 1 May 2008 01:31:18 +0000 (UTC)
Complaints-To: groups-abuse@google.com
Injection-Info: f63g2000hsf.googlegroups.com; posting-host=69.141.173.117; 
User-Agent: G2/1.0
X-HTTP-UserAgent: Mozilla/4.0 (compatible; MSIE 6.0; Windows NT 5.1; 
Xref: news.mathworks.com comp.soft-sys.matlab:466048


On Apr 30, 2:04=A0pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Apr 30, 12:01 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> > On Apr 30, 1:41 am, Yumnam Kirani Singh <kirani.si...@gmail.com>
> > wrote:
>
> >> Largest common factor is defined only for integers
> >> not for floating point numbers. See hep on gcd.
>
> > For my purposes I can truncate FP numbers to any
> > practical precision then multiply by the appropriate
> > power of 10 to get integers.

a is positive and monotonic increasing

function [g, h] =3D GCFgh(a,tol);
b =3D round(a/tol);
c =3D b(1)*ones(size(a));
d =3D gcd(b,c);
g =3D tol*min(d);
h =3D round(a/g);
return

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):

close all, clear, clc, disp(datestr(now))

tol =3D 1e-10;
format long e

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

% g00 =3D rand
g00 =3D 0.4567890123456789
g0  =3D tol*round(g00/tol)
dg  =3D g0-g00
a   =3D g0*h0;

% figure,hold on
% plot(g0*(1:128))
% plot(h0,a,'r.')

[g h]    =3D GCFgh(a,tol)
err(1,1) =3D maxabs(g-g0);
err(2,1) =3D maxabs(h-h0);
err(3,1) =3D maxabs(a-g*h);

[g h]    =3D GCFbl(a,tol)
err(1,2) =3D maxabs(g-g0);
err(2,2) =3D maxabs(h-h0);
err(3,2) =3D maxabs(a-g*h);

[g h]    =3D GCFrs(a,tol)
err(1,3) =3D maxabs(g-g0);
err(2,3) =3D maxabs(h-h0);
err(3,3) =3D maxabs(a-g*h);

disp(err)

%   0    1.831867990631508e-015    7.216449660063518e-016
%   0                         0                         0
%   0    2.344791028008331e-013    9.237055564881302e-014

Hope this helps.

Greg