A simple Blog for wyx I've been down the bottle hoping.
4815: [Cqoi2017]小Q的表格 数学+分块
发表于: | 分类: Oi | 评论:0 | 阅读:67

C球的这题好难啊QAQ我推了两个小时才推出一个能写的东西出来QAQ

我们来观察他给的这个东西

$b\times f(a,a+b) = (a+b)\times f(a,b)$

凭着我多年的做(口)题(胡)经验,这东西肯定和 $\gcd$ 有那么点关系

取这么几个值看看玩玩吧……

$\gcd(a,b)f(a,b) = b \times f(a,\gcd(a,b))$

$f(a,\gcd(a,b)) = f(\gcd(a,b),a)$

$\gcd(a,b)f(\gcd(a,b),a) = a \times f(\gcd(a,b),\gcd(a,b))$

看起来中间那个没什么太大的用处我们把他消掉,然后你会有惊人的发现!

$f(a,b) = \frac{ab}{\gcd^2(a,b)} f(\gcd(a,b), \gcd(a,b))$

我们把 $ab$ 看成是一个整体,容易发现这个东西只和 $\frac{f(\gcd(a,b),\gcd(a,b))}{\gcd^2(a,b)}$ 有那么一些关系最开始的时候他等于1,然后在不断修改的过程中这个东西在不断的改变,继而改变了整个值,为了方便我们让这个东西等于 $g(d), d = \gcd(a,b)$

我们现在思考现在求的东西是什么

$ans = \sum_{i=1}^k\sum_{j=1}^k f(i,j)$

$ans = \sum_{d=1}^k g(d) \sum_{i=1}^{k}\sum_{j=1}^k{ij[\gcd(i,j) == d]}$

$ans = \sum_{d=1}^k d^2 g(d) \sum_{i=1}^{\lfloor \frac{k}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{k}{d} \rfloor}ij [\gcd(i,j) == 1]$

把后面那个拿出来反演一下你就得到了一个可以分四块做的东西

这是一个傻逼的真实经历……

其实不用那个样子的……

$\sum_{i=1}^n\sum_{j=1}^n{ij}[\gcd(i,j) == 1]$

$\sum_{i=1}^n i \sum_{j=1}^n {j [\gcd(i,j)==1]}$

我们考虑先只算 $i>j$的部分另一半乘2就行了

那么假如 $j<i$ 就变成了 $i$ 以内和 $i$ 互质的数字的和!

这个东西我们在以前是知道怎么做的,这个东西就等于 $\frac{n\times \phi(n)}{2}$

然后这就是一道喜闻乐见的数据结构题了………………………………

然后树状数组无限TTTTTTT

我们需要一个查询更快的东西……所以分块吧

#include <cmath>
#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long LL;
const long long N = 4e6+5;
const long long mod = 1e9+7;
const long long M = 3e3+5;

inline LL read() {
    LL x=0,f=1;char ch = getchar();
    while(ch < '0' || ch > '9'){if(ch == '-')f=-1;ch = getchar();}
    while(ch >='0' && ch <='9'){x=(x<<1)+(x<<3)+ch-'0';ch = getchar();}
    return x*f;
}

bool F[N];
long long prime[N], cnt, phi[N];
LL sum[N];

inline void init() {
    phi[1] = 1;
    const long long Max = 4e6;
    for(long long i=2;i<=Max;++i) {
        if(!F[i]) {
            prime[++cnt] = i;
            phi[i] = i - 1;
        }
        for(long long j=1;prime[j]*i<=Max;++j) {
            F[i*prime[j]] = 1;
            if(i%prime[j]==0) {
                phi[i*prime[j]] = phi[i] * prime[j];
                break;
            }
            phi[i*prime[j]] = phi[i] * (prime[j]-1);
        }
    }
    for(long long i=1;i<=Max;++i) sum[i] = (sum[i-1] + (LL) i * i % mod * phi[i]%mod) % mod;
}

inline LL pow(LL a,long long b) {
    LL res = 1;
    for(;b;b>>=1,a=(LL)a*a%mod)
        if(b&1) res = res * a % mod;
    return res;
}

long long Q, n, a, b, k;
LL x;

struct Block {
    long long F[N], add[N], a[N];
    long long bl[N], sz, m;
    long long L[N], R[N];

    inline void init() {
        sz = (long long) sqrt(n);
        m = (n-1)/sz + 1;
        for(long long i=1;i<=n;++i) bl[i] = (i-1) / sz + 1, a[i] = (LL)i*i%mod, F[i] = a[i] + F[i-1] % mod;
        for(long long i=1;i<=m;++i) L[i] = (i-1)*sz + 1, R[i] = i * sz;
    }

    inline long long ask(long long x) {
        return (F[x] + add[bl[x]])  % mod;
    }

    inline void change(long long x,long long val) {
        long long tt = (val - a[x] + mod) % mod;
        long long r = R[bl[x]];
        a[x] = val;
        if(tt == 0) return;
        for(long long i=x;i<=r;++i) F[i] = (F[i] + tt) % mod;
        for(long long i=bl[x]+1;i<=m;++i) add[i] = (add[i] + tt)  % mod;
    }
} block;

inline void solve(long long n) {
    long long ans = 0,r , last = 0, now;
    for(long long i=1;i<=n;i=r+1,last=now) {
        r = n/(n/i);
        now = block.ask(r);
        ans = (ans + (LL) (now - last + mod)*sum[n/i] % mod) % mod;
    }
    printf("%d\n", (ans+mod)%mod);
}

int main() {
    Q = read() , n = read();
    init();
    block.init();
    for(long long i=1;i<=Q;++i) {
        a = read(), b = read();
        x = read(), k = read();
        x %= mod;
        long long d = __gcd(a, b);
        block.change(d, (LL)d*d%mod*x%mod*pow((LL)a*b%mod, mod-2)%mod);
        solve(k);
    }
}

Title - Artist
0:00

站点地图 网站地图
Copyright © 2015-2017 A simple Blog for wyx
Powered by Typecho自豪的采用Sgreen主题

TOP