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世界树

我们用一个栈维护当前构建虚树的最右链并将\(k\)个节点按照\(dfn\)排序,模拟\(dfs\)的过程依次插入。

对于每一个插入的节点\(x\),与栈顶元素取\(lca\),设\(lca(x,stk[top])=c\),那么依次取栈顶分情况讨论:

  1. \(dfn[c]=stk[top-1]\),即\(c\)为维护的栈中的元素
  2. \(dfn[c]>stk[top-1]\),即\(c\)在\(stk[top]\)和\(stk[top-1]\)之间
  3. \(dfn[c]<stk[top-1]\),即\(c\)在\(stk[top-1]\)之上

对于情况\(3\),直接把\(stk[top]\)退栈,并在\(stk[top]\)和\(stk[top-1]\)之间连边。

对于情况\(2\),把\(stk[top]\)退栈并在\(stk[top]\)和\(c\)之间连边,把\(c\)加入栈,退出。

对于情况\(1\),把\(stk[top]\)退栈并在\(stk[top]\)和\(stk[top-1]\)之间连边,退出。

每次能直接退栈的原因是该子树已经遍历完毕,不会对后来的建树产生影响。

实现

const int N = 10000 + 5;

int dfn[N];

bool cmp(int a, int b) {
  return dfn[a] < dfn[b];

}
int top;
int stk[N];

inline void init() {
  for (int i = 0; i < cnt; ++i) {
    b[i] = read();
  }
  sort(b, b + cnt, cmp);
  top = 0;
  stk[top++] = b[0];
  for (int i = 1; i < cnt; ++i) {
    if (top == 0) {
      stk[top++] = b[i];
      continue;
    }
    int c = lca(stk[top - 1], b[i]);
    while (top > 0 && dfn[c] < dfn[stk[top - 1]]) {
      if (top == 1 || dfn[c] >= dfn[stk[top - 2]]) {
        add(c, stk[top - 1], dep[stk[top - 1]] - dep[c]);
        --top;
        if (top == 0 || stk[top - 1] != c) {
          stk[top++] = c;
        }
        break;
      }
      add(stk[top - 2], stk[top - 1], dep[stk[top - 1]] - dep[stk[top - 2]]);
      --top;
    }
    stk[top++] = b[i];
  }
  while (top > 1) {
    add(stk[top - 2], stk[top - 1], dep[stk[top - 1]] - dep[stk[top - 2]]);
    --top;
  }
}

Description

给出一个字符串,求最大重复子串(重复次数最多,如果存在多个,求字典序最小的那一个)。

Solution

后缀数组的应用。直接下手不好解决,不妨枚举一下循环节的长度。我们发现,任何一个循环节为重复子串总会包含至少两个 字符。那么考虑枚举两个相邻的上述字符,可以通过后缀数组求出的长度,但是最长公共子串的开头并不一定是我们枚举的字符,所以还需要求出最长向前能匹配多少。这可以通过倒过来做一次后缀数组得到。那么我们现在有了一个极长区间,可以求得这个区间的循环节个数,也就可以求出一个区间满足开头落在这个区间内部的最大重复子串的循环节个数都为。只需要找字典序最小的一个。那么用表查一下这个区间内最小的的后缀就好了。时间复杂度。有一个优化,当求得一个极长重复子串之后,落在子串内的都不用再枚举了。

#include <map>
#include <cmath>
#include <vector>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>

using namespace std;

const int N = 100000 + 5;
const int LOG = 17;

int n;
int m;

char s[N];

struct suffix_array {
  char s[N];
  
  int c[N];
  int t0[N];
  int t1[N];
  int sa[N];

  int rk[N];
  int ht[N];

