MATLAB: Find greatest common denominator of two values to get integer

greatest common denominator

Hello,
I have two sampling rates:
fs_A = 24414;
fs_B = 92782;
And I want to find the greates value (x) that can be devided into both sampling rates to give integers. Thus, the two equations are:
24414 / x = an integer ; 92782 / x = another integer
Can anyone help me out?
I only get an integer for 'x' when I use the gcd function, but I would like to have non-integers solutions as well.
Thanks.

Best Answer

Ugh. This is not a pretty thing to do, when working in floating point arithmetic. You need to work with tolerances. Yes, GCD is perfect, when the numbers are integers. So what can you do?
Is there a solution when no integers are assumed? Hmm, I wonder. At first, you might think of trying to use the Euclidean algorithm, just not using integer arithmetic. In fact, that might even have a good chance to work, IF you were careful with the tolerances. But first, let me create an example.
format long g
>> A = 137*3*pi/2
A =
          645.597290312703
          
>> B = 137*7*pi
B =
          3012.78735479261
                   
gcd(A,B)
Error using gcd (line 31)
Inputs must be real integers.
A problem is that A and B are not known exactly as floating point numbers in MATLAB. In fact, I've constructed the example such that the largest value we can divide both A and B by, and get an integer result is simple: 137*pi/2. But suppose you just gave me those numbers, without telling me how they were created? Is there a simple trick? Well, yes. But you still need to use tolerances.
The idea is to use the function rat, with a tolerance.
[N,D] = rat(A/B,1e-13)
N =
     3
D =
    14
What does that tell us? What is the largest number we can divide both A and B by, and still get an integer? I'll put it in the variable div.
div = A/N
div =
          215.199096770901
I could have computed it equally well as
div = B/D
div =
          215.199096770901
As you can see, it is the same number. Is it the theoretically correct value? I said above that it should be...
137*pi/2
ans =
          215.199096770901
So rat actually returns the result, AFTER dividing A and B by the desired number div, and all I had to do was to extract it using N or D.
As you can see, this is indeed the largest real number you can divide into both A and B, and have them both yield an integer.
Can I try a more difficult case, where i don't know the answer in advance?
A = rand(1)*1e5
A =
          14188.6338627215
>> B = rand(1)*1e5
B =
          42176.1282626275
>> [N,D] = rat(A/B,1e-13)
N =
      218086
D =
      648267
>> div = A/N
div =
        0.0650598106376454
>> B/D
ans =
        0.0650598106376346
Here div is a pretty small number, but then what can you expect from two randomly chosen numbers? I might even be happy the divisor was that large given two entirely random numbers.
Ok, so how does this work, when applied to your problem? Remember, we need to use tolerances here. If you want a VERY accurate approximation, then you need to use a tiny tolerance.
fs_A = 24414;
fs_B = 92782;
>> [N,D] = rat(fs_A/fs_B,1e-13)
N =
       12207
D =
       46391
>> div = A/N
div =
          1.16233586161395
However, if I am willing to back down on the tolerance, then I can get a result that is less accurate, but perhaps closer to what you want. 
>> [N,D] = rat(fs_A/fs_B,1e-4)
N =
     5
D =
    19
    
>> fs_A/N
ans =
                    4882.8
>> fs_B/D
ans =
          4883.26315789474
As you see here, there is aa difference between how we might compute div, based on which fraction we try. But that suggests what may perhaps be a simple compromise. Just take the average.
>> div = (fs_A/N + fs_B/D)/2
div =
          4883.03157894737
>> fs_A/div
ans =
          4.99976287379712
          
>> fs_B/div
ans =
          19.0009010795709
So, to about 4 significant digits, that seems reasonable, and as good as we can get for that tolerance. A smaller tolerance yields a better solution:
>> [N,D] = rat(fs_A/fs_B,1e-6)
N =
   556
D =
        2113
>> div = (fs_A/N + fs_B/D)/2
div =
          43.9100761983882
>> fs_A/div
ans =
          555.999946110232
>> fs_B/div
ans =
           2113.0002048005
I'm afraid I can't do much better than a solution this simple. You will need to choose a tolerance, but then the rest is trivial.
How does this solution work when GCD would also have worked? In this next example, the GCD should be 12, and all are integers, so GCD will work, but so should RAT.
A = 144;
B = 60;
gcd(A,B)
ans =
    12
[N,D] = rat(A/B,1e-6);
div = (A/N + B/D)/2
div =
    12
RAT agrees with GCD.
Related Question