1+ #line 1 "cp-algo/math/laurent.hpp"
2+ #line 1 "cp-algo/util/big_alloc.hpp"
3+ #include < vector>
4+ #include < cstddef>
5+ #include < iostream>
6+ #if defined(__linux__) || defined(__unix__) || (defined(__APPLE__) && defined(__MACH__))
7+ # define CP_ALGO_USE_MMAP 1
8+ # include < sys/mman.h>
9+ #else
10+ # define CP_ALGO_USE_MMAP 0
11+ #endif
12+ namespace cp_algo {template <typename T,std::size_t Align=32 >
13+ class big_alloc {static_assert (Align>=alignof (void *)," Align must be at least pointer-size" );
14+ static_assert (std::popcount(Align)==1 ," Align must be a power of two" );
15+ public: using value_type=T;
16+ template <class U >struct rebind {using other=big_alloc<U,Align>;};
17+ constexpr bool operator ==(const big_alloc&)const =default ;
18+ constexpr bool operator !=(const big_alloc&)const =default ;
19+ big_alloc ()noexcept =default ;
20+ template <typename U,std::size_t A>
21+ big_alloc (const big_alloc<U,A>&)noexcept {}
22+ [[nodiscard]]T*allocate (std::size_t n){std::size_t padded=round_up (n*sizeof (T));
23+ std::size_t align=std::max<std::size_t >(alignof (T),Align);
24+ #if CP_ALGO_USE_MMAP
25+ if (padded>=MEGABYTE){void *raw=mmap (nullptr ,padded,
26+ PROT_READ|PROT_WRITE,
27+ MAP_PRIVATE|MAP_ANONYMOUS,-1 ,0 );
28+ madvise (raw,padded,MADV_HUGEPAGE);
29+ madvise (raw,padded,MADV_POPULATE_WRITE);
30+ return static_cast <T*>(raw);}
31+ #endif
32+ return static_cast <T*>(::operator new (padded,std::align_val_t (align)));}
33+ void deallocate (T*p,std::size_t n)noexcept {if (!p)return ;
34+ std::size_t padded=round_up (n*sizeof (T));
35+ std::size_t align=std::max<std::size_t >(alignof (T),Align);
36+ #if CP_ALGO_USE_MMAP
37+ if (padded>=MEGABYTE){munmap (p,padded);return ;}
38+ #endif ::operator delete (p,padded,std::align_val_t (align));}
39+ private: static constexpr std::size_t MEGABYTE=1 <<20 ;
40+ static constexpr std::size_t round_up (std::size_t x)noexcept {return (x+Align-1 )/Align*Align;}};
41+ template <typename T>
42+ using big_vector=std::vector<T,big_alloc<T>>;}
43+ #line 4 "cp-algo/math/laurent.hpp"
44+ #include < memory>
45+ #include < optional>
46+ #line 7 "cp-algo/math/laurent.hpp"
47+ #include < cassert>
48+ #include < algorithm>
49+ namespace cp_algo ::math{template <typename T>
50+ struct provider {mutable big_vector<T>cache;
51+ mutable int cache_offset=0 ;
52+ mutable bool initialized=false ;
53+ mutable bool all_cached=false ;
54+ virtual ~provider ()=default ;
55+ virtual std::optional<int >valuation ()const {return std::nullopt ;}
56+ virtual std::optional<int >degree ()const {return std::nullopt ;}
57+ virtual T coeff_lazy (int k)const =0;
58+ virtual void double_up ()const {int old_size=cache.size ();
59+ int new_size=old_size==0 ?1 :2 *old_size;
60+ if (auto deg=degree ()){int max_idx=*deg-cache_offset;
61+ if (max_idx<new_size){new_size=max_idx+1 ;
62+ all_cached=true ;}}
63+ cache.resize (new_size);
64+ for (int i=old_size;i<new_size;i++){cache[i]=coeff_lazy (cache_offset+i);}}
65+ virtual T coeff (int k)const {if (!initialized){auto v=valuation ();
66+ cache_offset=v?*v:0 ;
67+ initialized=true ;}
68+ int idx=k-cache_offset;
69+ if (idx<0 ){return T (0 );}
70+ if (all_cached&&idx>=(int )cache.size ()){return T (0 );}
71+ while (idx>=(int )cache.size ()&&!all_cached){double_up ();}
72+ if (idx<(int )cache.size ()){return cache[idx];}
73+ return T (0 );}
74+ T get (int k)const {return coeff (k);}};
75+ template <typename T>
76+ struct constant_provider :provider<T>{T value;
77+ int offset;
78+ constant_provider (T value,int offset=0 ):value(value),offset(offset){}
79+ std::optional<int >valuation ()const override {return value!=T (0 )?std::optional<int >(offset):std::nullopt ;}
80+ std::optional<int >degree ()const override {return value!=T (0 )?std::optional<int >(offset):std::nullopt ;}
81+ T coeff_lazy (int k)const override {return k==offset?value:T (0 );}
82+ T coeff (int k)const override {return coeff_lazy (k);}};
83+ template <typename T>
84+ struct polynomial_provider :provider<T>{polynomial_provider(big_vector<T>coeffs,int offset=0 ){auto non_zero=[](const T&x){return x!=T (0 );};
85+ auto first=std::ranges::find_if (coeffs,non_zero);
86+ auto last=std::ranges::find_if (coeffs|std::views::reverse,non_zero);
87+ if (first!=coeffs.end ()){int leading_zeros=first-coeffs.begin ();
88+ int trailing_zeros=last-coeffs.rbegin ();
89+ coeffs=big_vector<T>(first,coeffs.end ()-trailing_zeros);
90+ offset+=leading_zeros;}else {coeffs.clear ();}
91+ this ->cache =std::move (coeffs);
92+ this ->cache_offset =offset;
93+ this ->initialized =true ;
94+ this ->all_cached =true ;}
95+ std::optional<int >valuation ()const override {if (this ->cache .empty ())return std::nullopt ;
96+ return this ->cache_offset ;}
97+ std::optional<int >degree ()const override {if (this ->cache .empty ())return std::nullopt ;
98+ return this ->cache_offset +(int )this ->cache .size ()-1 ;}
99+ T coeff_lazy (int k)const override {int idx=k-this ->cache_offset ;
100+ if (idx<0 ||idx>=(int )this ->cache .size ()){return T (0 );}
101+ return this ->cache [idx];}
102+ T coeff (int k)const override {return coeff_lazy (k);}};
103+ template <typename T>
104+ struct unary_provider :provider<T>{std::shared_ptr<provider<T>>operand;
105+ unary_provider (std::shared_ptr<provider<T>>operand):operand(std::move(operand)){}
106+ virtual T transform (T const &a)const =0;
107+ std::optional<int >valuation ()const override {return operand->valuation ();}
108+ std::optional<int >degree ()const override {return operand->degree ();}
109+ T coeff_lazy (int k)const override {return transform (operand->coeff_lazy (k));}
110+ T coeff (int k)const {return transform (operand->coeff (k));}};
111+ template <typename T>
112+ struct binary_provider :provider<T>{enum class bound_type {valuation,degree};
113+ std::shared_ptr<provider<T>>lhs,rhs;
114+ binary_provider (std::shared_ptr<provider<T>>lhs,std::shared_ptr<provider<T>>rhs):lhs(std::move(lhs)),rhs(std::move(rhs)){}
115+ virtual T combine (T const &a,T const &b)const =0;
116+ std::optional<int >valuation ()const override {auto lv=lhs->valuation ();
117+ auto rv=rhs->valuation ();
118+ if (!lv||!rv)return std::nullopt ;
119+ if (*lv!=*rv){return std::min (*lv,*rv);}
120+ T val=combine (lhs->coeff_lazy (*lv),rhs->coeff_lazy (*lv));
121+ if (val!=T (0 )){return *lv;}
122+ return std::nullopt ;}
123+ std::optional<int >degree ()const override {auto ld=lhs->degree ();
124+ auto rd=rhs->degree ();
125+ if (!ld||!rd)return std::nullopt ;
126+ if (*ld!=*rd){return std::max (*ld,*rd);}
127+ T val=combine (lhs->coeff_lazy (*ld),rhs->coeff_lazy (*ld));
128+ if (val!=T (0 )){return *ld;}
129+ return std::nullopt ;}
130+ T coeff_lazy (int k)const override {return combine (lhs->coeff_lazy (k),rhs->coeff_lazy (k));}
131+ T coeff (int k)const {return combine (lhs->coeff (k),rhs->coeff (k));}};
132+ template <typename T>
133+ struct add_provider :binary_provider<T>{using binary_provider<T>::binary_provider;
134+ T combine (T const &a,T const &b)const override {return a+b;}};
135+ template <typename T>
136+ struct subtract_provider :binary_provider<T>{using binary_provider<T>::binary_provider;
137+ T combine (T const &a,T const &b)const override {return a-b;}};
138+ template <typename T>
139+ struct negate_provider :unary_provider<T>{using unary_provider<T>::unary_provider;
140+ T transform (T const &a)const override {return -a;}};
141+ template <typename T>
142+ struct scale_provider :unary_provider<T>{T scalar;
143+ scale_provider (std::shared_ptr<provider<T>>operand,T scalar):unary_provider<T>(std::move(operand)),scalar(scalar){}
144+ T transform (T const &a)const override {return a*scalar;}};
145+ template <typename T>
146+ struct multiply_provider :provider<T>{std::shared_ptr<provider<T>>lhs,rhs;
147+ multiply_provider (std::shared_ptr<provider<T>>lhs,std::shared_ptr<provider<T>>rhs):lhs(std::move(lhs)),rhs(std::move(rhs)){}
148+ std::optional<int >valuation ()const override {auto lv=lhs->valuation ();
149+ auto rv=rhs->valuation ();
150+ if (lv&&rv)return *lv+*rv;
151+ return std::nullopt ;}
152+ std::optional<int >degree ()const override {auto ld=lhs->degree ();
153+ auto rd=rhs->degree ();
154+ if (ld&&rd)return *ld+*rd;
155+ return std::nullopt ;}
156+ T coeff_lazy (int k)const override {T result=T (0 );
157+ for (int j=0 ;j<=k;j++){result+=lhs->coeff (j)*rhs->coeff (k-j);}
158+ return result;}};
159+ template <typename T>
160+ struct laurent {std::shared_ptr<provider<T>>impl;
161+ laurent ():impl(std::make_shared<constant_provider<T>>(T(0 ),0 )){}
162+ laurent (T value,int offset=0 ):impl(std::make_shared<constant_provider<T>>(value,offset)){}
163+ laurent (big_vector<T>coeffs,int offset=0 ):impl(std::make_shared<polynomial_provider<T>>(std::move(coeffs),offset)){}
164+ laurent (std::shared_ptr<provider<T>>impl):impl(std::move(impl)){}
165+ T operator [](int k)const {return impl->get (k);}
166+ laurent operator -()const {return std::make_shared<negate_provider<T>>(impl);}
167+ laurent operator +(const laurent&other)const {return std::make_shared<add_provider<T>>(impl,other.impl );}
168+ laurent operator -(const laurent&other)const {return std::make_shared<subtract_provider<T>>(impl,other.impl );}
169+ laurent operator *(const laurent&other)const {return std::make_shared<multiply_provider<T>>(impl,other.impl );}
170+ laurent&operator +=(const laurent&other){return *this =*this +other;}
171+ laurent&operator -=(const laurent&other){return *this =*this -other;}
172+ laurent&operator *=(const laurent&other){return *this =*this *other;}
173+ laurent operator *(T const &scalar)const {return std::make_shared<scale_provider<T>>(impl,scalar);}
174+ laurent&operator *=(T const &scalar){return *this =*this *scalar;}};
175+ template <typename T>
176+ laurent<T>operator *(T const &scalar,laurent<T>const &series){return series*scalar;}}
0 commit comments