Polynomial - 演算法筆記

Dense Polynomial ： array ，一格存一項。

索引值是次方，內容是係數。

₀ ₁ ₂ ₃ ₄ ₅ ₆ ₇ ₈ ₉. |8|4|0|2 ... Polynomial Polynomial 「多項式」。

變數的各種次方、乘上倍率、通通相加。

8x⁰-4x¹+0x²+2x³上升排列 2x³+0x²-4x¹+8x⁰下降排列 Polynomial資料結構 DensePolynomial：array，一格存一項。

索引值是次方，內容是係數。

₀₁₂₃₄₅₆₇₈₉ ___________________ |8|4|0|2||||||| ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾ SparsePolynomial：list或array，一份存一項。

內容是係數和次方。

係數為零的項，統統去掉，串成一串。

_________ ⤷|8|0|→|4|1|→|2|3| ‾‾‾‾‾‾‾‾‾ Polynomial類型 多項式視作數值：底數是x，位數是多項式係數。

多項式視作函數：輸入是x，輸出是多項式求值。

多項式可能無限項：有限多項式、無限多項式。

多項式必須搭配數域：整數多項式、有理數多項式、……。

Polynomial運算 加、減、乘、除、微、積、展開、分解、代入（求值）、求根（求解）。

addition(x²+3x+2)+(x+2)=x²+4x+4 subtraction(x²+3x+2)-(x+2)=x²+2x multiplication(x²+3x+2)×(x+2)=2x³+5x²+8x+4 division(x²+3x+2)÷(x+2)=x+1 expansion(x+1)²=x²+2x+1 factorizationx²+2x+1=(x+1)² differentiationd/dx(x²+3x+2)=2x+3 integration∫(x²+3x+2)dx=¹⁄₃x³+³⁄₂x²+2x+c evaluation(x²+3x+2)|x=1=6 resolution(x²+3x+2)=0⇒x=-1or-2 度數 最高次方。

f(x)=8x⁰-4x¹+0x²+2x³ deg(f)=3 加法、減法 各項一一配對、係數相加。

f(x)=8x⁰-4x¹+0x²+2x³ g(x)=1x⁰+2x¹ f(x)+g(x)=9x⁰-2x¹+0x²+2x³ f(x)-g(x)=7x⁰-6x¹+0x²+2x³ 乘法 各項交叉配對、係數相乘、次方相加。

f(x)=8x⁰-4x¹+0x²+2x³ g(x)=1x⁰+2x¹ f(x)g(x)=8x⁰1x⁰-4x¹1x⁰+0x²1x⁰+2x³1x⁰ +8x⁰2x¹-4x¹2x¹+0x²2x¹+2x³2x¹ 帶餘除法 帶餘除法的輸出是兩個多項式：商數、餘數。

長除法，一直補位、一直求餘數，直到餘數度數小於除數度數。

f(x)=8x⁰-4x¹+0x²+2x³+3x⁴ g(x)=1x⁰+2x¹+3x² f(x)=(¹⁄₃x⁰+0x¹+1x²)g(x)+(²⁄₃x⁰-¹⁰⁄₃x¹) 無窮除法 無窮除法的輸出是一個多項式：商數。

長除法，一直補位、一直求餘數，商數越來越長。

1-x ――――――――――――=1-3x+3x²+3x³+... 1+2x+3x² 知名範例是GeometricSeries。

a是任意數，例如1。

1/(1-ax)=(ax)⁰+(ax)¹+(ax)²+... 1/(1-x)=x⁰+x¹+x²+... 可能除得盡：餘數是零，商數是多項式。

可能除不盡：餘數不是零，商數是無窮多項式。

展開 多項式的連乘、多項式的次方，化作一個多項式。

f(x)=2x⁰(2x⁰+1x¹)(2x⁰-2x¹+x²) =8x⁰-4x¹+0x²+2x³ f(x)=(2x⁰+1x¹)² =4x⁰+4x¹+1x² 分解 展開的反方向。

f(x)=8x⁰-4x¹+0x²+2x³ =2x⁰(2x⁰+1x¹)(2x⁰-2x¹+1x²) 微分、積分 微分結果仍是多項式函數！dx竟然變不見！ 積分結果仍是多項式函數！常數項是任意數，通常設成0。

f(x)=8x⁰-4x¹+0x²+2x³ d/dxf(x)=-4x⁰+0x¹+6x² ∫f(x)dx=0x⁰+8x¹-2x²+0x³+½x⁴ 代入（求值） Horner'sRule。

從最高次方項開始，一乘一加，得到答案。

不需要次方運算。

f(x)=8x⁰-4x¹+0x²+2x³ =(((((2⋅x)+0)⋅x)-4)⋅x)+8 求根（求解） 用途不明。

f(x)=8x⁰-4x¹+0x²+2x³=0 x=-2or1+𝑖or1-𝑖 快速多項式運算 多項式運算，普通演算法是O(N²)。

幸運的是，多項式乘法，運用複數系統、餘數系統，運用快速傅立葉轉換、快速數論轉換，時間複雜度降為O(NlogN)，相當接近一次方時間O(N)。

以多項式乘法為根基，各種多項式運算皆可降至接近一次方時間。

多項式運算從此以後變得快速。

拓展數系的三個元件：對稱（自然數到整數）、數對（整數到有理數）、循環（有理數到複數）。

多項式乘法，同時運用三者，得到傅立葉矩陣，形成高速演算法，O(N²)變O(NlogN)。

箇中玄妙之處，令人拍案叫絕。

專著《AlgorithmesEfficacesenCalculFormel》。

專著《AlgorithmsforComputerAlgebra》。

專著《ModernComputerAlgebra》。

三本專著詳細介紹了快速多項式運算。

以下章節不再重複整理，只做簡單介紹。

