Lowest Common Multiple with doubles in C

Question

I am doing an assignment for a Coursera class that asks me to calculate the Lowest Common Multiple of two numbers, either of which are no larger than 2 * 10 ^ 9. I'm writing this in C and I'm running my code on a test case with the numbers 226553150 and 1023473145. The answer is 46374212988031350, but I'm getting 46374212988031344, which is off by 6!

I've written a correct solution in Python that uses essentially the same approach as the one I've posted below, but the details of numeric precision are obviously taken care of for me. I post this to SO to learn something about floating point precision in C, and because most of the questions I've seen on the internet and SO regarding LCM deal only with integers.

Here is my code, which I'm compiling with gcc -pipe -O2 -std=c11 lcm.c:

#include <stdio.h>
#include <math.h>

double gcd(double a, double b) {
    if (b == 0) {
        return a;
    }

    return gcd(b, fmod(a,b));
}

double lcm(double a, double b) {
    return (a * b) / gcd(a,b);
}

int main() {

    double a;
    double b;

    scanf("%lf %lf", &a, &b);

    printf("%.0lf\n", lcm(a,b));

    return 0;
}

Answer 1

I doubt why you use float not long. I change float to long as following, then it works fine.

#include <stdio.h>
#include <math.h>

long gcd(long a, long b) {
    if (b == 0) {
        return a;
    }

    return gcd(b, a%b);
}

long lcm(long a, long b) {
    return (a * b) / gcd(a,b);
}

int main() {

    long a;
    long b;

    scanf("%ld %ld", &a, &b);
    printf("%ld\n", lcm(a,b));

    return 0;
}

Answer 2

The closest number to 46374212988031350 that can be represented by a double is 46374212988031352 (off by 2). You can test that by using the most straight forward calculation.

#include <stdio.h>

int main() {

   // The gcd of 226553150 and 1023473145 is 5
   double a = 226553150;
   double b = 1023473145;
   double c = b/5;

   double lcm = a*c;

   printf("lcm: %.0lf\n", lcm);

   return 0;
}

Output:

lcm: 46374212988031352

You are making things worse by computing (a * b)/gcd(a, b). The closest number to 231871064940156750, (a*b), that can be represented by floating point numbers is 231871064940156736. In other words, you are losing more accuracy by computing (a*b) first.

You haven't posted the Python code you used to do the same computation. I am guessing Python is using integral types for the numbers. If I use:

a = 226553150;
b = 1023473145;
c = b/5;

lcm = a*c

print("lcm:", lcm)

I get the output:

('lcm:', 46374212988031350)

However, if I use floating point literals for a and b, I get a different output:

a = 226553150.0;
b = 1023473145.0;
c = b/5;

lcm = a*c

print("lcm:", "%18.0lf" % lcm)

Output:

('lcm:', ' 46374212988031352')

In summary, the difference you are seeing between the C program and the Python program is due to use of floating point types vs integral types. If you use long or long long instead of double, you should get the same output as the Python program, as shown in the answer by Deng Haijun.