欧拉项目 | 83题 | 最短路径

发表于 2019-03-10 分类于数学，欧拉项目 Disqus：

一个整数矩阵，从左上角到右下角，存在一条最短路径。走的方式没有要求，在某个点，可以往四个方向走；之前有两个题目，一个是限制只能往右或者往下，稍微复杂点的另一个题目的限制是不能往左。

我的想法很直接，构造一个矩阵，是从左上角到该点的最短路径，每个点的初始值都是int.MaxValue。
用一个queue来保存被更新的点，然后看看其周围的点是否可以被更新。
从左上角的点开始，赋值成给定矩阵左上角的值。该值被更新了，把(0, 0)点放入queue里，因为周围的点可能会被更新。
然后从queue拿出一个点，看看周围的值是否能被更新，如果是，也放到queue里。
循环往复，直到queue为空，没有需要被更新的点了。
此时，每个点上的值都是从左上角开始到该点的最短路径。

代码直接反应了上述的过程，data就是题目总给出的80 * 80的矩阵。

private const int N = 80;

public static int GetAnswer()
{
    int[,] sum = new int[N, N];
    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < N; j++)
        {
            sum[i, j] = int.MaxValue;
        }
    }
    sum[0, 0] = data[0, 0];
    Queue<Tuple<int, int>> pointsToBeUpdated = new Queue<Tuple<int, int>>();
    pointsToBeUpdated.Enqueue(new Tuple<int, int>(0, 0));
    while (pointsToBeUpdated.Count != 0)
    {
        var point = pointsToBeUpdated.Dequeue();
        if (CanMoveToLeft(point))
        {
            int newValue = sum[point.Item1, point.Item2] + data[point.Item1 - 1, point.Item2];
            if (newValue < sum[point.Item1 - 1, point.Item2])
            {
                sum[point.Item1 - 1, point.Item2] = newValue;
                pointsToBeUpdated.Enqueue(new Tuple<int, int>(point.Item1 - 1, point.Item2));
            }
        }

        if (CanMoveToRight(point))
        {
            int newValue = sum[point.Item1, point.Item2] + data[point.Item1 + 1, point.Item2];
            if (newValue < sum[point.Item1 + 1, point.Item2])
            {
                sum[point.Item1 + 1, point.Item2] = newValue;
                pointsToBeUpdated.Enqueue(new Tuple<int, int>(point.Item1 + 1, point.Item2));
            }
        }

        if (CanUp(point))
        {
            int newValue = sum[point.Item1, point.Item2] + data[point.Item1, point.Item2 - 1];
            if (newValue < sum[point.Item1, point.Item2 - 1])
            {
                sum[point.Item1, point.Item2 - 1] = newValue;
                pointsToBeUpdated.Enqueue(new Tuple<int, int>(point.Item1, point.Item2 - 1));
            }
        }

        if (CanDown(point))
        {
            int newValue = sum[point.Item1, point.Item2] + data[point.Item1, point.Item2 + 1];
            if (newValue < sum[point.Item1, point.Item2 + 1])
            {
                sum[point.Item1, point.Item2 + 1] = newValue;
                pointsToBeUpdated.Enqueue(new Tuple<int, int>(point.Item1, point.Item2 + 1));
            }
        }
    }

    return sum[N - 1, N - 1];
}

private static bool CanMoveToLeft(Tuple<int, intpoint)
{
    return point.Item1 > 0;
}

private static bool CanMoveToRight(Tuple<int, intpoint)
{
    return point.Item1 < N - 1;
}

private static bool CanUp(Tuple<int, int> point)
{
    return point.Item2 > 0;
}

private static bool CanDown(Tuple<int, int> point)
{
    return point.Item2 < N - 1;
}

欧拉项目 | 119题 | 某数各数字之和的若干次幂等于该数本

发表于 2019-03-03 分类于数学，欧拉项目 Disqus：

Problem 119

512是一个很有意思的数字，5+1+2=8，8^3=512
614656 = 28^4
an是满足这些条件的数列：
至少有两位数
各个数字之和的若干次幂等于该数
求a30

