筛法，透过肌肤的温度确实传达过来

一维积性函数前缀和

前置知识

powerful number

定义：\(\text{powerful number}\) 是没有 \(1\) 次质因子的数，也简称为 \(\text{PN}\)。

结论：\(n\) 以内的 \(\text{PN}\) 个数为 \(O(\sqrt n)\)。

证明：

所有 \(\text{PN}\) 必然可表示成 \(a^2b^3\) 的形式。

那么：

\[\sum_{a=1}^{\sqrt n}(\frac{n}{a^2})^{\frac{1}{3}}=\sqrt n \]

构造

构造 \(G(p^k)=F(p)^k,G\times H=F\)，可以注意到 \(G(P)\) 为完全积性函数。

那么：

\[\begin{aligned} F(p)&=G(1)H(p)+H(1)G(p)\\ &=G(p)+H(p) \end{aligned}\]

所以，\(H(p)=0\)。

\[\begin{aligned} \sum_{i=1}^n F(i)&=\sum_{i=1}^n\sum_{j|i}G(j)H(\frac{i}{j})\\ &=\sum_{i=1}^nH(i)\sum_{j=1}^{\frac{n}{i}}G(j)\\ &=\sum_{i=1}^n[i\in PN]H(i)SG_{\frac{n}{i}}\\ \end{aligned}\]

其中:

\[\begin{aligned} H(p)&=0\\ H(p^2)&=F(p^2)-F(p)^2\\ H(p^3)&=F(p^3)-F(p)F(p^2)\\ H(p^k)&=F(p^k)-F(p)F(p^{k-1}) \end{aligned}\]

所以，我们只需要求出 \(G\) 的块筛即可。

块筛

考虑两种方式。

法一

若 \(G\) 可杜教筛。

可以直接使用杜教筛。

此时，这个筛法就是常称的 \(\text{PN}\) 筛。

法二

若 \(F(p)\) 为一个低次多项式。

考虑 \(S(n,a)\)。

\[\begin{aligned} S(n,a)&=\sum_{i=1}^n[i\in P ~or~ min_i > P_a] G(i)\\ &=S(n,a-1)-g(p_a)(S(\frac{n}{P_a},a-1)-S(P_a-1,a-1))\\ S(n,a-1)&=S(n,a)+g(p_a)(S(\frac{n}{P_a},a-1)-S(P_a-1,a-1))\\ \end{aligned}\]

要求：\(G(x)\) 为完全积性函数。

若 \(\sum_{i=1}^n [i\in P]G(x)\) 很好求出，那么直接用第二个递推式即可。

否则，将多项式 \(G(x)\) 的每一项分别提出来，跑第一个递推式，合并后再用第二个递推式即可。

本质是 \(\text{min-25}\) 筛的前半部分。

完美利用了 \(G(x)\) 的性质。

个人更加推荐法二，这也可以看出这个筛法极强的拓展性。

速度

考虑 \(\text{PN}\) 的过程复杂度为 \(O(\sqrt n)\)。

块筛的过程要么是 \(O(n^{\frac{2}{3}})\) 或者 \(O(\frac{n^{\frac{3}{4}}}{\log n})\)（很可能不准）。

实际效率上比普通 \(\text{min-25}\) 筛快一点（板子题）。

\(10^{10}\)：\(0.2s\)。

\(10^{11}\)：\(0.6s\)。

\(10^{12}\)：\(3.0s\)。

代码难度个人认为比较简单，与杜教筛难度近似，感觉可以作为 \(\text{min-25}\) 筛的上位替代。

例题

P4213 【模板】杜教筛

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

using i64 = long long;

const int N = 1e5 + 10;

i64 n, ans1, ans2, w[N];
i64 s1[N], s2[N], p1[N], p2[N];
int m, ct, tt, p[N], v[N], d1[N], d2[N];

