# Submit Info #2649

Problem Lang User Status Time Memory
Sqrt of Formal Power Series python3 maspy AC 3638 ms 178.16 MiB

ケース詳細
Name Status Time Memory
all_zero_00 AC 475 ms 129.33 MiB
all_zero_01 AC 535 ms 141.63 MiB
example_00 AC 197 ms 67.49 MiB
example_01 AC 200 ms 67.49 MiB
lower_deg_zero_00 AC 3566 ms 168.37 MiB
lower_deg_zero_01 AC 1880 ms 155.11 MiB
lower_deg_zero_02 AC 259 ms 75.86 MiB
lower_deg_zero_03 AC 267 ms 82.64 MiB
lower_deg_zero_04 AC 245 ms 67.50 MiB
lower_deg_zero_05 AC 258 ms 77.70 MiB
lower_deg_zero_06 AC 257 ms 74.21 MiB
lower_deg_zero_07 AC 250 ms 67.47 MiB
max_random_00 AC 271 ms 84.96 MiB
max_random_01 AC 3577 ms 176.59 MiB
max_random_02 AC 3536 ms 178.16 MiB
monomial_00 AC 1902 ms 150.99 MiB
monomial_01 AC 241 ms 67.50 MiB
monomial_02 AC 3638 ms 169.73 MiB
monomial_03 AC 242 ms 67.58 MiB
random_00 AC 270 ms 84.98 MiB
random_01 AC 3551 ms 175.66 MiB
random_02 AC 3540 ms 176.84 MiB