UVa392126107191058610951463930104284981026810105 Polynomial(DataStructure) data constintN=100000; //degree,butnotsize inta[N+1]; //+1forsize structPoly { vectorx; intdeg; Poly(){} Poly(intdeg){x.resize(deg+1);} inlinevoidresize(intdeg){x.resize(deg+1);} friendinlineint&deg(Poly&p){returnp.deg;} inlineint&operator[](inti){returnx[i];} inlinevoidoperator=(intn){x[deg=0]=n;} }; read/write 截斷版本。

voidread(Poly&a) { for(inti=0;i<=deg(a);i++)cin>>a[i]; } voidwrite(Poly&a) { for(inti=0;i<=deg(a);i++)cout<>(istream&in,Poly&a) { staticcharbuf[N*23]; in>>buf; deg(a)=0; for(char*s=buf;*s;) { intsign=1; if(*s=='+')sign=+1,s++; elseif(*s=='-')sign=-1,s++; intvalue=0; if(!isdigit(*s))value=1; elsewhile(isdigit(*s))value=value*10+(*s++-'0'); intexp=0; if(*s!='x')exp=0; elseif(*++s!='^')exp=1; elsewhile(isdigit(*++s))exp=exp*10+(*s-'0'); deg(a)=max(deg(a),exp); a[exp]=(sign>0?value:-value); } returnin; } ostream&operator<=0;i--) { if(a[i]==0)continue; charsign=(a[i]>=0)?'+':'-'; intvalue=abs(a[i]); if(!leading)out<0;k>>=1) { if(k&1)mul(x,x,p); square(p,p); } b=x; } voidmyfunction() { ...... Polypow; ::pow(pow,a,k); //scoperesolutionoperator ...... } Polynomial(Arithmetic) addition+ 遞增法。

對應項相加。

O(N)。

voidadd(Poly&c,Poly&a,Poly&b,intC) { intA=min(deg(a),C); intB=min(deg(b),C); if(A=0;i--) { intk=c[i]=a[i+deg(b)]/b[deg(b)]; for(intj=deg(b);j>=0;j--) a[i+j]-=k*b[j]; } } voidrev(Poly&b,Poly&a) { if(&a==&b) for(inti=0;i<=deg(a)/2;++i) swap(a[i],a[deg(a)-i]); else for(inti=0;i<=deg(a);i++) b[deg(a)-i]=a[i]; deg(b)=deg(a); } voidcut(Poly&a,intA) { deg(a)=min(deg(a),A); } voiddiv(Poly&q,Poly&r,Poly&a,Poly&b) { intA=deg(a),B=deg(b); if(A=0;i--) { value=mul(value,x); value=add(value,a[i]); } returnvalue; } composition∘ 遞推法。

Horner'sRule。

O(𝓜(AB)A)。

分治法。

A對半分。

O(𝓜(AB)logA)。

加速技巧。

過程都在頻域。

O(𝓜(AB))。

順帶一提，進位轉換可以套用composition。

a₀x⁰+a₁x¹+a₂x²+a₃x³=(a₀x⁰+a₁x¹)+(a₂x⁰+a₃x¹)x² wherex=b₀x⁰+b₁x¹+b₂x² http://blog.csdn.net/liutian429073576/article/details/53644937 令計算結果是截斷多項式。

遞推法。

隨時模xᴺ。

O(𝓜(C)A)。

分治法。

A對半分。

O(𝓜(C)(AB/C)(1+log(min(A,ceil(C/B))))。

簡單寫成O(𝓜(C)A)或O(𝓜(C)BlogA)，後者需C≥A。

加速技巧。

Brent–KungAlgorithm。

B分成低次項k項、高次項B-k項。

泰勒展開，只取低階C/k項。

先求最低階，再以遞推法求得各階導數、冪數、分母。

O(𝓜(C)(klogA+C/k))。

令klogA=C/k讓時間複雜度達到最低。

O(𝓜(C)sqrt(ClogA))。

把B變小（長度k），分治法開頭變快（長度到達C之前），代價是泰勒展開（長度C/k）。

令k=sqrt(C/logA)讓兩者平衡。

voidcom(Poly&c,Poly&a,Poly&b,intC) { //divide-and-conquermethod:depth intL=0; while(1>=1,l++) { //deg(c)=(deg(b)<1;n>>=1,l++) if(n&1) { mul(c,c,p[l],C,m); add(c,c,s[l],C,m); } } //Brent–KungAlgorithm voidcom2(Poly&c,Poly&a,Poly&b,intC) { //balancefactorakamagicnumber constdoublealpha=4; //b_lowdegree intK=sqrt((C+1)/log2(deg(a)+1))*alpha; //Taylorserieslength intL=ceil((double)(C+1-1)/(K+1)); //fastFouriertransformlength intN=1; while(N<=C*2)N<<=1; //Taylorseries:zeroorder if(deg(b)<=K){com(c,a,b,C,m);return;} Polybl(C); //b_low copy(bl,b,K); com(c,a,bl,C+1,m); //C+1fora(b_low)′ //forsummation for(inti=deg(c)+1;i<=C;++i)c[i]=0; deg(c)=C; //derivative staticintder_abl[1< luoguP5373 compositeinverse:inversefunction⁻¹ 遞推法。

SeriesReversion。

O(N²+𝓜(N)N)。

公式解。

LagrangeInversion。

O(N²+𝓜(N)sqrtN)。

luoguP5809 decomposition Kozen–LandauAlgorithm。

O(N²)。

Polynomial(Calculus) differentiationd/dx 逐項微分。

O(N)。

a=a₀x⁰+a₁x¹+a₂x²+a₃x³ d/dxa=1a₁x¹⁻¹+2a₂x²⁻¹+3a₃x³⁻¹ voiddiff(Poly&b,Poly&a,intB) { intA=min(deg(a),B+1); for(inti=1;i<=A;i++)b[i-1]=mul(a[i],i); // b[A]=0; deg(b)=max(A-1,0); } integration∫dx 逐項積分。

O(N)。

a=a₀x⁰+a₁x¹+a₂x²+a₃x³ ∫adx=¹⁄₁a₀x⁰⁺¹+¹⁄₂a₁x¹⁺¹+¹⁄₃a₂x²⁺¹+¹⁄₄a₃x³⁺¹ voidinte(Poly&b,Poly&a,intB) { intA=min(deg(a),B-1); for(inti=A;i>=0;i--)b[i+1]=div(a[i],i+1); b[0]=0; deg(b)=A+1; } naturallogarithmln() 公式解。

多項式乘法視作數列卷積，拆出最後一項。

O(N²)。

公式解。

當輸入的常數項為1，則輸出的常數項為0——得以使用積分。

當輸入的常數項不為1，則提出倍率使之為1——倍率是常數項。

O(𝓜(N))。

whena>0 ln(a)=(a-1)¹/1-(a-1)²/2+(a-1)³/3-... whena>0 ln(a)=b ln′(a)a′=b′ (1/a)a′=b′ a′=ab′ naₙxⁿ⁻¹=(aₙxⁿ)(nbₙxⁿ⁻¹) naₙ=sum{iaₙ₋ᵢbₙ} 1≤i≤n aₙ=sum{iaₙ₋ᵢbₙ}/n 1≤i≤n aₙ-sum{iaₙ₋ᵢbₙ}/n=a₀bₙ 1≤i≤n-1 (aₙ-sum{iaₙ₋ᵢbₙ}/n)/a₀=bₙ 1≤i≤n-1 whena>0anda₀=1 ln(a)=b ln′(a)a′=b′ (1/a)a′=b′ ∫((1/a)a′)=b(ifa₀=1thenb₀=0) whena>0anda₀≠1 ln(a)=ln((a/a₀)a₀)=ln(a/a₀)+ln(a₀) voidlog(Poly&b,Poly&a,intB) { deg(b)=B; b[0]=log(a[0]); for(intn=1;n<=B;++n) { b[n]=0; for(inti=max(1,n-deg(a));i voidlog(Poly&b,Poly&a,intB) { staticPolyx(N); Poly&norm_a=b; //reusememory for(inti=0;i<=deg(a);i++) norm_a[i]=div(a[i],a[0]); Poly&rec_a=x; //reusememory rec(rec_a,norm_a,B); Poly&diff_a=b; //reusememory diff(diff_a,norm_a,B-1); Poly&diff_b=rec_a; //reusememory mul(diff_b,rec_a,diff_a,B-1); inte(b,diff_b,B); b[0]+=log(a[0]); } luoguP4725 naturalexponentiationexp() 公式解。

自然對數公式解，最後不做移項。

O(N²)。

遞推法。

HenselLifting。

O(𝓜(N))。

bₙ=sum{ibₙ₋ᵢaₙ}/n 1≤i≤n b≡exp(a)(modxⁿ) ln(b)≡a(modxⁿ) ln(b)-a≡0(modxⁿ) letf(b)=ln(b)-a,f′(b)≡1/b bꜛ≡b-f(b)/f′(b)(modx²ⁿ) bꜛ≡b-(ln(b)-a)/(1/b)(modx²ⁿ) bꜛ≡b-(ln(b)-a)b(modx²ⁿ) bꜛ≡b(1-ln(b)+a)(modx²ⁿ) voidexp(Poly&b,Poly&a,intB) { deg(b)=B; b[0]=exp(a[0]); for(intn=1;n<=B;++n) { b[n]=0; for(inti=1;i<=min(n,deg(a));++i) b[n]+=i*b[n-i]*a[i]; b[n]/=n; } } voidexp(Poly&b,Poly&a,intB) { deg(y)=0;y[0]=exp(a[0]); for(intn=1;n<=B<<1;n<<=1) { intX=min(n-1,deg(a)); intY=min(n-1,B); staticPolyx(N),y(N); log(x,y,X); sub(x,a,x,X); ++x[0]; mul(y,y,x,Y); } b=y; } luoguP4726P6613 exponentiationpow() 分治法。

次方對半分。

O(𝓜(N)logK)。

公式解。

O(𝓜(N))。

b=aᵏ ln(b)=kln(a) b=exp(kln(a)) //countoftrailingzero intctz(Poly&a) { inti; for(i=0;i<=deg(a);++i) if(a[i]!=0) returni; returni; } voidpow(Poly&b,Poly&a,intk,intB) { if(ctz(a)*k>B){b=0;return;} staticPolyx(N),p(N); for(x=1,p=a;k>0;k>>=1) { if(k&1)mul(x,x,p,B); square(p,p,B); } b=x; } voidpow(Poly&b,Poly&a,chark[],intB) { intz=ctz(a); intzk=z*k; if(zk>B){b=0;return;} intA=min(deg(a),B); deg(b)=A; //提項 intrec_az=rec(a[z]); for(inti=0;i<=A;i++) b[i]=(i+z<=A)?mul(a[i+z],rec_az):0; deg(b)-=z; //a₀=1 log(b,b,B-zk); for(inti=0;i<=B-zk;i++) b[i]=mul(b[i],k1); exp(b,b,B-zk); //乘回去 intpow_az=pow(a[z],k); //pow(a[z],k%phi(k)) for(inti=B;i>=0;i--) b[i]=(i-zk>=0)?mul(b[i-zk],pow_az):0; deg(b)+=zk; } luoguP5245P5273 squarerootsqrt() 遞推法。

HenselLifting。

O(𝓜(N))。

公式解。

O(𝓜(N))。

b≡sqrt(a)(modxⁿ) b²≡a(modxⁿ) b²-a≡0(modxⁿ) letf(b)=b²-a,f′(b)≡2b bꜛ≡b-f(b)/f′(b)(modx²ⁿ) bꜛ≡b-(b²-a)/(2b)(modx²ⁿ) bꜛ≡b-b/2+a/b/2(modx²ⁿ) bꜛ≡b/2+a/b/2(modx²ⁿ) bꜛ≡(b+a/b)/2(modx²ⁿ) b=a½ ln(b)=½ln(a) b=exp(½ln(a)) constintrec2=rec(2); voidsqrt(Poly&b,Poly&a,intB) { y[0]=sqrt(a[0]); for(intn=1;n<=B<<1;n<<=1) { intX=min(n-1,deg(a)); intY=min(n-1,B); staticPolyrec_y(N),x(N),y(N); rec(rec_y,y,Y); copy(x,a,X); mul(x,x,rec_y,Y); deg(y)=Y; for(inti=0;i<=Y;i++) y[i]=mul(x[i]+y[i],rec2); } b=y; } luoguP5205P5277 trigonometricfunctionssin()cos()tan() 公式解。

O(𝓜(N))。

exp(𝑖x)=cos(x)+𝑖sin(x)—→cos(x)=(exp(𝑖x)+exp(-𝑖x))/2 exp(-𝑖x)=cos(x)-𝑖sin(x)sin(x)=(exp(𝑖x)-exp(-𝑖x))/2𝑖 arcsin(x)=∫(x′/sqrt(1-x²))dx arccos(x)=-∫(x′/sqrt(1-x²))dx arctan(x)=∫(x′/(1+x²))dx constintm=998244353; constintg=3; constintI=pow(g,div(m-1,4)); constintrec_2I=rec(I*2); constintrec_2=rec(2); voidsin(Poly&b,Poly&a,intB) { staticPolyx(N),y(N); for(inti=0;i<=B;i++)b[i]=mul(a[i],I); exp(x,b,B); for(inti=0;i<=B;i++)b[i]=-b[i]; exp(y,b,B); for(inti=0;i<=B;i++)b[i]=mul(x[i]-y[i],rec_2I); } voidcos(Poly&b,Poly&a,intB) { staticPolyx(N),y(N); for(inti=0;i<=B;i++)b[i]=mul(a[i],I); exp(x,b,B); for(inti=0;i<=B;i++)b[i]=-b[i]; exp(y,b,B); for(inti=0;i<=B;i++)b[i]=mul(x[i]+y[i],rec_2); } voidtan(Poly&b,Poly&a,intB) { staticPolyx(N),y(N); sin(x,a,B); cos(y,a,B); rec(b,y,B); mul(b,x,b,B); } voidcot(Poly&b,Poly&a,intB) { staticPolyx(N),y(N); sin(x,a,B); cos(y,a,B); rec(b,x,B); mul(b,y,b,B); } voidasin(Poly&b,Poly&a,intB) { Poly&diff_a=b; diff(diff_a,a,B); staticPolyx(N); mul(x,a,a,B); for(inti=0;i<=B;++i)x[i]=-x[i]; ++x[0]; sqrt(x,x,B); rec(x,x,B); mul(x,diff_a,x,B-1); inte(b,x,B); } voidacos(Poly&b,Poly&a,intB) { asin(b,a,B); for(inti=0;i<=B;++i)b[i]=-b[i]; } voidatan(Poly&b,Poly&a,intB) { Poly&diff_a=b; diff(diff_a,a,B); staticPolyx(N); mul(x,a,a,B); ++x[0]; rec(x,x,B); mul(x,diff_a,x,B-1); inte(b,x,B); } luoguP5264P5265 high-orderfinitedifferenceΔᵏ 重複K次。

O(NK)。

乘以(1-x)ᵏ。

次方運算。

O(𝓜(N)logK)。

乘以(1-x)ᵏ。

自然對數指數，避開次方運算。

O(𝓜(N))。

乘以(1-x)ᵏ。

二項式展開，遞推法求係數。

O(𝓜(N))。

finitedifference:|cumulativesum: b=a(1-x)ᵏ|b=a/(1-x)ᵏ b=aexp(kln(1-x))|b=a/exp(kln(1-x)) voidfin(Poly&b,Poly&a,intk,intB) { deg(b)=B; b[0]=1; for(inti=1;i<=B;++i) b[i]=div(mul(b[i-1],-k+i-1),i); mul(b,a,b,B); } voidcum(Poly&b,Poly&a,intk,intB) { deg(b)=B; b[0]=1; for(inti=1;i<=B;++i) b[i]=div(mul(b[i-1],k+i-1),i); mul(b,a,b,B); } luoguP5488 Polynomial(Evaluation) multipointevaluation 遞推法。

Horner'sRule。

O(NK)。

遞增法。

求值a(xᵢ)化作求餘數amod(x-xᵢ)。

O(𝓜(N)K)。

分治法。

項數對半分、形成二元樹。

O(𝓜(N)+𝓜(K)logK)。

一、一次多項式連乘M(x₁,x₂,...)=(x-x₁)(x-x₂)...(x-xK)。

項數對半分M(all)=M(low)M(high)。

形成二元樹，從葉往根乘。

二、求餘數，一回合擴展成logK回合。

第一回合，模數是全部連乘M(all)。

第二回合，模數是前半部M(low)、後半部M(high)。

最後回合，模數是(x-xᵢ)。

形成二元樹，從根往葉模。

加速技巧。

Tellegen'sPrinciple。

O(𝓜(N)+𝓜(K)logK)。

constintK=100; //size voideval(inty[],Poly&a,intx[],intN) { staticPolyt[K<<1]; intTN=(N<<1)-1; for(inti=TN;i>=1;i--) if(i>=N) { t[i].resize(1); deg(t[i])=1; t[i][0]=neg(x[i-N]); t[i][1]=1; } else { intl=i<<1; //leftchild intr=l|1; //rightchild t[i].resize(deg(t[l])+deg(t[r])); mul(t[i],t[l],t[r],INT_MAX); } t[0]=a; for(inti=1;i<=TN;i++) { intp=i>>1; //parent mod(t[i],t[p],t[i]); if(i>=N)y[i-N]=t[i][0]; } } luoguP5050 interpolation 基礎知識請見本站文件「PolynomialInterpolation」。

遞增法。

LagrangeInterpolation。

O(𝓜(N)N)。

一、一次多項式連乘M(x₁,x₂,...)=(x-x₁)(x-x₂)...(x-xK)。

二、內插P(x₁,x₂,...)化作求商數M(x₁,x₂,...)/(x-xᵢ)。

分治法。

項數對半分、形成二元樹。

O(𝓜(N)logN)。

一、一次多項式連乘：M(all)=M(low)M(high)。

二、係數：L'Hôpital'sRule。

三、內插：P(all)=P(low)M(high)+P(high)M(low)。

voidinte(Poly&a,intx[],inty[],intN) //size { //multiplicationtree staticPolyt[::N<<1]; intTN=(N<<1)-1; for(inti=TN;i>=1;i--) if(i>=N) { t[i].resize(deg(t[i])=1); t[i][0]=neg(x[i-N]); t[i][1]=1; } else { intl=i<<1; intr=l|1; t[i].resize(deg(t[l])+deg(t[r])); mul(t[i],t[l],t[r],INT_MAX); } //callstack staticPolyr[lg2N+2],p[lg2N+2],b(N); intL=bit_width((unsignedint)TN); // intL=32-__builtin_clz(TN); //differentiationofmultiplication r[0]=t[1]; diff(r[0],r[0],deg(r[0])-1); //interpolation for(intn=0;n>1)+n,level=L; if(index>TN)index-=N,level-=1; //multipointevaluation intz=__builtin_ctz(index); inti,l; for(i=index>>z,l=level-z;i<=TN;i<<=1,l++) { r[l].resize(deg(t[i])-1); mod(r[l],r[l-1],t[i]); } // y[index-N]=r[level][0]; //evaluationtest //divide-and-conquermethod a=div(y[index-N],r[level][0]); for(i=index,l=level;(i&1)&&(i>1);i>>=1,l--) { mul(a,t[i^1],a,INT_MAX); mul(b,t[i],p[l],INT_MAX); add(a,a,b,INT_MAX); } p[l].resize(deg(a)); //savememory copy(p[l],a,deg(a)); } } luoguP5158P4781 moniclinearcongruencesf(x)≡rᵢ(x)(modmᵢ(x)) 中國餘數定理本質就是內插。

分治法。

式子對半分。

O(𝓜(N)+𝓜(K)logK)。

arithmetic 四則運算（除法不能有餘數）、循環版本四則運算（循環除法不會有餘數），擁有容易理解的解釋方式。

多點求值、點對點四則運算、內插。

O(𝓜(N)logN)。

多點求值用複數單位根，內插用快速傅立葉轉換。

O(𝓜(N))。

Polynomial(Resolution) root 遞推法。

Bairstow'sMethod。

複數共軛根形成實數二次因式。

原式除以二次因式得到一次餘式。

令一次餘式是零。

牛頓法求得二次因式的兩個係數。

時間複雜度O(N²T)，收斂速度是二次方。

遞推法。

Laguerre'sMethod。

原式是一次因式連乘。

一次導數是一次因式導數連加。

二次導數是一次因式平方連減。

以此推導遞推公式。

時間複雜度O(N²T)，收斂速度是三次方。

遞推法。

CompanionMatrix。

多項式求根＝CompanionMatrix特徵多項式求根＝CompanionMatrix求特徵值。

請見最後章節。

時間複雜度O(N²T)，不必猜測初始值。

MATLAB採用此方法。

equation 多項式方程式＝多項式求根。

遞迴加權總和【數學家稱作Ideal】 遞迴加權總和：給定集合，窮舉各種權重組合，得到各種加權總和。

A={a₁,a₂,a₃} a₄=w₁a₁+w₂a₂+w₃a₃ a₅=v₁a₁+v₂a₂+v₃a₃+v₄a₄ a₆=u₁a₁+u₂a₂+u₃a₃+u₄a₄+u₅a₅ :: 數學式子展開之後，其實等同於沒有遞迴的加權總和。

span(A)=k₁a₁+k₂a₂+k₃a₃ 多項式的遞迴加權總和：沒有差別。

A={x³-2xy,2x²y-y²+x,x²+1} a₄=2y(x³-2xy)+x(2x²y-y²+x)+0(x²+1)=-3xy²-x² a₅=y(x³-2xy)+0(2x²y-y²+x)+xy(x²+1)+(x-y+1)(-3xy²-x²) :: 遞迴餘數【尚無正式名稱】 遞迴餘數：給定模數集合，窮舉各種模數組合，得到各種餘數。

M={m₁,m₂,m₃} a⤳b(modm₁) b⤳c(modm₃)—→a⤳r(rmodM) c⤳d(modm₁) :: q⤳r(modm₂) 遞迴餘數的性質，宛如一般餘數的性質。

a₁⤳r₁(rmodM) a₂⤳r₂(rmodM) a₁+a₂⤳r₁+r₂(rmodM) 多項式的遞迴餘數：不斷消滅首項，度數越來越小，最終無法再模。

項有多種排序方式。

固定使用同一種，哪一種都可以。

字典順序：先以字母排序，再以度數排序。

（部分變數視作常數） x²y²+x²y+x²+xy²+xy+x+y²+y+1 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ x²x1 度數順序：先以度數排序，再以字母排序。

x²y²+x²y+xy²+x²+xy+y²+x+y+1 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 43210 Gröbnerbasis 給定被除數集合，實施遞迴加權總和；找到除數集合，可以遞迴整除每一個加權總和。

這樣的除數集合，稱作GröbnerBasis。

MisGröbnerBasisofA: aᵢ⤳0(rmodM)forallaᵢinspan(A) 遞迴加權總和，總共無限多個，無法一一檢查是否整除。

必須取巧。

遞推法。

Buchberger'sAlgorithm。

被除數集合納入除數集合。

窮舉兩兩除數，計算S-polynomial，一一檢查是否整除。

如果都被整除：當前除數集合就是GröbnerBasis。

如果無法整除：餘數納入除數集合，重頭檢查。

時間複雜度是指數時間。

遞推法。

M4GB。

目前最快的演算法。

lcm(lm(a),lm(b))lcm(lm(a),lm(b)) s(a,b)=————————————————a-————————————————b lm(a)lm(b) lm():leadingmonomial lcm():leastcommonmultiple a=x³-2xy b=2x²y-y²+x lm(a)=x³ lm(b)=2x²y lcm(lm(a),lm(b))=2x³y s(a,b)=(2x³y/x³)a-(2x³y/2x²y)b=-3xy²-x² equations GröbnerBasis。

我沒有研究。

Polynomial(Divisor) greatestcommondivisor 多項式搭配不同數域，多項式最大公因數產生不同定義。

一、有理數多項式：輾轉相除法，讓商數的首項係數保持是一。

「首一多項式MonicPolynomial」：首項係數為一。

二、整數多項式：係數事先提出最大公因數，以保證唯一解。

輾轉相除法，讓被除數乘上足夠倍率，以保證商數是整數。

「本原多項式PrimitivePolynomial」：係數最大公因數為一。

a=72x³-12x²-6x¹-12 b=-6x³+22x²-24x¹+8 |rationalpolynomial|content|primitivepart ---------|------------------------|--------|----------------- a|72(x-⅔)(x²+½x+¼)|1|72(x-⅔)(x²+½x+¼) b|-6(x-⅔)(x-1)(x-2)|1|-6(x-⅔)(x-1)(x-2) gcd(a,b)|(x-⅔)|1|(x-⅔) |integerpolynomial|content|primitivepart ---------|------------------------|--------|----------------- a|(2)(3)(3x-2)(4x²+2x+1)|(2)(3)|(3x-2)(4x²+2x+1) b|(-1)(2)(3x-2)(x-1)(x-2)|(-1)(2)|(3x-2)(x-1)(x-2) gcd(a,b)|(2)(3x-2)|(2)|(3x-2) a=a₀x⁰+a₁x¹+a₂x²+...polynomial a=cont(a)pp(a)integerpolynomial cont(a)=gcd(a₀,a₁,a₂,...)content pp(a)=a/cont(a)primitivepart cont(gcd(a,b))=gcd(cont(a),cont(b))commutativity pp(gcd(a,b))=gcd(pp(a),pp(b))commutativity gcd(a,b)=cont(gcd(a,b))pp(gcd(a,b))definition gcd(a,b)=gcd(cont(a),cont(b))gcd(pp(a),pp(b))algorithm alc(b)^(deg(a)-deg(b)+1)=bq+rpseudodivision pquo(a,b)=qpseudoquotient prem(a,b)=rpseudoremainder 遞歸法。

輾轉相除法。

O(N²)。

加速技巧。

Lehmer'sGCDAlgorithm。

O(N²)。

分治法。

Half-GCDAlgorithm。

位數對半分。

O(𝓜(N)logN)。

tuple(int,Poly)uf(Polya) { intcont=0; for(inti=0;i<=deg(a);i++)cont=gcd(cont,a[i]); Polypp=a; for(inti=0;i<=deg(a);i++)pp[i]=a[i]/cont; return{cont,pp}; } Polyprem(Polya,Polyb) { intd=pow(b[deg(b)],deg(a)-deg(b)+1); returnrem(a*d,b); } Polypgcd(Polya,Polyb) { Polyr; while(b!=0) { r=pp(prem(a,b)); a=b; b=r; } returnr; } Polygcd(Polya,Polyb) { auto[cont_a,pp_a]=uf(a); auto[cont_b,pp_b]=uf(b); returngcd(cont_a,cont_b)*pgcd(pp_a,pp_b); } factorization 回溯法。

Kronecker'sAlgorithm。

兩個多項式整除，所有對應點也會整除。

原式的函數點的某個因數，是因式的函數點。

遞迴窮舉因式的函數點，以「NewtonInterpolation」將函數點變成多項式；試除原式，如果整除，就是因式。

時間複雜度是指數時間。

貪心法。

LLLAlgorithm。

我沒有研究。

時間複雜度是多項式時間。

「不可約多項式IrreduciblePolynomial」：除了自己和1以外，沒有其他因式。

質式。

luoguP4506 square-freefactorization 擷取指數是一的因式連乘積。

換句話說，捨棄指數是平方、立方、……的因式。

注意到，只求得因式連乘積，不求得各個因式。

a(x)=(1+x)¹(1-2x+2x²)¹(1-x)²(3-2x)²(2-x+x²)²x³(2x-1)³ ^^^^^^^^^^^^^^^^^ square-freefactor 甚至更進一步，一邊擷取指數是一的因式連乘積，一邊實施微分降低指數，求得各種指數的因式連乘積。

換句話說，所有因式，依照指數，進行分類。

注意到，只求得因式連乘積暨指數，不求得各個因式。

a(x)=(1+x)¹(1-2x+2x²)¹(1-x)²(3-2x)²(2-x+x²)²x³(2x-1)³ ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 遞推法。

gcd(a,a′)可以降低指數，而且不影響係數。

時間複雜度是多項式時間。

a=a₁¹a₂²a₃³...各種指數的因式連乘積 a′=(1a₁¹⁻¹a₁′)(a₂²a₃³...)微分 +(2a₂²⁻¹a₂′)(a₁¹a₃³...) +(3a₃³⁻¹a₃′)(a₁¹a₂²...) +...... b=gcd(a,a′)=a₁¹⁻¹a₂²⁻¹a₃³⁻¹...次方減一 c=a/b=a₁a₂a₃...次方變一 d=gcd(b,c)=a₂a₃...得到指數大於一的因式連乘積 f=c/d=a₁得到指數為一的因式連乘積！ b/d=gcd(b,b′)=a₂²⁻²a₃³⁻²...開始第二回合，次方再減一 ::: 遞推法。

Yun'sAlgorithm。

過程有一點類似正規正交化。

時間複雜度是多項式時間。

a=a₁¹a₂²a₃³...各種指數的因式連乘積 a′=(1a₁¹⁻¹a₁′)(a₂²a₃³...)微分 +(2a₂²⁻¹a₂′)(a₁¹a₃³...) +(3a₃³⁻¹a₃′)(a₁¹a₂²...) +...... b=gcd(a,a′)=a₁¹⁻¹a₂²⁻¹a₃³⁻¹...次方減一 c=a/b=a₁a₂a₃...次方消失 p=a′/b=(1a₁′)(a₂a₃...) +(2a₂′)(a₁a₃...) +(3a₃′)(a₁a₂...) +...... q=(a/b)′=(a₁′)(a₂a₃...) +(a₂′)(a₁a₃...) +(a₃′)(a₁a₂...) +...... r=p-q=(1a₂′)(a₁a₃...) +(2a₃′)(a₁a₂...) +...... gcd(c,r)=a₁ algorithm=Gram–Schmidtorthonormalization square-free=linear dotproduct=gcd derivative=shearing basis=factorwithincrementalexponent 「無平方多項式Square-freePolynomial」：不存在因式，其指數是二次（以上）。

注意到，是指數（重複次數），不是度數（最高次方值）。

primitiverootofunity 「分圓多項式CyclotomicPolynomial」：單位原根連乘積。

遞推法。

利用N的因數。

O(𝓜(N)σ₀(N))≈O(𝓜(N)sqrtN)。

(xⁿ-1)=proΦd(x) d|n Φₙ(x)=(xⁿ-1)/proΦd(x) d|n d≠n 遞推法。

利用N的質因數。

O(𝓜(N)ω(N))≈O(𝓜(N)logN)。

n=p₁e₁×p₂e₂×...×pₖeₖ f₀(x)=(x-1) fᵢ(x)=fᵢ₋₁(xpᵢ)/fᵢ₋₁(x) Φₙ(x)=fₖ(xn/(pᵢ×...×pₖ)) vectordivisor(intn) { vectorv; intd; for(d=1;d*dcyclotomic_polynomial(intn) //deg { vectordivisor=::divisor(n); intN=divisor.size(); vectorcyc(N); Polya(n); for(inti=0;i luoguP1520 Polynomial(Recurrence) evaluation 已知線性遞迴函數共N項，求數列第K項。

「特徵多項式CharacteristicPolynomial」：線性遞迴函數，移項到等號同側，再變成生成函數。

動態規劃。

依序計算每一項。

O(NK)。

公式解。

Fiduccia'sAlgorithm。

不斷展開最高項＝長除法除以特徵多項式。

xᴷ模特徵多項式，求餘數。

O(𝓜(K))。

加速技巧。

一邊倍增、一邊模。

O(𝓜(N)logK)。

加速技巧。

數列第K項化作分式第K項。

以對稱函數進行擴分，消滅奇數項，傅立葉轉換少做一半長度。

O(𝓜(N)logK)。

linearrecursivefunction(degreeN) c₀f(n)+c₁f(n+1)+c₂f(n+2)+...+cN-1f(n+(N-1))=f(n+N) generatingfunction c₀xⁿ+c₁xⁿ⁺¹+c₂xⁿ⁺²+...+cN-1xn+(N-1)=xⁿ⁺ᴺ characteristicpolynomial -c₀x⁰+-c₁x¹+-c₂x²+...+-cN-1xN-1+xᴺ expansionofleadingmonomial←—→longdivisionbycharacteristicpolynomial F₅=F₄+F₃ F₅=(F₃+F₂)+F₃=2F₃+F₂←—→x⁵mod(x²-x¹-x⁰) F₅=2(F₂+F₁)+F₂=3F₂+2F₁ intreval(inta[],intA,ints[],intk) { if(k>=1) { if(k&1) { //r=r*p%c; mul(r,r,p,INT_MAX); mod(r,r,c,rec_rc); } //p=p*p%c; square(p,p,INT_MAX); mod(p,p,c,rec_rc); } //evaluation inty=0; for(inti=0;i<=deg(r);i++) y+=mul(r[i],s[i]); returny; } voidmod(Poly&r,Poly&a,Poly&b,Poly&rec_rb) { staticPolyra(N),q(N); intA=deg(a),B=deg(b); if(Adeg(r(x)) voidrinte(inta[],int&A,ints[],intK) { staticPolyc(::K),b(::K),t(::K); c=1; //reversecharacteristicpolynomial b=1; //prev //Newtoninterpolation intL=0,shift=1,prev_diff=1; for(intk=0;kk){shift++;continue;} //speedup L=k+1-L;shift=1;prev_diff=diff; b=t; } //recurrencecoefficients A=deg(c); for(inti=0;i luoguP5487UVa1539 Polynomial(Integer) 整數運算→多項式運算 整數可以視作多項式。

有些整數運算可以套用多項式運算。

比方來說，整數乘法套用多項式乘法。

因為3×37=111，所以3(3x+7)=9x+21，底數x=10。

記得認真考慮整數運算的時間複雜度。

快速傅立葉轉換內含整數乘法，遞迴使用整數乘法，時間複雜度成為O(NlogNloglogN)。

順帶一提，整數乘法的時間複雜度下限仍是懸案。

2019年出現新演算法，時間複雜度降至O(NlogN)，很可能就是下限了。

luoguP1601P1919P5432 整數運算↛多項式運算 有些整數運算則不可以套用多項式運算。

比方來說，整數分解套用多項式分解，雖然111=3×37，但是x²+x+1≠3(3x+7)。

我們不知道底數x應該是多少，也不知道是否有減法。

順帶一提，整數分解的時間複雜度類別仍是懸案，目前還不確定它究竟是P問題或是NP-complete問題，相當特別。

Polynomial(Matrix) 多項式運算→矩陣運算 有些多項式運算可以套用矩陣運算。

convolution乘法（卷積）⇒Toeplitzmatrix常對角矩陣求值 circularconvolution循環乘法　　⇔circulantmatrix循環矩陣求值 multipointevaluation多點求值　　⇔Vandermondematrix范德蒙矩陣求值 multipointevaluation多點求值特例⇔Fouriermatrix傅立葉矩陣求值 coprime互質　　　　⇔Sylvestermatrix結式矩陣行列式 root根　　　　　⇔companionmatrix同伴矩陣特徵值 linearrecurrence線性遞迴函數⇔companionmatrix同伴矩陣轉置 然而矩陣運算時間複雜度較高，缺乏實用價值。

唯一例外：多項式求根的演算法，利用了矩陣運算，目前沒有更快的演算法。

Multiplication⇨ToeplitzMatrix 「常對角矩陣」。

每一條左上右下斜線是相同元素的矩陣。

[1532] [2153] [7215] [0721] [8072] [9807] 多項式乘法（數列卷積）化作常對角矩陣求值。

反之則不然。

[0] [876000000][0] [087600000][1] (12345)[008760000][2] ∗———>[000876000][3] (678)[000087600][4] [000008760][5] [000000876][0] [0] convolutionToeplitzmatrixevaluation [60000] [76000][1] (12345)[87600][2] ∗———>[08760][3] (678)[00876][4] [00087][5] [00008] convolutionToeplitzmatrixevaluation 加法。

常對角矩陣只有2N-1種數字。

時間複雜度O(N)。

求值。

填充元素成為循環矩陣。

時間複雜度O(NlogN)。

[15372][1] [153][1][21537][2] [215][2]———>[72153][3] [721][3][37215][0] [53721][0] Toeplitzmatrixcirculantmatrix 求解。

時間複雜度O(N²)甚至O(Nlog²N)。

乘法。

時間複雜度O(N²)。

反矩陣。

時間複雜度O(N²)。

CircularMultiplication⇨CirculantMatrix 「循環矩陣」。

直向橫向皆循環的方陣。

[59876] [65987] [76598] [87659] [98765] 多項式循環乘法（數列循環卷積）＝循環矩陣求值。

[59876][1] (12345)[65987][2] ⊛[76598][3] (56789)[87659][4] [98765][5] circularconvolutioncirculantmatrixevaluation 求值、求解、乘法、反矩陣，都是O(NlogN)。

演算法是傅立葉轉換、數論轉換。

MultipointEvaluation⇨VandermondeMatrix 「范德蒙矩陣」。

橫條是任意數的0次方到N-1次方。

[9⁰9¹9²9³9⁴] [0⁰0¹0²0³0⁴] [3⁰3¹3²3³3⁴] [5⁰5¹5²5³5⁴] [2⁰2¹2²2³2⁴] 多項式多點求值＝范德蒙矩陣求值。

時間複雜度O(NK)。

多項式內插＝范德蒙矩陣求解。

時間複雜度O(N³)。

f(x)=c₀x⁰+c₁x¹+c₂x²+c₃x³+c₄x⁴ [x₀⁰x₀¹x₀²x₀³x₀⁴][c₀][f(x₀)] [x₁⁰x₁¹x₁²x₁³x₁⁴][c₁][f(x₁)] [x₂⁰x₂¹x₂²x₂³x₂⁴][c₂]=[f(x₂)] [x₃⁰x₃¹x₃²x₃³x₃⁴][c₃][f(x₃)] [x₄⁰x₄¹x₄²x₄³x₄⁴][c₄][f(x₄)] MultipointEvaluation⇨FourierMatrix 「傅立葉矩陣」可以視作范德蒙矩陣的特例：輸入們恰是單位根（複數平面單位圓的等分點），而且是方陣K=N。

f(x)=c₀x⁰+c₁x¹+c₂x²+c₃x³+c₄x⁴ [x₀⁰x₀¹x₀²x₀³x₀⁴][c₀][f(x₀)] [x₁⁰x₁¹x₁²x₁³x₁⁴][c₁][f(x₁)] [x₂⁰x₂¹x₂²x₂³x₂⁴][c₂]=[f(x₂)] [x₃⁰x₃¹x₃²x₃³x₃⁴][c₃][f(x₃)] [x₄⁰x₄¹x₄²x₄³x₄⁴][c₄][f(x₄)] wherexᵢ=e-𝑖(2π/N)⋅i Coprime⇨SylvesterMatrix 「西爾維斯特矩陣」或「結式矩陣」。

多項式係數逐步向右移位，其餘元素是0，方陣。

多項式互質＝結式矩陣行列式非零。

時間複雜度O(N³)。

a(x)=a₀x⁰+a₁x¹+...+aₘxᵐ b(x)=b₀x⁰+b₁x¹+...+bₙxⁿ [a₀...aₘ]\ [a₀...aₘ]nrows [:]| S=[a₀...aₘ]/ [b₀...bₙ]\ [b₀...bₙ]mrows [:]| [b₀...bₙ]/ 0.defineresultantres(a,b)=det(S) 1.gcd(a,b)isnotconstantpolynomialiffres(a,b)=det(S)=0 2.deg(gcd(a,b))=m+n-rank(S) 3.existxandysuchthatax+by=res(a,b)=det(S) wheredeg(x)