#{{{ header {.hints:off warnings:off optimization:speed.} import algorithm, sequtils, tables, macros, math, sets, strutils, strformat, sugar when defined(MYDEBUG): import header import streams proc scanf(formatstr: cstring){.header: "", varargs.} #proc getchar(): char {.header: "", varargs.} proc nextInt(): int = scanf("%lld",addr result) proc nextFloat(): float = scanf("%lf",addr result) proc nextString[F](f:F): string = var get = false result = "" while true: # let c = getchar() let c = f.readChar if c.int > ' '.int: get = true result.add(c) elif get: return proc nextInt[F](f:F): int = parseInt(f.nextString) proc nextFloat[F](f:F): float = parseFloat(f.nextString) proc nextString():string = stdin.nextString() template `max=`*(x,y:typed):void = x = max(x,y) template `min=`*(x,y:typed):void = x = min(x,y) template inf(T): untyped = when T is SomeFloat: T(Inf) elif T is SomeInteger: ((T(1) shl T(sizeof(T)*8-2)) - (T(1) shl T(sizeof(T)*4-1))) else: assert(false) proc discardableId[T](x: T): T {.discardable.} = return x macro `:=`(x, y: untyped): untyped = var strBody = "" if x.kind == nnkPar: for i,xi in x: strBody &= fmt""" {xi.repr} := {y[i].repr} """ else: strBody &= fmt""" when declaredInScope({x.repr}): {x.repr} = {y.repr} else: var {x.repr} = {y.repr} """ strBody &= fmt"discardableId({x.repr})" parseStmt(strBody) proc toStr[T](v:T):string = proc `$`[T](v:seq[T]):string = v.mapIt($it).join(" ") return $v proc print0(x: varargs[string, toStr]; sep:string):string{.discardable.} = result = "" for i,v in x: if i != 0: addSep(result, sep = sep) add(result, v) result.add("\n") stdout.write result var print:proc(x: varargs[string, toStr]) print = proc(x: varargs[string, toStr]) = discard print0(@x, sep = " ") template makeSeq(x:int; init):auto = when init is typedesc: newSeq[init](x) else: newSeqWith(x, init) macro Seq(lens: varargs[int]; init):untyped = var a = fmt"{init.repr}" for i in countdown(lens.len - 1, 0): a = fmt"makeSeq({lens[i].repr}, {a})" parseStmt(fmt""" block: {a}""") template makeArray(x:int; init):auto = var v:array[x, init.type] when init isnot typedesc: for a in v.mitems: a = init v macro Array(lens: varargs[typed], init):untyped = var a = fmt"{init.repr}" for i in countdown(lens.len - 1, 0): a = fmt"makeArray({lens[i].repr}, {a})" parseStmt(fmt""" block: {a}""") # }}} const Mod = 1000000007 # ModInt {{{ # ModInt[Mod] {{{ type ModInt[Mod: static[int]] = object v:int32 proc initModInt(a:SomeInteger, Mod:static[int]):ModInt[Mod] = var a = a.int a = a mod Mod if a < 0: a += Mod result.v = a.int32 proc getMod[Mod:static[int]](self: ModInt[Mod]):static int32 = self.Mod proc getMod[Mod:static[int]](self: typedesc[ModInt[Mod]]):static int32 = self.Mod macro declareModInt(Mod:static[int], t: untyped):untyped = var strBody = "" strBody &= fmt""" type {t.repr} = ModInt[{Mod.repr}] converter to{t.repr}(a:SomeInteger):{t.repr} = initModInt(a, {Mod.repr}) proc init{t.repr}(a:SomeInteger):{t.repr} = initModInt(a, {Mod.repr}) proc `$`(a:{t.repr}):string = $(a.v) """ parseStmt(strBody) when declared(Mod): declareModInt(Mod, Mint) ##}}} # DynamicModInt {{{ type DMint = object v:int32 proc setModSub(self:typedesc[not ModInt], m:int = -1, update = false):int32 = {.noSideEffect.}: var DMOD {.global.}:int32 if update: DMOD = m.int32 return DMOD proc fastMod(a:int,m:uint32):uint32{.inline.} = var minus = false a = a if a < 0: minus = true a = -a elif a < m.int: return a.uint32 var xh = (a shr 32).uint32 xl = a.uint32 d:uint32 asm """ "divl %4; \n\t" : "=a" (`d`), "=d" (`result`) : "d" (`xh`), "a" (`xl`), "r" (`m`) """ if minus and result > 0'u32: result = m - result proc initDMint(a:SomeInteger, Mod:int):DMint = result.v = fastMod(a.int, Mod.uint32).int32 proc getMod[T:not ModInt](self: T):int32 = T.type.setModSub() proc getMod(self: typedesc[not ModInt]):int32 = self.setModSub() proc setMod(self: typedesc[not ModInt], m:int) = discard self.setModSub(m, update = true) #}}} # Operations {{{ type ModIntC = concept x, type T x.v # x.v is int32 # x.getMod() is int32 # when T isnot ModInt: setMod(T, int) type SomeIntC = concept x x is SomeInteger or x is ModIntC proc Identity(self:ModIntC):auto = result = self;result.v = 1 proc init[Mod:static[int]](self:ModInt[Mod], a:SomeIntC):ModInt[Mod] = when a is SomeInteger: initModInt(a, Mod) else: a proc init(self:ModIntC and not ModInt, a:SomeIntC):auto = when a is SomeInteger: var r = self.type.default r.v = fastMod(a.int, self.getMod().uint32).int32 r else: a macro declareDMintConverter(t:untyped) = parseStmt(fmt""" converter to{t.repr}(a:SomeInteger):{t.repr} = let Mod = {t.repr}.getMod() if Mod > 0: result.v = fastMod(a.int, Mod.uint32).int32 else: result.v = a.int32 return result """) declareDMintConverter(DMint) macro declareDMint(t:untyped) = parseStmt(fmt""" type {t.repr} {{.borrow: `.`.}} = distinct DMint declareDMintConverter({t.repr}) """) proc `*=`(self:var ModIntC, a:SomeIntC) = when self is ModInt: self.v = (self.v.int * self.init(a).v.int mod self.getMod().int).int32 else: self.v = fastMod(self.v.int * self.init(a).v.int, self.getMod().uint32).int32 proc `==`(a:ModIntC, b:SomeIntC):bool = a.v == a.init(b).v proc `!=`(a:ModIntC, b:SomeIntC):bool = a.v != a.init(b).v proc `-`(self:ModIntC):auto = if self.v == 0: return self else: return self.init(self.getMod() - self.v) proc `$`(a:ModIntC):string = return $(a.v) proc `+=`(self:var ModIntC; a:SomeIntC) = self.v += self.init(a).v if self.v >= self.getMod(): self.v -= self.getMod() proc `-=`(self:var ModIntC, a:SomeIntC) = self.v -= self.init(a).v if self.v < 0: self.v += self.getMod() proc `^=`(self:var ModIntC, n:SomeInteger) = var (x,n,a) = (self,n,self.Identity) while n > 0: if (n and 1) > 0: a *= x x *= x n = (n shr 1) swap(self, a) proc inverse(self: ModIntC):auto = var a = self.v.int b = self.getMod().int u = 1 v = 0 while b > 0: let t = a div b a -= t * b;swap(a, b) u -= t * v;swap(u, v) return self.init(u) proc `/=`(a:var ModIntC,b:SomeIntC) = a *= a.init(b).inverse() proc `+`(a:ModIntC,b:SomeIntC):auto = result = a;result += b proc `-`(a:ModIntC,b:SomeIntC):auto = result = a;result -= b proc `*`(a:ModIntC,b:SomeIntC):auto = result = a;result *= b proc `/`(a:ModIntC,b:SomeIntC):auto = result = a;result /= b proc `^`(a:ModIntC,b:SomeInteger):auto = result = a;result ^= b # }}} # }}} #{{{ FastFourierTransform # clongdouble {{{ proc `+`(a, b:clongdouble):clongdouble {.importcpp: "(#) + (@)", nodecl.} proc `-`(a, b:clongdouble):clongdouble {.importcpp: "(#) - (@)", nodecl.} proc `*`(a, b:clongdouble):clongdouble {.importcpp: "(#) * (@)", nodecl.} proc `/`(a, b:clongdouble):clongdouble {.importcpp: "(#) / (@)", nodecl.} proc `-`(a:clongdouble):clongdouble {.importcpp: "-(#)", nodecl.} proc sqrt(a:clongdouble):clongdouble {.header: "", importcpp: "sqrtl(#)", nodecl.} proc exp(a:clongdouble):clongdouble {.header: "", importcpp: "expl(#)", nodecl.} proc sin(a:clongdouble):clongdouble {.header: "", importcpp: "sinl(#)", nodecl.} proc acos(a:clongdouble):clongdouble {.header: "", importcpp: "acosl(#)", nodecl.} proc cos(a:clongdouble):clongdouble {.header: "", importcpp: "cosl(#)", nodecl.} proc llround(a:clongdouble):int {.header: "", importcpp: "std::llround(#)", nodecl.} # }}} import math, sequtils, bitops #type Real = float type Real = clongdouble type Complex = tuple[x, y:Real] proc initComplex[S,T](x:S, y:T):Complex = (Real(x), Real(y)) proc `+`(a,b:Complex):Complex = initComplex(a.x + b.x, a.y + b.y) proc `-`(a,b:Complex):Complex = initComplex(a.x - b.x, a.y - b.y) proc `*`(a,b:Complex):Complex = initComplex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x) proc conj(a:Complex):Complex = initComplex(a.x, -a.y) type SeqComplex = object real, imag: seq[Real] proc initSeqComplex(n:int):SeqComplex = SeqComplex(real: newSeqWith(n, Real(0)), imag: newSeqWith(n, Real(0))) proc setLen(self: var SeqComplex, n:int) = self.real.setLen(n) self.imag.setLen(n) proc swap(self: var SeqComplex, i, j:int) = swap(self.real[i], self.real[j]) swap(self.imag[i], self.imag[j]) type FastFourierTransform = object of RootObj base:int rts: SeqComplex rev:seq[int] proc getC(self: SeqComplex, i:int):Complex = (self.real[i], self.imag[i]) proc `[]`(self: SeqComplex, i:int):Complex = self.getC(i) proc `[]=`(self: var SeqComplex, i:int, x:Complex) = self.real[i] = x.x self.imag[i] = x.y proc initFastFourierTransform():FastFourierTransform = return FastFourierTransform(base:1, rts: SeqComplex(real: @[Real(0), Real(1)], imag: @[Real(0), Real(0)]), rev: @[0, 1]) #proc init(self:typedesc[FastFourierTransform]):auto = initFastFourierTransform() proc ensureBase(self:var FastFourierTransform; nbase:int) = if nbase <= self.base: return let L = 1 shl nbase self.rev.setlen(1 shl nbase) self.rts.setlen(1 shl nbase) for i in 0..<(1 shl nbase): self.rev[i] = (self.rev[i shr 1] shr 1) + ((i and 1) shl (nbase - 1)) while self.base < nbase: let angle = acos(Real(-1)) * Real(2) / Real(1 shl (self.base + 1)) for i in (1 shl (self.base - 1))..<(1 shl self.base): self.rts[i shl 1] = self.rts[i] let angle_i = angle * Real(2 * i + 1 - (1 shl self.base)) self.rts[(i shl 1) + 1] = initComplex(cos(angle_i), sin(angle_i)) self.base.inc proc fft(self:var FastFourierTransform; a:var SeqComplex, n:int) = assert((n and (n - 1)) == 0) let zeros = countTrailingZeroBits(n) self.ensureBase(zeros) let shift = self.base - zeros for i in 0..0: fa[i shr 1].y else: fa[i shr 1].x) return ret var fft_t = initFastFourierTransform() #}}} const ArbitraryMod = true #{{{ ArbitraryModFFT type ArbitraryModFFT[ModInt] = object proc init[ModInt](t:typedesc[ArbitraryModFFT[ModInt]]):auto = ArbitraryModFFT[ModInt]() proc ceil_log2(n:int):int = result = 0 while (1 shl result) < n: result.inc proc fft[ModInt](self: var ArbitraryModFFT[ModInt], a:seq[ModInt]):SeqComplex = doAssert((a.len and (a.len - 1)) == 0) let l = ceil_log2(a.len) fft_t.ensureBase(l) result = initSeqComplex(a.len) for i in 0.. 0 and self[^1] == 0: discard self.pop() # operators +=, -=, *=, mod=, -, /= {{{ proc `+=`(self: var FormalPowerSeries, r:FormalPowerSeries) = if r.len > self.len: self.setlen(r.len) for i in 0.. self.len: self.setlen(r.len) for i in 0..= deg: continue if k notin r: r[k] = T(0) r[k] += c0 * c1 return toSeq(r.pairs) proc `*=`[T](self: var FormalPowerSeries[T], r: FormalPowerSeries[T]) = if self.len == 0 or r.len == 0: self.setlen(0) else: when FastMult: var fft = self.getFFT() self = fft[].multiply(self, r) else: assert(false) proc `mod=`[T](self: var FormalPowerSeries[T], r:FormalPowerSeries[T]) = self -= (self div r) * r proc `-`[T](self: FormalPowerSeries[T]):FormalPowerSeries[T] = var ret = self ret.applyIt(-it) return ret proc `/=`[T](self: var FormalPowerSeries[T], v:T) = self.applyIt(it / v) #}}} # operators +, -, *, div, mod {{{ macro declareOp(op) = fmt"""proc `{op}`[T](self:FormalPowerSeries[T];r:FormalPowerSeries[T] or T):FormalPowerSeries[T] = result = self;result {op}= r proc `{op}`[T](self: not FormalPowerSeries, r:FormalPowerSeries[T]):FormalPowerSeries[T] = result = initFormalPowerSeries[T](@[self]);result {op}= r""".parseStmt declareOp(`+`) declareOp(`-`) declareOp(`*`) declareOp(`/`) proc `div`[T](self, r:FormalPowerSeries[T]):FormalPowerSeries[T] = result = self;result.`div=` (r) proc `mod`[T](self, r:FormalPowerSeries[T]):FormalPowerSeries[T] = result = self;result.`mod=` (r) # }}} proc rev[T](self: FormalPowerSeries[T], deg = -1):auto = result = self if deg != -1: result.setlen(deg) result.reverse proc pre[T](self: FormalPowerSeries[T], sz:int):auto = result = self result.setlen(min(self.len, sz)) proc `div=`[T](self: var FormalPowerSeries[T], r: FormalPowerSeries[T]) = if self.len < r.len: self.setlen(0) else: let n = self.len - r.len + 1 self = (self.rev().pre(n) * r.rev().inv(n)).pre(n).rev(n) proc dot[T](self:FormalPowerSeries[T], r: FormalPowerSeries[T]):auto = result = initFormalPowerSeries[T](min(self.len, r.len)) for i in 0..= 1: result.delete(0, sz - 1) proc `shl`[T](self: FormalPowerSeries[T], sz:int):auto = result = initFormalPowerSeries[T](sz) result = result & self proc diff[T](self: FormalPowerSeries[T]):auto = let n = self.len result = initFormalPowerSeries[T](max(0, n - 1)) for i in 1.. 0: return FormalPowerSeries[T].none if deg - i div 2 <= 0: break var opt = (self shr i).sqrt(deg - i div 2) if not opt.isSome: return FormalPowerSeries[T].none var ret = opt.get shl (i div 2) if ret.len < deg: ret.setlen(deg) return ret.some return initFormalPowerSeries[T](deg).some var ret:FormalPowerSeries[T] if self.isSetSqrt: let opt = self.getSqrt()(self[0]) if not opt.isSome: return FormalPowerSeries[T].none ret = initFormalPowerSeries[T](@[T(opt.get)]) else: doAssert(self[0] == 1) ret = initFormalPowerSeries[T](@[T(1)]) let inv2 = T(1) / T(2); var i = 1 while i < deg: ret = (ret + self.pre(i shl 1) * ret.inv(i shl 1)) * inv2 i = i shl 1 return ret.pre(deg).some import typetraits # F(0) must be 0 proc exp[T](self: FormalPowerSeries[T], deg = -1):auto = doAssert self[0] == 0 deg.revise(self.len) when UseFFT: proc onlineConvolutionExp[T](self, conv_coeff:FormalPowerSeries[T]):auto = var fft = self.getFFT() let n = conv_coeff.len doAssert((n and (n - 1)) == 0) type FFTType = fft[].fft(initFormalPowerSeries[T](0)).type var conv_ntt_coeff = newSeq[FFTType]() i = n while (i shr 1) > 0: var g = conv_coeff.pre(i) conv_ntt_coeff.add(fft[].fft(g)) i = i shr 1 var conv_arg, conv_ret = initFormalPowerSeries[T](n) proc rec(l,r,d:int) = if r - l <= 16: for i in l.. deg: return initFormalPowerSeries[T](deg) result = (result shl (i * k)).pre(deg) if result.len < deg: result.setlen(deg) return return self proc eval[T](self: FormalPowerSeries[T], x:T):T = var r = T(0) w = T(1) for v in self: r += w * v w *= x return r proc powMod[T](self: FormalPowerSeries[T], n:int, M:FormalPowerSeries[T]):auto = let modinv = M.rev().inv() proc getDiv(base:FormalPowerSeries[T]):FormalPowerSeries[T] = if base.len < M.len: return initFormalPowerSeries[T](0) let n = base.len - M.len + 1 return (base.rev().pre(n) * modinv.pre(n)).pre(n).rev(n) var n = n x = self ret = initFormalPowerSeries[T](@[T(1)]) while n > 0: if (n and 1) > 0: ret *= x ret -= getDiv(ret) * M x *= x x -= getDiv(x) * M n = n shr 1 return ret # }}} # coef_of_generating_function {{{ #proc coef_of_generating_function(P,Q:FormalPowerSeries[FieldElem],N:int):auto = # assert Q[0] == 1 and Q.len == P.len + 1 # var (P, Q, N) = (P, Q, N) # while N > 0: # var Q1 = Q # for i in countup(0, Q.len - 1, 2): Q1[i] = Q[i] # for i in countup(1, Q.len - 1, 2): Q1[i] = -Q[i] # block: # var PQ1 = P * Q1 # P.setLen(0) # for i in countup(N mod 2, PQ1.len - 1, 2): P.add(PQ1[i]) # var QQ1 = Q * Q1 # Q.setLen(0) # for i in countup(0, QQ1.len - 1, 2): Q.add(QQ1[i]) # N = N div 2 # return P[0] proc coef_of_generating_function(P,Q:FormalPowerSeries[FieldElem],N:int):auto = assert Q[0] == 1 and Q.len == P.len + 1 let ix = (-Q) shr 1 let pm = powMod(ix, N, Q) let p = pm * P return (p mod Q) * Q.inv() # }}} #proc test() = # var # P = initFormalPowerSeries[Mint](@[2, 7, 1, 8]) # Q = initFormalPowerSeries[Mint](@[3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5]) # M = initFormalPowerSeries[Mint](@[1]).pre(Q.len - 1) # for i in 0..<10000: # dump(i) # M.setLen(Q.len - 1) # var M2 = powMod(P, i, Q) # M2.setLen(Q.len - 1) # for j in 0..= N: ans += p[i] * b[j] else: var P = @[Mint(1)] P.setLen(Q.len - 1) let ix = b shr 1 let base = N - b.len + 1 let r = coef_of_generating_function(P, Q, base) for i in 0..= N: continue for j in 0..= N: ans += r[i] * b[j] print ans