欧拉项目 549题 | 阶乘的整除性问题

题目的原始链接

定义s(n) = m，给定一个整数n，m是最小的整数，使得m!能够整除n。
举两个简单的例子，s(10)=5，s(25)=10。
定义S(n) = ∑s(i) 2 ≤ i ≤ n
题目中给了一个当n = 100的值用于测试：S(100)=2012

求S(10^8)

首先要意识到，不可能把阶乘本身算出来，比如1517，由质数37和41相乘得到，那么对应的m应该是41，但是41的阶乘非常之大。
给定一个数n，求满足题意的m，因为这里肯定需要遍历所有的n，10 ^ 8这么大了，所以求m的函数的时间复杂度最好是常量时间，nlg(n)也很快，n * sqrt(n)就太慢了，这都到10 ^ 12次了，复杂度已经高的吓人了，计算机真心无法短时间算出来了，所以我们下面的目标就是在合理的时间复杂搞定s(n)这个函数：一定要控制在nlg(n)以内。
首先，如果这个数n是质数，那么m肯定就是n本身。
对于一个非质数，一个可行的思路就是把n分解质因数，质因数出现的次数也是需要保留的，要不然25这个数，是5*5，如果不保留质因数的次数，那么计算得到的m就成了5了，显然是不对的。然后再如何处理呢？
假设分解质因数结果是p1 ^ n1 * p2 ^ n2 * … * pq^ nq
想办法凑够n1个p1，p1包含1个p1，2*p1又包含一个（如果p1等于2，这里还会多包含一个），3*p1又包含一个p1，依次类推。
再想办法凑够n2个p2。
依此类推。
共计q个结果，最大的那个数字就是我们要求的m。

以上想法计算S(100)没啥问题。但是：虽然我把之前分解质因数的函数性能提高了好几倍，但是离nlg(n)还是差很远，所以导致短时间无法运算完毕。。。

如果我能过滤掉很多值不用去计算了，那么时间就能大大降低了。
首先10 ^ 8以内的质数占了5.7%的样子，所以这5.7%是不用计算的。
其次，考虑某个质数p，乘以小于它的某个数m，得到的这个值n，S(n)也是p。写一个两层for循环就能得出这些n值，这占了大概66.7%，这一下子就节省了大多数的时间。
再次，在上一条的基础上再进一步，某个质数p，乘以两个小于它的数m1和m2，得到的n值，S(n)也是p。三层for循环能搞定，这大概占了22.3%，又能节约挺多时间的了。
这里可以做一个小的优化，m2可以从m1+1和p/m1+1两者比较大的一方开始，因为如果m2*m1小于p，那么这个数字肯定在上一步的时候就被标记了。优化之后，三层for循环，用时不到2s。
再再次，还能再排除3.4%的数据，这里也可以使用上面提到的优化思想。
最后，只剩约1.9%的数据需要按照之前说的逻辑进行判定了。

生成1亿以内的质数用时约3.9s，中间若干步耗时约5s，最后处理1.9%的数据耗时10s。。。
其实沿着中间那一段优化的想法，可能可以在线性时间复杂度做完，但是想了半天，没有太好的想法，暂时就这样子吧。

下面是我写的代码：
话说这次代码又丑又长，特别是中间的几个for循环嵌套。。。

public static long GetAnswer()
{
    var primes = Utils.GenPrimes(N);

    var result = new long[N + 1];
    foreach (var p in primes)
    {
        result[p] = p;
    }

    for (int index = 0; primes[index] < N / 2; index++)
    {
        for (int j = 2; j < primes[index]; j++)
        {
            var tmp = j * primes[index];
            if (tmp > N)
            {
                break;
            }
            else
            {
                result[tmp] = primes[index];
            }
        }
    }

    for (int index = 0; primes[index] < N / 6; index++)
    {
        for (int j = 2; j < primes[index] - 1; j++)
        {
            long k = Math.Max(j + 1, primes[index] / j);
            var test = j * k * primes[index];
            if (test < 0 || test > N)
            {
                break;
            }

            for (; k < primes[index]; k++)
            {
                var tmp = j * k * primes[index];
                if (tmp < 0 || tmp > N)
                {
                    break;
                }
                else
                {
                    result[tmp] = primes[index];
                }
            }

        }
    }

    for (int index = 0; primes[index] < N / 24; index++)
    {
        for (int j = 2; j < primes[index] - 2; j++)
        {
            long k = Math.Max(j + 1, primes[index] / j);
            var test = j * k * (k + 1) * primes[index];
            if (test < 0 || test > N)
            {
                break;
            }

            for (; k < primes[index] - 1; k++)
            {
                long m = Math.Max(k + 1, primes[index] * primes[index] / (j * k));
                var test2 = j * k * m * primes[index];
                if (test2 < 0 || test2 > N)
                {
                    break;
                }

                for (; m < primes[index]; m++)
                {
                    var tmp = j * k * m * primes[index];
                    if (tmp < 0 || tmp > N)
                    {
                        break;
                    }
                    else
                    {
                        result[tmp] = primes[index];
                    }
                }
            }
        }
    }

    for (int i = 2; i <= N; i++)
    {
        if (result[i] == 0)
        {
            var facs = Utils.TrialDivisioFac(i, primes);
            var maps = GetFactorsCounts(facs);
            long m = long.MinValue;
            foreach (var map in maps)
            {
                var temp = GetMaxMByPrime(map);
                if (temp > m)
                {
                    m = temp;
                }
            }

            result[i] = m;
        }
    }

    return result.Sum();
}

private static long GetMaxMByPrime(KeyValuePair<long, int> map)
{
    long result = 0;
    long count = map.Value;
    while (count > 0)
    {
        result += map.Key;
        long cur = result;
        while (cur > 0)
        {
            if (cur % map.Key == 0)
            {
                cur /= map.Key;
                count--;
            }
            else
            {
                break;
            }
        }
    }

    return result;
}

private static Dictionary<long, int> GetFactorsCounts(List<long> facs)
{
    var maps = new Dictionary<long, int>();
    foreach (var fac in facs)
    {
        if (maps.ContainsKey(fac))
        {
            maps[fac]++;
        }
        else
        {
            maps[fac] = 1;
        }
    }

    return maps;
}