// Find the continued fraction representation of P/Q.autofraction(intp,intq){std::vector<int>a;while(q){a.push_back(p/q);std::tie(p,q)=std::make_pair(q,p%q);}returna;}
1234567
# Find the continued fraction representation of P/Q.deffraction(p,q):a=[]whileq:a.append(p//q)p,q=q,p%qreturna
// Find the convergents of a continued fraction A.// Numerators and denominators stored separately in P and Q.autoconvergents(std::vector<int>a){std::vector<int>p={0,1};std::vector<int>q={1,0};for(autoit:a){p.push_back(p.back()*it+p.end()[-2]);q.push_back(q.back()*it+q.end()[-2]);}returnstd::make_pair(p,q);}
123456789
# Find the convergents of a continued fraction A.# Numerators and denominators stored separately in P and Q.defconvergents(a):p=[0,1]q=[1,0]foritina:p.append(p[-1]*it+p[-2])q.append(q[-1]*it+q[-2])returnp,q
// Return (x,y) such that Ax+By=C.// Assume that such (x,y) exists.autodio(intA,intB,intC){std::vector<int>p,q;std::tie(p,q)=convergents(fraction(A,B));C/=A/p.back();intt=p.size()%2?-1:1;returnstd::make_pair(t*C*q.end()[-2],-t*C*p.end()[-2]);}
1234567
# Return (x, y) such that Ax+By=C.# Assume that such (x, y) exists.defdio(A,B,C):p,q=convergents(fraction(A,B))C//=A//p[-1]# divide by gcd(A, B)t=(-1)iflen(p)%2else1returnt*C*q[-2],-t*C*p[-2]
// Expand [..., n] to [..., n-1, 1] if needed.voidexpand(std::vector<int>&a){if(a.back()>1){--a.back();a.push_back(1);}}// Check if a is smaller than b.boolless_than(std::vector<int>a,std::vector<int>b){expand(a);expand(b);for(inti=0;i<a.size()-1||i<b.size()-1;++i){intd=(i<a.size()-1?a[i]:0)-(i<b.size()-1?b[i]:0);if(i&1)d=-d;if(d<0){returntrue;}elseif(d>0){returnfalse;}}returnfalse;}
1 2 3 4 5 6 7 8 9101112131415
# Expand [..., n] to [..., n-1, 1] if needed.defexpand(a):ifa[-1]!=1:a[-1]-=1a.append(1)returna# Check if a is smaller than b.defless_than(a,b):a=expand(a)b=expand(b)a=[(-1)**i*a[i]foriinrange(len(a))]b=[(-1)**i*b[i]foriinrange(len(b))]returna<b
#include<algorithm>#include<iostream>#include<tuple>#include<vector>constexprintM=1e9+7;// FLTs. Essentially 2x2 matrix.structFracLinearTrans{intmat[4];FracLinearTrans():mat{}{}FracLinearTrans(intx):mat{x,1,1,0}{}FracLinearTrans(inta,intb,intc,intd):mat{a,b,c,d}{}FracLinearTransoperator*(constFracLinearTrans&rhs)const{returnFracLinearTrans(((longlong)mat[0]*rhs.mat[0]+(longlong)mat[1]*rhs.mat[2])%M,((longlong)mat[0]*rhs.mat[1]+(longlong)mat[1]*rhs.mat[3])%M,((longlong)mat[2]*rhs.mat[0]+(longlong)mat[3]*rhs.mat[2])%M,((longlong)mat[2]*rhs.mat[1]+(longlong)mat[3]*rhs.mat[3])%M);}FracLinearTransinv()const{returnFracLinearTrans(mat[3],M-mat[1],M-mat[2],mat[0]);}};intmain(){intn,q;std::cin>>n>>q;// Get prefix sum of FLTs.std::vector<FracLinearTrans>ps(1,{1,0,0,1});ps.reserve(n+1);for(inti=1;i<=n;++i){inta;std::cin>>a;ps[i]=ps[i-1]*FracLinearTrans(a);}// Query.for(;q;--q){intl,r;std::cin>>l>>r;// Difference.autores=ps[l-1].inv()*ps[r];intu=res.mat[0],d=res.mat[2];// Correct signs.if(!(l&1)){if(u)u=M-u;if(d)d=M-d;}printf("%d %d\n",u,d);}return0;}
# PYTHON IS TOO SLOW TO PASS THIS PROBLEM.# JUST FOR REFERENCE.M=10**9+7defmul(a,b):return((a[0]*b[0]+a[1]*b[2])%M,(a[0]*b[1]+a[1]*b[3])%M,(a[2]*b[0]+a[3]*b[2])%M,(a[2]*b[1]+a[3]*b[3])%M,)definv(a):return(a[3],M-a[1],M-a[2],a[0])n,q=map(int,input().split())ps=[(1,0,0,1)]# Get presum.forainmap(int,input().split()):ps.append(mul(ps[-1],(a,1,1,0)))for_inrange(q):l,r=map(int,input().split())res=mul(inv(ps[l-1]),ps[r])u,d=res[0],res[2]ifl%2==0:ifu:u=M-uifd:d=M-dprint(u,d)
// Return the continued fraction and minimal positive period// of a quadratic irrational (x + y * sqrt(n)) / z.autoquadratic_irrational(intx,inty,intz,intn){intp=x*z;intd=n*y*y*z*z;intq=z*z;intdd=(int)sqrt(n)*y*z;inti=0;std::vector<int>a;std::unordered_map<size_t,int>used;while(!used.count(((1LL<<32)*p)|q)){a.emplace_back((p+dd)/q);used[((1LL<<32)*p)|q]=i++;p=a.back()*q-p;q=(d-p*p)/q;}returnstd::make_pair(a,i-used[((1LL<<32)*p)|q]);}
1 2 3 4 5 6 7 8 91011121314151617
# Return the continued fraction and minimal positive period# of a quadratic irrational (x + y * sqrt(n)) / z.defquadratic_irrational(x,y,z,n):p=x*zd=n*y*y*z*zq=z*zdd=floor(sqrt(n))*y*zi=0a=[]used=dict()while(p,q)notinused:a.append((p+dd)//q)used[p,q]=ip=a[-1]*q-pq=(d-p*p)//qi+=1returna,i-used[p,q]
#include<algorithm>#include<cmath>#include<iostream>#include<tuple>#include<unordered_map>#include<vector>// Return the continued fraction and minimal positive period// of a quadratic irrational (x + y * sqrt(n)) / z.autoquadratic_irrational(intx,inty,intz,intn){intp=x*z;intd=n*y*y*z*z;intq=z*z;intdd=(int)sqrt(n)*y*z;inti=0;std::vector<int>a;std::unordered_map<size_t,int>used;while(!used.count(((1LL<<32)*p)|q)){a.emplace_back((p+dd)/q);used[((1LL<<32)*p)|q]=i++;p=a.back()*q-p;q=(d-p*p)/q;}returnstd::make_pair(a,i-used[((1LL<<32)*p)|q]);}// Fractional Linear Transformation.structFracLinearTrans{staticconstexprintM=1e9+7;intmat[4];FracLinearTrans(inta,intb,intc,intd):mat{a,b,c,d}{}FracLinearTransoperator*(constFracLinearTrans&rhs)const{returnFracLinearTrans(((longlong)mat[0]*rhs.mat[0]+(longlong)mat[1]*rhs.mat[2])%M,((longlong)mat[0]*rhs.mat[1]+(longlong)mat[1]*rhs.mat[3])%M,((longlong)mat[2]*rhs.mat[0]+(longlong)mat[3]*rhs.mat[2])%M,((longlong)mat[2]*rhs.mat[1]+(longlong)mat[3]*rhs.mat[3])%M);}};intmain(){intx,k,L;std::cin>>x>>k;std::vector<int>a;std::tie(a,L)=quadratic_irrational(0,1,1,x);// L==a.size()-1 for sqrt(x)FracLinearTranscyc(1,0,0,1);for(inti=a.size()-1;i;--i){cyc=FracLinearTrans(a[i],1,1,0)*cyc;}// 1/0=Inf.FracLinearTransres(0,1,0,0);// Tail terms.for(inti=k%L;i;--i){res=FracLinearTrans(a[i],1,1,0)*res;}// Binary exponentiation.for(intb=k/L;b;b>>=1){if(b&1)res=cyc*res;cyc=cyc*cyc;}// First term.res=FracLinearTrans(a[0],1,1,0)*res;printf("%d/%d",res.mat[1],res.mat[3]);return0;}
// Find [ah, ph, qh] such that points r[i]=(ph[i], qh[i]) constitute// upper convex hull of lattice points on 0 <= x <= N and 0 <= y <= r * x,// where r = [a0, a1, a2, ...] and there are ah[i]-1 integer points on the// segment between r[i] and r[i+1].autohull(std::vector<int>a,intN){std::vector<int>p,q;std::tie(p,q)=convergents(a);intt=N/q.back();std::vector<int>ah={t};std::vector<int>ph={0,t*p.back()};std::vector<int>qh={0,t*q.back()};for(inti=q.size()-1;i;--i){if(i%2){while(qh.back()+q[i-1]<=N){t=(N-qh.back()-q[i-1])/q[i];intdp=p[i-1]+t*p[i];intdq=q[i-1]+t*q[i];intk=(N-qh.back())/dq;ah.push_back(k);ph.push_back(ph.back()+k*dp);qh.push_back(qh.back()+k*dq);}}}returnmake_tuple(ah,ph,qh);}
1 2 3 4 5 6 7 8 9101112131415161718192021
# Find [ah, ph, qh] such that points r[i]=(ph[i], qh[i]) constitute# upper convex hull of lattice points on 0 <= x <= N and 0 <= y <= r * x,# where r = [a0, a1, a2, ...] and there are ah[i]-1 integer points on the# segment between r[i] and r[i+1].defhull(a,N):p,q=convergents(a)t=N//q[-1]ah=[t]ph=[0,t*p[-1]]qh=[0,t*q[-1]]foriinreversed(range(len(q))):ifi%2==1:whileqh[-1]+q[i-1]<=N:t=(N-qh[-1]-q[i-1])//q[i]dp=p[i-1]+t*p[i]dq=q[i-1]+t*q[i]k=(N-qh[-1])//dqah.append(k)ph.append(ph[-1]+k*dp)qh.append(qh[-1]+k*dq)returnah,ph,qh
# Find (x, y) such that y = (A*x+B) // C,# such that Cy - Ax is max and 0 <= x <= N.defclosest(A,B,C,N):# y <= (A*x + B)/C <=> diff(x, y) <= Bdefdiff(x,y):returnC*y-A*xp,q=convergents(fraction(A,C))qh,ph=0,B//Cforiinrange(2,len(q)-1):ifi%2==0:whilediff(qh+q[i+1],ph+p[i+1])<=B:t=1+(diff(qh+q[i-1],ph+p[i-1])-B-1)//(-diff(q[i],p[i]))dp=p[i-1]+t*p[i]dq=q[i-1]+t*q[i]k=(N-qh)//dqifk==0:returnqh,phifdiff(dq,dp)!=0:k=min(k,(B-diff(qh,ph))//diff(dq,dp))qh,ph=qh+k*dq,ph+k*dpreturnqh,phdefsolve(A,B,N):x,y=closest(A,N%A,B,N//A)returnN//A-x,y
// Find sum of floor(k * x) for k in [1, N] and x = [a0; a1, a2, ...]intsum_floor(std::vector<int>a,intN){N++;std::vector<int>ah,ph,qh;std::tie(ah,ph,qh)=hull(a,N);// The number of lattice points within a vertical right trapezoid// on points (0; 0) - (0; y1) - (dx; y2) - (dx; 0) that has// a+1 integer points on the segment (0; y1) - (dx; y2).autopicks=[](inty1,inty2,intdx,inta)->int{intb=y1+y2+a+dx;intA=(y1+y2)*dx;return(A+b)/2-y2;// = (A - b + 2) // 2 + b - (y2 + 1)};intans=0;for(size_ti=1;i<qh.size();i++){ans+=picks(ph[i-1],ph[i],qh[i]-qh[i-1],ah[i-1]);}returnans-N;}
1 2 3 4 5 6 7 8 9101112131415161718
# Find sum of floor(k * x) for k in [1, N] and x = [a0; a1, a2, ...].defsum_floor(a,N):N+=1ah,ph,qh=hull(a,N)# The number of lattice points within a vertical right trapezoid# on points (0; 0) - (0; y1) - (dx; y2) - (dx; 0) that has# a+1 integer points on the segment (0; y1) - (dx; y2) but with# the number of points on the vertical right line excluded.defpicks(y1,y2,dx,a):b=y1+y2+a+dxA=(y1+y2)*dxreturn(A+b)//2-y2# = (A - b + 2) // 2 + b - (y2 + 1)ans=0foriinrange(1,len(qh)):ans+=picks(ph[i-1],ph[i],qh[i]-qh[i-1],ah[i-1])returnans-N
# Find convex hull of lattice (x, y) such that C*y <= A*x+B.defhull(A,B,C,N):defdiff(x,y):returnC*y-A*xa=fraction(A,C)p,q=convergents(a)ah=[]ph=[B//C]qh=[0]definsert(dq,dp):k=(N-qh[-1])//dqifdiff(dq,dp)>0:k=min(k,(B-diff(qh[-1],ph[-1]))//diff(dq,dp))ah.append(k)qh.append(qh[-1]+k*dq)ph.append(ph[-1]+k*dp)foriinrange(1,len(q)-1):ifi%2==0:whilediff(qh[-1]+q[i+1],ph[-1]+p[i+1])<=B:t=(B-diff(qh[-1]+q[i+1],ph[-1]+p[i+1]))//(-diff(q[i],p[i]))dp=p[i+1]-t*p[i]dq=q[i+1]-t*q[i]ifdq<0orqh[-1]+dq>N:breakinsert(dq,dp)insert(q[-1],p[-1])foriinreversed(range(len(q))):ifi%2==1:whileqh[-1]+q[i-1]<=N:t=(N-qh[-1]-q[i-1])//q[i]dp=p[i-1]+t*p[i]dq=q[i-1]+t*q[i]insert(dq,dp)returnah,ph,qhdefsolve(N,M,A,B):ah,ph,qh=hull(A,B,M,N)# The number of lattice points within a vertical right trapezoid# on points (0; 0) - (0; y1) - (dx; y2) - (dx; 0) that has# a+1 integer points on the segment (0; y1) - (dx; y2) but with# the number of points on the vertical right line excluded.defpicks(y1,y2,dx,a):b=y1+y2+a+dxA=(y1+y2)*dxreturn(A+b)//2-y2# = (A - b + 2) // 2 + b - (y2 + 1)ans=0foriinrange(1,len(qh)):ans+=picks(ph[i-1],ph[i],qh[i]-qh[i-1],ah[i-1])returnans-N
// Find Q that minimizes Q*r mod m for 1 <= k <= n < m.intmod_min(intr,intn,intm){autoa=fraction(r,m);std::vector<int>p,q;std::tie(p,q)=convergents(a);for(inti=2;i<q.size();++i){if(i%2==1&&(i+1==q.size()||q[i+1]>n)){intt=(n-q[i-1])/q[i];returnq[i-1]+t*q[i];}}return0;}
123456789
# Find Q that minimizes Q*r mod m for 1 <= k <= n < m.defmod_min(r,n,m):a=fraction(r,m)p,q=convergents(a)foriinrange(2,len(q)):ifi%2==1and(i+1==len(q)orq[i+1]>n):t=(n-q[i-1])//q[i]returnq[i-1]+t*q[i]return0