Factor Counting Algorithm

[RAWBERRY]

Recently I've gotten into the project Euler Problems, and have been doing pretty well with coming up with minute solutions. However, I'm having difficulties with problem 012. Brute forcing it is taking much too long. Can anyone help me come up with a good algorithm for factoring numbers quickly?

Here is the problem's page: http://projecteuler.net/index.php?section=problems&id=12

And here is my current code:

#include <iostream>
using namespace std;

int main()
{
    unsigned int n = 0;
    long long int trin = 0;
    int divisors;

    do
    {
        divisors = 0; //Reset the number of divisors
        //Calculate the triangle number
        n++;
        trin = (n*(n+1))/2;

        //count the divisors
        for (long long int count = 1; count <= trin; count++)
        {
            if (trin % count == 0)
            {
                divisors++;
            }
        }
        cout << n << ", " << trin << ", " << divisors << endl;
    } while (divisors < 500);

    cout << fixed << "\n\n" << trin << endl << divisors;
    return 0;
}

[Athar]

You only need check for divisors up to sqrt(n), because every divisor of n has a counterpart divisor on the "other side" of the square root.
With that, I get the result in just 0.7s.

[Browni3141]

checking the output like that will also kill your speed. If you already know that your program works then get rid of it.

[simeonz]

I did not solve any of these problems, so the approach may turn out inefficient, but if you are not doing it for the points, you could sweat harder. :)

