我需要一个更快的解决方案来解决以下问题

Question

您将获得N 一对数字n、k。对于每一对，计算P 的除数。

P = k^n * (1 + 2 + 3 + ... + k)

从nr_div_huge.in 读取N，然后读取N 对n and k。在输出文件nr_div_huge.out 中，每一行都写着每对P 的除数。因为它可能很大，所以会以 1.000.000.007 为模显示。

• 1≤N≤15.000
• 1≤k，n≤1.000.000.000
• 内存限制：0.1MB
• 时间限制：0.2s

例子：

nr_div_huge.in

2
2 3
4 4

nr_div_huge.out

8
20

说明：

对于第一对 P=54 和 54 有 8 个除数。 第二对 P=2560 和 2560 有 20 个除数。

这是我的代码：

#include <iostream>
#include <fstream>
#define MOD 1000000007
using namespace std;

ifstream fin ("nr_div_huge.in");
ofstream fout("nr_div_huge.out");

bool sieve[30000];
long long int prime[3000];
int N, k, n, r, P;

//making the Sieve of Eratosthenes have in the `prime` array the prime numbers up to 30000
void Erat()
{
    int r=0;
    sieve[0]=sieve[1]=1;

    for(int i=2; i*i<=30000; i++)
        if(sieve[i]==0)
            for(int j=2; j<=30000/i; j++)sieve[i*j]=1;

    for(int i=1; i<=30000; ++i)
        if(!sieve[i])prime[++r]=i;
}

//finding the numbers of divisors with Euler formula
int divisors(unsigned int n, unsigned int power)
{
    int r=0, d, p, nr=1;
    d=prime[++r];

    while(n>1)
    {
        p=0;
        while(n%d==0)
        {
            p++;
            n/=d;
        }

        if(d==2 && p>0)
            nr=(1LL*nr*(p*power))%MOD;
        else if(p>0)
            nr=(1LL*nr*(p*power+1))%MOD;

        if(prime[r+1]==0)d++;
        else d=prime[++r];

        if(d*d>n)d=n;
    }

    return nr;
}
/*
    I've transformed the formula P = k^n * (1 + 2 + 3 + ... +k) into P= k^(n+1)*(k+1)/2

    Also, I've used this property:
        If we have 2 numbers a and b, which are co-prime numbers, a will have n divisors and b m divisors,
        the (a*b) number will have m*n divisors.

*/
int NrP(long long int n, long long int k)
{
    unsigned long long int power=1,nrdk1,nrdkn,res=0;

    nrdk1=divisors(k+1,1);
    nrdkn=divisors(k,n+1);
    res=(1LL*nrdk1*nrdkn)%MOD;
    return res;
}


int main()
{
    Erat();
    fin>>N;

    for(int i=1; i<=N; i++)
    {
        fin>>n>>k;
        fout<<NrP(n,k)<<'\n';
    }
    return 0;
}

算法说明：

我一开始就对 Eratosthenes 进行筛选，以使 prime 数组中的素数达到 30000。为了更快地计算1 + 2 + 3 + ... + k 的总和，我使用公式(k*(k+1))/2 并且P 变为(k^(n+1) * (k+1))/2。由于k 和k+1 是互质的，那么k^(n+1) 和k+1 也将是互质的，其中一个是偶数。为了消除/2，我将偶数相除。因为数字是互质数，所以除数是n*m，其中n是k^(n+1)的除数，m是k+1的除数。

为了计算除数的个数，我使用了欧拉公式：

n = p1^e1 * p2^e2 * … * pk^ek - prime factorization

number of divisors = (e1 + 1) * (e2 + 1) * ... * (ek + 1)

如何加快解决速度？或者有没有更快的解决方案？

Answer 1

素因数分解k。假设k 是2*2*2*3*5 = 120。
所以k 的任何除数必须是2^x * 3^y * 5^z，其中x、y、z 可以分别是[0,1,2]、[0,1]、[0,1]。 同样，k^n 的除数必须与2^x*3^y*5^z 相同，其中x、y、z 可以分别为[0,1 .. 3*n]、[0, .. n]、[0, .. n]。 那么，您可以使用这些元素形成的独特组合的数量是所有选择的乘积：(3*n+1)*(n+1)*(n+1)

