Project Euler - Problem 512

Sums of totients of powers

$\varphi(n)$ をオイラーの $\varphi$ 関数とする．
関数 $f(n)$ を $f(n) = ( \sum_{i = 1}^{n} \varphi(n^{i}) ) \mod (n + 1)$ によって定める．
さらに関数 $g(n)$ を $g(n) = \sum_{i = 1}^{n} f(i)$ によって定める．
$g(100) = 2007$ である．

$g(5 \times 10^{8})$ はいくつか？

2⁹ 問目ということでパワーのある問題

オイラーの $\varphi$ 関数ですが 1 から n までの自然数のうち n と互いに素なものの個数を $\varphi(n)$ で表します．オイラーのトーシェント関数と読みます．トーシェント

12 であれば 1, 5, 7, 11 と互いに素なので $\varphi(12) = 4$ です．

寄付くれとか言ってクソデカい広告を掲示している Wikipedia さんにも書いてありますが

$\varphi(p^{e}) = p^{e} - p^{e - 1} = p^{e - 1}(p - 1) = p^{e}(1 - 1/p) \ \ (p : prime)$

$\varphi(mn) = \varphi(m) \varphi(n)$

なんかが成立します．n の素因数分解が

$n = \prod_{k = 1}^{d} {p_{k}}^{e_{k}} = {p_{1}}^{e_{1}} {p_{2}}^{e_{2}} \cdots {p_{d}}^{e_{d}}$

であるとすると，

$\displaystyle \varphi(n) = \prod_{k = 1}^{d} \varphi({p_{k}}^{e_{k}}) = \prod_{k = 1}^{d} {p_{k}}^{e_{k}} (1 - \frac{1}{ p_{k} }) = n \prod_{k = 1}^{d} (1 - \frac{1}{ p_{k} })$

となり，指数 e に関係なく素因数 p のみに依存することがわかります．

上記性質より(ふつうに計算してもわかりますが)，

$\displaystyle \varphi(n^{i}) = n^{i} \prod_{k = 1}^{d} (1 - \frac{1}{ p_{k} })$

ですね．f はこれの 1 ~ n までの足し合わせなので，

$\displaystyle f(n) = \sum_{i = 1}^{n} \varphi(n^{i}) = \sum_{i = 1}^{n} n^{i} \prod_{k = 1}^{d} (1 - \frac{1}{ p_{k} }) = (1 + n + \cdots + n^{n - 1}) \cdot n \prod_{k = 1}^{d} (1 - \frac{1}{ p_{k} }) = (1 + n + \cdots + n^{n - 1}) \varphi(n)$

といい感じになりました．mod はあとでちゃんととります．

等比級数の部分をまとめるか否かうんうん悩みましたが

$n^{2} \equiv n^{2} + 2(n + 1) = (n + 1)^{2} + 1 \equiv 1 \ \ \mod (n + 1)$

$n^{2k - 1} \equiv n$
$n^{2k} \equiv 1$

に気付いたらあとはすんなりいきました．

n が偶数の場合は

\begin{align} f(2m) &= \underbrace{(1 + 2m + 1 + 2m + \cdots + 1 + 2m)}_{2m} \varphi(2m) \\\\ &= \underbrace{(2m + 1)}_{\equiv \ 0} \cdot m \cdot \varphi(2m) \\\\ &\equiv 0 \ \ \mod (2m + 1) \end{align}

n が奇数の場合は

\begin{align} f(2m - 1) &= \underbrace{ \{1 + (2m - 1) + 1 + (2m - 1) + \cdots + 1 + (2m - 1) + 1\}}_{2m - 1}\varphi(2m - 1) \\\\ &= \{ \underbrace{2m}_{\equiv \ 0} (m - 1) + 1 \} \varphi(2m - 1) \\\\ &\equiv \varphi(2m - 1) \ \ \mod 2m \end{align}

よって n が奇数の場合のみ考慮すればよいことがわかる．

あとはエラトステネスの篩っぽい方法(コード参照)で $f(n) = \varphi(n)$ をバーっと求めてだーっと足して $g(n)$ を求めればいい．

