**Write-up** [Project Euler](https://projecteuler.net)
---

**Problem 169: Sums of Powers of Two** - 99,99999%
---
* Đề [Problem 169](https://projecteuler.net/problem=169)


* Sau khi chạy trâu các số từ 1 -> 20 và thử submit lên [oeis.org](https://oeis.org/) thì được kết quả là dãy [Stern's Diatomic Series](https://mathworld.wolfram.com/SternsDiatomicSeries.html). Do đó có công thức truy hồi số cách biểu diễn của số n là:
* $a_{n}$ = $\begin{cases} a_{n / 2} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ n \ chẵn \\ a_{(n + 1) / 2} + a_{(n - 1) / 2} \ \ \ \ \ n \ lẻ \end{cases}$
* Tuy nhiên vì dãy Stern's Diatomic Series bắt đầu đếm từ 0, còn dãy ta cần tính bắt đầu từ 1 nên khi tính, để ra kết quả chính xác thì ta tính $n+1$ sẽ được kết quả của số $n$
* Nhìn vào có thể hoàn toàn cài đặt bằng đệ quy, tuy nhiên số cần tính là 10^25^, vượt quá kiểu long long trong C++ nhưng vẫn có thể tính toán trong Python 1 cách bình thường

num = {0: 1, 1: 1}

def cal(n):
    if n in num:
        return num[n]
    if n % 2 == 1:
        num[n] = cal((n + 1) // 2) + cal((n - 1) // 2)
        return num[n]
    else:
        num[n] = cal(n // 2)
        return num[n]


print(cal(10 ** 25 + 1))

=> 178653872807

**Problem 185: Number Mind** - 100%
---
* Đề: [Problem 185](https://projecteuler.net/problem=185)
* Đề chỉ yêu cầu đơn giản là tìm ra số n có 16 số thỏa các điều kiện các chữ số đúng với các số cho sẵn ban đầu
* Vì số có 16 chữ số là 10^16, nếu sinh 10^16 số và duyệt thì đợi rất lâu, vậy nên ý tưởng duyệt trâu từng điều kiện, giả sử ở số s hiện tại k số đúng là s~1~, s~2~,.. s~k~ thì ở các số tiếp theo sẽ cài đặt sao cho n - k số còn lại đều không đúng
* Thử chạy trâu khi chưa xác định được độ phức tạp thì vô tình có kết quả trong khoảng 10s
* Hoàn toàn cài theo yêu cầu của đề, k số đúng thì n - k số còn lại không đúng, nếu những số đã đúng trùng với số đang xét thì số đang xét cũng đúng - hoàn toàn quay lui không có gì thêm

#include<bits/stdc++.h>
using namespace std;

vector<pair<int, string>> num;
string b;
int vis[20][10], a;
int patch[20];
int n_patch[20];
bool ret = false;

bool check()
{
    for(int i = 0; i < 16; i++)
    {
        if(!patch[i]) continue;
        if(vis[i][n_patch[i]]) return false;
    }
    return true;
}

bool f_ans()
{
    for(int i = 0; i < 16; i++)
    {
        if(patch[i] && vis[i][n_patch[i]]) return false;
        bool tmp = false;
        for(int j = 0 + (i == 0); j < 10; j++)
            tmp |= !(vis[i][j]);
        if(!tmp) return false;
    }
    for(int i = 0; i < 16; i++)
    {
        if(patch[i]) cout << n_patch[i];
        else
            for(int j = 0 + (i == 0); j < 10; j++)
                if(!vis[i][j])
                {
                    cout << j;
                    break;
                }
    }
    return true;
}

void solve(int k)
{
    if (ret) return;
    if(k == num.size())
    {
        ret = f_ans();
        return;
    }
    string tmp = num[k].second;
    int u = num[k].first;
    if(u == 0)
    {
        for(int i = 0; i < tmp.size(); i++)
            vis[i][tmp[i] - '0']++;
        solve(k + 1);
        return;
    }
    if(u == 1)
    {
        for(int i = 0; i < tmp.size(); i++)
        {
            if(ret) return;
            if(vis[i][tmp[i] - '0']) continue;
            if(patch[i] && n_patch[i] != (tmp[i] - '0')) continue;
            patch[i]++;
            n_patch[i] = tmp[i] - '0';

            bool atmp = true;
            for(int j = 0; j < tmp.size(); j++)
            {
                if(i == j) continue;
                if(patch[j] && n_patch[j] == (tmp[j] - '0'))
                {
                    atmp = false;
                    break;
                }
            }

            if(atmp)
            {
                for(int j = 0; j < tmp.size(); j++)
                {
                    if(j == i) continue;
                    vis[j][tmp[j] - '0']++;
                }

                solve(k + 1);

                for(int j = 0; j < tmp.size(); j++)
                {
                    if(j == i) continue;
                    vis[j][tmp[j] - '0']--;
                }
            }
            patch[i]--;
        }
        return;
    }
    if(u == 2)
    {
        for(int i = 0; i < tmp.size(); i++)
            for(int x = i + 1; x < tmp.size(); x++)
            {
                if(ret) return;

                if(vis[i][tmp[i] - '0'] || vis[x][tmp[x] - '0']) continue;
                if((patch[i] && n_patch[i] != (tmp[i] - '0')) || (patch[x] && n_patch[x] != (tmp[x] - '0'))) continue;
                patch[i]++;
                patch[x]++;
                n_patch[i] = tmp[i] - '0';
                n_patch[x] = tmp[x] - '0';

                bool atmp = true;
                for(int j = 0; j < tmp.size(); j++)
                {
                    if(i == j) continue;
                    if(x == j) continue;
                    if(patch[j] && n_patch[j] == (tmp[j] - '0'))
                    {
                        atmp = false;
                        break;
                    }
                }

                if(atmp)
                {
                    for(int j = 0; j < tmp.size(); j++)
                    {
                        if(j == i) continue;
                        if(j == x) continue;
                        vis[j][tmp[j] - '0']++;
                    }

                    solve(k + 1);

                    for(int j = 0; j < tmp.size(); j++)
                    {
                        if(j == i) continue;
                        if(j == x) continue;
                        vis[j][tmp[j] - '0']--;
                    }
                }
                patch[i]--;
                patch[x]--;
            }
        return;
    }

    if(u == 3)
    {
        for(int i = 0; i < tmp.size(); i++)
            for(int x = i + 1; x < tmp.size(); x++)
                for(int y = x + 1; y < tmp.size(); y++)
                {
                    if(ret) return;

                    if(vis[i][tmp[i] - '0'] || vis[x][tmp[x] - '0'] || vis[y][tmp[y] - '0']) continue;
                    if((patch[i] && n_patch[i] != (tmp[i] - '0')) || (patch[x] && n_patch[x] != (tmp[x] - '0')) || (patch[y] && n_patch[y] != (tmp[y] - '0'))) continue;
                    patch[i]++;
                    patch[x]++;
                    patch[y]++;
                    n_patch[i] = tmp[i] - '0';
                    n_patch[x] = tmp[x] - '0';
                    n_patch[y] = tmp[y] - '0';

                    bool atmp = true;
                    for(int j = 0; j < tmp.size(); j++)
                    {
                        if(i == j) continue;
                        if(x == j) continue;
                        if(y == j) continue;
                        if(patch[j] && n_patch[j] == (tmp[j] - '0'))
                        {
                            atmp = false;
                            break;
                        }
                    }

                    if(atmp)
                    {
                        for(int j = 0; j < tmp.size(); j++)
                        {
                            if(j == i) continue;
                            if(j == x) continue;
                            if(j == y) continue;
                            vis[j][tmp[j] - '0']++;
                        }

                        solve(k + 1);

                        for(int j = 0; j < tmp.size(); j++)
                        {
                            if(j == i) continue;
                            if(j == x) continue;
                            if(j == y) continue;
                            vis[j][tmp[j] - '0']--;
                        }
                    }
                    patch[i]--;
                    patch[x]--;
                    patch[y]--;
                }
            return;
    }
}

int main()
{
    ios_base::sync_with_stdio(false); cin.tie(0);

    freopen("text.txt", "r", stdin);

    while(cin >> a >> b)
        num.push_back({a, b});
    sort(num.begin(), num.end());
    solve(0);
}


**Problem 231 : Prime factorisation of binomial coefficients** - 100%
---
* Đề: [Problem 231](https://projecteuler.net/problem=231)


- Phân tích tổ hợp chập nCk thành tích các thừa số nguyên tố có dạng 2^a^x3^b^x5^c^..., sau đó tính tổng của 2a+3b+5c+...
* Vì trong công thức tính nCk là các giai thừa nên nó sẽ là tích tất cả các số <= n, vậy nên toàn bộ ước nguyên tố lẫn lũy thừa của ước nguyên tố đều sẽ được tính
* Do đó có thể sử dụng phân tích n! thành tích nguyên tố để đếm trong độ phức tạp O(log~e~n), dùng sàng nguyên tố để prebuild các số nguyên tố trong khoảng từ 2 -> n sau đó tính toán để giảm độ phức tạp
* Hoàn toàn có thể chạy trong thời gian nhỏ vì độ phức tạp không quá lớn: O(nlog~e~n) với n = 2e7

#include<bits/stdc++.h>
using namespace std;
#define N (long long)2e7
#define M (long long)15e6

bool prime[N + 5];
vector <long long> num;
void build()
{
    for(int i = 2; i * i <= N; i++)
    {
        if(prime[i]) continue;
        for(int j = i * i; j <= N; j += i)
            prime[j] = true;
    }
    for(int i = 2; i <= N; i++)
        if(!prime[i])
            num.push_back(i);
}

long long cal(long long k)
{
    long long ret = 0;
    long long c_p, tmp;
    for(int i = 0; i < num.size(); i++)
    {
        tmp = num[i];
        while(tmp <= k)
        {
            ret += num[i] * (k / tmp);
            c_p++;
            tmp *= num[i];
        }
    }
    return ret;
}

long long n, m;
int main()
{
    ios_base::sync_with_stdio(false); cin.tie(0);
    build();
    cout << cal(N) - (cal(M) + cal(N - M));
}


**Problem 237: Tours on a 4xN Playing Board** - 100%
---
* Đề: [Problem 237](https://projecteuler.net/problem=237)

* Biết được bài toán có dạng quy hoạch động nên tìm công thức quy hoạch động bằng cách:
* Chạy trâu quay lui tìm đường đi trên bảng, sinh 10 giá trị đầu của T(n)

#include <bits/stdc++.h>
using namespace std;
#define N 20

int num[5][N + 5];
int ans;
bool vis[5][N + 5];
int x[5]{0, 0, 1, -1};
int y[5]{1, -1, 0, 0};

void atmp(int n, int i, int j, int s)
{
    if(s == n * 4 - 1 && i == 4 && j == 1)
    {
        ans++;
        return;
    }
    for(int k = 0; k < 4; k++)
    {
        int u = i + x[k];
        int v = j + y[k];
        if(u < 1 || u > 4) continue;
        if(v < 1 || v > n) continue;
        if(vis[u][v]) continue;
        vis[u][v] = true;
        atmp(n, u, v, s + 1);
        vis[u][v] = false;
    }
}

int main()
{
    vis[1][1] = true;
    for(int i = 1; i <= 10; i++)
    {
        ans = 0;
        atmp(i, 1, 1, 0);
        cout << ans << " ";
    }
    return 0;
}

* => Sau đó đưa dãy số vào [oeis.org](https://oeis.org/ ) tìm được công thức truy hồi:
* $a_{n}$ = $a_{n-4}$ - $2a_{n-3}$ + $2a_{n-2}$ + $2a_{n-1}$
* Tuy nhiên n lớn (10^12^) nên quy đổi bài toán sang nhân ma trận để tính lũy thừa trong độ phức tạp log, giảm thời gian chờ kết quả:
* Ma trận ban đầu $\begin{bmatrix} a_{1} \ \ \ a_{2} \ \ \ a_{3} \ \ \ a_{4}\end{bmatrix}$ = $\begin{bmatrix} 1 \ \ \ 1 \ \ \ 4 \ \ \ 8 \end{bmatrix}$
* Ma trận cơ sở: $\begin{bmatrix} 0 \ \ \ 0 \ \ \ 0 \ \ \ \ 1 \\ 1 \ \ \ 0 \ \ \ 0 \ {-2} \\ 0 \ \ \ 1 \ \ \ 0 \ \ \ 2 \\ 0 \ \ \ 0 \ \ \ 1 \ \ \ 2 \end{bmatrix}$
* Phép nhân ma trận:
$\begin{bmatrix} a_{n - 3} \ \ \ a_{n - 2} \ \ \ a_{n - 1} \ \ \ a_{n}\end{bmatrix}$ * $\begin{bmatrix} 0 \ \ \ 0 \ \ \ 0 \ \ \ \ 1 \\ 1 \ \ \ 0 \ \ \ 0 \ {-2} \\ 0 \ \ \ 1 \ \ \ 0 \ \ \ 2 \\ 0 \ \ \ 0 \ \ \ 1 \ \ \ 2 \end{bmatrix}$ = $\begin{bmatrix} a_{n - 2} \ \ \ a_{n - 1} \ \ \ a_{n} \ \ \ a_{n + 1}\end{bmatrix}$

* => $\begin{bmatrix} a_{n - 3} \ \ \ a_{n - 2} \ \ \ a_{n - 1} \ \ \ a_{n}\end{bmatrix}$ * $\begin{bmatrix} 0 \ \ \ 0 \ \ \ 0 \ \ \ \ 1 \\ 1 \ \ \ 0 \ \ \ 0 \ {-2} \\ 0 \ \ \ 1 \ \ \ 0 \ \ \ 2 \\ 0 \ \ \ 0 \ \ \ 1 \ \ \ 2 \end{bmatrix} ^k$ = $\begin{bmatrix} a_{n - 3 + k} \ \ \ a_{n - 2 + k} \ \ \ a_{n - 1 + k} \ \ \ a_{n + k}\end{bmatrix}$

Kết quả cuối cùng sẽ là lấy ma trận ban đầu nhân cột 4 của ma trận $\begin{bmatrix} 0 \ \ \ 0 \ \ \ 0 \ \ \ \ 1 \\ 1 \ \ \ 0 \ \ \ 0 \ {-2} \\ 0 \ \ \ 1 \ \ \ 0 \ \ \ 2 \\ 0 \ \ \ 0 \ \ \ 1 \ \ \ 2 \end{bmatrix} ^{10^{12}-4}$
* Độ phức tạp nhỏ: nhân ma trận O(4x4x4) x lũy thừa tốn O(log~2~(N)) (N = 10^12^) = tổng thể độ phức tạp sẽ là O(12xlog~2~(N))

#include <bits/stdc++.h>
using namespace std;
#define N (long long)1e12
#define M (long long)1e8

long long a[5]{0, 1, 1, 4, 8};

struct mat{
    long long num[5][5]{{},{},{},{},{}};
};

mat Base;

mat mul(mat a, mat b)
{
    mat ret;
    for(int i = 1; i <= 4; i++)
        for(int j = 1; j <= 4; j++)
            for(int k = 1; k <= 4; k++)
            {
                ret.num[i][j] += a.num[i][k] * b.num[k][j];
                ret.num[i][j] %= M;
            }
    return ret;
}

mat m_pow(mat a, long long b)
{
    if(b == 1) return a;
    mat c = m_pow(a, b / 2);
    mat k = mul(c, c);
    if ((b % 2))
        k = mul(k, a);
    return k;
}

long long ans()
{
    long long ret = 0;
    for(int i = 1; i <= 4; i++)
    {
        ret += a[i] * Base.num[i][4];
        ret %= M;
    }
    return ret % M;
}

int main()
{
    Base.num[1][4] = Base.num[2][1] = Base.num[3][2] = Base.num[4][3] = 1;
    Base.num[3][4] = Base.num[4][4] = 2;
    Base.num[2][4] = -2;
    Base = m_pow(Base, N - 4);
    cout << ans();
}


**Problem 247: Squares Under a Hyperbola**
---
* Đề: [Prblem 247](https://projecteuler.net/problem=247)