Project Euler

[Cogman]

http://www.projecteuler.net

I'm sure we all know about this.. however, I thought it would be fun to discuss what solutions we come up with here.

I'm currently working on #12

#include <iostream>
#include <vector>
#include <sstream>
#include <cmath>
#include <algorithm>

using namespace std;

int getNumDivisors(long long val)
{
    int numDivisors = 0;
    for (int i = 2; i <= val; ++i)
    {
        if (val % i == 0)
        {
            ++numDivisors;
        }
    }
    return numDivisors;
}

long long triNum(long long n)
{
    return ((n * n + n) >> 1);
}

int main()
{
    int i = 1;
    int numDiv;
    numDiv = getNumDivisors(triNum(i));
    while (numDiv <= 500)
    {
        i += 1;
        numDiv = getNumDivisors(triNum(i));
    }
    cout << i << endl;
}

I have a general brute force solution, but not the actual fast solution. The big slowdown is in the getNumDivisors function. I'm not sure how to speed it up (I have an idea with using primes, but am not sure if it is the correct one.)

 

[Train]

one simple trick, a divisor can't be greater than half of the value. So you can avoid checking every number up to the value, and just up to its halfway point.

i.e. 100 can't be divided by anything larger then 50.

 

[Cogman]

one simple trick, a divisor can't be greater than half of the value. So you can avoid checking every number up to the value, and just up to its halfway point.

i.e. 100 can't be divided by anything larger then 50.

Ah, good point.

 

[eLiu]

Heh some of these problems are really trivial... like #15 is a classic result of probability. And others are like... do something simple with primes or compute the GCM or something.

WARNING: the suggestion given here is wrong (as explained in edit3). Maybe it'll give you some ideas? *shrug*
Couldn't you build up a table of computed results. For example, say I wanted to find the number of divisors of 1002. Then I start dividing 1002... so 1002%2 == 0, so I know numdivisors(1002) = 2*numdivisors(501). So look in my table: do I know the # of divisors of 501? If yes, return that value.
If no, then start dividing: 501%2 != 0, 501%3 == 0, so numdivisors(501) = 2*numdivisors(167), etc. Eventually you reach a base case where the input to numdivisors is prime, in which case you return 2 (divisors: 1 and itself).

Note it's 2* b/c if I have a number alpha*N, then in addition to everything that divides N (1,...,N) , I have the new divisors: (alpha*1,...,alpha*N).

So you could preload this table with numdivisors of 0 through N, and then add additional values as you run across them. Choose N so that for very large numbers, you end up dividing by only a few small numbers before you end up back in the domain of your table.

As for Train's comment, 100 divides 100. According to that euler website, you do need to count that. I mean algorithmically this changes nothing, but wouldn't want to get the wrong answer b/c you're off by 1, lol.

edit: also, wasn't implying that all the problems are easy. Most of them don't appear obvious at all, lol.

edit2: probably you want to at least preload your table with a shitload of primes, come to think about it. My guess is that realizing a base case that some # is prime will be pretty costly, b/c that's the only situation where you'll be doing a large number of modulo ops.


edit3: BALLS. It occured to me that my method is full of FAIL! 🙁 B/c like... 28 has factors:
1,2,4,7,14,28
14 has factors:
1,2,7,14
Can't just double b/c 2*1 is already counted, and 2*7 is already counted! AGH!

The most obvious remedy would be to actually track all divisors (not just how many). But that starts to get a lot more cumbersome... it may be cheaper to track only the prime factorizations & compute the # of divisors from that: if N = x^a * y^b * z^c * ... then the number of divisors is (a+1)*(b+1)*(c+1)*...

So my new suggestion would be to write a prime factorization code that computes factorizations by storing known results. Like 1002 has one more power of 2 than 501. (Like use a seive or something?) Then you can easily compute the # of factors based on the prime factorization. Or I bet there are smart prime factorize-ers out there.

 

[Cogman]

Ok, so I went through a couple of ideas. The first was to store off values of factorization and then look up the number of divisors for a prime factorization... That was wrong.