inline void init_all(int m = N - 10) {
  for (int i = 2; i <= m; i++) {
    if (!v[i]) p[++ct] = i, p1[ct] = 1 + p1[ct - 1], p2[ct] = i + p2[ct - 1];
    for (int j = 1; j <= ct && p[j] * i <= m; j++) {
      v[p[j] * i] = 1; if (i % p[j] == 0) break;
    }
  }
}
inline int get(i64 x) { return (x <= m ? d1[x] : d2[n / x]); }
inline void init() {
  m = sqrt(n) + 100, ct = tt = ans1 = ans2 = 0;
  for (i64 l = 1, r; l <= n; l = r + 1) {
    r = (n / (n / l)), w[++tt] = n / l;
    w[tt] <= m ? d1[w[tt]] = tt : d2[r] = tt;
    s1[tt] = w[tt], s2[tt] = (w[tt] + 1) * w[tt] / 2;
    s1[tt]--, s2[tt]--;
  }
  while (1ll * p[ct + 1] * p[ct + 1] <= n) ct++;
  for (int i = 1; i <= ct; i++) {
    const i64 k1 = p[i];
    const i64 k2 = k1 * k1;
    int lim = 0;
    while (lim < tt && k2 <= w[lim + 1]) lim++;
    for (int j = 1; j <= lim; j++) {
      if (k2 > w[j]) break;
      int id = get(w[j] / k1);
      s1[j] = (s1[j] - (s1[id] - p1[i - 1]));
      s2[j] = (s2[j] - k1 * (s2[id] - p2[i - 1]));
    }
  }
  for (int i = 1; i <= tt; i++) s1[i] = -s1[i], s2[i] = s2[i] + s1[i];
  for (int i = 1; i <= tt; i++) p1[i] = -p1[i], p2[i] = p2[i] + p1[i];
  for (int i = ct; i >= 1; i--) {
    const i64 k1 = p[i];
    const i64 k2 = k1 * k1;
    int lim = 0;
    while (lim < tt && k2 <= w[lim + 1]) lim++;
    for (int j = lim; j >= 1; j--) {
      int id = get(w[j] / k1);
      s1[j] = (s1[j] + -1 * (s1[id] - p1[i - 1]));
      s2[j] = (s2[j] + (k1 - 1) * (s2[id] - p2[i - 1]));
    }
  }
  for (int i = 1; i <= tt; i++) s1[i]++, s2[i]++;
  for (int i = 1; i <= tt; i++) p2[i] = p2[i] - p1[i], p1[i] = -p1[i];
}
inline void dfs(int now, i64 cur, i64 h1, i64 h2, i64&ans1, i64&ans2) {
  ans1 = (ans1 + h1 * s1[get(n / cur)]);
  ans2 = (ans2 + h2 * s2[get(n / cur)]);
  for (int i = now; i <= ct; i++) {
    const i64 k = p[i];
    if (k * k > n / cur) break;
    for (i64 x = k * k, y = 1, z = 2; x * cur <= n; x = x * k, y = y * k, z++) {
      dfs(i + 1, cur * x, h1 * (z == 2 ? -1 : 0), h2 * (k - 1) * y, ans1, ans2);
    }
  }
}

int main() {
  int t;
  cin >> t, init_all();
  while (t--) {
    cin >> n, init();
    dfs(1, 1, 1, 1, ans1, ans2);
    cout << ans2 << " ";
    cout << ans1 << "\n";
    for (int i = 1; i <= tt; i++) s1[i] = s2[i] = d1[i] = d2[i] = w[i] = 0;
  }
  return 0;
}

P5325 【模板】Min_25 筛

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

using i64 = long long;

const int N = 1e6 + 10;
const int mod = 1e9 + 7;

i64 n, ans, w[N];
i64 s1[N], s2[N], p1[N], p2[N];
int m, ct, tt, p[N], v[N], d1[N], d2[N];

