How to efficiently compute weird numbers

user_13576827 0
Viewed 247

I am trying to print n weird numbers where n is really big number (eg: 10000).

I found this site to check the algorithm for n 600 if I have some errors: http://www.numbersaplenty.com/set/weird_number/more.php

However, my algorithm is really slow in bigger numbers:

import java.util.ArrayList;
import java.util.List;

public class Test {

    public static void main(String[] args) {
        int n = 2;

        for ( int count = 1 ; count <= 15000 ; n += 2 ) {
            if (n % 6 == 0) {
                continue;
            }

            List<Integer> properDivisors = getProperDivisors(n);
            int divisorSum = properDivisors.stream().mapToInt(i -> i.intValue()).sum();

            if ( isDeficient(divisorSum, n) ) {
                continue;
            }

            if ( isWeird(n, properDivisors, divisorSum) ) {
                System.out.printf("w(%d) = %d%n", count, n);
                count++;
            }
        }
    }

    private static boolean isWeird(int n, List<Integer> divisors, int divisorSum) {
        return isAbundant(divisorSum, n) && ! isSemiPerfect(divisors, n);
    }

    private static boolean isDeficient(int divisorSum, int n) {
        return divisorSum < n;
    }

    private static boolean isAbundant(int divisorSum, int n) {
        return divisorSum > n;
    }

    private static boolean isSemiPerfect(List<Integer> divisors, int sum) {
        int size = divisors.size();

        //  The value of subset[i][j] will be true if there is a subset of divisors[0..j-1] with sum equal to i 
        boolean subset[][] = new boolean[sum+1][size+1];

        // If sum is 0, then answer is true 
        for (int i = 0; i <= size; i++) {
            subset[0][i] = true; 
        }

        //  If sum is not 0 and set is empty, then answer is false 
        for (int i = 1; i <= sum; i++) {
            subset[i][0] = false; 
        }

        // Fill the subset table in bottom up manner 
        for ( int i = 1 ; i <= sum ; i++ ) {
            for ( int j = 1 ; j <= size ; j++ ) {
                subset[i][j] = subset[i][j-1];
                int test = divisors.get(j-1);
                if ( i >= test ) {
                    subset[i][j] = subset[i][j] || subset[i - test][j-1]; 
                }
            } 
        } 

        return subset[sum][size];
    }

    private static final List<Integer> getProperDivisors(int number) {
        List<Integer> divisors = new ArrayList<Integer>();
        long sqrt = (long) Math.sqrt(number);
        for ( int i = 1 ; i <= sqrt ; i++ ) {
            if ( number % i == 0 ) {
                divisors.add(i);
                int div = number / i;
                if ( div != i && div != number ) {
                    divisors.add(div);
                }
            }
        }
        return divisors;
    }

}

I have three easy breakouts:

  1. If a number is divisable by 6 it is semiperfect which means it cannot be weird

  2. If a number is deficient this means it cannot be weird

The above points are based on https://mathworld.wolfram.com/DeficientNumber.html

  1. If a a number is odd it cannot be weird at least for 10^21 numbers (which is good for the numbers I am trying to obtain).

The other optimization that I used is the optimization for finding all the dividers of a number. Instead of looping to n, we loop to SQRT(n).

However, I still need to optimize: 1. isSemiPerfect because it is really slow 2. If I can optimize further getProperDivisors it will be good too.

Any suggestions are welcome, since I cannot find any more optimizations to find 10000 weird numbers in reasonable time.

PS: Any code in Java, C#, PHP and JavaScript are OK for me.

EDIT: I found this topic and modified isSemiPerfect to look like this. However, it looks like it does not optimize but slow down the calculations:

private static boolean isSemiPerfect(List<Integer> divisors, int n) {
        BigInteger combinations = BigInteger.valueOf(2).pow(divisors.size());
        for (BigInteger i = BigInteger.ZERO; i.compareTo(combinations) < 0; i = i.add(BigInteger.ONE)) {
          int sum = 0;
          for (int j = 0; j < i.bitLength(); j++) {
            sum += i.testBit(j) ? divisors.get(j) : 0;
          }

          if (sum == n) {
            return true;
          }
        }

        return false;
      }