对于这种题目，欧拉项目有个套路，你不太可能遍历每个数字去检查它是否满足条件的时候，就要换种思路，快速构造满足题目的集合，然后按照题目找到答案。
回到这个题目，我个人感觉a30不会超过long能表示的范围，所以我只要看long以内的就可以了。
第一步，找到数字之和s的范围。最小值是2，因为至少有两位数；最大值是9*19，因为long.MaxValue有19位数字。其实这个范围比实际范围大一些，但是无所谓了，多的不多，不会影响程序效率。
第二步，求s的若干次幂，从2次幂开始，3次幂，4次幂等等，直到超出long的范围。这些都是可能的值n。满足这一步的n不到1700个。
第三步，求n的各个数字之和sum，如果sum等于s，那么这个n就是满足题意的。
第四步，把所有的n排序，找出第三十个即可。很幸运，long以内有34个满足题目的数。

下面是代码

var candidates = new List<long>();
for (int i = 2; i < 9 * long.MaxValue.ToString().Length; i++)
{
    int power = 2;
    while (true)
    {
        var ret = BigInteger.Pow(i, power);
        if (ret > long.MaxValue)
        {
            break;
        }

        long number = (long)ret;
        long candidate = number;

        long sum = 0;
        while (number != 0)
        {
            sum += number % 10;
            number /= 10;
        }

        if (sum == i)
        {
            candidates.Add(candidate);
        }

        power++;
    }
}

candidates.Sort();

return candidates[29];

欧拉项目 | 104题 | 头尾是不重复数字组成的斐波那契数

发表于 2018-12-18 分类于数学，欧拉项目 Disqus：

104题链接

这里就不赘述什么是斐波那契数列了。
f(541)包含113个数字，最后九个数字恰好是由数字1-9组成的；f(2749)包含575个数字，前九个数字恰好是由数字1-9组成的。
求k，f(k)的前九个数字和最后九个数字都是由数字1-9组成的。

我一开始的思路是使用BigInteger来存放斐波那契数，逐个计算上去，判断前九个数字和最后九个数字是否是由1-9组成的。
不过我发现程序很慢，用Visual Studio分析了一下程序的性能，发现BigInteger.Log10和BigInteger.Pow非常耗时。我用BigInteger.Log10计算出数字的长度L，然后除以BigInteger.Pow(10, L - 9)得到最前面的九个数字。
稍微改变了一下策略，先判断最后九个数字，后判断前九个数字。速度提高很多，因为判断最后九个数字，只要模10的九次方就行了，比计算长度再除快非常多，而需要进行后一个判断的几率很小。于是乎，大约30s能解决。
能不能再快一点呢？我发现BigInteger的模运算也挺慢的。其实我可以用一个长度为九的数组进行最后九位的加法，然后判断是否是由1-9九个数字组成就好。

private static int[] Add(int[] fn_1, int[] fn_2)
{
    int[] ret = new int[9];
    int n1 = fn_1.Length - 1;
    int n2 = fn_2.Length - 1;
    int nr = ret.Length - 1;
    int carry = 0;
    for (; nr >= 0; n1--, n2--, nr--)
    {
        int a = n1 >= 0 ? fn_1[n1] : 0;
        int b = n2 >= 0 ? fn_2[n2] : 0;
        int sum = a + b + carry;
        ret[nr] = sum % 10;
        carry = sum / 10;
    }

    return ret;
}

判断是否是由数字1-9组成也很简单。

private static bool IsPandigitalFibonacci(int[] fn)
{
    int[] digitalCounts = new int[10];
    for (int i = 0; i < 9; i++)
    {
        digitalCounts[fn[i]]++;
    }

    for (int i = 1; i < 10; i++)
    {
        if (digitalCounts[i] != 1)
        {
            return false;
        }
    }

    return true;
}

使用这个方式来代替原来的想法，得到了完整的算法。

int[] fn_2 = new int[] { 1 };
int[] fn_1 = new int[] { 1 };
BigInteger f1 = 1;
BigInteger f2 = 1;
int i = 3;
while (true)
{
    BigInteger fnBig = f1 + f2;
    int[] fn = Add(fn_1, fn_2);
    if (IsPandigitalFibonacci(fn))
    {
        int digits = (int)Math.Floor(BigInteger.Log10(fnBig) + 1);
        var dividend = BigInteger.Pow(10, digits - 9);
        var firstNineDigits = fnBig / dividend;
        if (Utils.IsPandigital(firstNineDigits.ToString(), false))
        {
            return i;
        }
    }

    fn_2 = fn_1;
    fn_1 = fn;
    f2 = f1;
    f1 = fnBig;
    i++;
}