import sys read = sys.stdin.buffer.read readline = sys.stdin.buffer.readline readlines = sys.stdin.buffer.readlines import random import numpy as np MOD = 998244353 def modSqrt(a,p): # sqrt(a) mod p の1つを返す。存在しない場合は -1 を返す k = (p - 1) // 2 if p == 2: return a & 1 if a == 0: return 0 if pow(a,k,p) != 1: return -1 while True: # b^2 - a が非剰余となる元 b を乱択 b = random.randint(2,p-1) D = (b * b - a) % p if not D: return b if pow(D,k,p) == p-1: break k += 1 # (b + \sqrt{D})^k を計算する, # T^2 = D のもとで、(b+T)^k を計算する。 f0,f1 = b,1 # (b+T)^{2^n}を管理 g0,g1 = 1,0 while k: if k & 1: g0, g1 = f0 * g0 + D * f1 * g1, f1 * g0 + f0 * g1 f0, f1 = f0 * f0 + D * f1 * f1, 2 * f0 * f1 f0 %= p; f1 %= p; g0 %= p; g1 %= p k >>= 1 return g0 def cumprod(A, MOD = MOD): L = len(A); Lsq = int(L**.5+1) A = np.resize(A, Lsq**2).reshape(Lsq,Lsq) for n in range(1,Lsq): A[:,n] *= A[:,n-1]; A[:,n] %= MOD for n in range(1,Lsq): A[n] *= A[n-1,-1]; A[n] %= MOD return A.ravel()[:L] def make_fact(U, MOD = MOD): x = np.arange(U, dtype = np.int64); x[0] = 1 fact = cumprod(x, MOD) x = np.arange(U, 0, -1, dtype=np.int64); x[0] = pow(int(fact[-1]), MOD-2, MOD) fact_inv = cumprod(x, MOD)[::-1] return fact,fact_inv class Polynomial: MOD = MOD ModFact, ModFactInv = make_fact(10**6,MOD) ModInv = np.zeros(10**6, np.int64) ModInv[1:] = ModFactInv[1:] * ModFact[:-1] % MOD def __init__(self, array): assert array.dtype == np.int64 self.v = array @classmethod def zeros(cls, N): return cls(np.zeros(N,np.int64)) @classmethod def ones(cls, N): return cls(np.ones(N,np.int64)) @classmethod def arange(cls,N): return cls(np.arange(N,dtype=np.int64)) def __getitem__(self,item): if type(item) == slice: return Polynomial(self.v[item]) else: return self.v[item] def __repr__(self): return 'Polynomial modulo {}\n'.format(Polynomial.MOD) + self.v.__repr__() def print(self, sep = ' '): print(sep.join(self.v.astype(str))) def __eq__(self,other): return np.all(self.v == other.v) def __len__(self): return len(self.v) def __add__(self, other): Ls = len(self.v); Lo = len(other.v) if Ls < Lo: Ls, Lo = Lo, Ls; self,other = other,self array = self.v.copy(); array[:Lo] += other.v; array[:Lo] %= Polynomial.MOD return Polynomial(array) def __iadd__(self, other): Ls = len(self.v); Lo = len(other.v) if Ls < Lo: self.resize(Lo) self.v[:Lo] += other.v; self.v %= Polynomial.MOD; return self def __sub__(self, other): Ls = len(self); Lo = len(other) if Ls < Lo: Ls, Lo = Lo, Ls; self,other = other,self array = self.v.copy(); array[:Lo] -= other.v; array[:Lo] %= Polynomial.MOD return Polynomial(array) def __isub__(self, other): Ls = len(self.v); Lo = len(other.v) if Ls < Lo: self.resize(Lo) self.v[:Lo] -= other.v; self.v %= Polynomial.MOD; return self def __mul__(self, other): f = self.v; g = other.v Lf = len(f); Lg = len(g); L = Lf + Lg - 1 if not Lf: return Polynomial.zeros(1) if not Lg: return Polynomial.zeros(1) fl = f & (1 << 15) - 1; fh = f >> 15 gl = g & (1 << 15) - 1; gh = g >> 15 if L < (1<<8): conv = lambda f,g: \ np.convolve(f,g) % Polynomial.MOD else: fft = np.fft.rfft; ifft = np.fft.irfft; fft_len = 1 << L.bit_length() conv = lambda f,g: \ np.rint(ifft(fft(f,fft_len) * fft(g,fft_len))[:L]).astype(np.int64) % Polynomial.MOD x = conv(fl, gl); z = conv(fh, gh); y = conv(fl+fh, gl+gh) return Polynomial((x + ((y - x - z) << 15) + (z << 30)) % Polynomial.MOD) def __truediv__(self,other): L = len(other) return (self * other.inv())[:L] def resize(self,N): L = len(self.v); self.v = np.resize(self.v,N); self.v[L:] = 0; return self def inv(self): F = self; L = len(F); x = int(F[0]); assert x != 0 i = (L-1).bit_length(); N = 1 << i; F.resize(N) x = pow(x,MOD-2,MOD); R = Polynomial(np.int64([x])) n = 1 for _ in range(i): RF = (R * F[:n+n])[:n+n] RRF = (R * RF)[:n+n] R += R; R -= RRF n += n F.resize(L) return R[:L] def differentiate(self): L = len(self) if L == 1: return Polynomial.zeros(1) return Polynomial(self.v[1:] * np.arange(1,L,dtype=np.int64) % MOD) def integrate(self): L = len(self) A = np.zeros(L+1,np.int64) A[1:] = self.v * Polynomial.ModInv[1:L+1] % MOD return Polynomial(A) def exp(self): # http://www.csd.uwo.ca/~eschost/publications/BoSc09-final.pdf h = self; L = len(h) assert h[0] == 0 dh = h.differentiate() i = (L-1).bit_length(); N = 1 << i; h.resize(N) f = Polynomial.ones(1); g = Polynomial.ones(1) m = 1 for _ in range(i): fg = (f * g)[:m]; fgg = (fg * g)[:m] g+= g; g -= fgg q = dh[:m-1] df = f.differentiate() df -= (f * q)[:m+m-1] df *= g; w = df[:m+m-1] + q w = w.integrate() f = (h - w + Polynomial.ones(1)) * f; f = f[:m+m] m += m return f[:L] def log(self): return (self.differentiate() / self).integrate()[:len(self)] def sqrt_1(self): F = self; assert F[0] == 1 assert Polynomial.MOD != 2 L = len(F) i = (L-1).bit_length(); N = 1 << i; F.resize(N) n = 1 R = Polynomial.ones(1) for _ in range(i+1): Ri = R.inv() FRi = (F[:n+n] * Ri)[:n+n] FRi += R; R = FRi R.v *= Polynomial.ModInv[2]; R.v %= Polynomial.MOD n += n F.resize(L) return R[:L] def sqrt_any(self): f = self; L = len(f) if np.all(f.v == 0): return f D = np.where(f.v != 0)[0][0] if D & 1: raise ValueError y = f[D] x = modSqrt(int(y),Polynomial.MOD) if x == -1: raise ValueError invy = pow(int(y),MOD-2,MOD) g = Polynomial(f.v[D:] * invy % Polynomial.MOD) h = Polynomial.zeros(L) h.v[D//2:D//2+len(g)] = g.sqrt_1().v * x % Polynomial.MOD return h N = int(readline()) f = Polynomial(np.int64(read().split())) try: g = f.sqrt_any() g.print() except(ValueError): print(-1)