mathematical-optimization optimization algorithm
1 Answers

The issue is indeed in function isSemiPerfect. I transposed your code in C++, it was still quite slow.

Then I modified this function by using backtracking. I now obtain the first 15000 weird values in about 15s. My interpretation is that in about all the cases, the value is semiperfect, and the backtracking function converges rapidly.

Note also that in my backtracking implementation, I sort the divisors, which allow to reduce the number of cases to be examined.

Edit 1: an error was corrected in getProperDivisors. Final results did not seem to be modified !

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

//  return true if sum is obtained

bool test_sum (std::vector<int>& arr, int amount) {
    int n = arr.size();

    std::sort(arr.begin(), arr.end(), std::greater<int>());
    std::vector<int> bound (n);
    std::vector<int> select (n);
    bound[n-1] = arr[n-1];
    for (int i = n-2; i >= 0; --i) {
        bound[i] = bound[i+1] + arr[i];
    }

    int sum = 0;    // current sum
    int i = 0;      // index of the coin being examined
    bool up_down = true;

    while (true) {
        if (up_down) {
            if (i == n || sum + bound[i] < amount) {
                up_down = false;
                i--;
                continue;
            }
            sum += arr[i];
            select[i] = 1;
            if (sum == amount) return true;
            if (sum < amount) {
                i++;
                continue;
            }
            up_down = false;
            if (select[i] == 0) i--;
        } else {            // DOWN
            if (i < 0) break;
            if (select[i] == 0) {
                i--;
            } else {
                sum -= arr[i];
                select[i] = 0;
                i++;
                up_down = true;
            }
        }
    }
    return false;
}

bool isDeficient(int divisorSum, int n) {
    return divisorSum < n;
}

bool isAbundant(int divisorSum, int n) {
    return divisorSum > n;
}

bool isSemiPerfect(std::vector<int> &divisors, int sum) {
    int size = divisors.size();

    //  The value of subset[i][j] will be true if there is a subset of divisors[0..j-1] with sum equal to i 
    //bool subset[sum+1][size+1];
    std::vector<std::vector<bool>> subset(sum+1, std::vector<bool> (size+1));

    // If sum is 0, then answer is true 
    for (int i = 0; i <= size; i++) {
        subset[0][i] = true; 
    }

    //  If sum is not 0 and set is empty, then answer is false 
    for (int i = 1; i <= sum; i++) {
        subset[i][0] = false; 
    }

    // Fill the subset table in bottom up manner 
    for ( int i = 1 ; i <= sum ; i++ ) {
        for ( int j = 1 ; j <= size ; j++ ) {
            subset[i][j] = subset[i][j-1];
            int test = divisors[j-1];
            if ( i >= test ) {
                subset[i][j] = subset[i][j] || subset[i - test][j-1]; 
            }
        } 
    } 
    return subset[sum][size];
}

bool isWeird(int n, std::vector<int> &divisors, int divisorSum) {
    //return isAbundant(divisorSum, n) && !isSemiPerfect(divisors, n);
    return isAbundant(divisorSum, n) && !test_sum(divisors, n);
}

std::vector<int> getProperDivisors_old(int number) {
    std::vector<int> divisors;
    long sqrtn = sqrt(number);
    for ( int i = 1 ; i <= sqrtn ; i++ ) {
        if ( number % i == 0 ) {
            divisors.push_back(i);
            int div = number / i;
            if (div != i && div != number) {
                divisors.push_back(div);
            }
        }
    }
    return divisors;
}

    std::vector<int> getProperDivisors(int number) {
    std::vector<int> divisors;
    long sqrtn = sqrt(number);
    divisors.push_back(1);
    for ( int i = 2 ; i <= sqrtn ; i++ ) {
        if (number % i == 0) {
            divisors.push_back(i);
            int div = number/i;
            if (div != i) divisors.push_back(div);
        }
    }
    return divisors;
}


