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.