10
\$\begingroup\$

Here is a description of Project Euler problem 530:

Every divisor \$d\$ of a number \$n\$ has a complementary divisor \$n/d\$.

Let \$f(n)\$ be the sum of the greatest common divisor of \$d\$ and \$n/d\$ over all positive divisors \$d\$ of \$n\$, that is \$\displaystyle f(n)=\sum_{d/n}\gcd\left(d, \frac{n}{d}\right)\$.

Let \$F\$ be the summatory function of \$f\$, that is \$\displaystyle F(k) = \sum_{n=1}^{k}f(n)\$.

You are given that \$F(10)=32\$ and \$F(1000)=12776\$.

Find \$F(10^{15})\$.

My solution for Euler Problem 530 is as follows:

import java.math.BigInteger;

public class problem530 {

public static void main(String[] args) {

    double start = System.nanoTime();

    BigInteger result = BigInteger.ZERO;
    BigInteger one = BigInteger.ONE;

    for (long i = 1; i <= 1000; i++){
            if ( i == 1) {
                result = result.add(one);
            }
            if (PrimeTest.testLong(i) == true) {
                result = result.add(one).add(one);
            }
            else if (i != 1 && PrimeTest.testLong(i) == false) {
                result = result.add(one);
                for(long j = 1; j*2 <= i; j++) {
                    if (i%j == 0) {
                    long sub = i / j;
                    result = result.add(BigInteger.valueOf(GCD(sub, j)));
                }
            }
        }
    }

    System.out.println(result);

    double duration = (System.nanoTime() - start)/1000000000;

    System.out.println("Your code took " + duration + " seconds to execute.");

}
}

