如何有效地计算奇怪的数字

Question

我正在尝试打印 n weird numbers 其中 n 是非常大的数字（例如：10000）。

如果我有一些错误，我发现这个网站来检查 n 600 的算法： http://www.numbersaplenty.com/set/weird_number/more.php

但是，我的算法在大数字时真的很慢：

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;
    }

}

我有三个简单的突破：

1. 如果一个数字可以被 6 整除，那么它就是semiperfect，这意味着它不会很奇怪

2. 如果一个数字是deficient，这意味着它不会很奇怪

以上几点基于https://mathworld.wolfram.com/DeficientNumber.html

1. 如果 a 数字是奇数，至少对于 10^21 个数字来说它不会很奇怪（这对我要获取的数字很有用）。

我使用的另一个优化是查找数字的所有除数的优化。我们不是循环到 n，而是循环到 SQRT(n)。

但是，我仍然需要优化： 1. isSemiPerfect 因为它真的很慢 2.如果能进一步优化getProperDivisors就好了。

欢迎提出任何建议，因为我找不到更多优化方法来在合理的时间内找到 10000 个奇怪的数字。

PS：Java、C#、PHP 和 JavaScript 中的任何代码对我来说都可以。

编辑：我找到了this topic 并修改了 isSemiPerfect 看起来像这样。但是，它看起来并没有优化而是减慢了计算速度：

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;
      }

Answer 1

问题确实出在函数 isSemiPerfect 上。我用 C++ 转换了你的代码，它仍然很慢。

然后我使用回溯修改了这个函数。我现在在大约 15 秒内获得了前 15000 个奇怪的值。我的解释是，在几乎所有情况下，该值都是半完美的，并且回溯函数快速收敛。

还要注意，在我的回溯实现中，我对除数进行排序，这样可以减少要检查的案例数量。

编辑 1：在 getProperDivisors 中纠正了一个错误。最终结果似乎没有被修改！

#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 完全重新定义了除数的生成。它现在使用素数分解。更复杂，但全球时间除以 7.5。在我的电脑上生成奇怪的数字现在需要 2 秒。

#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";
}