这比以前又快了很多，大约5s就能搞定。
分析下程序的瓶颈，发现最耗时的是BigInteger.Add，这个没办法优化了，除非能有个算法能够直接得到前九个数字是啥，这样就不要BigInteger来保存结果，我觉得没有这种算法了；次耗时的是BigInteger.Pow，它的调用次数已经优化了，而且我也没有想到更好的构造被除数的方式。
其实五秒已经挺快的了，很满意。

欧拉项目 | 66题 | 丢番图方程

发表于 2018-12-15 分类于数学，欧拉项目 Disqus：

66题链接

二次丢番图方程如下

$x^2 – Dy^2 = 1$ 当D = 13时，x的最小解是

$649^2 – 13*180^2 = 1$ 当D是平方数时，无解。
题目中给出了几个示例，当D = {2, 3, 5, 6, 7}时，列出了x的最小解，其中最小解的最大值是x = 9，对应的D = 5
找到这样一个D，使得x最小解是最大值，其中D ≤ 1000。

如果我们能对于某个D，快速找到最小解x。这个题就很容易了。

List<int> D = Enumerable.Range(1, 1000)
                        .Where(i => Math.Pow((int)Math.Sqrt(i), 2) != i)
                        .ToList();
BigInteger max = 0;
long ret = 0;
foreach (var d in D)
{
    BigInteger cur = GetX(d);
    if (max < cur)
    {
        max = cur;
        ret = d;
    }
}

return ret;

现在问题就是如何GetX。论文An Algorithm to Solve a Pell Equation给出了一种快速的解法。

接下来，照猫画虎，把算法实现一下。

private static BigInteger GetX(int n)
{
    List<BigInteger> a = new List<BigInteger>() { 0, 1 };
    List<BigInteger> b = new List<BigInteger>() { 1, 0 };
    List<int> c = new List<int>() { n, 1 };

    int i = 1;
    do
    {
        int q = (int)((Math.Sqrt(n - c[i - 1] * c[i]) + Math.Sqrt(n)) / c[i]);
        c.Add((int)(2 * q * Math.Sqrt(n - c[i - 1] * c[i]) + c[i - 1] - q * q * c[i]));
        a.Add(q * a[i] + a[i - 1]);
        b.Add(q * b[i] + b[i - 1]);
        i++;
    } while (a.Last() * a.Last() - n * b.Last() * b.Last() != 1);

    return a.Last();
}

为什么要用BigInteger类型呢？虽然D不大，但是x和y的值可能会很大，甚至超过long的表示范围。

欧拉项目-613题-毕达哥拉斯蚂蚁

发表于 2018-10-27 分类于数学，欧拉项目 Disqus：

613题链接

题目是毕达哥拉斯蚂蚁，用我们的说法是勾股蚂蚁。
题目大概是说，有个小蚂蚁落到了一个宽30cm，长40cm的直角三角形上，它会随机的选择一个方向爬，问从斜边爬出去的概率是多少。

上图是一个示意图，x轴的长度是4，y轴长度是3。把题目的30，40缩放为3，4，对求概率是没有影响的。
假设落到点P(x, y)，那么它从斜边爬出去的概率是多少呢？
角PAC的大小是arctan(x/(3-y))，角PBC的大小是arctan(y/(4-x))，那么角APB的大小是PI/2+arctan(x/(3-y))+arctan(y/(4-x))，那么从点P开始从斜边爬出去的概率就是PI/2+arctan(x/(3-y))+arctan(y/(4-x)) / (2*PI)。对三角形中的每一点计算概率然后除以三角形面积就是要求的概率了。

1 / 12 * PI * (integrate PI/2+arctan(x/(3-y))+arctan(y/(4-x)) dy dx, x=0..4, y=0..-3*x/4+3)

PI/2是常量，可以很容易计算出来，公式化简为

1 / 12 * PI * (3 * PI + (integrate arctan(x/(3-y))+arctan(y/(4-x)) dy dx, x=0..4, y=0..-3*x/4+3))