inline int get(i64 x) { return (x <= m ? d1[x] : d2[n / x]); }
inline void init() {
  m = sqrt(n) + 100;
  for (i64 l = 1, r; l <= n; l = r + 1) {
    r = (n / (n / l)), w[++tt] = n / l;
    w[tt] <= m ? d1[w[tt]] = tt : d2[r] = tt;
    s1[tt] = (__int128)(w[tt] + 1) * w[tt] * 500000004 % mod;
    s2[tt] = (__int128)(2 * w[tt] + 1) * s1[tt] * 333333336 % mod;
    s1[tt]--, s2[tt]--;
  }
  for (int i = 2; i <= m; i++) {
    if (!v[i]) p[++ct] = i, p1[ct] = (i + p1[ct - 1]) % mod, p2[ct] = (1ll * i * i + p2[ct - 1]) % mod;
    for (int j = 1; j <= ct && p[j] * i <= m; j++) {
      v[p[j] * i] = 1;
      if (i % p[j] == 0) break;
    }
  }
  for (int i = 1; i <= ct; i++) {
    const i64 k1 = p[i];
    const i64 k2 = k1 * k1;
    int lim = 0;
    while (lim < tt && k2 <= w[lim + 1]) lim++;
    for (int j = 1; j <= lim; j++) {
      if (k2 > w[j]) break;
      int id = get(w[j] / k1);
      s1[j] = (s1[j] - k1 % mod * (s1[id] - p1[i - 1])) % mod;
      s2[j] = (s2[j] - k2 % mod * (s2[id] - p2[i - 1])) % mod;
    }
  }
  for (int i = 1; i <= tt; i++) s1[i] = (s2[i] - s1[i]) % mod;
  for (int i = 1; i <= ct; i++) p1[i] = (p2[i] - p1[i]) % mod;
  for (int i = ct; i >= 1; i--) {
    const i64 k1 = p[i];
    const i64 k2 = k1 * k1;
    int lim = 0;
    while (lim < tt && k2 <= w[lim + 1]) lim++;
    for (int j = lim; j >= 1; j--) {
      int id = get(w[j] / k1);
      s1[j] = (s1[j] + k1 * (k1 - 1) % mod * (s1[id] - p1[i - 1])) % mod;
    }
  }
  for (int i = 1; i <= tt; i++) s1[i]++;
}
inline void dfs(int now, i64 cur, i64 h, i64&ans) {
  ans = (ans + h * s1[get(n / cur)]) % mod;
  for (int i = now; i <= ct; i++) {
    const i64 k = p[i];
    if (k * k > n / cur) break;
    for (i64 x = k * k, y = k; x * cur <= n; x = x * k, y = y * k) {
      dfs(i + 1, cur * x, (__int128)h * x * (y + k - 2) % mod, ans);
    }
  }
}

int main() {
  cin >> n, init();
  dfs(1, 1, 1, ans);
  cout << (ans + mod) % mod << "\n";
  return 0;
}

二维积性函数前缀和

前置知识

一维积性函数前缀和。

定义

我们称一个函数是二维积性函数当且仅当满足：

\[\forall \gcd(ab,cd)=1,f(ac,bd)=f(a,b)\times f(c,d) \]

并且有：\(f(1,1)=1\)。

相似的，我们也可以对 \(x,y\) 分解质因数，有：

\[f(x,y)=\prod_i f(p_i^{a_{i}},p_i^{b_{i}}) \]

求解

考虑构造：\(g(p^a,p^b)=f(p^a,1)\times f(1,p^b)\)。

不难发现 \(g\) 也是一个二维积性函数。

考虑狄利克雷卷积，\(g*h=f\)：

\[f(a,b)=\sum_{u|a,v|b}g(u,v)h(\frac{a}{u},\frac{b}{v}) \]

就可以得到狄利克雷除法，\(h=f/g\)：

\[h(a,b)=f(a,b)-\sum_{u|a,v|b,(u\not=1||v\not=1)}g(u,v)h(\frac{a}{u},\frac{b}{v}) \]

观察，\(h(1,p)\) 处的值。

\[\begin{align} h(1,p)=h(p,1)&=f(1,p)-g(1,p)h(1,1)\\ &=f(1,p)-f(1,p)\\ &=0 \end{align} \]