I have PrimeTest and GCD written in libraries. The code works (it returns the correct answer for the test cases). It takes unimaginably long to get the answer for \$F(10^{15})\$ though (essentially, I won't get the answer). I'm looking for ways to optimise it.

At this link, @alan wrote something there about a sieve (presumably in C++?) which I don't quite understand. Does anyone have any ideas for optimisation?

  • I perform GCD operations using the Eulerian algorithm.
  • Prime Testing is performed for:

    n = 1; n*n <= (number being tested); n++
    

To give an idea of how long this code takes to execute:

  • It ran F(10)=32 in 0.008 seconds.
  • It ran F(1000)=12766 in 0.019 seconds.
  • It ran F(100,000)=2907546 in 29.724 seconds.

This was the optimisation suggested in another post on codereview.

It seems to me that the alternative summation suggested there is wrong! F(10) would give you 22, not 32. I'm not sure how the alternative summation was derived.

\$\endgroup\$
  • \$\begingroup\$ If you do double duration = (System.nanoTime()- start)/1000000000;, and removing the end statement it will be more accurate. \$\endgroup\$ – piepi Dec 24 '15 at 2:42
  • \$\begingroup\$ @piepi You're absolutely right. Thanks for that! Any ideas on the rest of it? \$\endgroup\$ – Thevesh Theva Dec 24 '15 at 2:42
  • \$\begingroup\$ @TheveshTheva Figuring out! \$\endgroup\$ – piepi Dec 24 '15 at 2:45
  • \$\begingroup\$ @TheveshTheva So, as you have further mentioned about the Sieve, I asked [that question](codereview.stackexchange.com/questions/114627/extending-sieve-of-eratosthenes-beyond-a-billion-follow-up) a few days ago. I still haven't found a way to extend it beyond 10^9. I will be waiting to know how it can be used to test primes. \$\endgroup\$ – piepi Dec 24 '15 at 2:55
  • 5
    \$\begingroup\$ Welcome to Code Review! I have rolled back the last edit. Please see what you may and may not do after receiving answers. \$\endgroup\$ – Quill Dec 24 '15 at 3:43
4
\$\begingroup\$

Special Cases and Branching

If you look at your prime checking, your code is:

if (prime(i) == true) {
    // ...
}
else if (prime(i) == false) {
    // ...
}

This ends up with you checking composite numbers for primality twice! Once in each branch! The second part could have just been else. Additionally, explicitly checking for == true is an anti-pattern. Prefer:

if (prime(i)) {
    // ...
} else {
    // ...
}

That said, why is being prime a special case? Your non-prime logic would've worked just fine on primes too! And for 1! So just drop all the branches completely and simply iterate.

Need a better algorithm

10^15 is simply too big to brute force the answer. Even dropping all the unnecessary extra work you're doing is likely insufficient. You'll have to try to come up with some optimal substructure this problem has so that you can use dynamic programming.

\$\endgroup\$
  • \$\begingroup\$ Thanks for that! I removed the code which checks for primality twice. The reason why being prime is a special case is that prime numbers only have themselves and 1 as divisors, and the GCD to add in both those cases is 1. I tried running the code with all branches dropped and it ran slower. :/ With primality testing was 4 seconds faster for case F(100,000). \$\endgroup\$ – Thevesh Theva Dec 24 '15 at 3:31
1
\$\begingroup\$

I don't usually use java, but as to the algorithm, to make the suggestion of the real optimization work, you need to consider 2 things.

First, if you want to calculate F(10), this algorithm will go from d=1 to 3 in the first loop. Then from 1 to n*d <= 10. However, with this solution you are doing the gcd of (2,5) but never doing the one for (5,2).

To solve this, you could just add the gcd result*2. However, in special cases like (1,1),(2,2),(3,3) you only want to add it once because they are square numbers. So you need to add gcd*1 if d == n, and gcd*2 if d != n

This gets us to the second consideration. Problem with this is that we will do gcd*2 for (1,2) but then we will do gcd*2 with (2,1) so we will add some cases 4 times instead of 2 as we wanted.

The fix for this is pretty easy, instead of going from n=1 to nd <=k, go from n=d to nd <= k. Why? Because if n < d, we have already passed over those pairs before ((2,1) pair was already seen when d=1 and n=2)

To make all this more clear, here is a short python implementation: (I used d as x and n as y)

def solve(n):
    solution = 0
    for x in xrange(1,int(math.sqrt(n))+1): # from 1 to sqrt(n)
        y = x # to avoid repetitions, from x to x*y <= n
        while x*y<=n :
            solution+= gcd(x,y)*(2 if x != y else 1) # this adds gcd*1 if x ==y and gcd*2 if x != y
            y+=1 
    print solution   

This runs pretty fast as we are not even doing prime checks (which would be faster with a sieve, but not necessary with this solution). Python is pretty slow so it would take a lot to solve 10^15, but you could try to code it in java and check the speed improvement

In python for n = 500000, it takes a bit less than 2 seconds. I did try the same in java, but it seems to take 2.8 seconds for n = 500000, even slower than python. Might be because I had never used BigInteger before, so there might be some improvements to make, here is the code:

    double start = System.nanoTime();

    BigInteger n = new BigInteger("500000");

    BigInteger x = BigInteger.ONE;

    BigInteger solution = BigInteger.ZERO;

    while(x.multiply(x).compareTo(n) != 1){


        BigInteger y = x;


        while(y.multiply(x).compareTo(n) != 1){

            BigInteger gc = gcd(x,y);

            if(x.compareTo(y) == 0){
                solution = solution.add(gc);
            }else{
                solution = solution.add(gc).add(gc);

            }

            y = y.add(BigInteger.ONE);

        }

        x = x.add(BigInteger.ONE);

    }

    double duration = (System.nanoTime() - start)/1000000000;

    System.out.println(solution);
    System.out.println("Your code took " + duration + " seconds to execute.");

For further optimizations, you could memorize some gcd calls as we are clearly repeating a lot of them

\$\endgroup\$
  • \$\begingroup\$ You're absolutely right! I adapted your ideas as follows: (see edit in my post). It definitely compiles faster than the previous version! However, still can't go up to F(10^15) unfortunately. \$\endgroup\$ – Thevesh Theva Dec 25 '15 at 3:06
  • \$\begingroup\$ @TheveshTheva yeah, even with a O(n) solution it would take over 200 hours to solve, so there must be a way to solve it in O(sqrt(n)), not sure how \$\endgroup\$ – juvian Dec 25 '15 at 19:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.