Fermat's factorisation in C++

Question

For fun, I've been implementing some maths stuff in C++, and I've been attempting to implement Fermats Factorisation Method, however, I don't know that I understand what it's supposed to return. This implementation I have, returns 105 for the example number 5959 given in the Wikipedia article.

The pseudocode in Wikipedia looks like this:

One tries various values of a, hoping that is a square.

FermatFactor(N): // N should be odd
    a → ceil(sqrt(N))
    b2 → a*a - N
    while b2 isn't a square:
        a → a + 1    // equivalently: b2 → b2 + 2*a + 1
        b2 → a*a - N //               a → a + 1
    endwhile
    return a - sqrt(b2) // or a + sqrt(b2)

My C++ implementation, look like this:

int FermatFactor(int oddNumber)
{
    double a = ceil(sqrt(static_cast<double>(oddNumber)));
    double b2 = a*a - oddNumber;
    std::cout << "B2: " << b2 << "a: " << a << std::endl;

    double tmp = sqrt(b2);
    tmp = round(tmp,1);
    while (compare_doubles(tmp*tmp, b2))  //does this line look correct?
    {
        a = a + 1;
        b2 = a*a - oddNumber;
        std::cout << "B2: " << b2 << "a: " << a << std::endl;
        tmp = sqrt(b2);
        tmp = round(tmp,1);
    }

    return static_cast<int>(a + sqrt(b2));
}

bool compare_doubles(double a, double b)
{
    int diff = std::fabs(a - b);
    return diff < std::numeric_limits<double>::epsilon();
}

What is it supposed to return? It seems to be just returning a + b, which is not the factors of 5959?

EDIT

double cint(double x){
    double tmp = 0.0;
    if (modf(x,&tmp)>=.5)
        return x>=0?ceil(x):floor(x);
    else
        return x<0?ceil(x):floor(x);
}

double round(double r,unsigned places){
    double off=pow(10,static_cast<double>(places));
    return cint(r*off)/off;
}

Answer 1

Do note that you should be doing all those calculations on integer types, not on floating point types. It would be much, much simpler (and possibly more correct).

---

Your compare_doubles function is wrong. diff should be a double.

And once you fix that, you'll need to fix your test line. compare_doubles will return true if its inputs are "nearly equal". You need to loop while they are "not nearly equal".

So:

bool compare_doubles(double a, double b)
{
    double diff = std::fabs(a - b);
    return diff < std::numeric_limits<double>::epsilon();
}

And:

while (!compare_doubles(tmp*tmp, b2))  // now it is
{

And you will get you the correct result (101) for this input.

You'll also need to call your round function with 0 as "places" as vhallac points out - you shouldn't be rounding to one digit after the decimal point.

The Wikipedia article you link has the equation that allows you to identify b from N and a-b.

Answer 2

There are two problems in your code:

1. compare_doubles return true when they are close enough. So, the while loop condition is inverted.
2. The round function requires number of digits after decimal point. So you should use round(x, 0).

As I've suggested, it is easier to use int for your datatypes. Here's working code implemented using integers.

Answer 3

The two factors are (a+b) and (a-b). It is returning one of those. You can get the other easily.

N = (a+b)*(a-b)
a-b = N/(a+b)