MATLAB Central Newsreader - Largest Common Factor

Largest Common Factor

Greg Heath — Wed, 30 Apr 2008 00:12:05 -0400

What is the best way to find the largest common
factor of an array? For example, 4.1 is the largest common factor of
[12.3 20.5 32.8]].
I've looped over whole fractions of 12.3 to
get an integer array but feel there must be a better way.

Thanks,

Greg

Re: Largest Common Factor

Roger Stafford — Wed, 30 Apr 2008 01:50:02 -0400

Greg Heath <heath@alumni.brown.edu> wrote in message
<7ebde538-62ed-4bba-
aed4-15843490c573@24g2000hsh.googlegroups.com>...
> What is the best way to find the largest common
> factor of an array? For example, 4.1 is the largest common factor of
> [12.3 20.5 32.8]].
> I've looped over whole fractions of 12.3 to
> get an integer array but feel there must be a better way.
>
> Thanks,
>
> Greg
---------
  Euclid's algorithm would still work for non-integer elements like 12.3, but it
would have to have a tolerance that allows for round-off errors. For example,
the binary floating point value matlab uses to approximate 12.3 is not an
exact integral multiple of the number with which it approximates 4.1. I don't
know whether the current version of 'gcd' makes such an allowance or
whether it even permits non-integers as input arguments. (My own version
does not.)

  Using Euclid's algorithm, the numbers 12.3 and 20.5 would go:

rem(20.5,12.3) = 8.2
rem(12.3, 8.2) = 4.1
rem( 8.2, 4.1) = 0.0

except that this last one wouldn't be exactly zero but only within some
tolerated distance from zero. If you started with, say, pi and sqrt(2), this
procedure would not stop until it had dropped far down into numbers in the
tolerated noise level.

  I can think of no better way of handling more than two elements than
repeating the above operation using the largest common factor found up to a
point in an array, together with the array's next element. Short of a matlab
function designed for this, it sounds like using a common garden variety for-
loop to accomplish this feat is your best bet. This algorithm does not appear
to lend itself to any kind of vectorizaton that I can think of.

Roger Stafford

Re: Largest Common Factor

Yumnam Kirani Singh — Wed, 30 Apr 2008 05:41:01 -0400

Largest common factor is defined only for integers not for floating point numbers. See hep on gcd.

Re: Largest Common Factor

Bruno Luong — Wed, 30 Apr 2008 06:58:02 -0400

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
wrote in message <fv8j8a$2j8$1@fred.mathworks.com>...

>
> I can think of no better way of handling more than two
elements than
> repeating the above operation using the largest common
factor found up to a
> point in an array, together with the array's next element.

I also think euclide's algorithm is a great idea. But we
need to define how to apply on the array. What about

1) find the minimum element "p" (absolute value) of the array
2) Replace all other elements by the remaining with "p".
3) if all the replacements are "zero" (within a tolerance)
then stop
4) Otherwise loop on 1)

When break the loop, p is the answer.

Bruno

Re: Largest Common Factor

Bruno Luong — Wed, 30 Apr 2008 10:42:02 -0400

"Bruno Luong" <b.luong@fogale.fr> wrote in message
<fv959q$34q$1@fred.mathworks.com>...
> "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
> wrote in message <fv8j8a$2j8$1@fred.mathworks.com>...
>
> >
> > I can think of no better way of handling more than two
> elements than
> > repeating the above operation using the largest common
> factor found up to a
> > point in an array, together with the array's next element.
>
> I also think euclide's algorithm is a great idea. But we
> need to define how to apply on the array. What about
>
> 1) find the minimum element "p" (absolute value) of the array
> 2) Replace all other elements by the remaining with "p".
> 3) if all the replacements are "zero" (within a tolerance)
> then stop
> 4) Otherwise loop on 1)
>
> When break the loop, p is the answer.
>
> Bruno

Here is a code of the above idea:

asize=100;
maxa=10;
myeps=1e-10; % small number
% generate asize array, multiple of 4.1
a=4.1*(1+round(rand(1,asize)*(maxa-1)));

% The engine!

aworking=a;
while 1
    [p imin]=min(aworking);
    irest=[1:imin-1 imin+1:length(aworking)];
    arest=mod(aworking(irest),p);
    arest(abs(arest)<myeps)=nan; % lock
    if any(abs(arest))
        aworking(irest)=arest;
    else
        break
    end
end
disp(p);

Re: Largest Common Factor

Roger Stafford — Wed, 30 Apr 2008 11:12:03 -0400

