51Nod 1238 最小公倍数之和 V3

Reading time ~2 minutes

Description

求: \[ \sum_{i=1}^n \sum_{j=1}^n lcm(i,j) \] $n \leq 10^{10}$

Solution

我们要求的就是下面这个式子: \[ \sum_{i=1}^n \sum_{j=1}^n \frac{ij}{\gcd(i,j)} \] 显然可以反演,但是反演式并不好快速求和,复杂度难以接受。

化成$\mu$不行但是可以试试化成$\varphi$。我们设: \[ f(x)=x\sum_{i=1}^x \frac{i}{\gcd(x,i)} \] 枚举$x$的因子作为$\gcd$,原式化成: \[ f(x)=x\sum_{d|x} j \cdot \varepsilon(\gcd(\frac{i}{d},j)) \] 注意到后面的部分就是对小于等于$\frac{i}{d}$且与$\frac{i}{d}$互质的数之和,也就是$\frac{\frac{i}{d}\varphi(\frac{i}{d})+\varepsilon(\frac{i}{d})}{2}$。

我们的答案等于: \[ ANS=2\sum_{i=1}^{n} f(i) - \frac{n(n+1)}{2} \] 带入$f(x)$并化简,得到: \[ ANS=\sum_{k=1}^{n} k\sum_{d=1}^{\lfloor\frac{n}{d}\rfloor} d^2\varphi(d) \] 设: \[ S(n)=\sum_{i=1}^{n} i^2\varphi(i) \] 问题转变为怎么快速求$S(n)$。

显然有:$id^2 \cdot[\varphi \times1]=id^3$,我们化简这个式子,得到: \[ id^3(n)=\sum_{k=1}^{n} k^2 S(\lfloor \frac{n}{k}\rfloor) \]

直接上杜教筛就行了。预处理前$n^{\frac{2}{3}}$的前缀和,注意到我们求出$S(n)$的时候就已经求出了所有的$S(\lfloor \frac{n}{d}\rfloor)$,复杂度跟求一次前缀和相同,为$O(n^{\frac{2}{3}})$。

#include <bits/stdc++.h>

using namespace std;

const int P = 1e9 + 7;
const int LIM = 4641600;
const int N = 4641600 + 5;

long long n;

inline int add(int a, int b) {
  a += b;
  return a >= P ? a - P : a;
}

inline int sub(int a, int b) {
  a -= b;
  return a < 0 ? a + P : a;
}

inline int mul(int a, int b)  {
  unsigned long long x = (long long) a * b;
  unsigned xh = (unsigned) (x >> 32), xl = (unsigned) x, d, m;
  asm(
    "divl %4; \n\t"
    : "=a" (d), "=d" (m)
    : "d" (xh), "a" (xl), "r" (P)
  );
  return m;
}

int sum(long long x) {
  int k = x % P;
  return mul(mul(k, k + 1), 500000004);
}

int sum2(long long x) {
  int k = x % P;
  return mul(mul(mul(k, k + 1), (2 * k + 1) % P), 166666668);
}

struct LinearSieve {
  int p[N];
  int tot;
  int f[N];
  int sum[N];
  int notP[N];

  void solve(int n) {
    f[1] = 1;
    for (int i = 2; i <= n; ++i) {
      if (!notP[i]) {
        p[tot++] = i;
        f[i] = mul(mul(i, i), i - 1);
      }
      for (int j = 0; j < tot; ++j) {
        int k = p[j] * i;
        if (k > n) {
          break;
        }
        notP[k] = 1;
        if (i % p[j] == 0) {
          f[k] = mul(f[i], mul(mul(p[j], p[j]), p[j]));
          break;
        } else {
          f[k] = mul(f[i], mul(mul(p[j], p[j]), p[j] - 1));
        } 
      }
    }

    for (int i = 1; i <= n; ++i) {
      sum[i] = add(sum[i - 1], f[i]);
    }
  }
} s;

int mem[N];

int calc(long long x) {
  if (x <= LIM) {
    return s.sum[x];
  } else if (mem[n / x]) {
    return mem[n / x];
  }
  int ans = mul(sum(x), sum(x));
  for (long long i = 2; i <= x; ++i) {
    long long vi = x / i;
    long long nxt = x / vi;
    ans = sub(ans, mul(sub(sum2(nxt), sum2(i - 1)), calc(vi)));
    i = nxt;
  }
  return mem[n / x] = ans;
}

int main() {
  s.solve(LIM);
  scanf("%lld", &n);
  int ans = 0;
  for (long long i = 1; i <= n; ++i) {
    long long vi = n / i;
    long long nxt = n / vi;
    ans = add(ans, mul(sub(sum(nxt), sum(i - 1)), calc(vi)));
    i = nxt;
  }
  printf("%d\n", ans);
  return 0;
}

虚树学习笔记

算法思想 实现算法思想当问题的求解只涉及到树中的\(k\)个节点时,为了确保复杂度只与\(k\)相关,可选用的做法是把这\(k\)个节点提出来新建一棵树,我们管这颗新建的树叫虚树。资料:虚树详解+例子分析+模板BZOJ3572 Hnoi2014世界树我们用一个栈维护当...… Continue reading

POJ 3693 Maximum repetition substring

Published on November 24, 2017

NOIP2017 爆炸记

Published on November 24, 2017