博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
BZOJ4816 SDOI2017 数字表格 莫比乌斯反演
阅读量:4654 次
发布时间:2019-06-09

本文共 2940 字,大约阅读时间需要 9 分钟。


做莫比乌斯反演题显著提高了我的\(\LaTeX\)水平

推式子(默认\(N \leq M\),分数下取整,会省略大部分过程)

\(\begin{align*} \prod\limits_{i=1}^N \prod\limits_{j=1}^M f[gcd(i,j)] & = \prod\limits_{d=1}^N f[d]^{\sum\limits_{i=1}^\frac{N}{d} \sum\limits_{j=1}^\frac{M}{d}[gcd(i,j)==1]} \\ & = \prod\limits_{d = 1}^N f[d]^{\sum\limits_{p=1}^\frac{N}{d} \mu(p) \frac{N}{dp} \frac{M}{dp}} \end{align*}\)

推到这里可以\(O(TN)\)地通过两个数论分块得出答案,可以获得70pts

当然这还不够,按照老套路枚举\(dp\)继续推式子

\(\begin{align*} \prod\limits_{d = 1}^N f[d]^{\sum\limits_{p=1}^\frac{N}{d} \mu(p) \frac{N}{dp} \frac{M}{dp}} & = \prod\limits_{T=1} ^ N \prod\limits_{d | T} f[d] ^ {\mu (\frac{T}{d}) \frac{N}{T} \frac{M}{T}} \\ & = \prod\limits_{T=1} ^ N \prod\limits_{d | T} (f[d] ^ {\mu (\frac{T}{d})} ) ^ { \frac{N}{T} \frac{M}{T}} \end{align*}\)

\(\frac{N}{T} \frac{M}{T}\)可以数论分块,所以要求\(f[d]^{\mu(\frac{T}{d})}\)的前缀积

当你以为要线性筛的时候……\(N \leq 10^6\)直接枚举倍数乘进去就行了……

总复杂度\(O(NlogN + T\sqrt{N})\)

#include
//This code is written by Itstusing namespace std;inline int read(){ int a = 0; char c = getchar(); bool f = 0; while(!isdigit(c)){ if(c == '-') f = 1; c = getchar(); } while(isdigit(c)){ a = (a << 3) + (a << 1) + (c ^ '0'); c = getchar(); } return f ? -a : a;}const int MOD = 1e9 + 7 , MAXN = 1e6 + 7;bool nprime[MAXN];int fib[MAXN] , inv[MAXN] , mu[MAXN] , times[MAXN] , prime[MAXN];int N , M , cnt;inline int poww(long long a , int b){ int times = 1; while(b){ if(b & 1) times = times * a % MOD; a = a * a % MOD; b >>= 1; } return times;}void init(){ fib[1] = inv[1] = times[1] = times[0] = mu[1] = 1; for(int i = 2 ; i <= N ; ++i){ times[i] = 1; fib[i] = (fib[i - 1] + fib[i - 2]) % MOD; inv[i] = poww(fib[i] , MOD - 2); } for(int i = 2 ; i <= N ; ++i){ if(!nprime[i]){ prime[++cnt] = i; mu[i] = -1; } for(int j = 1 ; j <= cnt && i * prime[j] <= N ; ++j){ nprime[i * prime[j]] = 1; if(i % prime[j] == 0) break; mu[i * prime[j]] = -1 * mu[i]; } } for(int i = 1 ; i <= N ; ++i) for(int j = 1 ; j * i <= N ; ++j) if(mu[j]) times[j * i] = 1ll * times[j * i] * (mu[j] > 0 ? fib[i] : inv[i]) % MOD; for(int i = 2 ; i <= N ; ++i) times[i] = 1ll * times[i - 1] * times[i] % MOD; for(int i = 0 ; i <= N ; ++i) inv[i] = poww(times[i] , MOD - 2);}int main(){ N = 1e6; init(); for(int T = read() ; T ; --T){ N = read(); M = read(); if(N > M) swap(N , M); int ans = 1; for(int i = 1 , pi ; i <= N ; i = pi + 1){ pi = min(N / (N / i) , M / (M / i)); ans = 1ll * ans * poww(1ll * times[pi] * inv[i - 1] % MOD , 1ll * (N / i) * (M / i) % (MOD - 1)) % MOD; } printf("%d\n" , ans); } return 0;}

转载于:https://www.cnblogs.com/Itst/p/10354853.html

你可能感兴趣的文章
杭电多校 Harvest of Apples 莫队
查看>>
C/C++心得-结构体
查看>>
函数名作为参数传递
查看>>
apt-get for ubuntu 工具简介
查看>>
数值计算算法-多项式插值算法的实现与分析
查看>>
day8-异常处理与网络编程
查看>>
Python基础-time and datetime
查看>>
Linux epoll 笔记(高并发事件处理机制)
查看>>
shell脚本练习01
查看>>
WPF图标拾取器
查看>>
通过取父级for循环的i来理解闭包,iife,匿名函数
查看>>
HDU 3374 String Problem
查看>>
数据集
查看>>
[Leetcode] unique paths ii 独特路径
查看>>
HDU 1217 Arbitrage (Floyd + SPFA判环)
查看>>
IntelliJ idea学习资源
查看>>
Django Rest Framework -解析器
查看>>
ExtJs 分组表格控件----监听
查看>>
Hibernate二级缓存配置
查看>>
LoadRunner常用术语
查看>>