两个反正切，我不会算。。。借助了一下wolframalphal

有了解析解之后，继续化简公式

1 / 12 * PI * (3 * PI + 3 * PI + ln((3^4*2^16*sqrt(3/5))/(5^12)))

1/2 + PI/12*ln((3^4*4^8*sqrt(3/5))/(5^12))

上式包含了3，4，5，同时分子和分母幂次相同，3的4次方乘以4的8次方，分母是5的12次方，是不是很神奇哈哈。其实想想就知道结果和3，4，5是相关的。

有了解析式，计算出答案也就很容易了。

欧拉项目 | 96题 | 数独

发表于 2018-10-20 分类于数学，欧拉项目 Disqus：

96题链接

数独，一个很经典的填数字游戏。规则也很简单，每一行、每一列、九个九宫格的九个数字不重复。
设计很好的数独游戏只有一个解，题目中也假设了这一点。
题目最后要求五十个数独的解，然后计算左上角三位数之和。
这个题对我来说还算比较简单，只需要把我人肉做数独游戏的方法写成程序就好了。

由于整体思路非常清晰，我采用了自上而下的设计方式，几乎一气呵成了整个程序，除了debug了一个小地方fix了一个bug外。
我需要一个函数把题目给的txt下载下来并一行一行的读到一系列二位数组中。
对于每个二维数组，就是一个数独题目，然后求解，得到左上角的三位数。求和。

private static readonly int Length = 9;

public static int GetAnswer()
{
    int sum = 0;
    var lines = GetContent();
    for (int i = 0; i < 50; i++)
    {
        int[,] soDoku = new int[Length, Length];
        for (int j = 0; j < Length; j++)
        {
            var line = lines[10 * i + j + 1];
            for (int k = 0; k < Length; k++)
            {
                soDoku[j, k] = Convert.ToInt32(line[k].ToString());
            }
        }

        sum += GetTopLeft(SolveSuDoku(soDoku));
    }

    return sum;
}

GetContent()函数很简单，下载文件，切分行。如果愿意的话，甚至能写成一行。

private static string[] GetContent()
{
    WebClient client = new WebClient();
    string content = client.DownloadStri("https://projecteuler.net/project/resources/p096_sudoktxt");
    return content.Split('\n');
}

GetTopLeft(int[,] soDoku)函数相对也比较简单，拿到二维数组的左上角三个数字，拼成string然后转成int即可。或者是第一位乘以100加第二位乘以10加第三位。

private static int GetTopLeft(int[,] soDoku)
{
    return int.Parse($"{soDoku[0, 0]}{soDoku[0, 1]}{soDoku[0, 2]}");
}

最复杂的部分就是解数独了。
解数独就是在不停的猜测和测试，整个过程显然是递归的。
找个一个空，可能有几个情况，选一个填进去，依旧是未完成的数组，再来一轮，一轮，直到完成。整个过程是一个树，分了很多叉。题目告诉我们只有一个解，那么只有一个叶子节点是完整的解，其他叶子节点都不合法，有矛盾的地方，比如某个点所在的行和列出现了重复数字。
首先，我们判断数独是否完成，如果是，则返回该解。
如果没有，我们需要找一个点，该点没有被填写数字，且可能填的数字尽可能的少，这种情况下该节点的子树才最少，效率能尽可能的高。FindMinimalPoint返回参数包含三项，前两者描述的该节点的位置，第三项是可能的值。如果Item3是空的话，那么说明当前的数独有矛盾，返回null。试图尝试所有的可能性，如果有解，则返回，如果所有的可能性都无解的话，走到最后一行返回null。

private static int[,] SolveSuDoku(int[,] soDoku)
{
    if (IsCompleted(soDoku))
    {
        return soDoku;
    }

    var point = FindMinimalPoint(soDoku);
    foreach (var v in point.Item3)
    {
        int[,] newSoDoku = soDoku.Clone() as int[,];
        newSoDoku[point.Item1, point.Item2] = v;
        var ret = SolveSuDoku(newSoDoku);
        if (ret != null)
        {
            return ret;
        }
    }

    return null;
}

