博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
洛谷 3784(bzoj 4913) [SDOI2017]遗忘的集合——多项式求ln+MTT
阅读量:4480 次
发布时间:2019-06-08

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

题目:

   

和洛谷3489“付公主的背包”一样的套路。

要设 a[ i ] 表示第 i 个值有没有出现。

然后就有 \( \prod\limits_i(\frac{1}{1-x^i})^{a_i} = f(x) \)

因为有 \( \prod \) ,所以两边取 ln 。

\( \sum\limits_{i}a_{i}ln(\frac{1}{1-x^i}) = ln(f(x)) \)

现在想求一个 \( ln(\frac{1}{1-x^i}) \) 的更优美的形式(一般是形如 \( \sum \) 的),来更简单地刻画 a[ i ] 和 f[ i ] 的关系。(f[ i ] 是 ln( f(x) ) 的第 i 项系数)

因为有 \( ln \) ,所以先求导再积分来化式子。

并且 \( \frac{f'(x)}{f(x)} \) 了之后,把 \( f'(x) \) 写成 \( \sum \) 的形式,用 \( f(x) \) 和 \( \int \) 化出一个更好看的 \( \sum \) 的式子。

\( \int (1-x^i)\sum\limits_{j=1}i*j*x^{i*j-1} \) // j 从 1 开始

\( = \int \sum\limits_{j=1}i*j*x^{i*j-1} - \sum\limits_{j=1}i*j*x^{i*(j+1)-1} \)

\( = \int \sum\limits_{j=1}i*x^{i*j-1} \)

\( = \sum\limits_{j=0}\frac{1}{j}*x^{i*j} \)

所以 \( \sum\limits_{i=1}a_i\sum\limits_{j=0}\frac{1}{j}x^{i*j} = ln(f(x)) \)

  \( \sum\limits_{i=1}\sum\limits_{j=0}a_i*\frac{1}{j} = f[i*j] \)

\( f[i]=\sum\limits_{j\|i}a_j*\frac{j}{i} \)

把分母的 i 乘到左边,然后莫比乌斯反演一下就知道 \( a_i *i= \sum\limits_{j\|i}f[j]*j*u(i/j) \)

实现的时候要写 MTT 。写拆系数 FFT 的话需要 long double 。自己写的三模数 NTT 还没调出来,不知是哪里出错。

有许许多多的细节需要注意。

#include
#include
#include
#include
#define db long double#define ll long longusing namespace std;int rdn(){ int ret=0;bool fx=1;char ch=getchar(); while(ch>'9'||ch<'0'){
if(ch=='-')fx=0;ch=getchar();} while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar(); return fx?ret:-ret;}const int N=(1<<19)+5;int n,p,f[N],g[N],u[N],pri[N]; bool vis[N];int upt(int x){
if(x>=p)x-=p;if(x<0)x+=p;return x;}int pw(int x,int k){
int ret=1;while(k){
if(k&1)ret=(ll)ret*x%p;x=(ll)x*x%p;k>>=1;}return ret;}namespace poly{ const db pi=acos(-1); struct cpl{ db x,y; cpl(db x=0,db y=0):x(x),y(y) {} cpl operator+ (const cpl &b)const {
return cpl(x+b.x,y+b.y);} cpl operator- (const cpl &b)const {
return cpl(x-b.x,y-b.y);} cpl operator* (const cpl &b)const {
return cpl(x*b.x-y*b.y,x*b.y+y*b.x);} cpl operator/ (const int &b)const {
return cpl(x/b,y/b);} }; cpl conj(cpl a){
return cpl(a.x,-a.y);} int len,r[N],inv[N]; cpl Wn[N]; int bs,pbs,bs2; cpl pa[N],pb[N],pc[N],pd[N]; int A[N],B[N],tp[N]; void init() { int tmp=sqrt(p); for(bs=0,pbs=1;pbs<=tmp;bs++,pbs<<=1); bs2=bs<<1; pbs--; } void fft_pre() { for(int i=0,j=len>>1;i
>1]>>1)+((i&1)?j:0); for(int R=2,m=1;R<=len;m=R,R<<=1) Wn[R]=cpl( cos(pi/m),sin(pi/m) ); } void fft(cpl *a,bool fx) { for(int i=0;i
>1;i
>15,a[i]&32767); //for(int i=0;i
>15,b[i]&32767); for(int i=0;i
>bs,a[i]&pbs); for(int i=0;i
>bs,b[i]&pbs); for(int i=n1;i
>1;i
拆系数FFT
#include
#include
#include
#define ll long longusing namespace std;int rdn(){ int ret=0;bool fx=1;char ch=getchar(); while(ch>'9'||ch<'0'){
if(ch=='-')fx=0;ch=getchar();} while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar(); return fx?ret:-ret;}int upt(int x,int mod){
while(x>=mod)x-=mod;while(x<0)x+=mod;return x;}int pw(int x,int k,int mod){
int ret=1;while(k){
if(k&1)ret=(ll)ret*x%mod;x=(ll)x*x%mod;k>>=1;}return ret;}const int N=(1<<19)+5; int p;namespace poly{ const double eps=1e-6; int m[3]={
998244353,1004535809,469762049}; ll M=(ll)m[0]*m[1], A[N],B[N],C[3][N]; int len,r[N],Wn[N][2],inv[N]; int tp[N],ta[N],tb[N]; ll mul(ll a,ll b,ll mod) { a=(a%mod+mod)%mod; b=(b%mod+mod)%mod;///// ll ret=(a*b- (ll)((long double)a/mod*b+eps) *mod)%mod; if(ret<0)ret+=mod; return ret; } void ntt_pre(int len,int mod) { for(int R=2;R<=len;R<<=1) Wn[R][0]=pw( 3,(mod-1)/R,mod ), Wn[R][1]=pw( 3,(mod-1)-(mod-1)/R,mod ); } void ntt(ll *a,bool fx,int mod) { for(int i=0;i
>1;i
>1;i
>1]>>1)+((i&1)?j:0); for(int i=0;i<3;i++) { mod=m[i]; for(int j=0;j
三模数NTT(TLE+WA)

 

转载于:https://www.cnblogs.com/Narh/p/10405656.html

你可能感兴趣的文章
Redis客户端基本命令
查看>>
第6章 XHTML:Web重构
查看>>
left join多表操作
查看>>
在pcDuino上刷了AndDroid,Ubuntu,XBMC
查看>>
1.C#基础学习笔记3---C#字符串(转义符和内存存储无关)
查看>>
delete了,析构函数却没有调用
查看>>
POJ 3264 Balanced Lineup
查看>>
removeFromSuperview 添加动画-by小雨
查看>>
java 操作格子问题(线段树)
查看>>
ARM32 Linux kernel virtual address space
查看>>
用C语言扩展Python的功能
查看>>
[MySQL 5.6] 初识5.6的optimizer trace
查看>>
LabVIEW(数据库连接)
查看>>
【CSDN博客之星】2013年CSDN博客之星正在评选,希望大家支持投票,非常感谢!...
查看>>
更新内容:关于Windows Azure技术内容搜索的新页面可用
查看>>
Web Service 的工作原理
查看>>
LaTeX技巧003:实现一个章标题
查看>>
第一次冲刺--站立会议07
查看>>
linux下pip的安装
查看>>
自我介绍
查看>>