  inline void build() {
    int *x = t0;
    int *y = t1;
    m = 256;
    s[n++] = 0;
    for (int i = 0; i < m; ++i) {
      c[i] = 0;
    }
    for (int i = 0; i < n; ++i) {
      x[i] = s[i];
      ++c[x[i]];
    }
    for (int i = 1; i < m; ++i) {
      c[i] += c[i - 1];
    }
    for (int i = n - 1; i >= 0; --i) {
      int t = --c[x[i]];
      sa[t] = i;
    }
    for (int k = 1; k <= n; k <<= 1) {
      int p = 0;
      for (int i = n - k; i < n; ++i) {
        y[p++] = i;
      }
      for (int i = 0; i < n; ++i) {
        if (sa[i] >= k) {
          y[p++] = sa[i] - k;
        }
      }
      for (int i = 0; i < m; ++i) {
        c[i] = 0;
      }
      for (int i = 0; i < n; ++i) {
        ++c[x[y[i]]];
      }
      for (int i = 1; i < m; ++i) {
        c[i] += c[i - 1];
      }
      for (int i = n - 1; i >= 0; --i) {
        int t = --c[x[y[i]]];
        sa[t] = y[i];
      }
      p = 1;
      swap(x, y);
      x[sa[0]] = 0;
      for (int i = 1; i < n; ++i) {
        if (y[sa[i]] == y[sa[i - 1]] && y[sa[i] + k] == y[sa[i - 1] + k]) {
          x[sa[i]] = p - 1;
        } else {
          x[sa[i]] = p++;
        }
      }
      if (p >= n) {
        break;
      }
      m = p;
    }
    --n;
    for (int i = 0; i < n; ++i) {
      sa[i] = sa[i + 1];
    }
    for (int i = 0; i < n; ++i) {
      rk[sa[i]] = i;
    }
    int k = 0;
    for (int i = 0; i < n; ++i) {
      if (k) {
        --k;
      }
      if (rk[i] == 0) {
        k = 0;
        continue;
      }
      int j = sa[rk[i] - 1];

      while (s[i + k] == s[j + k]) {
        ++k;
      }
      ht[rk[i]] = k;
    }
  }

  int mn[N][LOG];

  inline void init_st() {
    for (int i = 0; i < n; ++i) {
      mn[i][0] = ht[i];
    }
    for (int k = 1; k < LOG; ++k) {
      for (int i = 0; i < n; ++i) {
        mn[i][k] = mn[i][k - 1];
        if (i + (1 << (k - 1)) < n) {
          mn[i][k] = min(mn[i][k], mn[i + (1 << (k - 1))][k - 1]);
        }
      }
    }
  }

  inline int query(int l, int r) {
    if (l > r) {
      swap(l, r);
    }
    ++l;
    int k = log(r - l + 1) / log(2);
    return min(mn[l][k], mn[r - (1 << k) + 1][k]);
  }

};

suffix_array pre;
suffix_array suf;

int mx = 0;
int cur;

int mn[N][LOG];

inline void init_rk() {
  for (int i = 0; i < n; ++i) {
    mn[i][0] = suf.rk[i];
  }
  for (int k = 1; k < LOG; ++k) {
    for (int i = 0; i < n; ++i) {
      mn[i][k] = mn[i][k - 1];
      if (i + (1 << (k - 1)) < n) {
        mn[i][k] = min(mn[i][k], mn[i + (1 << (k - 1))][k - 1]);
      }
    }
  }
}

inline int query_rk(int l, int r) {
  int k = log(r - l + 1) / log(2);
  return min(mn[l][k], mn[r - (1 << k) + 1][k]);
}

int main() {
  int kase = 0;
  while (true) {
    scanf("%s", s);
    if (s[0] == '#') {
      return 0;
    }
    n = strlen(s);
    
    memcpy(suf.s, s, sizeof(s));
    suf.build();
    suf.init_st();
    reverse(s, s + n);
    memcpy(pre.s, s, sizeof(s));
    pre.build();
    pre.init_st();
    reverse(s, s + n);

    init_rk();

    mx = 0;
    int ans_len = 0;
    for (int l = 1; l <= n; ++l) {
      for (int i = 0; (i + 1) * l < n; ++i) {
        int p = i * l;
        int q = (i + 1) * l;
        if (s[p] != s[q]) {
          continue;
        }
        int suf_max = suf.query(suf.rk[p], suf.rk[q]);
        int pre_max = pre.query(pre.rk[n - p - 1], pre.rk[n - q - 1]);
        int low = p - pre_max + 1;
        int t = (suf_max + pre_max - 1) / l + 1;
        int t0 = l * t - l - 1;
        int high = p + suf_max - 1 - t0;
        if (high < low) {
          continue;
        }
        int mn_rk = suf.sa[query_rk(low, high)];

        if (t == mx) {
          if (suf.rk[cur] > suf.rk[mn_rk]) {
            cur = mn_rk;
            ans_len = t * l;
          }
        } else if (t > mx) {
          mx = t;
          cur = mn_rk;
          ans_len = t * l;
        }
        i = (p + pre_max - 1) / l;
      }
     
    }
    printf("Case %d: ", ++kase);
    for (int i = cur; i < cur + ans_len; ++i) {
      putchar(s[i]);
    }
    puts("");
  }
  return 0;
}