那么我们可以枚举 \(h\) 的点值。

有：

\[\begin{align} \sum\sum f(i,j)&=\sum\sum h(i,j)\sum^\frac{n}{i}\sum^\frac{m}{j}g(l,r)\\ &=\sum\sum h(i,j)(\sum^\frac{n}{i}f(l,1))(\sum^\frac{m}{j}f(1,r)) \end{align} \]

那么我们只需要求出 \(f\) 的块筛即可，大多数情况 \(n,m\) 比较小可以直接前缀和，\(n,m\) 大的时候可以用前面的筛法。

复杂度证明比较困难，是 \(O((nm)^{\frac{1}{2}})\) 的，这样就有了一个二维积性函数前缀和的通用解法。

更高维

若是想要拓展到更高维度的情况，也可以使用这个做法的到开根的复杂度。

但是常数极大。

三维就有了几百倍的常数了，一般只有 \(h\) 的性质更强的时候才会使用。

二维的时候可以看作线性。

实现细节

包括前面的一维积性函数前缀和，\(h\) 怎么求一直是一个问题。

在前面的二维积性函数前缀和中，我们使用了手动推式子得到了 \(h\) 的通项。

但其实有一个不要脑子的做法。

就是我们直接进行狄利克雷除法，因为只要统计质数的幂次处的点值，所以跑得是很快的，当然要数据范围比较小的时候才可以，\(10^8,10^9\) 就不可以了。

所以更多的时候，二维积性函数前缀和才可以用这个做法。

代码

这个做法有一点太强了。

根本不需要推式子，直接秒掉。

P4619 [SDOI2018] 旧试题

与正解复杂度相同，跑的飞快，这就是三维积性函数前缀和，并且 \(h\) 的性质更强的情况，时间复杂度：\(O(n\sqrt n)\)。

/*
  ! 如果没有天赋，那就一直重复
  ! Created: 2024/08/13 19:29:18
*/
#include <bits/stdc++.h>
using namespace std;

#define x first
#define y second
// #define int long long
#define mp(x, y) make_pair(x, y)
#define eb(...) emplace_back(__VA_ARGS__)
#define fro(i, x, y) for (int i = (x); i <= (y); i++)
#define pre(i, x, y) for (int i = (x); i >= (y); i--)
inline void JYFILE19();

using i64 = long long;
using pii = pair<int, int>;

bool ST;
const int N = 1e5 + 10;
const int mod = 1e9 + 7;

int n, m, k, ct, ans;
int p[N];
int h[N][2][2][2];
int vs[N], pr[N], pd[N], dp[N];