My next attempt came after doing to google research. I found that getting the prime factorization of the number, and then multiplying the exponents + 1 of those prime factorizations yielded the correct result for the number of divsors.

The following code resulted (and gave the correct answer). BTW, the sqrt code is newtons method for finding sqrts... Yeah, the sqrt function is probably faster but whatever, it works well enough.

**Project euler doesn't really like the giving of answers.. However, if you aren't going to attempt the problem, or just want to see some clever code, here it is.

Spoiler

#include <iostream>
#include <vector>
#include <sstream>
#include <cmath>
#include <algorithm>
#include <map>

using namespace std;

long long primers[] = {2, 3, 5, 7};
vector<long long> primeList(primers, primers + 4);

map<long long, short> numDivisors;

int sqrtVal(long long val)
{
    long long guess = (val >> 2) + 1;
    long long lastGuess = val;
    while ((lastGuess < guess - 1) & (lastGuess > guess + 1))
    {
        guess = guess - (guess >> 1) + (val >> 1) / guess;
    }
    return guess + 1;
}

void getMorePrimes()
{
    primeList.reserve(primeList.size() + 100);
    for (int i = 0; i < 100; ++i)
    {
        bool prime = false;
        long long guess = *(primeList.end() - 1);
        while (!prime)
        {
            prime = true;
            guess += 2;
            long long max = sqrtVal(guess);
            for (long long i = 0; (i < primeList.size()) & (primeList[i] < max); ++i)
            {
                if (guess &#37; primeList[i] == 0)
                {
                    prime = false;
                    break;
                }
            }
        }
        primeList.push_back(guess);
    }
}

short getNumDivisors(long long val)
{
    int numDivisors = 1;
    long long sqrtV = sqrtVal(val);
    int numPrimes = primeList.size();
    for (int i = 0; val != 1; ++i)
    {
        if (i >= numPrimes)
        {
            getMorePrimes();
            numPrimes += 100;
        }
        int multiple = 0;
        while (val % primeList[i] == 0)
        {
            ++multiple;
            val /= primeList[i];
        }
        if (multiple != 0)
        {
            numDivisors *= (multiple + 1);
        }
    }
    return numDivisors;
}

long long triNum(long long n)
{
    return ((n * n + n) >> 1);
}

int main()
{
    getMorePrimes();
    int largest = 0;
    int i = 2;
    int numDiv;
    numDiv = getNumDivisors(triNum(i));
    while (numDiv <= 500)
    {
        i += 1;
        numDiv = getNumDivisors(triNum(i));
        if (numDiv > largest)
        {
            cout << numDiv << endl;
            largest = numDiv;
        }
    }
   cout << "index: " << i << endl << "triVal: " << triNum(i) << endl;
}

 

[eLiu]

Cool, if memoizing prime factorizations & using the (a+1)*(b+1)*(c+1)*... technique didn't work, I was out of ideas, lol. Seems like I posted that (my "edit3") after you'd already figured it out.

 

[iCyborg]

Btw, it's enough to check for divisors from 1 to sqrt(n), and add them in pairs, e.g. for 28, since 2 is a divisor, then add two divisors when you run into 2: 2 and 28/2=14. Just be careful about sqrt(n) since its pair is the same number, so don't count it twice. You could also use Erathostenes sieve to quickly build up primes list.
Also, recomputing (n*n+n)/2 every time seems a waste given the recursive formula. Just keep it through the loop and add n to get n'th sum.

A better idea for the main loop, at the expense of O(n) space is the following:
1. Compute D[n] - num of divisors for n
2. If n is even then num of divisors for n'th triangle number is D[n-1]*D[n/2], if odd D[(n-1)/2]*D[n]

The reason is that since triangle numbers are of the form n(n+1)/2 and n,n+1 are always coprime, the number of divisors for n(n+1) is the product of the number of divisors for D[n] and D[n+1]. One of them is halved by the 2 in denominator, and unfortunately this is why we need to keep track of num of divisors for everybody between n/2 and n.

The code would be:

int getNumDivisors(int n) 
{
	if (n == 1) return 1;
	int num = 2, i;
	for (i=2; i<=sqrt(double(n)); ++i)
		if (n % i == 0) num+=2;
	if ((i-1)*(i-1) == n) --num;

	return num;
}


int main()
{

	int D[25000];
	D[1]=1;
	int n=2;
	int divs = 0;

	while (divs<=500) {
		D[n] = getNumDivisors(n);
		(n%2 == 0) ? divs = D[n-1]*D[n/2] : divs = D[(n-1)/2]*D[n];
		++n;
	}
	n-=2; // n-1, n used for n-1'th triangle number, and decrement once more for the last loop inc
	cout << "index: " << n << endl << "triVal: " << (n*n+n)/2 << " has " << divs << " divisors" << endl
	return 0;
}

 

[Cogman]

Btw, it's enough to check for divisors from 1 to sqrt(n), and add them in pairs, e.g. for 28, since 2 is a divisor, then add two divisors when you run into 2: 2 and 28/2=14. Just be careful about sqrt(n) since its pair is the same number, so don't count it twice. You could also use Erathostenes sieve to quickly build up primes list.

Yeah, Erathostenes sieve would have been faster than what I implemented, I just never liked the memory wasted from it 🙂. As for the sqrt, the prime searcher did use it. However, you'll notice that the divisor counter decreased the size of the number every time it found a matching prime. Recomputing the sqrt (or even computing it in the first place) is wasteful as all numbers will eventually decompose to 1. The only case where this isn't true is prime numbers which are fairly rare when it comes to triangle numbers. The majority of cases benefit from a fast != 1 comparison.

Also, recomputing (n*n+n)/2 every time seems a waste given the recursive formula. Just keep it through the loop and add n to get n'th sum.

This is true 😀 I should have thought that adding is faster than a multiplication and bit shift.

A better idea for the main loop, at the expense of O(n) space is the following:
1. Compute D[n] - num of divisors for n
2. If n is even then num of divisors for n'th triangle number is D[n-1]*D[n/2], if odd D[(n-1)/2]*D[n]

The reason is that since triangle numbers are of the form n(n+1)/2 and n,n+1 are always coprime, the number of divisors for n(n+1) is the product of the number of divisors for D[n] and D[n+1]. One of them is halved by the 2 in denominator, and unfortunately this is why we need to keep track of num of divisors for everybody between n/2 and n.

Interestingly enough, I tried is approach right after the first attempt and surprisingly it is VERY slow compared to the division by primes trick.

The code would be:

int getNumDivisors(int n) 
{
	if (n == 1) return 1;
	int num = 2, i;
	for (i=2; i<=sqrt(double(n)); ++i)
		if (n % i == 0) num+=2;
	if ((i-1)*(i-1) == n) --num;

	return num;
}


int main()
{

	int D[25000];
	D[1]=1;
	int n=2;
	int divs = 0;

	while (divs<=500) {
		D[n] = getNumDivisors(n);
		(n%2 == 0) ? divs = D[n-1]*D[n/2] : divs = D[(n-1)/2]*D[n];
		++n;
	}
	n-=2; // n-1, n used for n-1'th triangle number, and decrement once more for the last loop inc
	cout << "index: " << n << endl << "triVal: " << (n*n+n)/2 << " has " << divs << " divisors" << endl
	return 0;
}

 

[iCyborg]

Yeah, Erathostenes sieve would have been faster than what I implemented, I just never liked the memory wasted from it 🙂. As for the sqrt, the prime searcher did use it.

I was replying to an earlier post, yep, no point with sqrt when you're computing prime factor representation.

Interestingly enough, I tried is approach right after the first attempt and surprisingly it is VERY slow compared to the division by primes trick.

my getNumDivisors is simple and stupid/slow, that's probably the difference. computing divisors for n instead of (n*n+n)/2 plus additional mult should be somewhat faster though. Probably not by much though...