结束了…似乎按照惯例要写一篇游记?

Day 0

上午做个3个多小时的高铁到秦皇岛。陪妹子玩尼尔真有意思,不过似乎在火车打游戏上引来很多关注?刚下车看到不少学校没有妹子学OI,真可怜。还是跟去年一样住在7天,跟一间屋子很资瓷啊。晚上跟学长一起吃火锅,我虐狗了吗? 被灌一点酒,真是不要命了。没敢吃太多,毕竟明天考试是不是?时间很快就过去,晚上被叫去开会,听了听大家犯过的脑残错误,听了听学长传授经验,什么二中校歌,我咋不知道。 看了看以前做过的题收了收心,我还做过这种题!10点了还在补番,我就去睡觉了,意外的睡的不错。

Day 1

一觉醒来发现6点50了,穿好衣服急匆匆感到楼下吃早餐5分钟噎完饭又跑回房间上厕所,md果然还是吃多了,擦着点集合去燕山大学,路上到是有说有笑,不知道到了考场上还能不能笑的出来呢,看的出来大家都是很放松的。到了考场突然意识到没带水,又满楼转悠找水,最后还是没找到。看了一眼座位号6-4,正好坐在大屏幕底下。开题开题!先看T1。woc画风不对啊,模拟呢?先一发xjb推了个式子,拍了拍第一组就拍出错了。。。看了一眼表刚过5分钟,淦!不虚,改!突然发现式子有点问题,改完就拍过了。我有100了不虚,看T2。大模拟?T1T2是不是放反了?想了想细节开始码,一次过大样例。不过不好拍就放了。看T3。分层图?大概一下就好了?等等怎么拓扑?看了一眼表还有3h,时间怎么说也够了把(flag)。看部分分,70分都没有0权边?先把70分码了。 然后淦正解。。。。如果不拓扑用强跑,大样例,盲目改了一会儿没有用,看来必须拓扑,突然意识到只拓扑0权边就好了?看了一眼表还有20min,陷入沉思,开始怀疑自己的算法事实证明是正确的。。。想了想270了,一狠心md30分不要了查了查文件名就坐等收题了。

出考场发现全世界T1都是找规律。T2似乎爆炸了?似乎很多人T1都没做出来?妹子似乎没考好?下午陪妹子去海边,开导了半天。晚上回来继续开会。看了看以前做过的题收了收心,我还做过这种题! ,就去睡觉了,意外的睡的不错。

Day 2

起的稍微早了一点,早早到了考场。md还坐在大屏幕底下。开题!这T1?想了想浮点数的精度误差把删了强行平方。10min敲完鏼过大样例。突然发现可能爆。改完看T2。? 先状压一发直接用表示中的点都被用的最小代价,转移直接建出树来。突然发现没法处理权值相等的情况,算了算时间复杂度md这啥啊看来写的不是正解,不过鏼过了大样例就不虚了,应该会被卡掉?谁知道呢。再看T3。等等吉司机的数据结构题?开始考数据结构了?先把30暴力敲完,看了看的,应该线段树也能做,啪啪啪敲完,有60了然而忘了开long long。正解鏼不出来弃疗,看了一眼表还有2h。怎么这么长时间,于是又开始想T3的正解,发现能平衡树暴搞?很久都没有敲过平衡树了于是开始现调。。。。最后自然没有调出来。致死都没有发现要开long long

出分::100+100+60。被老爷机卡掉一个点。:100+95+50。才发现T3要开。良心出题人D2T2才卡我5分?总分505,貌似是河北省除集训队选手外(ORZ yzy)rk1?我。。我可是暴力选手啊。。不过貌似全国排名跪了,还要多多增加姿势水平!北校的也都挂了?高分被衡水刷屏了啊。妹子也考挂了。。sad。。

尼尔真好玩

机房的人一下子少了不少呢。。。。