inline void add(int&x, int y) { x = x + y; if (x >= mod) x -= mod; }
inline void del(int&x, int y) { x = x - y; if (x < 0) x += mod; }
inline void add(int&x, int y, int z) { x = (x + 1ll * y * z) % mod; }
inline void del(int&x, int y, int z) { x = (x - 1ll * y * z % mod); if (x < 0) x += mod; }
inline void init(int n) {
  dp[1] = 1;
  fro(i, 2, n) {
    if (!vs[i]) {
      pr[++ct] = i, dp[i] = 2, pd[i] = i;
      for (int j = i; j <= n / i; j = j * i) {
        dp[j * i] = dp[j] + 1;
      }
    }
    for (int j = 1; j <= ct && pr[j] * i <= n; j++) {
      vs[pr[j] * i] = 1;
      if (i % pr[j] == 0) {
        pd[pr[j] * i] = pd[i] * pr[j];
        break;
      }
      pd[pr[j] * i] = pr[j];
    }
  }
  fro(i, 2, n) if (vs[i]) dp[i] = dp[pd[i]] * dp[i / pd[i]];
  fro(i, 2, n) add(dp[i], dp[i - 1]);
  fro(i, 1, ct) {
    int r = pr[i];
    for (int o1 = 0; o1 <= 1; o1++)
    for (int o2 = 0; o2 <= 1; o2++)
    for (int o3 = 0; o3 <= 1; o3++) {
      int f = (o1 + o2 + o3 + 1);
      for (int q1 = o1, w1 = 0; w1 <= o1; q1--, w1++)
      for (int q2 = o2, w2 = 0; w2 <= o2; q2--, w2++)
      for (int q3 = o3, w3 = 0; w3 <= o3; q3--, w3++) {
        if (w1 == o1 && w2 == o2 && w3 == o3) break;
        del(f, (q1 + 1) * (q2 + 1) * (q3 + 1), h[i][w1][w2][w3]);
      }
      h[i][o1][o2][o3] = f;
    }
  }
}
inline void dfs(int d, int p1, int p2, int p3, int w) {
  ans = (ans + 1ll * dp[n / p1] * dp[m / p2] % mod * dp[k / p3] % mod * w) % mod;
  for (int i = d; i <= ct; i++) {
    i64 p = pr[i], flg = 0;
    if (p1 * p <= n && p2 * p <= m) dfs(i + 1, p1 * p, p2 * p, p3, 1ll * w * h[i][1][1][0] % mod), flg = 1;
    if (p1 * p <= n && p3 * p <= k) dfs(i + 1, p1 * p, p2, p3 * p, 1ll * w * h[i][1][0][1] % mod), flg = 1;
    if (p2 * p <= m && p3 * p <= k) dfs(i + 1, p1, p2 * p, p3 * p, 1ll * w * h[i][0][1][1] % mod), flg = 1;
    if (p1 * p <= n && p2 * p <= m && p3 * p <= k) dfs(i + 1, p1 * p, p2 * p, p3 * p, 1ll * w * h[i][1][1][1] % mod), flg = 1;
    if (flg == 0) break;
  }
}


signed main() {
  JYFILE19();
  int t;
  cin >> t, init(100000);
  while (t--) {
    cin >> n >> m >> k, ans = 0;
    dfs(1, 1, 1, 1, 1);
    cout << ans << "\n";
  }
  return 0;
}

bool ED;
inline void JYFILE19() {
  // freopen("", "r", stdin);
  // freopen("", "w", stdout);
  srand(random_device{}());
  ios::sync_with_stdio(0), cin.tie(0);
  double MIB = fabs((&ED - &ST) / 1048576.), LIM = 1024;
  cerr << "MEMORY: " << MIB << endl, assert(MIB <= LIM);
}

P8570 [JRKSJ R6] 牵连的世界

直接爆掉标算了，时间复杂度：\(O(n)\)，标算复杂度：\(O(n\log n\log\log n)\)。

/*
  ! 如果没有天赋，那就一直重复
  ! Created: 2024/08/13 21:13:39
*/
#include <bits/stdc++.h>
using namespace std;

#define x first
#define y second
// #define int long long
#define mp(x, y) make_pair(x, y)
#define eb(...) emplace_back(__VA_ARGS__)
#define fro(i, x, y) for (int i = (x); i <= (y); i++)
#define pre(i, x, y) for (int i = (x); i >= (y); i--)
inline void JYFILE19();

using i64 = long long;
using pii = pair<int, int>;

bool ST;
const int N = 3e6 + 10;
const int M = 3e5 + 10;
const int K = 3e3 + 10;
const int mod = 1e9 + 7;

int n, m, k, ct, ans;
int vs[N], pr[N], pd[N];
i64 sp[K][25][25];
i64 ps[M][2][2];
i64 dp[N];