如果实现FindMinimalPoint呢？遍历非0点，找到可能的值（FindVaulesAt），记录下可能的值最少的点坐标和对应的值并返回。对于FindVaulesAt函数，很可能返回0个解，就是说该数独有矛盾。minList在初始化的时候写了十个值，目的是比九大，其实这里写九个也没问题，因为最小值肯定比九要小。

private static (int, int, List<int>) FindMinimalPoint(int[,] soDoku)
{
    int minRow = -1;
    int minCol = -1;
    var minList = new List<int>() { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
    for (int i = 0; i < Length; i++)
    {
        for (int j = 0; j < Length; j++)
        {
            if (soDoku[i, j] == 0)
            {
                var values = FindVaulesAt(i, j, soDoku);
                if (values.Count < minList.Count)
                {
                    minList = values;
                    minRow = i;
                    minCol = j;
                }
            }
        }
    }

    return (minRow, minCol, minList);
}

FindVaulesAt想法很直接。每个空可能填入一到九九个数字，排除所在行、所在列和所在的九宫格已经出现的数字即可。
文章开头提到fix了一个小bug就是在这个函数，一开始row的初始值写成了i % 3，row < i - i % 3 + 3，进而没能得到正确的可能的值。

private static List<int> FindVaulesAt(int i, int j, int[,] soDoku)
{
    List<int> values = new List<int> { 1, 2, 3, 4, 5, 6, 7, 8, 9 };
    for (int k = 0; k < Length; k++)
    {
        values.Remove(soDoku[i, k]);
    }

    for (int k = 0; k < Length; k++)
    {
        values.Remove(soDoku[k, j]);
    }

    for (int row = i - i % 3; row < i - i % 3 + 3; row++)
    {
        for (int col = j - j % 3; col < j - j % 3 + 3; col++)
        {
            values.Remove(soDoku[row, col]);
        }
    }

    return values;
}

最后，给出整个拼图的最后一块IsCompleted。如果所有空都不是0，说明数独已经求解完毕了。

private static bool IsCompleted(int[,] soDoku)
{
    foreach (var cell in soDoku)
    {
        if (cell == 0)
        {
            return false;
        }
    }

    return true;
}

欧拉项目 | 601题 | Divisibility streaks

发表于 2018-07-28 分类于数学，欧拉项目 Disqus：

601题链接

对于每个整数n我们可以定义一个函数streak(n)=k，其中k是最小的n+k不能被k+1整除的数。
比如
13 可以被 1 整除
14 可以被 2 整除
15 可以被 3 整除
16 可以被 4 整除
17 不可以被 5 整除
所以streak(13)=4
同样的
120 可以被 1 整除
121 不可以被 2 整除
所以streak(120)=1

P(s,N)是满足下面两个条件的n的和

1 < n < N
streak(n)=s

题目给了两个值P(3,14)=1 和 P(6,106)=14286帮助验证自己的想法。

求

$\sum_{i=1}^{31} P(i, 4^i)$

首先看下N的数量级，4^31=2^62，都接近long的数量级了，遍历小于N的每个数显然是不合适的。

整除k+1，其实就是被2，3，4等整除。手算了i=1，2，3，4的情况，发现有规律的：
比如i=4，那么到了2 3 4 5的最小公倍数往后的情况和之前的情况是一样的，循环起来了。
那么如果我能计算最小公倍数以内的情况，那么就能很容易得到全部的情况。
不过事与愿违，从2到32这31个数的最小公倍数是一个长达15的数字，虽然比2^62小好几个数量级，但是还是不能遍历。
这条路走不通。

各种想了2个小时之后，找到了一种直接计算的方式。

什么样的n满足streak(n)=s呢？
n+1 % 2 == 0，否则streak(n)=1
同理，n+2 % 3 == 0
类推到 n+s-1 % s == 0
最后 n+s % s+1 != 0
对于前面s-1个等式被整除的数都减去k，k=2,3…s，得到一个统一的式子
n-1 % k == 0
也就是说，n-1是2,3…s最小公倍数的若干倍。利用这一点，我们可以容易的找到可能的候选n，n的数量不会很多，对于i等于31的情况，粗略估计也只有四位数个吧。
最后利用n+s % s+1 !=0得到n，然后求和。

下面是P函数的代码

static long P(int s, long N)
{
    long lcm = Enumerable.Range(1, s)
        .Select(Convert.ToInt64)
        .Aggregate((cur, next) => Utils.GetLcm(cur, next));

    return Enumerable.Range(1, (int)((N - 2) / lcm))
        .Select(i => i * lcm + 1)
        .Count(c => (c + s) % (s + 1) != 0);
}

Enumerable.Range(1, s)，多了个1，对于最小公倍数没有影响。
N-2的原因是

1 < n < N

调用P函数得到和

1	Enumerable.Range(1, 31).Select(i => P(i, 1L << 2 * i)).Sum();

欧拉项目 | 516题 | 5-smooth totients

发表于 2018-07-27 分类于数学，欧拉项目 Disqus：

516题链接

5-smooth数是指质因数不超过5的数。换句话说，质因数只能是2，3，5。
S(L)是指不超过的L的n求和，其中n满足欧拉函数，Euler’s totient functionφ(n)是一个5-smooth数。
S(100)=3728. 题目给出这个提示是进行一个简单的验证，检查自己的算法/程序是否有问题。
求S(10^12) % 2^32

10^12次方很大，电脑是GHz级别，10^9，所以遍历小于L的每一个数是不现实的。我们要想办法构造出所有满足条件的n，然后求和。
φ(n)有很多种写法，下面这种对于我们分析这个题目很有帮助。

$\phi(n)=\prod_{i=1}^k p_i^{r_k-1}(p_i-1)$

$n=\prod_{i=1}^k p_i^{r_k}$

如果质因数只有2，3，5，减一之后是1，2，4，那么这个n是满足条件的。

质因数能大于5吗？
可以。比如7，7-1之后是6，质因数是2和3，如果这种质因数的次数r只能是1，r-1之后是0，那么也还是满足要求的。如果多余1次，φ(n)就会包含比5大的因数。再比如11，41等都是满足条件的。
这种质数是可以任意相乘的，组合之后仍然满足题意。他们和质因数只有235的数相乘，依旧满足题意。顺便说一下，只有这些质数作为质因数的数也是满足题意的。

如何构造这种质数呢？
找到所有质因数是235的数，加一，再判断是不是质数。

让我们先写个程序找到所有小于10^12且质因数都是235的数。

for (int i = 0; i < Math.Log(max) / Math.Log(2); i++)
{
    for (int j = 0; j < Math.Log(max) / Math.Log(3); j++)
    {
        for (int k = 0; k < Math.Log(max) / Math.Log(5); k++)
        {
            if (Math.Pow(2, i) * Math.Pow(3, j) * Math.Pow(5, k) <= max)
            {
                smooths.Add((long)(Math.Pow(2, i) * Math.Pow(3, j) * Math.Pow(5, k)));
            }
        }
    }
}

smooths.Sort();

有个微不足道的优化，j和k的上限可以限制的更小一些，但是没必要了，数不大。排序的目的后面会提到。

接下来，生成所有大于5且满足上述分析的质数。

var candidates =
    smooths
    .Select(l => l + 1)
    .Where(Utils.IsPrime)
    .Where(p => p > 5)
    .ToArray();

那么如果对这些质数进行组合呢？普通的组合算法是行不通的。因为candidates有500多个数据，从中任选5个就达到了10^12这个量级，所以需要合适的裁剪。
我的做法是分层，第一层都是1个数组成的组合，第二层都是2个数组成的组合，依此类推。同时，后一层是从前一层的组合生成的。对于i层的某一个组合，i+1层就是从smooths里面选一个数，然后加上去。smooths排序的好处有两个，一是找到i层组合的最后一个数，只需要往上添加比它大的数即可，排序可以利用二分查找；二是如果第j个加上去之后，乘积大于10^12，后面的数就不用看了。

List<List<List<long>>> tiers = new List<List<List<long>>>();
var firstTier = candidates.Select(p => new List<long>() { p }).ToList();
tiers.Add(firstTier);

while (tiers.Last().Count != 0)
{
    var lastTier = tiers.Last();
    var currentTier = new List<List<long>>();
    tiers.Add(currentTier);

    foreach (var combination in lastTier)
    {
        double product = combination.Aggregate((total, next) => total * next);

        int index = Array.BinarySearch(candidates, combination.Last());
        for (int i = index + 1; i < candidates.Length; i++)
        {
            if (product * candidates[i] <= max)
            {
                currentTier.Add(new List<long>(combination)
                {
                    candidates[i]
                });
            }
            else
            {
                break;
            }
        }
    }
}

所有的组合都找到了，很容易得到乘积。

var factors = new List<long>();
foreach (var tier in tiers)
{
    foreach (var combination in tier)
    {
        factors.Add(combination.Aggregate((total, next) => total * next));
    }
}

我们已经把这种情况得到了。

这种质数是可以任意相乘的，组合之后仍然满足题意。

接着处理下一种情况

他们和质因数只有235的数相乘，依旧满足题意。

factors.Sort();

foreach (var factor in factors)
{
    foreach (var smooth in smooths)
    {
        if ((double)smooth * factor <= max)
        {
            sum += smooth * factor;
        }
        else
        {
            break;
        }
    }
}

sum += smooths.Sum();

排序的好处就是可以提前终止。上面代码的最后一句是把基本case包含上了。

如果质因数只有2，3，5，减一之后是1，2，4，那么这个n是满足条件的。

smooths里面包含1（ijk都是0的情况），所以循环也包含了前面提到的这种情况。

顺便说一下，只有这些质数作为质因数的数也是满足题意的。

最后，返回结果。

1	return (uint)sum;

sum以开始定义为long，uint截断后就是模2^32的结果。

One Percenter

发表于 2018-05-10 分类于杂 Disqus：

在欧拉项目完成了113个题目，看了下历史，16，17，18年，一共完成了12个。。。要加油了，每年12题的目标的呢
不过好消息是，发现获得了一个成就，One Percerter，也就是所有玩这个的人里面排进了前1%

新的起点，新的动力！

最起码不能再掉出前百分之一 ::>_<::

欧拉项目 | 429题 | 元因数的平方和

发表于 2018-05-10 分类于数学，欧拉项目 Disqus：

429题链接

整数n的元因数（unitary divisor）d满足gcd(d, n/d) = 1
4! = 24的元因数有1，3，8，24
它们的平方和是1^2 + 3^2 + 8^2 + 24^2 = 650

S(n)表示n的元因数的平方和
求S(100 000 000!)模1 000 000 009

将n分解质因数，可以写作
n = p1^k1 + p2^k2 + ... + pm^km
那么元因数就是完全占有某些质数，比如p1^k1和p2^k2 + ... + pm^km就是一对元因数
有多少元因数呢？2^m个，一个很典型的排列组合问题。
我们下面看看平方和有什么规律。先看几个简单的。
将pi^ki记作qi
有两个质因数，那么平方和

1 2	S = 1 + q1^2 + q2^2 + (q1q2)^2 = (1+q1^2)(1+q2^2)

有三个质因数呢？

1 2	S = 1 + q1^2 + q2^2 + q3^2 + (q1q2)^2 + (q1q3)^2 + (q2q3)^2 + (q1q2q3)^2 = (1+q1^2)(1+q2^2)(1+q3^2)

我在纸上推理了四个质因数的情况，进而猜测到S = (1+q1^2)(1+q2^2)...(1+qm^2)
更通用的公式在这里

假定我们已经有了质因数的分解，那么很容易就能求出平方和了。

long m = 1000000009;
long ret = 1;
foreach (var pair in factors)
{
    var p = pair.Key;
    var k = pair.Value * 2;
    long r = 1;
    for (int i = 0; i < k; i++)
    {
        r *= p;
        r %= m;
    }
    ret *= (r + 1);
    ret %= m;
}

return ret;

现在的问题是如何得到100 000 000!的质因数，我们需要知道有多少个2，多少个3，多少个5等等。
Legendre’s formula给出了n!分解质因数的方法，这个方法简单易懂。
得到小于100 000 000的所有质数（筛选法），然后循环计算指数就可以了。

var factors = new Dictionary<long, int>();
foreach (var p in primes)
{
    if (p != 0)
    {
        int i = 1;
        int power = 0;
        while (true)
        {
            var cur = n / (long)Math.Pow(p, i);
            if (cur == 0)
            {
                break;
            }

            power += (int)cur;
            i++;
        }

        factors[p] = power;
    }
}