#include #include #include #include #include #include #include #ifndef CPPLIB_SRC_ALGORITHM_MATH_INTEGER_BIG_INTEGER_HPP_INCLUDED #define CPPLIB_SRC_ALGORITHM_MATH_INTEGER_BIG_INTEGER_HPP_INCLUDED #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #if defined(__x86_64__) && (defined(__GNUC__) || defined(__clang__)) #include #endif #if !defined(__SIZEOF_INT128__) #error "This optimized ExactInteger requires unsigned __int128." #endif namespace exact_integer_detail{templateinline constexpr bool native_integer=std::is_integral_v>||std::same_as,__int128_t>||std::same_as,__uint128_t>;templateconcept NativeInteger=native_integer;templatestruct MakeUnsigned{using type=std::make_unsigned_t;};template<>struct MakeUnsigned<__int128_t>{using type=__uint128_t;};template<>struct MakeUnsigned<__uint128_t>{using type=__uint128_t;};templateusing MakeUnsignedT=typename MakeUnsigned>::type;}class BigInteger;namespace big_integer_detail{struct AddExpression{const BigInteger&left;const BigInteger&right;};struct SubtractExpression{const BigInteger&left;const BigInteger&right;};struct MultiplyExpression{const BigInteger&left;const BigInteger&right;};struct DivideExpression{const BigInteger&left;const BigInteger&right;};struct ModuloExpression{const BigInteger&left;const BigInteger&right;};templateinline constexpr bool is_expression=std::same_as,AddExpression>||std::same_as,SubtractExpression>||std::same_as,MultiplyExpression>||std::same_as,DivideExpression>||std::same_as,ModuloExpression>;templateconcept Expression=is_expression;}class ExactInteger{using Limb=std::uint64_t;using Wide=__uint128_t;static constexpr std::size_t inline_limb_capacity=4;static constexpr Limb decimal_base=10'000'000'000'000'000'000ULL;template[[gnu::always_inline]]static inline Limb add_equal_fixed(Limb*b,const Limb*c,const Limb*e)noexcept{Wide f=0;for(std::size_t g=0;g(c[g])+e[g]+f;b[g]=static_cast(k);f=k>>64;}return static_cast(f);}template[[gnu::always_inline]]static inline Limb subtract_equal_fixed(Limb*l,const Limb*o,const Limb*p)noexcept{Wide q=0;for(std::size_t s=0;s(p[s])+q;const Wide x=o[s];l[s]=static_cast(x-w);q=x(q);}[[gnu::always_inline]]static inline Limb add_equal_length(Limb*y,const Limb*z,const Limb*A,std::size_t B)noexcept{switch(B){case 2:return add_equal_fixed<2>(y,z,A);case 4:return add_equal_fixed<4>(y,z,A);case 8:return add_equal_fixed<8>(y,z,A);default:break;} #if defined(__x86_64__) && (defined(__GNUC__) || defined(__clang__)) unsigned char D;Limb*E=y;const Limb*F=z;const Limb*G=A;std::size_t H=B;__asm__ volatile("clc\n\t""1:\n\t""movq (%[left]), %%rax\n\t""adcq (%[right]), %%rax\n\t""movq %%rax, (%[output])\n\t""leaq 8(%[left]), %[left]\n\t""leaq 8(%[right]), %[right]\n\t""leaq 8(%[output]), %[output]\n\t""decq %[remaining]\n\t""jnz 1b\n\t""setc %[carry]":[output]"+r"(E),[left]"+r"(F),[right]"+r"(G),[remaining]"+r"(H),[carry]"=qm"(D): :"rax","cc","memory");return D; #else Wide I=0;for(std::size_t J=0;J(z[J])+A[J]+I;y[J]=static_cast(K);I=K>>64;}return static_cast(I); #endif }[[gnu::always_inline]]static inline Limb subtract_equal_length(Limb*L,const Limb*M,const Limb*N,std::size_t O)noexcept{switch(O){case 2:return subtract_equal_fixed<2>(L,M,N);case 4:return subtract_equal_fixed<4>(L,M,N);case 8:return subtract_equal_fixed<8>(L,M,N);default:break;} #if defined(__x86_64__) && (defined(__GNUC__) || defined(__clang__)) unsigned char P;Limb*Q=L;const Limb*R=M;const Limb*S=N;std::size_t T=O;__asm__ volatile("clc\n\t""1:\n\t""movq (%[left]), %%rax\n\t""sbbq (%[right]), %%rax\n\t""movq %%rax, (%[output])\n\t""leaq 8(%[left]), %[left]\n\t""leaq 8(%[right]), %[right]\n\t""leaq 8(%[output]), %[output]\n\t""decq %[remaining]\n\t""jnz 1b\n\t""setc %[borrow]":[output]"+r"(Q),[left]"+r"(R),[right]"+r"(S),[remaining]"+r"(T),[borrow]"=qm"(P): :"rax","cc","memory");return P; #else Wide U=0;for(std::size_t V=0;V(N[V])+U;const Wide X=M[V];L[V]=static_cast(X-W);U=X(U); #endif }[[gnu::always_inline]]static inline Limb divide_two_by_one(Limb Y,Limb Z,Limb aa,Limb&ab)noexcept{ #if defined(__x86_64__) && (defined(__GNUC__) || defined(__clang__)) Limb ac;Limb ad;__asm__("divq %[divisor]":"=a"(ac),"=d"(ad):"a"(Z),"d"(Y),[divisor]"r"(aa):"cc");ab=ad;return ac; #else const Wide ae=(static_cast(Y)<<64)|Z;const Limb af=static_cast(ae/aa);ab=static_cast(ae%aa);return af; #endif }class LimbStorage{std::arrayag{};std::unique_ptrah;std::size_t ai=0;std::size_t aj=inline_limb_capacity;Limb*pointer()noexcept{return ah?ah.get():ag.data();}const Limb*pointer()const noexcept{return ah?ah.get():ag.data();}public:LimbStorage()=default;LimbStorage(const LimbStorage&ak){resize(ak.ai);std::copy_n(ak.data(),ai,data());}LimbStorage(LimbStorage&&al)noexcept{if(al.ah){ah=std::move(al.ah);ai=al.ai;aj=al.aj;al.ai=0;al.aj=inline_limb_capacity;}else{ai=al.ai;std::copy_n(al.ag.data(),ai,ag.data());al.ai=0;}}LimbStorage&operator=(const LimbStorage&am){if(this==&am)return*this;resize(am.ai);std::copy_n(am.data(),ai,data());return*this;}LimbStorage&operator=(LimbStorage&&an)noexcept{if(this==&an)return*this;clear();ah.reset();aj=inline_limb_capacity;if(an.ah){ah=std::move(an.ah);ai=an.ai;aj=an.aj;an.ai=0;an.aj=inline_limb_capacity;}else{ai=an.ai;std::copy_n(an.ag.data(),ai,ag.data());an.ai=0;}return*this;}Limb*data()noexcept{return pointer();}const Limb*data()const noexcept{return pointer();}std::size_t size()const noexcept{return ai;}std::size_t capacity()const noexcept{return aj;}bool empty()const noexcept{return ai==0;}Limb&operator[](std::size_t ao)noexcept{return pointer()[ao];}const Limb&operator[](std::size_t ap)const noexcept{return pointer()[ap];}Limb&back()noexcept{return pointer()[ai-1];}const Limb&back()const noexcept{return pointer()[ai-1];}Limb*begin()noexcept{return data();}Limb*end()noexcept{return data()+ai;}const Limb*begin()const noexcept{return data();}const Limb*end()const noexcept{return data()+ai;}void reserve(std::size_t aq){if(aq<=aj)return;std::size_t ar=aj;while(ar(std::numeric_limits::max)()/2){ar=aq;break;}ar*=2;}auto as=std::make_unique_for_overwrite(ar);std::copy_n(data(),ai,as.get());ah=std::move(as);aj=ar;}void resize(std::size_t at,Limb au=0){reserve(at);if(at>ai){std::fill(data()+ai,data()+at,au);}ai=at;}void assign(std::size_t av,Limb aw){resize(av);std::fill(data(),data()+ai,aw);}void clear()noexcept{ai=0;}void push_back(Limb ax){reserve(ai+1);data()[ai++]=ax;}void pop_back()noexcept{--ai;}};LimbStorage ay;bool az=false;static std::spantrim_span(std::spanaA){while(!aA.empty()&&aA.back()==0){aA=aA.first(aA.size()-1);}return aA;}void normalize()noexcept{while(!ay.empty()&&ay.back()==0)ay.pop_back();if(ay.empty())az=false;}templatevoid assign_integral(Integer aB){ay.clear();az=false;using Value=std::remove_cv_t;if constexpr(std::same_as){if(aB)ay.push_back(1);return;}using Unsigned=exact_integer_detail::MakeUnsignedT;Unsigned aC=static_cast(aB);if constexpr(std::numeric_limits::is_signed){if(aB<0){az=true;aC=Unsigned{0}-aC;}}if constexpr(std::numeric_limits::digits<=64){if(aC!=0)ay.push_back(static_cast(aC));}else{while(aC!=0){ay.push_back(static_cast(aC));aC>>=64;}}}static int compare_magnitude(const ExactInteger&aD,const ExactInteger&aE){if(aD.ay.size()!=aE.ay.size()){return aD.ay.size()0;){if(aD.ay[aF]!=aE.ay[aF]){return aD.ay[aF](ay[aK])+aG.ay[aK]+aJ;ay[aK]=static_cast(aL);aJ=aL>>64;}while(aJ!=0&&aK(ay[aK])+aJ;ay[aK]=static_cast(aM);aJ=aM>>64;++aK;}if(aJ!=0)ay.push_back(static_cast(aJ));}void subtract_magnitude(const ExactInteger&aN){Wide aO=0;for(std::size_t aP=0;aP(aN.ay[aP]):0)+aO;const Wide aR=ay[aP];ay[aP]=static_cast(aR-aQ);aO=aR(ay[aU])+aT;ay[aU]=static_cast(aV);aT=aV>>64;}if(aT!=0)ay.push_back(static_cast(aT));}void multiply_magnitude_limb(Limb aW){if(aW==0||ay.empty()){ay.clear();az=false;return;}if(aW==1)return;Wide aX=0;for(Limb&aY:ay){const Wide aZ=static_cast(aY)*aW+aX;aY=static_cast(aZ);aX=aZ>>64;}if(aX!=0)ay.push_back(static_cast(aX));}Limb divide_magnitude_limb_inplace(Limb ba){Limb bb=0;for(std::size_t bc=ay.size();bc-->0;){ay[bc]=divide_two_by_one(bb,ay[bc],ba,bb);}normalize();return bb;}static std::vectoradd_vectors(std::spanbd,std::spanbe){const std::size_t bf=std::max(bd.size(),be.size());std::vectorbg(bf+1,0);Wide bh=0;for(std::size_t bi=0;bi(bj);bh=bj>>64;}bg[bf]=static_cast(bh);while(!bg.empty()&&bg.back()==0)bg.pop_back();return bg;}static void subtract_vector_inplace(std::vector&bk,std::spanbl){Wide bm=0;for(std::size_t bn=0;bn(bl[bn]):0)+bm;const Wide bp=bk[bn];bk[bn]=static_cast(bp-bo);bm=bp[[gnu::always_inline]]static inline void schoolbook_multiply_fixed(const Limb*bq,const Limb*br,Limb*bs)noexcept{std::fill_n(bs,LeftSize+RightSize,Limb{0});for(std::size_t bt=0;bt(bq[bw])*bv+bs[bx]+bu;bs[bx]=static_cast(bz);bu=bz>>64;}bs[bt+LeftSize]=static_cast(bu);}}static void schoolbook_multiply_raw_bmi2(const Limb*bA,std::size_t bB,const Limb*bC,std::size_t bD,Limb*bE){if(bB==bD){switch(bB){case 2:schoolbook_multiply_fixed<2,2>(bA,bC,bE);return;case 3:schoolbook_multiply_fixed<3,3>(bA,bC,bE);return;case 4:schoolbook_multiply_fixed<4,4>(bA,bC,bE);return;default:break;}}std::fill(bE,bE+bB+bD,0);if(bB(bA[bI]),static_cast(bG),&bK);unsigned long long bM;unsigned char bN=_addcarry_u64(0,bL,static_cast(bE[bJ]),&bM);bK+=bN;unsigned long long bO;bN=_addcarry_u64(0,bM,static_cast(bH),&bO);bK+=bN;bE[bJ]=static_cast(bO);bH=static_cast(bK);}bE[bF+bB]=bH; #else Wide bP=0;for(std::size_t bQ=0;bQ(bA[bQ])*bG+bE[bR]+bP;bE[bR]=static_cast(bS);bP=bS>>64;}bE[bF+bB]=static_cast(bP); #endif }}static void schoolbook_multiply_raw(const Limb*bT,std::size_t bU,const Limb*bV,std::size_t bW,Limb*bX){if(bU==bW){switch(bU){case 2:schoolbook_multiply_fixed<2,2>(bT,bV,bX);return;case 3:schoolbook_multiply_fixed<3,3>(bT,bV,bX);return;case 4:schoolbook_multiply_fixed<4,4>(bT,bV,bX);return;default:break;}}std::fill(bX,bX+bU+bW,0);if(bU(bT[cb])*ca+bX[cd]+bZ;bX[cd]=static_cast(ce);bZ=ce>>64;}bX[bY+bU]=static_cast(bZ);}}static int compare_equal_arrays(const Limb*cf,const Limb*cg,std::size_t ch)noexcept{for(std::size_t ci=ch;ci-->0;){if(cf[ci]!=cg[ci]){return cf[ci](cn[cu])+cp[ct]+cs;cn[cu]=static_cast(cv);cs=cv>>64;}std::size_t cw=cr+ct;while(cs!=0&&cw(cn[cw])+cs;cn[cw]=static_cast(cx);cs=cx>>64;++cw;}}static void subtract_array_at(Limb*cy,std::size_t cz,const Limb*cA,std::size_t cB,std::size_t cC)noexcept{Wide cD=0;std::size_t cE=0;for(;cE(cA[cE])+cD;const Wide cH=cy[cF];cy[cF]=static_cast(cH-cG);cD=cH=0){subtract_equal_arrays(cU,cL+cR,cL,cR);}else{subtract_equal_arrays(cU,cL,cL+cR,cR);}const int cZ=compare_equal_arrays(cM,cM+cR,cR);if(cZ>=0){subtract_equal_arrays(cV,cM,cM+cR,cR);}else{subtract_equal_arrays(cV,cM+cR,cM,cR);}karatsuba_equal_raw(cU,cV,cR,cW,cX);std::fill(cO,cO+cN*2,0);std::copy_n(cS,cN,cO);std::copy_n(cT,cN,cO+cN);add_array_at(cO,cN*2,cS,cN,cR);add_array_at(cO,cN*2,cT,cN,cR);if((cY<0)!=(cZ<0)){subtract_array_at(cO,cN*2,cW,cN,cR);}else{add_array_at(cO,cN*2,cW,cN,cR);}}static bool can_use_power_two_karatsuba(std::size_t da,std::size_t db)noexcept{const std::size_t dc=std::min(da,db);const std::size_t dd=std::max(da,db);if(dc<64||dd>dc*2)return false;const std::size_t de=std::bit_ceil(dd);return dc*4>=de*3;}static void power_two_karatsuba_into(std::spandf,std::spandg,Limb*dh,std::size_t di){struct KaratsubaScratch{std::vectordj;std::vectordk;std::vectordl;};static thread_local KaratsubaScratch scratch;scratch.dj.resize(karatsuba_scratch_size(di));const Limb*dm=df.data();const Limb*dn=dg.data();if(df.size()!=di){scratch.dk.assign(di,0);std::copy(df.begin(),df.end(),scratch.dk.begin());dm=scratch.dk.data();}if(dg.size()!=di){scratch.dl.assign(di,0);std::copy(dg.begin(),dg.end(),scratch.dl.begin());dn=scratch.dl.data();}karatsuba_equal_raw(dm,dn,di,dh,scratch.dj.data());}static std::vectorschoolbook_multiply(std::spandp,std::spandq){dp=trim_span(dp);dq=trim_span(dq);if(dp.empty()||dq.empty())return{};std::vectordr(dp.size()+dq.size(),0);schoolbook_multiply_raw(dp.data(),dp.size(),dq.data(),dq.size(),dr.data());while(!dr.empty()&&dr.back()==0)dr.pop_back();return dr;}static void add_shifted(std::vector&ds,std::spandt,std::size_t du){if(dt.empty())return;const std::size_t dv=du+dt.size();if(ds.size()(ds[dy])+dt[dx]+dw;ds[dy]=static_cast(dz);dw=dz>>64;}std::size_t dA=du+dx;while(dw!=0){const Wide dB=static_cast(ds[dA])+dw;ds[dA]=static_cast(dB);dw=dB>>64;++dA;if(dA==ds.size()&&dw!=0){ds.push_back(0);}}}static std::vectormultiply_recursive(std::spandC,std::spandD){dC=trim_span(dC);dD=trim_span(dD);if(dC.empty()||dD.empty())return{};if(dC.size()dD.size()*2){return schoolbook_multiply(dC,dD);}const std::size_t dF=dC.size()/2;const auto dG=dC.first(dF);const auto dH=dC.subspan(dF);const auto dI=dD.first(std::min(dF,dD.size()));const auto dJ=dD.subspan(std::min(dF,dD.size()));auto dK=multiply_recursive(dG,dI);auto dL=multiply_recursive(dH,dJ);auto dM=add_vectors(dG,dH);auto dN=add_vectors(dI,dJ);auto dO=multiply_recursive(dM,dN);subtract_vector_inplace(dO,dK);subtract_vector_inplace(dO,dL);std::vectordP(dC.size()+dD.size()+1,0);add_shifted(dP,dK,0);add_shifted(dP,dO,dF);add_shifted(dP,dL,dF*2);while(!dP.empty()&&dP.back()==0)dP.pop_back();return dP;}static ExactInteger add_values(const ExactInteger&dQ,const ExactInteger&dR){if(dQ.is_zero())return dR;if(dR.is_zero())return dQ;ExactInteger dS;if(dQ.az==dR.az){const std::size_t dT=std::max(dQ.ay.size(),dR.ay.size());dS.ay.resize(dT+1,0);Wide dU=0;std::size_t dV=0;const std::size_t dW=std::min(dQ.ay.size(),dR.ay.size());for(;dV(dQ.ay[dV])+dR.ay[dV]+dU;dS.ay[dV]=static_cast(dX);dU=dX>>64;}const ExactInteger&dY=dQ.ay.size()>=dR.ay.size()?dQ:dR;for(;dV(dY.ay[dV])+dU;dS.ay[dV]=static_cast(dZ);dU=dZ>>64;}dS.ay[dT]=static_cast(dU);dS.az=dQ.az;dS.normalize();return dS;}const int ea=compare_magnitude(dQ,dR);if(ea==0)return dS;const ExactInteger&eb=ea>0?dQ:dR;const ExactInteger&ec=ea>0?dR:dQ;dS.ay.resize(eb.ay.size(),0);Wide ed=0;std::size_t ee=0;for(;ee(ec.ay[ee])+ed;const Wide eg=eb.ay[ee];dS.ay[ee]=static_cast(eg-ef);ed=eg(eh-ed);ed=eh(ei.ay[eo])+ej.ay[eo]+en;el.ay[eo]=static_cast(eq);en=eq>>64;}const ExactInteger&er=ei.ay.size()>=ej.ay.size()?ei:ej;for(;eo(er.ay[eo])+en;el.ay[eo]=static_cast(es);en=es>>64;}el.ay[em]=static_cast(en);el.az=ei.az;el.normalize();return el;}const int et=compare_magnitude(ei,ej);if(et==0)return el;const ExactInteger&eu=et>0?ei:ej;const ExactInteger&ev=et>0?ej:ei;el.ay.resize(eu.ay.size(),0);Wide ew=0;std::size_t ex=0;for(;ex(ev.ay[ex])+ew;const Wide ez=eu.ay[ex];el.ay[ex]=static_cast(ez-ey);ew=ez(eA-ew);ew=eA0?ei.az:!ei.az;el.normalize();return el;}static ExactInteger schoolbook_multiply_value(std::spaneB,std::spaneC,bool eD){eB=trim_span(eB);eC=trim_span(eC);ExactInteger eE;if(eB.empty()||eC.empty())return eE;eE.ay.assign(eB.size()+eC.size(),0);schoolbook_multiply_raw(eB.data(),eB.size(),eC.data(),eC.size(),eE.ay.data());eE.az=eD;eE.normalize();return eE;}static ExactInteger multiply_values(const ExactInteger&eF,const ExactInteger&eG){if(eF.is_zero()||eG.is_zero())return ExactInteger{};const bool eH=eF.az!=eG.az;if(eG.ay.size()==1){ExactInteger eI=eF;eI.multiply_magnitude_limb(eG.ay[0]);eI.az=eH;return eI;}if(eF.ay.size()==1){ExactInteger eJ=eG;eJ.multiply_magnitude_limb(eF.ay[0]);eJ.az=eH;return eJ;}constexpr std::size_t eK=36;const std::size_t eL=std::min(eF.ay.size(),eG.ay.size());const std::size_t eM=std::max(eF.ay.size(),eG.ay.size());if(can_use_power_two_karatsuba(eF.ay.size(),eG.ay.size())){const std::size_t eN=std::bit_ceil(eM);ExactInteger eO;eO.ay.assign(eN*2,0);power_two_karatsuba_into(std::span(eF.ay.data(),eF.ay.size()),std::span(eG.ay.data(),eG.ay.size()),eO.ay.data(),eN);eO.az=eH;eO.normalize();return eO;}if(eL<=eK||eM>eL*2){return schoolbook_multiply_value(std::span(eF.ay.data(),eF.ay.size()),std::span(eG.ay.data(),eG.ay.size()),eH);}auto eP=multiply_recursive(std::span(eF.ay.data(),eF.ay.size()),std::span(eG.ay.data(),eG.ay.size()));return from_vector(std::move(eP),eH);}static ExactInteger from_vector(std::vector&&eQ,bool eR=false){ExactInteger eS;eS.ay.resize(eQ.size());std::copy(eQ.begin(),eQ.end(),eS.ay.begin());eS.az=eR&&!eQ.empty();eS.normalize();return eS;}static std::vectornormalized_left_shift(std::spaneT,unsigned eU,bool eV){std::vectoreW(eT.size()+static_cast(eV),0);if(eU==0){std::copy(eT.begin(),eT.end(),eW.begin());return eW;}Limb eX=0;for(std::size_t eY=0;eY>(64-eU);}if(eV)eW[eT.size()]=eX;return eW;}static std::vectornormalized_right_shift(std::spanfa,unsigned fb){std::vectorfc(fa.size(),0);if(fb==0){std::copy(fa.begin(),fa.end(),fc.begin());}else{Limb fd=0;for(std::size_t fe=fa.size();fe-->0;){const Limb ff=fa[fe];fc[fe]=(ff>>fb)|fd;fd=ff<<(64-fb);}}while(!fc.empty()&&fc.back()==0)fc.pop_back();return fc;}static void divide_magnitudes_into(const ExactInteger&fg,const ExactInteger&fh,ExactInteger*fi,ExactInteger*fj){const int fk=compare_magnitude(fg,fh);if(fk<0){if(fi){fi->ay.clear();fi->az=false;}if(fj){*fj=fg;fj->az=false;}return;}if(fk==0){if(fi)*fi=1;if(fj){fj->ay.clear();fj->az=false;}return;}if(fh.ay.size()==1){Limb fl=0;if(fi){fi->ay.resize(fg.ay.size());std::copy(fg.ay.begin(),fg.ay.end(),fi->ay.begin());fi->az=false;fl=fi->divide_magnitude_limb_inplace(fh.ay[0]);}else{for(std::size_t fm=fg.ay.size();fm-->0;){(void)divide_two_by_one(fl,fg.ay[fm],fh.ay[0],fl);}}if(fj)*fj=fl;return;}const std::size_t fn=fh.ay.size();const std::size_t fo=fg.ay.size();const std::size_t fp=fo-fn;const unsigned fq=std::countl_zero(fh.ay.back());struct DivisionScratch{std::vectorfr;std::vectorfs;};static thread_local DivisionScratch scratch;auto&ft=scratch.fr;auto&fu=scratch.fs;ft.resize(fn);fu.resize(fo+1);if(fq==0){std::copy_n(fh.ay.data(),fn,ft.data());std::copy_n(fg.ay.data(),fo,fu.data());fu[fo]=0;}else{Limb fv=0;for(std::size_t fw=0;fw>(64-fq);}fv=0;for(std::size_t fy=0;fy>(64-fq);}fu[fo]=fv;}if(fi){fi->ay.resize(fp+1,0);fi->az=false;}constexpr Wide fA=Wide{1}<<64;const Limb fB=ft[fn-1];const Limb fC=ft[fn-2];for(std::size_t fD=fp+1;fD-->0;){const std::size_t fE=fD;Limb fF;Wide fG;if(fu[fE+fn]==fB){fF=(std::numeric_limits::max)();fG=static_cast(fu[fE+fn-1])+fB;}else{Limb fH=0;fF=divide_two_by_one(fu[fE+fn],fu[fE+fn-1],fB,fH);fG=fH;}while(fG(fF)*fC>(fG<<64)+fu[fE+fn-2]){--fF;fG+=fB;} #if defined(__x86_64__) && defined(__BMI2__) \ && (defined(__GNUC__) || defined(__clang__)) Limb fI=0;unsigned char fJ=0;for(std::size_t fK=0;fK(fF),static_cast(ft[fK]),&fL);unsigned long long fN;const unsigned char fO=_addcarry_u64(0,fM,static_cast(fI),&fN);fL+=fO;unsigned long long fP;fJ=_subborrow_u64(fJ,static_cast(fu[fE+fK]),fN,&fP);fu[fE+fK]=static_cast(fP);fI=static_cast(fL);}unsigned long long fQ;unsigned char fR=_subborrow_u64(0,static_cast(fu[fE+fn]),static_cast(fI),&fQ);unsigned long long fS;fR=_subborrow_u64(fR,fQ,static_cast(fJ),&fS);const bool fT=fR!=0;fu[fE+fn]=static_cast(fS);if(fT){--fF;unsigned char fU=0;for(std::size_t fV=0;fV(fu[fE+fV]),static_cast(ft[fV]),&fW);fu[fE+fV]=static_cast(fW);}fu[fE+fn]+=static_cast(fU);} #else Wide fX=0;Wide fY=0;for(std::size_t fZ=0;fZ(fF)*ft[fZ]+fX;fX=ga>>64;const Wide gb=static_cast(ga)+fY;const Wide gc=fu[fE+fZ];fu[fE+fZ]=static_cast(gc-gb);fY=gc(fu[fE+fn])(static_cast(fu[fE+fn])-gd);if(ge){--fF;Wide gf=0;for(std::size_t gg=0;gg(fu[fE+gg])+ft[gg]+gf;fu[fE+gg]=static_cast(gh);gf=gh>>64;}fu[fE+fn]+=static_cast(gf);} #endif if(fi)fi->ay[fE]=fF;}if(fi)fi->normalize();if(fj){fj->ay.resize(fn,0);fj->az=false;if(fq==0){std::copy_n(fu.data(),fn,fj->ay.data());}else{Limb gi=0;for(std::size_t gj=fn;gj-->0;){const Limb gk=fu[gj];fj->ay[gj]=(gk>>fq)|gi;gi=gk<<(64-fq);}}fj->normalize();}}void add_magnitude_one(){add_magnitude_limb(1);}templateUnsigned magnitude_to_unsigned()const{static_assert(!std::numeric_limits::is_signed);Unsigned gl=0;if constexpr(std::numeric_limits::digits<=64){if(!ay.empty())gl=static_cast(ay[0]);}else{for(std::size_t gm=ay.size();gm-->0;){gl=static_cast(gl<<64);gl=static_cast(gl|ay[gm]);}}return gl;}public:ExactInteger()=default;templateExactInteger(Integer gn){assign_integral(gn);}templateExactInteger&operator=(Integer go){assign_integral(go);return*this;}bool is_zero()const noexcept{return ay.empty();}bool is_negative()const noexcept{return az;}std::size_t bit_length()const noexcept{if(ay.empty())return 0;return(ay.size()-1)*64+static_cast(64-std::countl_zero(ay.back()));}ExactInteger absolute()const{ExactInteger gp=*this;gp.az=false;return gp;}templateInteger checked_to()const{using Value=std::remove_cv_t;if constexpr(std::same_as){if(*this==0)return false;if(*this==1)return true;throw std::overflow_error("ExactInteger does not fit target integer type");}else{const ExactInteger gq=std::numeric_limits::is_signed?ExactInteger((std::numeric_limits::min)()):ExactInteger(0);const ExactInteger gr((std::numeric_limits::max)());if(*thisgr){throw std::overflow_error("ExactInteger does not fit target integer type");}using Unsigned=exact_integer_detail::MakeUnsignedT;const Unsigned gs=magnitude_to_unsigned();if constexpr(!std::numeric_limits::is_signed){return static_cast(gs);}else{if(!az)return static_cast(gs);const Unsigned gt=static_cast((std::numeric_limits::max)())+Unsigned{1};if(gs==gt){return(std::numeric_limits::min)();}return static_cast(-static_cast(gs));}}}std::pairdivmod(Limb gu)const{if(gu==0){throw std::domain_error("ExactInteger division by zero");}ExactInteger gv=absolute();const Limb gw=gv.divide_magnitude_limb_inplace(gu);gv.az=az&&!gv.is_zero();return{std::move(gv),gw};}static std::pairdivmod(const ExactInteger&gx,const ExactInteger&gy){if(gy.is_zero()){throw std::domain_error("ExactInteger division by zero");}std::pairgz;divide_magnitudes_into(gx,gy,&gz.first,&gz.second);gz.first.az=!gz.first.is_zero()&&gx.az!=gy.az;gz.second.az=!gz.second.is_zero()&&gx.az;return gz;}void assign_quotient(const ExactInteger&gA,const ExactInteger&gB){if(gB.is_zero()){throw std::domain_error("ExactInteger division by zero");}if(this==&gA||this==&gB){ExactInteger gC;divide_magnitudes_into(gA,gB,&gC,nullptr);gC.az=!gC.is_zero()&&gA.az!=gB.az;*this=std::move(gC);return;}divide_magnitudes_into(gA,gB,this,nullptr);az=!is_zero()&&gA.az!=gB.az;}void assign_remainder(const ExactInteger&gD,const ExactInteger&gE){if(gE.is_zero()){throw std::domain_error("ExactInteger division by zero");}if(this==&gD||this==&gE){ExactInteger gF;divide_magnitudes_into(gD,gE,nullptr,&gF);gF.az=!gF.is_zero()&&gD.az;*this=std::move(gF);return;}divide_magnitudes_into(gD,gE,nullptr,this);az=!is_zero()&&gD.az;}[[gnu::always_inline]]inline void assign_sum(const ExactInteger&gG,const ExactInteger&gH){if(this==&gG){*this+=gH;return;}if(this==&gH){*this+=gG;return;}if(gG.is_zero()){*this=gH;return;}if(gH.is_zero()){*this=gG;return;}if(gG.az!=gH.az){const int gI=compare_magnitude(gG,gH);if(gI==0){ay.clear();az=false;return;}const ExactInteger&gJ=gI>0?gG:gH;const ExactInteger&gK=gI>0?gH:gG;ay.resize(gJ.ay.size(),0);if(gJ.ay.size()==gK.ay.size()){(void)subtract_equal_length(ay.data(),gJ.ay.data(),gK.ay.data(),gJ.ay.size());az=gJ.az;normalize();return;}Wide gL=0;std::size_t gM=0;for(;gM(gK.ay[gM])+gL;const Wide gO=gJ.ay[gM];ay[gM]=static_cast(gO-gN);gL=gO(gP-gL);gL=gP(gG.ay[gS])+gH.ay[gS]+gR;ay[gS]=static_cast(gU);gR=gU>>64;}const ExactInteger&gV=gG.ay.size()>=gH.ay.size()?gG:gH;for(;gS(gV.ay[gS])+gR;ay[gS]=static_cast(gW);gR=gW>>64;}ay[gQ]=static_cast(gR);az=gG.az;normalize();}[[gnu::always_inline]]inline void assign_difference(const ExactInteger&gX,const ExactInteger&gY){if(this==&gX){*this-=gY;return;}if(this==&gY){ExactInteger gZ=subtract_values(gX,gY);*this=std::move(gZ);return;}if(gY.is_zero()){*this=gX;return;}if(gX.is_zero()){*this=gY;if(!is_zero())az=!az;return;}if(gX.az!=gY.az){const std::size_t ha=std::max(gX.ay.size(),gY.ay.size());ay.resize(ha+1,0);if(gX.ay.size()==gY.ay.size()){ay[ha]=add_equal_length(ay.data(),gX.ay.data(),gY.ay.data(),ha);az=gX.az;normalize();return;}Wide hb=0;std::size_t hc=0;const std::size_t hd=std::min(gX.ay.size(),gY.ay.size());for(;hc(gX.ay[hc])+gY.ay[hc]+hb;ay[hc]=static_cast(he);hb=he>>64;}const ExactInteger&hf=gX.ay.size()>=gY.ay.size()?gX:gY;for(;hc(hf.ay[hc])+hb;ay[hc]=static_cast(hg);hb=hg>>64;}ay[ha]=static_cast(hb);az=gX.az;normalize();return;}const int hh=compare_magnitude(gX,gY);if(hh==0){ay.clear();az=false;return;}const ExactInteger&hi=hh>0?gX:gY;const ExactInteger&hj=hh>0?gY:gX;ay.resize(hi.ay.size(),0);if(hi.ay.size()==hj.ay.size()){(void)subtract_equal_length(ay.data(),hi.ay.data(),hj.ay.data(),hi.ay.size());az=hh>0?gX.az:!gX.az;normalize();return;}Wide hk=0;std::size_t hl=0;for(;hl(hj.ay[hl])+hk;const Wide hn=hi.ay[hl];ay[hl]=static_cast(hn-hm);hk=hn(ho-hk);hk=ho0?gX.az:!gX.az;normalize();}void assign_product(const ExactInteger&hp,const ExactInteger&hq){if(this==&hp||this==&hq){ExactInteger hr=multiply_values(hp,hq);*this=std::move(hr);return;}if(hp.is_zero()||hq.is_zero()){ay.clear();az=false;return;}const bool hs=hp.az!=hq.az;if(hq.ay.size()==1){*this=hp;multiply_magnitude_limb(hq.ay[0]);az=hs;return;}if(hp.ay.size()==1){*this=hq;multiply_magnitude_limb(hp.ay[0]);az=hs;return;}constexpr std::size_t ht=36;const std::size_t hu=std::min(hp.ay.size(),hq.ay.size());const std::size_t hv=std::max(hp.ay.size(),hq.ay.size());if(can_use_power_two_karatsuba(hp.ay.size(),hq.ay.size())){const std::size_t hw=std::bit_ceil(hv);ay.assign(hw*2,0);power_two_karatsuba_into(std::span(hp.ay.data(),hp.ay.size()),std::span(hq.ay.data(),hq.ay.size()),ay.data(),hw);az=hs;normalize();return;}if(hu<=ht||hv>hu*2){const ExactInteger*hx=&hp;const ExactInteger*hy=&hq;if(hp.ay.size()ay.size()+hy->ay.size(),0);schoolbook_multiply_raw(hx->ay.data(),hx->ay.size(),hy->ay.data(),hy->ay.size(),ay.data());az=hs;normalize();return;}auto hz=multiply_recursive(std::span(hp.ay.data(),hp.ay.size()),std::span(hq.ay.data(),hq.ay.size()));ay.resize(hz.size());std::copy(hz.begin(),hz.end(),ay.begin());az=hs;normalize();}void multiply_add(Limb hA,Limb hB){if(is_zero()){if(hB!=0)ay.push_back(hB);return;}if(az){throw std::logic_error("multiply_add requires a nonnegative ExactInteger");}Wide hC=hB;for(Limb&hD:ay){const Wide hE=static_cast(hD)*hA+hC;hD=static_cast(hE);hC=hE>>64;}if(hC!=0)ay.push_back(static_cast(hC));}ExactInteger operator-()const{ExactInteger hF=*this;if(!hF.is_zero())hF.az=!hF.az;return hF;}ExactInteger&operator+=(const ExactInteger&hG){if(hG.is_zero())return*this;if(is_zero()){*this=hG;return*this;}if(az==hG.az){add_magnitude(hG);return*this;}const int hH=compare_magnitude(*this,hG);if(hH==0){ay.clear();az=false;}else if(hH>0){subtract_magnitude(hG);}else{ExactInteger hI=hG;hI.subtract_magnitude(*this);*this=std::move(hI);}return*this;}templateExactInteger&operator+=(Integer hJ){return*this+=ExactInteger(hJ);}ExactInteger&operator-=(const ExactInteger&hK){if(hK.is_zero())return*this;if(this==&hK){ay.clear();az=false;return*this;}ExactInteger hL=hK;hL.az=!hL.az;return*this+=hL;}templateExactInteger&operator-=(Integer hM){return*this-=ExactInteger(hM);}ExactInteger&operator*=(const ExactInteger&hN){*this=multiply_values(*this,hN);return*this;}templateExactInteger&operator*=(Integer hO){using Value=std::remove_cv_t;if constexpr(std::numeric_limits::digits<=64){if constexpr(std::numeric_limits::is_signed){if(hO<0){const auto hP=static_cast(exact_integer_detail::MakeUnsignedT{0}-static_cast>(hO));const bool hQ=!az;multiply_magnitude_limb(hP);az=hQ&&!is_zero();return*this;}}multiply_magnitude_limb(static_cast(hO));return*this;}else{return*this*=ExactInteger(hO);}}ExactInteger&operator<<=(std::size_t hR){if(is_zero()||hR==0)return*this;const std::size_t hS=hR/64;const unsigned hT=static_cast(hR%64);const std::size_t hU=ay.size();ay.resize(hU+hS+1,0);if(hS!=0){for(std::size_t hV=hU;hV-->0;){ay[hV+hS]=ay[hV];}std::fill(ay.begin(),ay.begin()+static_cast(hS),0);}if(hT!=0){Limb hW=0;for(std::size_t hX=hS;hX>(64-hT);}ay[hS+hU]=hW;}else{ay[hS+hU]=0;}normalize();return*this;}ExactInteger&operator>>=(std::size_t hZ){if(is_zero()||hZ==0)return*this;const bool ia=az;const std::size_t ib=hZ/64;const unsigned ic=static_cast(hZ%64);bool id=false;const std::size_t ie=std::min(ib,ay.size());for(std::size_t ig=0;ig=ay.size()){ay.clear();}else{const std::size_t ii=ay.size()-ib;if(ib!=0){for(std::size_t ij=0;ij0;){const Limb im=ay[il];ay[il]=(im>>ic)|ik;ik=im<<(64-ic);}}}normalize();if(ia&&id){add_magnitude_one();az=true;}else if(!is_zero()){az=ia;}return*this;}std::string to_string()const{if(is_zero())return"0";ExactInteger in=absolute();std::vectorio;io.reserve((bit_length()/63)+1);while(!in.is_zero()){io.push_back(in.divide_magnitude_limb_inplace(decimal_base));}std::string ip;ip.reserve(io.size()*19+static_cast(az));if(az)ip.push_back('-');char iq[32];auto[ir,is]=std::to_chars(iq,iq+sizeof(iq),io.back());if(is!=std::errc{})throw std::runtime_error("ExactInteger decimal conversion failed");ip.append(iq,ir);for(std::size_t it=io.size()-1;it-->0;){auto[iu,iv]=std::to_chars(iq,iq+sizeof(iq),io[it]);if(iv!=std::errc{})throw std::runtime_error("ExactInteger decimal conversion failed");const std::size_t iw=static_cast(iu-iq);ip.append(19-iw,'0');ip.append(iq,iu);}return ip;}friend ExactInteger abs(const ExactInteger&ix){return ix.absolute();}friend bool operator==(const ExactInteger&iy,const ExactInteger&iz){if(iy.az!=iz.az||iy.ay.size()!=iz.ay.size()){return false;}return std::equal(iy.ay.begin(),iy.ay.end(),iz.ay.begin());}friend std::strong_ordering operator<=>(const ExactInteger&iA,const ExactInteger&iB){if(iA.az!=iB.az){return iA.az?std::strong_ordering::less:std::strong_ordering::greater;}const int iC=compare_magnitude(iA,iB);if(iC==0)return std::strong_ordering::equal;const bool iD=iA.az?iC>0:iC<0;return iD?std::strong_ordering::less:std::strong_ordering::greater;}friend ExactInteger operator+(const ExactInteger&iE,const ExactInteger&iF){return add_values(iE,iF);}friend ExactInteger operator-(const ExactInteger&iG,const ExactInteger&iH){return subtract_values(iG,iH);}friend ExactInteger operator*(const ExactInteger&iI,const ExactInteger&iJ){return multiply_values(iI,iJ);}friend ExactInteger operator<<(ExactInteger iK,std::size_t iL){iK<<=iL;return iK;}friend ExactInteger operator>>(ExactInteger iM,std::size_t iN){iM>>=iN;return iM;}friend std::ostream&operator<<(std::ostream&iO,const ExactInteger&iP){return iO<'9'){throw std::invalid_argument("invalid BigInteger decimal digit");}}while(iT(iZ);const char*jd=jc+static_cast(ja);const auto[je,jf]=std::from_chars(jc,jd,jb);if(jf!=std::errc{}||je!=jd){throw std::invalid_argument("invalid BigInteger decimal digit");}return jb;};iV=iY(iT,iX);iT+=iX;constexpr std::uint64_t jg=10'000'000'000'000'000'000ULL;while(iTBigInteger(Integer jh):iQ(jh){}explicit BigInteger(std::string_view ji):iQ(parse_decimal(ji)){}explicit BigInteger(const ExactInteger&jk):iQ(jk){}explicit BigInteger(ExactInteger&&jl):iQ(std::move(jl)){}BigInteger(const big_integer_detail::AddExpression&jm);BigInteger(const big_integer_detail::SubtractExpression&jn);BigInteger(const big_integer_detail::MultiplyExpression&jo);BigInteger(const big_integer_detail::DivideExpression&jp);BigInteger(const big_integer_detail::ModuloExpression&jq);BigInteger&operator=(const big_integer_detail::AddExpression&jr);BigInteger&operator=(const big_integer_detail::SubtractExpression&js);BigInteger&operator=(const big_integer_detail::MultiplyExpression&jt);BigInteger&operator=(const big_integer_detail::DivideExpression&ju);BigInteger&operator=(const big_integer_detail::ModuloExpression&jv);templateBigInteger&operator=(Integer jw){iQ=jw;return*this;}BigInteger&assign(std::string_view jx){iQ=parse_decimal(jx);return*this;}bool is_zero()const noexcept{return iQ.is_zero();}bool is_negative()const noexcept{return iQ.is_negative();}std::size_t bit_length()const noexcept{return iQ.bit_length();}BigInteger absolute()const{return BigInteger(iQ.absolute());}std::string to_string()const{return iQ.to_string();}templateInteger checked_to()const{return iQ.template checked_to();}const ExactInteger&exact_integer()const noexcept{return iQ;}static std::pairdivmod(const BigInteger&jy,const BigInteger&jz){auto jA=ExactInteger::divmod(jy.iQ,jz.iQ);return{BigInteger(std::move(jA.first)),BigInteger(std::move(jA.second))};}BigInteger operator-()const{return BigInteger(-iQ);}BigInteger&operator+=(const BigInteger&jB){iQ+=jB.iQ;return*this;}BigInteger&operator-=(const BigInteger&jC){iQ-=jC.iQ;return*this;}BigInteger&operator*=(const BigInteger&jD){iQ*=jD.iQ;return*this;}BigInteger&operator/=(const BigInteger&jE){iQ.assign_quotient(iQ,jE.iQ);return*this;}BigInteger&operator%=(const BigInteger&jF){iQ.assign_remainder(iQ,jF.iQ);return*this;}BigInteger&operator<<=(std::size_t jG){iQ<<=jG;return*this;}BigInteger&operator>>=(std::size_t jH){iQ>>=jH;return*this;}BigInteger&operator++(){iQ+=1;return*this;}BigInteger operator++(int){BigInteger jI=*this;++*this;return jI;}BigInteger&operator--(){iQ-=1;return*this;}BigInteger operator--(int){BigInteger jJ=*this;--*this;return jJ;}friend BigInteger abs(const BigInteger&jK){return jK.absolute();}friend bool operator==(const BigInteger&,const BigInteger&)=default;friend std::strong_ordering operator<=>(const BigInteger&jL,const BigInteger&jM){return jL.iQ<=>jM.iQ;}friend big_integer_detail::AddExpression operator+(const BigInteger&jN,const BigInteger&jO)noexcept{return{jN,jO};}friend big_integer_detail::SubtractExpression operator-(const BigInteger&jP,const BigInteger&jQ)noexcept{return{jP,jQ};}friend big_integer_detail::MultiplyExpression operator*(const BigInteger&jR,const BigInteger&jS)noexcept{return{jR,jS};}friend big_integer_detail::DivideExpression operator/(const BigInteger&jT,const BigInteger&jU)noexcept{return{jT,jU};}friend big_integer_detail::ModuloExpression operator%(const BigInteger&jV,const BigInteger&jW)noexcept{return{jV,jW};}friend BigInteger operator<<(BigInteger jX,std::size_t jY){jX<<=jY;return jX;}friend BigInteger operator>>(BigInteger jZ,std::size_t ka){jZ>>=ka;return jZ;}friend std::ostream&operator<<(std::ostream&kb,const BigInteger&kc){return kb<>(std::istream&kd,BigInteger&ke){std::string kf;if(!(kd>>kf))return kd;try{ke.assign(kf);}catch(const std::invalid_argument&){kd.setstate(std::ios::failbit);}return kd;}};inline BigInteger::BigInteger(const big_integer_detail::AddExpression&kg){iQ.assign_sum(kg.left.iQ,kg.right.iQ);}inline BigInteger::BigInteger(const big_integer_detail::SubtractExpression&kh){iQ.assign_difference(kh.left.iQ,kh.right.iQ);}inline BigInteger::BigInteger(const big_integer_detail::MultiplyExpression&ki){iQ.assign_product(ki.left.iQ,ki.right.iQ);}inline BigInteger::BigInteger(const big_integer_detail::DivideExpression&kj){iQ.assign_quotient(kj.left.iQ,kj.right.iQ);}inline BigInteger::BigInteger(const big_integer_detail::ModuloExpression&kk){iQ.assign_remainder(kk.left.iQ,kk.right.iQ);}inline BigInteger&BigInteger::operator=(const big_integer_detail::AddExpression&kl){iQ.assign_sum(kl.left.iQ,kl.right.iQ);return*this;}inline BigInteger&BigInteger::operator=(const big_integer_detail::SubtractExpression&km){iQ.assign_difference(km.left.iQ,km.right.iQ);return*this;}inline BigInteger&BigInteger::operator=(const big_integer_detail::MultiplyExpression&kn){iQ.assign_product(kn.left.iQ,kn.right.iQ);return*this;}inline BigInteger&BigInteger::operator=(const big_integer_detail::DivideExpression&ko){iQ.assign_quotient(ko.left.iQ,ko.right.iQ);return*this;}inline BigInteger&BigInteger::operator=(const big_integer_detail::ModuloExpression&kp){iQ.assign_remainder(kp.left.iQ,kp.right.iQ);return*this;}inline BigInteger operator+(BigInteger&&left,const BigInteger&right){left+=right;return std::move(left);}inline BigInteger operator+(const BigInteger&left,BigInteger&&right){right+=left;return std::move(right);}inline BigInteger operator+(BigInteger&&left,BigInteger&&right){left+=right;return std::move(left);}inline BigInteger operator-(BigInteger&&left,const BigInteger&right){left-=right;return std::move(left);}inline BigInteger operator-(const BigInteger&left,BigInteger&&right){BigInteger result(left);result-=right;return result;}inline BigInteger operator-(BigInteger&&left,BigInteger&&right){left-=right;return std::move(left);}inline BigInteger operator*(BigInteger&&left,const BigInteger&right){left*=right;return std::move(left);}inline BigInteger operator*(const BigInteger&left,BigInteger&&right){right*=left;return std::move(right);}inline BigInteger operator*(BigInteger&&left,BigInteger&&right){left*=right;return std::move(left);}inline BigInteger operator/(BigInteger&&left,const BigInteger&right){left/=right;return std::move(left);}inline BigInteger operator/(BigInteger&&left,BigInteger&&right){left/=right;return std::move(left);}inline BigInteger operator%(BigInteger&&left,const BigInteger&right){left%=right;return std::move(left);}inline BigInteger operator%(BigInteger&&left,BigInteger&&right){left%=right;return std::move(left);}templateinline BigInteger operator+(const Expression&left,const BigInteger&right){BigInteger result(left);result+=right;return result;}templateinline BigInteger operator+(const BigInteger&left,const Expression&right){BigInteger result(right);result+=left;return result;}templateinline BigInteger operator+(const LeftExpression&left,const RightExpression&right){BigInteger result(left);result+=BigInteger(right);return result;}templateinline BigInteger operator-(const Expression&left,const BigInteger&right){BigInteger result(left);result-=right;return result;}templateinline BigInteger operator-(const BigInteger&left,const Expression&right){BigInteger result(left);result-=BigInteger(right);return result;}templateinline BigInteger operator-(const LeftExpression&left,const RightExpression&right){BigInteger result(left);result-=BigInteger(right);return result;}templateinline BigInteger operator*(const Expression&left,const BigInteger&right){BigInteger result(left);result*=right;return result;}templateinline BigInteger operator*(const BigInteger&left,const Expression&right){BigInteger result(right);result*=left;return result;}templateinline BigInteger operator*(const LeftExpression&left,const RightExpression&right){BigInteger result(left);result*=BigInteger(right);return result;}templateinline BigInteger operator/(const Expression&left,const BigInteger&right){BigInteger result(left);result/=right;return result;}templateinline BigInteger operator/(const BigInteger&left,const Expression&right){BigInteger result(left);result/=BigInteger(right);return result;}templateinline BigInteger operator/(const LeftExpression&left,const RightExpression&right){BigInteger result(left);result/=BigInteger(right);return result;}templateinline BigInteger operator%(const Expression&left,const BigInteger&right){BigInteger result(left);result%=right;return result;}templateinline BigInteger operator%(const BigInteger&left,const Expression&right){BigInteger result(left);result%=BigInteger(right);return result;}templateinline BigInteger operator%(const LeftExpression&left,const RightExpression&right){BigInteger result(left);result%=BigInteger(right);return result;}templateinline BigInteger operator<<(const Expression&value,std::size_t shift){BigInteger result(value);result<<=shift;return result;}templateinline BigInteger operator>>(const Expression&value,std::size_t shift){BigInteger result(value);result>>=shift;return result;}templateinline bool operator==(const Expression&left,const BigInteger&right){return BigInteger(left)==right;}templateinline bool operator==(const BigInteger&left,const Expression&right){return left==BigInteger(right);}templateinline std::ostream&operator<<(std::ostream&stream,const Expression&value){return stream< rhs.distance; } return lhs.vertex > rhs.vertex; } }; } // namespace int main() { std::ios::sync_with_stdio(false); std::cin.tie(nullptr); int n, m; std::cin >> n >> m; std::vector rawEdges; rawEdges.reserve(m); // maximumExponent[p] is the largest exponent of p occurring in any b_i. std::vector maximumExponent(MAX_AB + 1, 0); for (int i = 0; i < m; ++i) { int u, v, a, b; std::cin >> u >> v >> a >> b; --u; --v; rawEdges.push_back({u, v, a, b}); int value = b; for (int prime = 2; prime * prime <= value; ++prime) { if (value % prime != 0) { continue; } int exponent = 0; while (value % prime == 0) { value /= prime; ++exponent; } maximumExponent[prime] = std::max(maximumExponent[prime], exponent); } if (value > 1) { maximumExponent[value] = std::max(maximumExponent[value], 1); } } // Every b_i divides this common denominator. cpp_int commonDenominator = 1; std::vector> primePowers; for (int prime = 2; prime <= MAX_AB; ++prime) { if (maximumExponent[prime] == 0) { continue; } primePowers.push_back({prime, maximumExponent[prime]}); for (int exponent = 0; exponent < maximumExponent[prime]; ++exponent) { commonDenominator *= prime; } } std::vector> graph(n); for (const RawEdge& raw : rawEdges) { cpp_int scaledWeight = commonDenominator; scaledWeight /= raw.denominator; scaledWeight *= raw.numerator; // One copy and one move are necessary because the graph is undirected. graph[raw.u].push_back({raw.v, scaledWeight}); graph[raw.v].push_back({raw.u, std::move(scaledWeight)}); } std::vector distance(n); std::vector reached(n, false); std::priority_queue, StateGreater> queue; reached[0] = true; distance[0] = 0; queue.push({cpp_int(0), 0}); while (!queue.empty()) { State current = queue.top(); queue.pop(); if (!reached[current.vertex] || current.distance != distance[current.vertex]) { continue; } for (const Edge& edge : graph[current.vertex]) { // Evaluate the arbitrary-precision addition exactly once. cpp_int nextDistance = current.distance + edge.weight; if (!reached[edge.to] || nextDistance < distance[edge.to]) { reached[edge.to] = true; distance[edge.to] = nextDistance; queue.push({std::move(nextDistance), edge.to}); } } } // The official output limit is approximately 8 MiB, so buffering the // complete answer avoids tens of thousands of formatted stream writes. std::string output; output.reserve(8U * 1024U * 1024U); for (int vertex = 1; vertex < n; ++vertex) { cpp_int numerator = distance[vertex]; cpp_int denominator = commonDenominator; // The complete prime factorization of the denominator is already // known. Reduce by small primes instead of running a general cpp_int // Euclidean gcd for every output vertex. for (const auto& [prime, exponent] : primePowers) { for (int count = 0; count < exponent; ++count) { if (numerator % prime != 0) { break; } numerator /= prime; denominator /= prime; } } output += numerator.to_string(); output.push_back(' '); output += denominator.to_string(); output.push_back('\n'); } std::cout << output; return 0; }