int main() {        

    int n = 2, count;
    std::vector<int> weird;
    int Nweird = 15000;

    for (count = 0; count < Nweird; n += 2) {
        if (n % 6 == 0) continue;

        auto properDivisors = getProperDivisors(n);
        int divisorSum = std::accumulate (properDivisors.begin(), properDivisors.end(), 0);

        if (isDeficient(divisorSum, n) ) {
            continue;
        }

        if (isWeird(n, properDivisors, divisorSum)) {
            //std::cout << count << " " << n << "\n";
            weird.push_back (n);
            count++;
        }
    }
    for (int i = Nweird - 10; i < Nweird; ++i) {
        std::cout << weird.at(i) << " ";
    }
    std::cout << "\n";
}

EDIT 2 The generation of Divisors were completely redefined. It uses now prime decomposition. Much more complex, but global time divided by 7.5. Generation of weird numbers take now 2s on my PC.

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

template <typename T>
struct factor {T val = 0; T mult = 0;};

template <typename T>
class decompo {
    private:
        std::vector<T> memory = {2, 3, 5, 7, 11, 13, 17, 19, 23, 31, 37, 39, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97};
        T index = 0;
    public:
        decompo () {};
        void reset () {index = 0;};
        T pop () {index = memory.size() - 1; return memory[index];};
        T get_next ();
        std::vector<T>  find_all_primes (T n);
        std::vector<factor<T>> decomp (T n);
        std::vector<T> GetDivisors (T n);
        void  complete (T n);
};

template <typename T>
T decompo<T>::get_next () {
    ++index;
    if (index <= memory.size()) {
        return memory[index-1];
    }
    T n = memory.size();
    T candidate = memory[n-1] + 2;
    while (1) {
        bool found = true;
        for (T i = 1; memory[i] * memory[i] <= candidate; ++i) {
            if (candidate % memory[i] == 0) {
                found = false;
                break;
            }
        }
        if (found) {
            memory.push_back (candidate);
            return candidate;
        }
        candidate += 2;
    }
}

template <typename T>
std::vector<T> decompo<T>::find_all_primes (T n) {
    reset();
    std::vector<T> result;
    while (1) {
        T candidate = get_next();
        if (candidate <= n) {
            result.push_back (candidate);
        } else {
            return result;
        }
    }
}

template <typename T>
void decompo<T>::complete (T n) {
    T last = pop();
    while (last < n) {
        last = get_next();
    }
    return;
}

template <typename T>
std::vector<factor<T>> decompo<T>::decomp (T n) {
    reset();
    std::vector<factor<T>> result;
    if (n < 2) return result;
    T candidate = get_next();
    T last_prime = 0;
    while (candidate*candidate <= n) {
        if (n % candidate == 0) {
            if (candidate == last_prime) {
                result[result.size()-1].mult ++;
            } else {
                result.push_back ({candidate, 1});
                last_prime = candidate;
            }
            n /= candidate;
        } else {
            candidate = get_next();
        }
    }
    if (n > 1) {
        if (n != last_prime) result.push_back ({n, 1});
        else result[result.size()-1].mult ++;
    }
    return result;
}

template <typename T>
std::vector<T> decompo<T>::GetDivisors (T n) {
    std::vector<T>  div;
    auto primes = decomp (n);
    int n_primes = primes.size();
    std::vector<int> exponent (n_primes, 0);
    div.push_back(1);
    int current_index = 0;
    int product = 1;
    std::vector<int> product_partial(n_primes, 1);;

    while (true) {
        current_index = 0;
        while (current_index < n_primes && exponent[current_index] == primes[current_index].mult) current_index++;
        if (current_index == n_primes) break;
        for (int index = 0; index < current_index; ++index) {
            exponent[index] = 0;
            product /= product_partial[index];
            product_partial[index] = 1; 
        }
        exponent[current_index]++;
        product *= primes[current_index].val;
        product_partial[current_index] *= primes[current_index].val;
        if (product != n && product != 1) div.push_back (product);      
    }

    return div;
}