template<typename T> inline void add(T&x, int y) { x = x + y; if (x >= mod) x -= mod; }
template<typename T> inline void del(T&x, int y) { x = x - y; if (x < 0) x += mod; }
template<typename T> inline void add(T&x, int y, int z) { x = (x + 1ll * y * z) % mod; }
template<typename T> inline void del(T&x, int y, int z) { x = (x - 1ll * y * z % mod); if (x < 0) x += mod; }
inline int f(i64 x, int y) {
  if (x == 1) return 1;
  return x / y * (y - 1) % mod;
}
inline void init(int n) {
  dp[1] = 1;
  fro(i, 2, n) {
    if (!vs[i]) {
      pr[++ct] = i, pd[i] = i, dp[i] = 2;
      for (int j = i; j <= n / i; j = j * i) {
        dp[j * i] = dp[j] + 1;
      }
      dp[i] = dp[i] * (i - 1) % mod;
      for (int j = i; j <= n / i; j = j * i) {
        dp[j * i] = dp[j * i] * (i - 1) * j % mod;
      }
    }
    for (int j = 1; j <= ct && pr[j] * i <= n; j++) {
      vs[i * pr[j]] = 1;
      if (i % pr[j] == 0) {
        pd[i * pr[j]] = pd[i] * pr[j];
        break;
      }
      pd[i * pr[j]] = pr[j];
    }
  }
  fro(i, 2, n) {
    if (pd[i] != i) {
       dp[i] = dp[i / pd[i]] * dp[pd[i]] % mod;
    }
  }
  fro(i, 2, n) add(dp[i], dp[i - 1]);
  fro(i, 1, ct) {
    i64 r = pr[i];
    if (r * r > n) {
      ps[i][0][0] = 1;
      ps[i][1][1] = (5 * r - r * r - 4) % mod;
      if (ps[i][1][1] < 0) ps[i][1][1] += mod;
    } else {
      for (i64 p1 = 1, o1 = 0; p1 <= n; p1 = p1 * r, o1++) {
      for (i64 p2 = 1, o2 = 0; p2 <= n; p2 = p2 * r, o2++) {
        int d = (o1 + o2 + 1) * f(p1 * p2, r) % mod;
        for (i64 q1 = p1, w1 = 0; w1 <= o1; q1 = q1 / r, w1++)
        for (i64 q2 = p2, w2 = 0; w2 <= o2; q2 = q2 / r, w2++) {
          if (w1 == o1 && w2 == o2) break;
          del(d, (o1 - w1 + 1) * (o2 - w2 + 1) % mod * f(q1, r) % mod * f(q2, r) % mod * sp[i][w1][w2] % mod);
        }
        sp[i][o1][o2] = d;
      }
      }
    }
  }
}
inline void dfs(int d, int p1, int p2, int w) {
  ans = (ans + dp[n / p1] * dp[m / p2] % mod * w) % mod;
  fro(i, d, ct) {
    i64 p = pr[i];
    if (p1 * p > n || p2 * p > m) break;
    if (p * p > max(n, m)) {
      dfs(i + 1, p1 * p, p2 * p, w * ps[i][1][1] % mod);
    } else {
      for (i64 q1 = p1 * p, o1 = 1; q1 <= n; q1 = q1 * p, o1++)
      for (i64 q2 = p2 * p, o2 = 1; q2 <= m; q2 = q2 * p, o2++)
        dfs(i + 1, q1, q2, w * sp[i][o1][o2] % mod);
    }
  }
}

signed main() {
  JYFILE19();
  cin >> n >> m, init(max(n, m)), dfs(1, 1, 1, 1), cout << ans << "\n";
  return 0;
}

bool ED;
inline void JYFILE19() {
  // freopen("", "r", stdin);
  // freopen("", "w", stdout);
  srand(random_device{}());
  ios::sync_with_stdio(0), cin.tie(0);
  double MIB = fabs((&ED - &ST) / 1048576.), LIM = 256;
  cerr << "MEMORY: " << MIB << endl, assert(MIB <= LIM);
}

#885. 【UR #27】红场阅兵

一切的开始，根号只有 \(70\) 分，但后续优化太难了，根本不会，一般根号也够了，时间复杂度：\(O(n)\)。

/*
  ! 如果没有天赋，那就一直重复
  ! Created: 2024/08/13 08:47:10
*/
#include <bits/stdc++.h>
using namespace std;

#define x first
#define y second
// #define int long long
#define mp(x, y) make_pair(x, y)
#define eb(...) emplace_back(__VA_ARGS__)
#define fro(i, x, y) for (int i = (x); i <= (y); i++)
#define pre(i, x, y) for (int i = (x); i >= (y); i--)
inline void JYFILE19();