#define _CRT_SECURE_NO_WARNINGS
#include <iostream>
#include <cmath>
#include <string>
#include <vector>
#include <algorithm>
#include <queue>
#include <map>
#include <functional>
#include <set>
#include <numeric>
#include <stack>
#include <utility>
#include <time.h>
//#include "util.h"

using namespace std;
typedef unsigned uint;
typedef long long ll;
typedef unsigned long long ull;



//呪文
template <typename _KTy, typename _Ty> ostream& operator << (ostream& ostr, const pair<_KTy, _Ty>& m) { cout << "{" << m.first << ", " << m.second << "}"; return ostr; }
template <typename _KTy, typename _Ty> ostream& operator << (ostream& ostr, const map<_KTy, _Ty>& m) { if (m.empty()) { cout << "{ }"; return ostr;    } cout << "{" << *m.begin(); for (auto itr = ++m.begin(); itr != m.end(); itr++) { cout << ", " << *itr; } cout << "}"; return ostr; }
template <typename _Ty> ostream& operator << (ostream& ostr, const vector<_Ty>& v) { if (v.empty()) { cout << "{ }"; return ostr; } cout << "{" << v.front(); for (auto itr = ++v.begin(); itr != v.end(); itr++) { cout << ", " << *itr; }    cout << "}"; return ostr; }
template <typename _Ty> ostream& operator << (ostream& ostr, const set<_Ty>& s) { if (s.empty()) { cout << "{ }"; return ostr; } cout << "{" << *(s.begin()); for (auto itr = ++s.begin(); itr != s.end(); itr++) { cout << ", " << *itr; }    cout << "}"; return ostr; }
#define PI 3.14159265358979323846
#define EPS 1e-6
#define MIN(a,b) ((a)<(b)?(a):(b))
#define MAX(a,b) ((a)>(b)?(a):(b))
#define all(x) (x).begin(), (x).end()



#define N 500000000

//エラトステネスの篩用ブール型配列
bool p[N + 1];
//トーシェント関数の配列(n 奇数の場合のみ必要なのでメモリは半分)
int tot[N / 2];

int PE0512() 
{
    //N 以下の素数を求める
    for (int i = 2; i <= N; i++) p[i] = true;
    for (int i = 2; i * i <= N; i++)
        if (p[i])
            for (int j = i * i; j <= N; j += i)
                p[j] = false;

    //トーシェント関数を 1 で初期化
    for (int i = 0; i < N / 2; i++) tot[i] = 1;

    //奇素数 i に対して
    for (ll i = 3; i <= N; i += 2) {
        if (p[i]) {

            //奇素数 i の奇数倍 j に対して
            for (ll j = i; j <= N; j += 2 * i) {
                // φ(i) = i - 1 を利用
                // 必要なのは奇数項目のみなので，メモリ節約のため添字を圧縮(j / 2 の部分)
                tot[j / 2] *= i - 1;
            }
            
            //奇素数 i の累乗 j に対して
            for (ll j = i * i; j <= N; j *= i) {
                // j の奇数倍に対して
                for (ll k = j; k <= N; k += 2 * j) {
                    // φ(i^e) = i^(e - 1) * (i - 1) を利用する．(配列に順次 i をかけていく)
                    tot[k / 2] *= i;
                }
            }

        }
    }

    ll ans = 0;
    for (ll i = 0; i < N / 2; i++) ans += tot[i];

    cout << ans << endl;

    return 0;
}

int main()
{
    clock_t start, end;
    start = clock();

    PE0512();

    end = clock();
    printf("%d msec.\n", end - start);

    return 0;
}

環境にもよりますがおれの VS2015 なんかだと 2³¹ B = 2 GB を超えるとメモリが死ぬっぽくて，偶数番目が 0 で不要なことと $\varphi(n) < n$ より tot[n] が int に収まることなどを利用して辛うじてプログラムが動いた．ジャスト 60 秒だった．

だいぶカツカツだったのであとで他の解説も読んでみよう……

絶滅

どうでもいい

Project Euler - Problem 512

Sums of totients of powers