"Bruno Luong" <b.luong@fogale.fr> wrote in message <fv959q$34q
$1@fred.mathworks.com>...
> I also think euclide's algorithm is a great idea. But we
> need to define how to apply on the array. What about
>
> 1) find the minimum element "p" (absolute value) of the array
> 2) Replace all other elements by the remaining with "p".
> 3) if all the replacements are "zero" (within a tolerance)
> then stop
> 4) Otherwise loop on 1)
>
> When break the loop, p is the answer.
>
> Bruno
-----------
  I would envision implementing Euclid's algorithm for arrays of non-integers
with something like the following. Let x be a vector of values, not necessarily
integers and not necessarily of like sign, for which we wish to find the largest
common factor, f. That is, we want the largest magnitude, f, such that all
elements of x are integer multiples of f (ignoring round-off errors.) For this
to work properly, the vector x should not contain any zero or nearly zero
values - the multiples should be non-zero multiples.

tol = 1e-8;
f = x(1);
for k = 2:length(x)
  r = f;
  f = x(k);
  while abs(r) > tol
   t = f;
   f = r;
   r = rem(t,f);
  end
end
f = abs(f); % This is the largest common factor of elements in x

  The value for 'tol' here is not necessarily optimum. The best choice would
depend on the magnitudes of elements in x and on the length of x, but I
haven't worked out a good criterion as yet. It does have to be chosen
carefully. If it is too small, it will miss cases where a common factor exists
but gets overlooked because of accumulating round-off errors. If it is too
large, it will allow false, overly-large common factors to be admitted.

Roger Stafford

Re: Largest Common Factor

Greg Heath — Wed, 30 Apr 2008 16:01:39 -0400

On Apr 30, 1:41=A0am, 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.

Thanks.

Greg

Re: Largest Common Factor

Greg Heath — Wed, 30 Apr 2008 16:04:56 -0400

On Apr 30, 7:12=A0am, "Roger Stafford"
<ellieandrogerxy...@mindspring.com.invalid> wrote:
> "Bruno Luong" <b.lu...@fogale.fr> wrote in message <fv959q$34q
>
> $...@fred.mathworks.com>...> I also think euclide's algorithm is a great i=
dea. But we
> > need to define how to apply on the array. What about
>
> > 1) find the minimum element "p" (absolute value) of the array
> > 2) Replace all other elements by the remaining with "p".
> > 3) if all the replacements are "zero" (within a tolerance)
> > then stop
> > 4) Otherwise loop on 1)
>
> > When break the loop, p is the answer.
>
> > Bruno
>
> -----------
> =A0 I would envision implementing Euclid's algorithm for arrays of non-int=
egers
> with something like the following. =A0Let x be a vector of values, not nec=
essarily
> integers and not necessarily of like sign, for which we wish to find the l=
argest
> common factor, f. =A0That is, we want the largest magnitude, f, such that =
all
> elements of x are integer multiples of f (ignoring round-off errors.) =A0F=
or this
> to work properly, the vector x should not contain any zero or nearly zero
> values - the multiples should be non-zero multiples.
>
> =A0tol =3D 1e-8;
> =A0f =3D x(1);
> =A0for k =3D 2:length(x)
> =A0 r =3D f;
> =A0 f =3D x(k);
> =A0 while abs(r) > tol
> =A0 =A0t =3D f;
> =A0 =A0f =3D r;
> =A0 =A0r =3D rem(t,f);
> =A0 end
> =A0end
> =A0f =3D abs(f); % This is the largest common factor of elements in x
>
> =A0 The value for 'tol' here is not necessarily optimum. =A0The best choic=
e would
> depend on the magnitudes of elements in x and on the length of x, but I
> haven't worked out a good criterion as yet. =A0It does have to be chosen
> carefully. =A0If it is too small, it will miss cases where a common factor=
exists
> but gets overlooked because of accumulating round-off errors. =A0If it is =
too
> large, it will allow false, overly-large common factors to be admitted.
>
> Roger Stafford

Roger and Bruno,

Thanks. The application will be to overlay a uniform grid on a
nonuniform grid for the purposes of using the fft instead of a dft.

Greg

Re: Largest Common Factor

Greg Heath — Wed, 30 Apr 2008 18:04:26 -0400

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.

Example:

a = 2*pi*sort(rand(128,1)); % realistic data
a = [12.3 20.5 32.8] % demo data
tol = 4;

% [g h] = GCFgh(a,tol)
% function [g h] = GCFgh(a,tol)

mult = 10^tol;
b = round(mult*a);
c = b(1)*ones(size(a));
d = gcd(b,c);
g = min(d)/mult
h = a/g

a = 12.3 20.5 32.8
g = 4.1
h = 3 5 8

Greg

Re: Largest Common Factor

Greg Heath — Thu, 01 May 2008 01:31:17 -0400

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

Re: Largest Common Factor

Bruno Luong — Thu, 01 May 2008 09:10:06 -0400

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