'Quite Good' (but not Perfect) numbers, with a timer question ...

[castroboy]

Homework for Uni, so I'm not asking anyone to solve it for me but any tips would be really appreciated.

I have to find all 'quite good' (not perfect) numbers between 1-1000000. A 'quite good' number is one which is close to perfect, with a 'badness' of anywhere between 1-10 (i.e. 6 is perfect because its factors 1+2+3 = 6; 15 is 'quite good', because its factors 1+3+5 = 9 so it has a 'badness of 6; 32 is 'quite good' because its factors 1+2+4+8+16 = 31 so it has a 'badness' of 1).

The problem I'm having is that the test environment has a time limit of ~10 seconds before it times out.

With a 'badness' of 0 I have no problems, because only the actual perfect numbers are required and I get them in about 100ms.

With a 'badness' of 1, I also have no problems, because only the actual perfect numbers _and_ powers of 2 are required, that takes about 200ms.

Anything beyond that, I'm stumped. Up to about 100,000 I seem to be fine but then the code just takes forever (easily more than a minute at 300,000). I can't think of any algorithm to speed it up, because the 'badness' can vary hugely. I can't make my for statement (i = 0; i <= arga; i+= 2), because some numbers in the desired solution are multiples of 3 and 5.

HELP!

void func2(int arga, int argb) {

    int a = 0;

    for (int i = 2; i < arga; i++) {
        for (int j = 1; j < i; j++) {

            if (i % j == 0) {
                a += j;
            }
        }
        if ((i > 1) && (a >= (i - argb)) && (a <= (i + argb))) {

            cout << i << " ";
        }
        a = 0;
    }
    cout << endl;
}

[jonnin]

printing is very slow. Redirect to a file might save you.

still thinking on the actual algorithm.

you should be able to re-arrange the 2 loops so that the if statement inside is more often true (I%j == 0)

if I > 1 is not possible to fail. remove it. (I starts at 2)
make arga and argb constants.

I%1 == 0, always. Is that exploitable?

if I is prime, the inner loop only hits on 1.
one stupid way to do this is to have a list of the first N primes and skip them (do the %1 increment though). Even if its just the first 10000 or something, the time save would be notable. If your check for prime is slow, you don't gain anything.

can the inner loop be up to sqrt I? I think this is true, and I think it is the key to getting this done.

[castroboy]

The problem is, not all the 'i's are square numbers.

For instance, with an arga of 100,000 (the range I need to search), and an argb of 10 (the 'badness'), one of the numbers I need to get is 1155.

The sqrt of 1155 is 33.98529.

[JLBorges]

See: http://planetmath.org/FormulaForSumOfDivisors

Tip: For fast repeated prime factorisation, first generate all prime numbers up to the square root of the upper bound, say, 1'000'000 (use a sieve algorithm)

[castroboy]

Hey guys-

Thanks for all your advice.

[lastchance]

The "square root thing" is the maximum value you need to check to determine whether it is PRIME. However, that's not what you are doing here. Your method is a brute-force search for proper divisors, and the largest proper divisor of any even number n is n/2, not sqrt(n).

You can, however, stop searching once your sum exceeds "maximum badness", and you will save a lot of computer resources not having to work out the pow() function.

[JLBorges]

Corrected (see lastchance's comment): About a second for 1,000,000 on coliru

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>

std::vector<int> prime_numbers_till( int n ) // invariant: n >= 0
{
     ++n ;
     std::vector<int> sieve( n, true ) ;
     for( int i = 2 ; i < n ; ++i ) if( sieve[i] )
              for( auto j = i*i ; j < n ; j += i ) sieve[j] = false ;

     std::vector<int> primes ;
     for( int i = 2 ; i < n ; ++i ) if( sieve[i] ) primes.push_back(i) ;
     return primes ;
}

int multiplicity( int number, int factor ) // invariant: factor > 1
{
    int n = 0 ;
    for( ; number%factor == 0 ; ++n ) number /= factor ;
    return n ;
}

int ipow( int n, int e ) // invariant: non-negative e
{
    if( e == 0 ) return 1 ;
    else return n * ipow( n, e-1 ) ;
}

// see: http://planetmath.org/FormulaForSumOfDivisors
int term( int prime_factor, int multiplicity ) // invariant: multiplicity >= 0
{ return ( ipow( prime_factor, multiplicity+1 ) - 1 ) / (prime_factor-1) ; }

int sum_proper_divisors( int number, const std::vector<int>& primes )
{
    int sum = 1 ;
    for( int prime_no : primes )
    {
        if( prime_no > (number/2) ) break ;
        const auto mult = multiplicity( number, prime_no ) ;
        if(mult) sum *= term( prime_no, mult ) ;
    }
    return sum - number ;
}

// invariant: non-negative max_badness
bool quite_good( int number, int max_badness, const std::vector<int>& primes ) 
{ return std::abs( number - sum_proper_divisors( number, primes ) ) <= max_badness ; }

int main()
{
    int N = 1'000'000 ;
    int badness = 10 ;

    const auto primes = prime_numbers_till( std::ceil( std::sqrt(N) ) + 2 ) ;
    
    int cnt = 0 ; 
    for( int number = 1 ; number <= N ; ++number )
    {
        if( quite_good( number, badness, primes ) ) 
        { 
            std::cout << std::setw(8) << number << ' ' ; 
            if( ++cnt%10 == 0 ) std::cout << '\n' ;
        }
    }
    std::cout << '\n' ;
}

http://coliru.stacked-crooked.com/a/ff08ef3451b04c4d

[jonnin]

divisors come in pairs though.

so take 100.
100/2 is 2 and 50 paired, so yes, n/2 is correct to pick up the divisors. But you already saw 50 when you checked 2, if set up your algorithm to look in pairs.

if you check up to its square root using that concept, I think you get additional lift.
the factors are 1, 2, 4,5,10 and their "others" 100,50,25,20,10

Does that cut the work down farther or am I misunderstanding the actual problem at hand? On even a very low value like 100, that saves 40 iterations!

[lastchance]

divisors come in pairs though

- @jonnin

Ah yes - I missed that: you have a good point @jonnin. As long as you take the divisors in pairs it would take less iterations. Perhaps to avoid the overhead of computing sqrt(n) you could look out for when the divisors "crossed" each other.

[jonnin]

sqrt is built into x86 FPU. It does not cost all that much. Its a good idea but I would bet the FPU can outperform any logic you spend avoiding it. The only time ive avoided a root and boosted performance successfully was by comparing the squared numbers instead of the roots of the numbers when checking distances. I don't see a way around it here. Even a lookup table of roots (and sin/cos) had worse performance than the FPU (yes, I tried it, lol).

[castroboy]

- @jonnin

You're correct as far as I can see. The algorithm I ended up using finds the lowest divisor, then works out the inverse (e.g. 36 can be divided by 2, which also gives the inverse that 18 is a factor of 36).

No need for me to calculate all the iterations up to i <= (n/2), since I already know the value of n/2 from my first calculation.

Oh, and thanks again for your advice, I was stuck in laborious for() loops that were killing my time limit!

[48T7M4Gy]

#include <iostream>
#include <vector>
#include <iomanip>

int main()
{
    const int limit = 1'000'000;
    std::vector<int> array;
    
    array.resize(limit, 1);
    
    // 'SIEVE'
    int start = 0;
    for(int i = 2; i < limit; i++){
        start = i + i;
        for(int k = start; k < limit; k += i)
            array[k] += i;
    }
    
    // DISPLAY RESULTS
    int column_no = 0;
    int no_columns = 7;
    
    int badness = 0;
    for(int i = 1; i < limit + 1; i++){
        badness = abs(array[i] - i);
        if(badness > 0 and badness < 11){
            std::cout << std::setw(9) << std::left << i;
            column_no++;
            if(column_no % no_columns == 0)
                std::cout << '\n';
        }
    }
    std::cout << '\n';
    
    return 0;
}

2        3        4        5        7        8        9        
10       11       12       14       15       16       18       
20       21       22       26       32       40       44       
50       52       56       64       68       70       88       
104      110      128      130      136      152      184      
196      256      315      368      464      512      592      
650      656      836      884      1012     1024     1155     
1696     1888     1952     2048     2144     2272     2336     
4030     4096     5830     8192     8384     8768     8925     
11096    16384    17816    18632    18904    32128    32445    
32768    32896    33664    45356    65536    70564    77744    
85936    91388    100804   116624   128768   130304   131072   
133376   254012   262144   388076   391612   442365   518656   
521728   522752   524288   527872   528896   
Program ended with exit code: 0

[lastchance]

@kemort - what a brilliant solution! Saves doing individual checks and any % or / operations.

Definitely one to remember.

[castroboy]

@kemort

This doesn't seem like a more elegant solution to me compared to what I'm using currently, sorry.

^^ **EDIT:** Apologies for the above, I've looked further at your code, and it does everything I need in a fraction of the time! Very nice.

[48T7M4Gy]

@lastchance thanks, but the real honour goes to Eratosthenes, all I did was adapt his sieve for primes - that bloke was smart!

( I had to delete my earlier post but the problem I had was eventually so simple to solve - I had one too many zeroes in the 1'0000'000 - hard for me to pick at the time. )

[JLBorges]

kemort!

This is a great solution: simple, elegant and fast.

Scales well too (linear space, sub-quadratic time): about 1.3 seconds for 10,000,000 on coliru
http://coliru.stacked-crooked.com/a/7af40f2f6e4be6f8

[48T7M4Gy]

:)