using i64 = long long;
using pii = pair<int, int>;

bool ST;
const int N = 3e7 + 10;
const int M = 2e6 + 10;
const int mod = 998244353;

int n, ct, w0, w1, w2, ans;
int b[N];
int p[M];
int h[M][61];
int pr[N];
int pd[N];
bool vs[N];

inline void add(int&x, int y) { x = x + y; if (x > mod) x -= mod; }
inline void add(int&x, int y, int z) { x = (x + 1ll * y * z) % mod; }
inline void del(int&x, int y) { x = x - y; if (x < 0) x += mod; }
inline void del(int&x, int y, int z) { x = (x - 1ll * y * z) % mod; if (x < 0) x += mod; }
inline i64 F(i64 x) {
  return (w0 + x * w1 + x * x % mod * w2) % mod;
}
inline i64 f(i64 x, i64 y) {
  return F(x % mod * y % mod);
}
inline i64 g(i64 x, i64 y) {
  if (x == 1 && y == 1) return 1;
  if (x == 1) return f(1, y);
  if (y == 1) return f(x, 1);
  return f(x, 1) * f(1, y) % mod;
}
inline void init() {
  fro(i, 2, n) {
    if (!vs[i]) {
      pr[++ct] = i, pd[i] = i, b[i] = F(i);
      for (int j = i; j <= n / i; j = j * i) b[j * i] = F(j * i);
    }
    for (int j = 1; j <= ct && pr[j] * i <= n; j++) {
      vs[pr[j] * i] = 1;
      if (i % pr[j] == 0) {
        pd[pr[j] * i] = pd[i] * pr[j];
        break;
      }
      pd[pr[j] * i] = pr[j];
    }
  }
  b[1] = 1;
  fro(i, 2, n) {
    if (pd[i] != i) {
      b[i] = 1ll * b[pd[i]] * b[i / pd[i]] % mod;
    }
  }
  fro(i, 1, n) add(b[i], b[i - 1]);
  fro(i, 1, ct) {
    for (i64 j = 1; j <= n / pr[i]; j = j * pr[i]) p[i]++;
    h[i][0] = 1;
    for (i64 j = 1, p1 = pr[i]; j <= p[i] + p[i]; p1 = p1 * pr[i], j++) {
      int d = f(p1, pr[i]); del(d, g(p1, pr[i]));
      for (i64 l = 1, q1 = p1 / pr[i]; l < j; q1 = q1 / pr[i], l++)
        del(d, g(q1, 1), h[i][l]);
      h[i][j] = d;
    }
  }
}
inline void dfs(int x, int p1, int p2, int w) {
  add(ans, 1ll * b[n / p1] * b[n / p2] % mod * w % mod);
  for (i64 i = x; i <= ct && p1 <= n / pr[i] && p2 <= n / pr[i]; i++) {
    i64 p = pr[i];
    if (max(p1, p2) > n / p) return;
    if (p * p > n) {
      dfs(i + 1, p1 * p, p2 * p, 1ll * w * h[i][1] % mod);
    } else {
      for (i64 l = 1, q1 = p * p1; q1 <= n; q1 = q1 * p, l++)
        for (i64 r = 1, q2 = p * p2; q2 <= n; q2 = q2 * p, r++)
          dfs(i + 1, q1, q2, 1ll * w * h[i][l + r - 1] % mod);
    }
  }
}

signed main() {
  JYFILE19();
  cin >> n >> w0 >> w1 >> w2, init(), dfs(1, 1, 1, 1), cout << ans << "\n";
  return 0;
}

bool ED;
inline void JYFILE19() {
  // freopen("", "r", stdin);
  // freopen("", "w", stdout);
  srand(random_device{}());
  ios::sync_with_stdio(0), cin.tie(0);
  double MIB = fabs((&ED - &ST) / 1048576.), LIM = 1024;
  cerr << "MEMORY: " << MIB << endl, assert(MIB <= LIM);
}