（如果不清楚，请考虑计算圣代冰淇淋独特组合的数量的类比，其中冰淇淋圣代由 3 种可能的冰淇淋口味中的 1 种和 4 种可能的冰淇淋配料中的 1 种制成......那么唯一圣代的总数是 3*4。在这里，我们的选择是使用多少个2s，使用多少个3s，以及使用多少个5s。每个独特的组合产生一个独特的产品。）

还将(k * (k+1)) / 2 的素数添加到素数计数中，它应该可以工作。

编辑：我看到您提出的解决方案已经在使用欧拉函数 (n = p1^e1 * p2^e2 * … * pk^ek - prime factorization)。

此处概述的解决方案也使用了它，但缺少的见解是k^2 的素数分解很容易从k 的素数分解推导出来：12 的素数分解是2^2 * 3， 12*12 = 144 的质因式分解是 2^2 * 3 * 2^2 * 3 = 2^4 * 3^2。因此，您不需要计算k^n。您只需对k 进行质因式分解，然后使用此数据推导出k^n 的质因数分解。

== 编辑 2 == 在那之后你仍然有错误，所以也许你的一些计算值ints 需要是long longs？另外，我看到你的筛子上升到 3000，但它不应该上升到 sqrt(10^9) 吗？

Answer 2

你必须先做一点数学才能使问题更简单。如果不这样做，您将结束计算不会包含在 long long 变量中的数字，这需要像 gmp 这样的多精度库并且花费太多时间。

但是一旦你观察到 1+2+...+k 就是 k*(k+1)/2，你会得到：

P = k^(n+1) * (k+1) / 2

而k 和k+1 是互质数（连续整数的平凡属性）。然后，您只需将k 和k+1 分解为素数。其中一个可以被 2 整除。但这让您甚至无需计算 P 的主要因素分解。

从那里，你有P = a0^b0 * ... * ax^bx，除数是(b0+1)*...*(bx+1)。

你的例子的应用：

• 2,3: P = 3^3 * 4 / 2 = 2^1*3^3 - 除数数 2 * 4 = 8
• 4,4: P = 4^5 * 5 / 2 = 5 * (2^2)^5 / 2 = 2^9 * 5^1 - 除数数 10 * 2 = 20

您仍然需要用 C++ 编写算法，但它会比大数字运算更有效率。

---

这是一个可能的 C++ 实现

#include <vector>
#include <iostream>
#include <fstream>

// Will search all primes up to (SQSIZE * SQSIZE) - 1
constexpr const unsigned long SQSIZE = 200;

std::vector<unsigned> make_sieve(unsigned sqsize) {
    unsigned size = sqsize * sqsize;
    std::vector<unsigned> sieve(size, 0);
    // Classic Eratosthene
    for (unsigned i=2; i<sqsize; i++) {
        if(sieve[i] == 0) {
            for(unsigned j = i*i; j<size; j+=i) {
                sieve[j] = 1;
            }
        }
    }
    unsigned len=0;
    // Pack prime number in the vector
    for (unsigned i=2; i<size; i++) {
        if (sieve[i] == 0) sieve[len++] = i;
    }
    sieve.resize(len);
    return sieve;
}

unsigned long n_divisors(unsigned long k, unsigned long n, const std::vector<unsigned> &sieve) {
    unsigned long n_div = 1;
    for(unsigned i: sieve) {
        unsigned exp = 0;
        while ((k % i) == 0) {
            exp += 1;
            k /= i;
        }
        // special processing if k divisible by 2: use k^n/2
        if ((i == 2) && (exp != 0)) n_div *= (exp * n - 1) + 1;
        else n_div *= exp * n + 1;
        if (k == 1) break;
    }
    if (n_div == 1) n_div = 2;  // k is prime and greater sieve max
    return n_div;
}

int main() {
    std::vector<unsigned> sieve = make_sieve(SQSIZE);
    unsigned long k, n, n_div;
    unsigned N;

    std::ifstream fin ("nr_div_huge.in");
    std::ofstream fout("nr_div_huge.out");

    fin >> N;
    for(unsigned i=0; i<N; i++) {
        fin >> n >> k;
        n_div = n_divisors(k, n+1, sieve) * n_divisors(k + 1, 1, sieve);
        fout << n_div << '\n';
    }
    return 0;
}

这里不需要long long。 n 和 k 小于 1000000000 并适合 32 位整数。所以使用unsigned long。

注意：上面的代码盲目假设没有 IO 问题并且从不测试它们。永远不要对生产级代码这样做...