You could check the divisors of n and n+1 separately. Then, the next triangle number can reuse the divisors of n+1 from the previous iteration and then compute only the divisors of n+2. Those are smaller quantities and you should be able to identify their divisors faster. Of course, use the square root trick anyways. (And don't recompute the square root for each consecutive number...) You already know that the divisors are not more than 500, so you can preallocate spaces as appropriate.

Also, I am fictionalizing on something else I have never tested - array of remainders. Well, not exactly remainders. Let's call them pseudo-remainders. Every time you see that the pseudo-remainder is 0, this would indicate that the corresponding divisor is exact. Now, my plan is to reload the pseudo-remainder with the value of the corresponding divisor when it turns out 0 on the previous step. For all divisors, no matter if they were factors or not, you decrement the pseudo-remainders with 1 before the next iteration. This works only when you process consecutive numbers, as in n, n+1,... Again the space will never exceed 500 elements or whatever the goal is. But you will process only the part that corresponds to divisors not exceeding the square root, of course.

If you want to use the decomposition that I suggested, then you can couple the pseudo-remainder arrays of n and n+1. You overwrite one with the reload+decrement result from the other.

Regards

EDIT: The remainders array approach requires more memory resources than I thought. I will probably try to do this after all for the experience sake.

EDIT: Another simplification is to search for the prime factors only instead of the all divisors. The technique then is to look for divisor as before, but when you find one, you divide by it as many times as is possible to do it exactly. The number of times you can divide this way gives you the multiplicity of the prime factor. Add 1 to every multiplicity and multiply them together. This should give you the number of divisors, because as wikipedia says: "Any positive divisor of n is a product of prime divisors of n raised to some power."

[simeonz]

Ok, I ended up using two classes.
- One supports the expansion of a number in positional numeral system with custom radix.
- One maintains the prime integer factorization of a number.

Both classes aim to support efficient increment-with-one operations. The former (the expansion) can start from arbitrary number, but the latter (the factorization) starts from 0 only. The purpose is to have incremental factorization of consecutive numbers. I didn't bother to implement booting from arbitrary point.

The expansion holds array of 'digits'. The factorization holds array of 'expansions' with radix-es for each prime number less than the current quantity.

The increment-with-one of the expansions returns the multiplicity of the radix as divisor. This equals the number of trailing zeros in the expansion and is easy to compute from (is actually equal to) the number of carries during the increment.

The increment operation of the factorization returns the total number of divisors for the new value. Computes them according to the rule in my last post (the second edit:), using the multiplicity returned from the expansions. It can detect prime numbers easily, because they have only trivial divisors. In this case a new expansion is created with this prime number as radix and is added to the list.

I use lists. I could rearrange stuff to use vectors for locality of access, although since the lists are allocated right after program start-up they are in consecutive addresses and considering that the memory consumption is low, probably everything gets cached.

The added memory consumption is roughly in the kilos, but the main processing loop does not work from registers anymore, as would be the case with the direct solution. I think there is some trade-off. On my box the result is computed in about 3/4 of a second in debug mode and under decisecond in release mode.

All this business works because I compute the divisors of the two consecutive factors of the triangle number, and consecutive numbers are coprime (they have no common divisors). This means that the number of divisors for the doubled triangle number is the product from the number of divisors of its two consecutive factors. I still needed a (moderately ugly) hack for that halving thing. Also, here is a nice article: http://en.wikipedia.org/wiki/Pronic_number.

Regards

[RAWBERRY]

The squareroot function in the math.h header returns a double, and I cannot use the modulus operator with a floating point number. Is there some sort of work around?

EDIT (Check Second Edit): Need to use google first before I start asking questions. Sorry. There's a function in math.h header which modulates two doubles. fmod(double x, double y). Also I assumed that Aether meant sqrt(trin), and not sqrt(n).

However, with that modification everything is returning 0 divisors, because the square roots of the triangle numbers are never an even 0.

#include <iostream>
#include <math.h>
using namespace std;

int main()
{
    unsigned int n = 0;
    double trin = 0;
    int divisors;

    do
    {
        divisors = 0; //Reset the number of divisors
        //Calculate the triangle number
        n++;
        trin = sqrt((n*(n+1))/2);


        //count the divisors
        for (double count = 1.0; count <= trin; count++)
        {
            if (fmod(trin,count) == 0)
            {
                divisors++;
            }
        }
        cout << fixed << n << ", " << trin << ", " << divisors << endl;
    } while (divisors < 500);

    cout << fixed << "\n\n" << trin * trin << endl << divisors;
    return 0;
}

BIG EDIT: I completely misunderstood, and I feel so foolish right now because I had learned this before. I needed to check only up to the square root of each triangle number, instead of the entire triangle number, so my for loop looked like:

for (long long int count = 1; count <= sqrt(trin); count++)

This follows the minute rule, however my answer of 842161320 is wrong.

EDIT 3: Forgot to double my divisors. Thanks for all the help, my solution is now correct and was calculated in 1.731 seconds.

[simeonz]

You achieved approximately the same performance as mine
- using much simpler brute-force approach (which is arguably more readable)
- using long long int (I used int)
- using trigonometric call (..err, inside your loop condition)
My thinking is too convoluted to be a programmer? LOL

Best Regards

[hamsterman]

a faster way to go up to sqrt(trin) is count*count <= trin

[RAWBERRY]

Indeed it is. Thanks hamster!

[simeonz]

@RAWBERRY
Don't use floating point counter. Use an integer one. Precompute floor(sqrt(n * (n+1) / 2)) before the loop and store it in variable of integral type.

This is the fastest I figured:

#include <iostream>
#include <ctime>
using namespace std;

int main()
{
    //**** START TIMING ****
    clock_t start = clock();

    unsigned int n_plus1 = 0;

    int divisors_n;
    int divisors_n_plus1 = 0;

    do
    {
        n_plus1++;

        divisors_n = divisors_n_plus1;
        divisors_n_plus1 = 1;

        int trial_number = 2; //trial number for prime factor of n_plus1
        int remaining_divisor = n_plus1; //part from n_plus1 not yet factored

        while (remaining_divisor > 1)
        {

            if ((remaining_divisor % trial_number) == 0)
            {
                int factor_multiplicity = 0;

                do
                {
                    remaining_divisor /= trial_number;
                    ++factor_multiplicity;
                }
                while((remaining_divisor % trial_number) == 0);

                divisors_n_plus1 *=
                        (int)(trial_number > 2) /*UGLY HACK for the halving thing*/
                        + factor_multiplicity;
            }

            ++trial_number;
        }

        if (divisors_n_plus1 == 1) { divisors_n_plus1 = 2; }
        //Handles primes

    } while (divisors_n * divisors_n_plus1 < 500);

    //**** END TIMING ****
    cout << "Time: " << (double)(clock() - start) / CLOCKS_PER_SEC << endl;

    cout << "Triangle number: " << (((long long)(n_plus1 - 1) * n_plus1) / 2) << endl;
    cout << "Divisors #: " << divisors_n * divisors_n_plus1 << endl;
    return 0;
}

I post it, because for me it illustrates the proper way to factor a number, and therefore to count its divisors. I used some of my earlier ideas (to count the divisors of n and n+1 separately), but I have to admit that when I tested it against my original idea with arrays of remainders, this current solution performed much better. The incremental approach that I thought was good turned out very bloated :(

This uses 0.07s to compute the number.

Regards