//  return true if sum is obtained

bool test_sum (std::vector<int>& arr, int amount) {
    int n = arr.size();

    std::sort(arr.begin(), arr.end(), std::greater<int>());
    std::vector<int> bound (n);
    std::vector<int> select (n);
    bound[n-1] = arr[n-1];
    for (int i = n-2; i >= 0; --i) {
        bound[i] = bound[i+1] + arr[i];
    }

    int sum = 0;    // current sum
    int i = 0;      // index of the coin being examined
    bool up_down = true;

    while (true) {
        if (up_down) {
            if (i == n || sum + bound[i] < amount) {
                up_down = false;
                i--;
                continue;
            }
            sum += arr[i];
            select[i] = 1;
            if (sum == amount) return true;
            if (sum < amount) {
                i++;
                continue;
            }
            up_down = false;
            if (select[i] == 0) i--;
        } else {            // DOWN
            if (i < 0) break;
            if (select[i] == 0) {
                i--;
            } else {
                sum -= arr[i];
                select[i] = 0;
                i++;
                up_down = true;
            }
        }
    }
    return false;
}

bool isDeficient(int divisorSum, int n) {
    return divisorSum < n;
}

bool isAbundant(int divisorSum, int n) {
    return divisorSum > n;
}

bool isSemiPerfect(std::vector<int> &divisors, int sum) {
    int size = divisors.size();

    //  The value of subset[i][j] will be true if there is a subset of divisors[0..j-1] with sum equal to i 
    //bool subset[sum+1][size+1];
    std::vector<std::vector<bool>> subset(sum+1, std::vector<bool> (size+1));

    // If sum is 0, then answer is true 
    for (int i = 0; i <= size; i++) {
        subset[0][i] = true; 
    }

    //  If sum is not 0 and set is empty, then answer is false 
    for (int i = 1; i <= sum; i++) {
        subset[i][0] = false; 
    }

    // Fill the subset table in bottom up manner 
    for ( int i = 1 ; i <= sum ; i++ ) {
        for ( int j = 1 ; j <= size ; j++ ) {
            subset[i][j] = subset[i][j-1];
            int test = divisors[j-1];
            if ( i >= test ) {
                subset[i][j] = subset[i][j] || subset[i - test][j-1]; 
            }
        } 
    } 
    return subset[sum][size];
}

bool isWeird(int n, std::vector<int> &divisors, int divisorSum) {
    //return isAbundant(divisorSum, n) && !isSemiPerfect(divisors, n);
    return isAbundant(divisorSum, n) && !test_sum(divisors, n);
}

std::vector<int> getProperDivisors(int number) {
    std::vector<int> divisors;
    long sqrtn = sqrt(number);
    divisors.push_back(1);
    for ( int i = 2 ; i <= sqrtn ; i++ ) {
        if (number % i == 0) {
            divisors.push_back(i);
            int div = number/i;
            if (div != i) divisors.push_back(div);
        }
    }
    return divisors;
}

int main() {        
    decompo <int> decomposition;
    decomposition.complete (1e3);       // not relly useful

    int n = 2, count;
    std::vector<int> weird;
    int Nweird = 15000;

    for (count = 0; count < Nweird; n += 2) {
        if (n % 6 == 0) continue;

        //auto properDivisors = getProperDivisors(n);
        auto properDivisors = decomposition.GetDivisors(n);
        int divisorSum = std::accumulate (properDivisors.begin(), properDivisors.end(), 0);

        if (isDeficient(divisorSum, n) ) {
            continue;
        }

        if (isWeird(n, properDivisors, divisorSum)) {
            //std::cout << count << " " << n << "\n";
            weird.push_back (n);
            count++;
        }
    }
    for (int i = Nweird - 10; i < Nweird; ++i) {
        std::cout << weird.at(i) << " ";
    }
    std::cout << "\n";
}
user_105533410
Related