import java.io.IOException; import java.io.InputStream; import java.io.OutputStream; import java.io.PrintStream; import java.io.PrintWriter; import java.lang.reflect.Array; import java.math.BigInteger; import java.util.ArrayDeque; import java.util.ArrayList; import java.util.Arrays; import java.util.Collection; import java.util.Collections; import java.util.Comparator; import java.util.HashMap; import java.util.HashSet; import java.util.Iterator; import java.util.List; import java.util.Map.Entry; import java.util.Map; import java.util.NoSuchElementException; import java.util.Objects; import java.util.Queue; import java.util.Random; import java.util.Set; import java.util.TreeMap; import java.util.function.BiFunction; import java.util.function.DoubleUnaryOperator; import java.util.function.IntBinaryOperator; import java.util.function.IntFunction; import java.util.function.IntToDoubleFunction; import java.util.function.IntToLongFunction; import java.util.function.IntUnaryOperator; import java.util.function.LongBinaryOperator; import java.util.function.LongToDoubleFunction; import java.util.function.LongUnaryOperator; import java.util.function.Predicate; import java.util.function.ToIntFunction; import java.util.random.RandomGenerator; import java.util.stream.IntStream; import java.util.stream.LongStream; import java.util.stream.Stream; public class Main { static MyPrintWriter pw = MyPrintWriter.getInstance(); static FastScanner sc = FastScanner.getInstance(); public static void main(String[] args) throws IOException { Thread.setDefaultUncaughtExceptionHandler((t, e) -> System.exit(1)); new Main().run(); pw.flush(); } void run() { long M = sc.nextLong(); int K = sc.nextInt(); long R = sc.nextLong(); long ans = f(M, K, R) - f(M, K, 0); ans = fp.reduce(ans); pw.println(ans); } Fp fp = Fp.MOD998244353; long f(long M, int K, long R) { var P = PolynomialFpDynamic.MOD998244353; long[] A = new long[K + 1]; A[0] = 1; for (int i = 0; i < K; i++) { A[i + 1] = fp.reduce(-fp.inv(K)); } var BmodA = P.geometricSumMod(M, A); BmodA[0] = fp.reduce(BmodA[0] - M); BmodA = P.mod(BmodA, A); long[] BleqR_modA = P.geometricSumMod(R + 1, A); BleqR_modA[0] = fp.reduce(BleqR_modA[0] - M); var ans = P.solveCyclicCongruencePoint(A, M, BmodA, BleqR_modA, R); return ans; } } class ArrayUtils { public static void swap(int i, int j, long[] A) { if (i == j) { return; } long tmp = A[i]; A[i] = A[j]; A[j] = tmp; } public static void swap(long[] A, long[] B) { if (A.length != B.length) { throw new AssertionError(); } for (int i = 0; i < A.length; i++) { long tmp = A[i]; A[i] = B[i]; B[i] = tmp; } } public static void reverse(long[] a) { int s = 0; int t = a.length - 1; while (s < t) { swap(s, t, a); ++s; --t; } } } interface CommutativeRingStrategy extends RingStrategy {} interface EuclideanDomainStrategy extends GCDDomainStrategy { T div(T a, T b); T mod(T a, T b); long norm(T a); default T canonicalUnit(T a) { return one(); } @Override default T gcd(T a, T b) { while (!equals(b, zero())) { a = mod(a, b); T t = a; a = b; b = t; } if (equals(a, zero())) { return a; } return div(a, canonicalUnit(a)); } record ExtGCDResult(T x, T y, T gcd) {} default ExtGCDResult extgcd(T a, T b) { T x0 = one(); T y0 = zero(); T g0 = a; T x1 = zero(); T y1 = one(); T g1 = b; while (!equals(g1, zero())) { T q = div(g0, g1); T nextG = sub(g0, mul(q, g1)); T nextX = sub(x0, mul(q, x1)); T nextY = sub(y0, mul(q, y1)); x0 = x1; y0 = y1; g0 = g1; x1 = nextX; y1 = nextY; g1 = nextG; } if (equals(g0, zero())) { return new ExtGCDResult<>(x0, y0, g0); } T u = canonicalUnit(g0); return new ExtGCDResult<>(div(x0, u), div(y0, u), div(g0, u)); } } interface ExactDivRingStrategy extends IntegralDomainStrategy { T exactDiv(T a, T b); } class FastScanner { private static FastScanner instance = null; private final InputStream in = System.in; private final byte[] buffer = new byte[1 << 16]; private int ptr = 0; private int buflen = 0; private FastScanner() { } public static FastScanner getInstance() { if (instance == null) { instance = new FastScanner(); } return instance; } private boolean hasNextByte() { if (ptr < buflen) { return true; } ptr = 0; try { buflen = in.read(buffer); } catch (IOException e) { e.printStackTrace(); } return buflen > 0; } private int readByte() { if (hasNextByte()) { return buffer[ptr++]; } else { return -1; } } private boolean isPrintableChar(int c) { return (33 <= c) && (c <= 126); } public boolean hasNext() { while (hasNextByte() && (!isPrintableChar(buffer[ptr]))) { ptr++; } return hasNextByte(); } public long nextLong() { if (!hasNext()) { throw new NoSuchElementException(); } long n = 0; boolean minus = false; int b = readByte(); if (b == '-') { minus = true; b = readByte(); } while ((b >= '0') && (b <= '9')) { n = ((n << 1) + (n << 3)) + (b - '0'); b = readByte(); } return minus ? -n : n; } public int nextInt() { return ((int) (nextLong())); } } class Fp extends Zn implements LongFieldStrategy { public static final Fp MOD998244353 = new Fp(998244353); public Fp(long mod) { super(mod); } int[] inv = new int[0]; public long inv(long n) { if (n < 0) { n = reduce(n); } return n < inv.length ? inv[((int) (n))] & 0xffffffffL : MathUtils.modInv(n, mod); } @Override public long div(long a, long b) { return mul(a, inv(b)); } @Override public long mod(long a, long b) { if (equals(b, zero())) { throw new ArithmeticException("Division by zero"); } return zero(); } @Override public long norm(long a) { return equals(a, zero()) ? 0 : 1; } @Override public long canonicalUnit(long a) { if (equals(a, zero())) { return one(); } return a; } public long pow(long a, long n) { if (n < 0) { a = inv(a); n = -n; } return MathUtils.modPow(a, n, mod); } } interface GCDDomainStrategy extends IntegralDomainStrategy { T gcd(T a, T b); } interface IntegralDomainStrategy extends CommutativeRingStrategy {} interface LongCommutativeRingStrategy extends LongRingStrategy {} interface LongEuclideanDomainStrategy extends LongGCDDomainStrategy { long div(long a, long b); long mod(long a, long b); long norm(long a); default long canonicalUnit(long a) { return one(); } @Override default long gcd(long a, long b) { while (!equals(b, zero())) { a = mod(a, b); long t = a; a = b; b = t; } if (equals(a, zero())) { return a; } return div(a, canonicalUnit(a)); } record ExtGCDResult(long x, long y, long gcd) {} default ExtGCDResult extgcd(long a, long b) { long x0 = one(); long y0 = zero(); long g0 = a; long x1 = zero(); long y1 = one(); long g1 = b; while (!equals(g1, zero())) { long q = div(g0, g1); long nextG = sub(g0, mul(q, g1)); long nextX = sub(x0, mul(q, x1)); long nextY = sub(y0, mul(q, y1)); x0 = x1; y0 = y1; g0 = g1; x1 = nextX; y1 = nextY; g1 = nextG; } if (equals(g0, zero())) { return new ExtGCDResult(x0, y0, g0); } long u = canonicalUnit(g0); return new ExtGCDResult(div(x0, u), div(y0, u), div(g0, u)); } } interface LongExactDivRingStrategy extends LongIntegralDomainStrategy { long divExact(long a, long b); } interface LongFieldStrategy extends LongEuclideanDomainStrategy , LongExactDivRingStrategy { @Override default long divExact(long a, long b) { return div(a, b); } long inv(long a); default long div(long a, long b) { return mul(a, inv(b)); } default long geometricSum(long a) { return inv(sub(one(), a)); } @Override default ExtGCDResult extgcd(long a, long b) { if (!equals(a, zero())) { return new ExtGCDResult(inv(a), zero(), one()); } else if (!equals(b, zero())) { return new ExtGCDResult(zero(), inv(b), one()); } else { return new ExtGCDResult(zero(), zero(), zero()); } } } interface LongGCDDomainStrategy extends LongIntegralDomainStrategy { long gcd(long a, long b); } interface LongIntegralDomainStrategy extends LongCommutativeRingStrategy {} interface LongRingStrategy extends LongSemiRingStrategy { long neg(long a); default long sub(long a, long b) { return add(a, neg(b)); } } interface LongSemiRingStrategy { long zero(); long one(); long add(long a, long b); long mul(long a, long b); boolean equals(long a, long b); } class MathUtils { public static long modPow(long a, long n, long mod) { if (n < 0) { long inv = MathUtils.modInv(a, mod); return MathUtils.modPow(inv, -n, mod); } if (n == 0) { return 1; } return (MathUtils.modPow((a * a) % mod, n / 2, mod) * ((n % 2) == 1 ? a : 1)) % mod; } public static long modInv(long a, long mod) { a = ((a % mod) + mod) % mod; long[] f0 = new long[]{ 1, 0, mod }; long[] f1 = new long[]{ 0, 1, a }; while (f1[2] != 0) { long q = f0[2] / f1[2]; for (int i = 0; i < 3; i++) { f0[i] -= q * f1[i]; } ArrayUtils.swap(f0, f1); } return f0[1] < 0 ? mod + f0[1] : f0[1]; } } class MyPrintWriter extends PrintWriter { private static MyPrintWriter instance = null; private MyPrintWriter() { super(System.out); } public static MyPrintWriter getInstance() { if (instance == null) { instance = new MyPrintWriter(); } return instance; } } class PolynomialFpDynamic implements EuclideanDomainStrategy , UFDStrategy , ExactDivRingStrategy { public final long mod; Fp fp; public final boolean isNTTFriendly; public final long primitiveRoot; public final int maxPow2; long[][] bitreversedRoots; long[][] bitreversedInvRoots; public static final int FFT_NAIVE_THRESHOLD = 128; public static final int FFT_MIN_LENGTH_THRESHOLD = 10; public static final PolynomialFpDynamic MOD998244353 = new PolynomialFpDynamic(998244353L, 3); public static final PolynomialFpDynamic MOD469762049 = new PolynomialFpDynamic(469762049L, 3); public static final PolynomialFpDynamic MOD167772161 = new PolynomialFpDynamic(167772161L, 3); public int countTerms(long[] a, int limit) { int count = 0; for (long v : a) { if (fp.reduce(v) != 0) { count++; if (count > limit) { return count; } } } return count; } public PolynomialFpDynamic(long mod, long primitiveRoot) { this.mod = mod; fp = new Fp(mod); this.isNTTFriendly = true; this.primitiveRoot = primitiveRoot; this.maxPow2 = Long.numberOfTrailingZeros(mod - 1); this.bitreversedRoots = new long[maxPow2 + 1][]; this.bitreversedInvRoots = new long[maxPow2 + 1][]; } long addMod(long a, long b) { long sum = a + b; return sum >= mod ? sum - mod : sum; } long subMod(long a, long b) { long diff = a - b; return diff < 0 ? diff + mod : diff; } void prepareRoots(int n) { if (Integer.bitCount(n) != 1) { throw new AssertionError(); } int sz = Integer.numberOfTrailingZeros(n); if (sz > maxPow2) { throw new AssertionError("NTT length exceeds mod - 1 power of two"); } if (bitreversedRoots[sz] != null) { return; } long root = MathUtils.modPow(primitiveRoot, (mod - 1) / n, mod); long iroot = MathUtils.modInv(root, mod); bitreversedRoots[sz] = new long[n]; bitreversedInvRoots[sz] = new long[n]; for (int n_ = n / 2; n_ >= 1; n_ /= 2 , root = (root * root) % mod , iroot = (iroot * iroot) % mod) { long w = 1; long iw = 1; for (int j = 0; j < n_; ++j) { bitreversedRoots[sz][n_ + j] = w; bitreversedInvRoots[sz][n_ + j] = iw; w = (w * root) % mod; iw = (iw * iroot) % mod; } int cur = 0; for (int j = 0; j < n_; ++j) { if (cur < j) { ArrayUtils.swap(n_ + cur, n_ + j, bitreversedRoots[sz]); ArrayUtils.swap(n_ + cur, n_ + j, bitreversedInvRoots[sz]); } for (int k = n_ / 2; k > (cur ^= k); k /= 2); } } } public void fftToBitReversed(long[] a) { int n = a.length; int sz = Integer.numberOfTrailingZeros(n); prepareRoots(n); for (int m = 1, t = n / 2; m <= (n / 2); m *= 2 , t /= 2) { for (int i = 0, k = 0; i < m; ++i , k += 2 * t) { long s = bitreversedRoots[sz][m + i]; for (int j = k; j < (k + t); ++j) { long u = a[j]; long v = (a[j + t] * s) % mod; a[j] = addMod(u, v); a[j + t] = subMod(u, v); } } } } public void ifftFromBitReversed(long[] a) { long invN = MathUtils.modInv(a.length, mod); int n = a.length; int sz = Integer.numberOfTrailingZeros(n); prepareRoots(n); for (int m = n / 2, t = 1; m >= 1; m /= 2 , t *= 2) { for (int i = 0, k = 0; i < m; ++i , k += 2 * t) { long s = bitreversedInvRoots[sz][m + i]; if (m == 1) { s = (s * invN) % mod; } for (int j = k; j < (k + t); ++j) { long u = a[j]; long v = a[j + t]; if (m == 1) { a[j] = ((u + v) * invN) % mod; } else { a[j] = addMod(u, v); } a[j + t] = (((u + mod) - v) * s) % mod; } } } } public long[] mulFFT(long[] a, long[] b) { if ((a.length == 0) || (b.length == 0)) { return new long[0]; } int n = 1; int len = (a.length + b.length) - 1; while (n < len) { n *= 2; } if (Integer.numberOfTrailingZeros(n) > maxPow2) { throw new AssertionError("NTT length exceeds mod - 1 power of two"); } long[] fa = new long[n]; long[] fb = new long[n]; for (int i = 0; i < a.length; i++) { fa[i] = fp.reduce(a[i]); } for (int i = 0; i < b.length; i++) { fb[i] = fp.reduce(b[i]); } prepareRoots(n); fftToBitReversed(fa); fftToBitReversed(fb); for (int i = 0; i < n; ++i) { fa[i] = (fa[i] * fb[i]) % mod; } ifftFromBitReversed(fa); return Arrays.copyOf(fa, len); } public long[] mulNaive(long[] a, long[] b) { long[] c = new long[(a.length + b.length) - 1]; for (int i = 0; i < a.length; i++) { if (a[i] == 0) { continue; } for (int j = 0; j < b.length; j++) { if (b[j] == 0) { continue; } c[i + j] = (c[i + j] + (a[i] * b[j])) % mod; } } return c; } private long[] mulCRT(long[] a, long[] b) { int n = a.length; int m = b.length; long m1 = 998244353L; long m2 = 469762049L; long m3 = 167772161L; double maxVal = (((double) (Math.min(n, m))) * (mod - 1)) * (mod - 1); int k_count = 3; if (maxVal < m1) { k_count = 1; } else if (maxVal < (((double) (m1)) * m2)) { k_count = 2; } long[] res1 = MOD998244353.mulFFT(a, b); if (k_count == 1) { for (int i = 0; i < res1.length; i++) { res1[i] %= mod; } return res1; } long[] res2 = MOD469762049.mulFFT(a, b); if (k_count == 2) { int len = res1.length; long[] res = new long[len]; long[] v = new long[2]; long[] ms = new long[]{ m1, m2 }; for (int i = 0; i < len; i++) { v[0] = res1[i]; v[1] = res2[i]; res[i] = Zn.crt(v, ms) % mod; } return res; } long[] res3 = MOD167772161.mulFFT(a, b); int len = res1.length; long[] res = new long[len]; long inv123 = MathUtils.modInv(((m1 % m3) * (m2 % m3)) % m3, m3); long m12m = ((m1 % mod) * (m2 % mod)) % mod; long[] v = new long[2]; long[] ms = new long[]{ m1, m2 }; for (int i = 0; i < len; i++) { v[0] = res1[i]; v[1] = res2[i]; long x12 = Zn.crt(v, ms); long k3 = ((((res3[i] - (x12 % m3)) + m3) % m3) * inv123) % m3; res[i] = ((x12 % mod) + (((k3 % mod) * m12m) % mod)) % mod; } return res; } @Override public long[] mul(long[] a, long[] b) { if ((a.length == 0) || (b.length == 0)) { return new long[0]; } if ((a.length == 1) && (b.length == 1)) { return new long[]{ (a[0] * b[0]) % mod }; } int n = a.length; int m = b.length; if (((n + m) - 1) <= 128) { return mulNaive(a, b); } int sparseThreshold = 50; int countB = countTerms(b, sparseThreshold); if (countB <= sparseThreshold) { return sparseMul(a, getTerms(b, countB), b.length); } int countA = countTerms(a, sparseThreshold); if (countA <= sparseThreshold) { return sparseMul(b, getTerms(a, countA), a.length); } if ((isNTTFriendly && (((n + m) - 1) > FFT_NAIVE_THRESHOLD)) && (Math.min(n, m) > FFT_MIN_LENGTH_THRESHOLD)) { return mulFFT(a, b); } else if ((!isNTTFriendly) && (((n + m) - 1) > 128)) { return mulCRT(a, b); } return mulNaive(a, b); } public long[] squared(long[] a) { if (a.length == 0) { return new long[0]; } if (a.length == 1) { return new long[]{ (a[0] * a[0]) % mod }; } int len = (2 * a.length) - 1; if (isNTTFriendly && (len > FFT_NAIVE_THRESHOLD)) { return squaredFFT(a); } else if ((!isNTTFriendly) && (len > 128)) { return mulCRT(a, a); } return squaredNaive(a); } private long[] squaredFFT(long[] a) { if (a.length == 0) { return new long[0]; } int n = 1; int len = (2 * a.length) - 1; while (n < len) { n *= 2; } if (Integer.numberOfTrailingZeros(n) > maxPow2) { throw new AssertionError("NTT length exceeds mod - 1 power of two"); } long[] fa = new long[n]; for (int i = 0; i < a.length; i++) { fa[i] = fp.reduce(a[i]); } prepareRoots(n); fftToBitReversed(fa); for (int i = 0; i < n; ++i) { fa[i] = (fa[i] * fa[i]) % mod; } ifftFromBitReversed(fa); return Arrays.copyOf(fa, len); } public long[] squaredNaive(long[] a) { int len = (2 * a.length) - 1; long[] ret = new long[len]; for (int i = 0; i < a.length; ++i) { if (a[i] == 0) { continue; } for (int j = i + 1; j < a.length; ++j) { if (a[j] == 0) { continue; } ret[i + j] = (ret[i + j] + ((2 * a[i]) * a[j])) % mod; } } for (int i = 0; i < a.length; ++i) { if (a[i] == 0) { continue; } ret[2 * i] = (ret[2 * i] + (a[i] * a[i])) % mod; } return ret; } @Override public long[] zero() { return new long[0]; } @Override public long[] one() { return new long[]{ 1 }; } @Override public long[] add(long[] a, long[] b) { long[] ret = new long[Math.max(a.length, b.length)]; for (int i = 0; i < ret.length; ++i) { ret[i] = (i < a.length ? a[i] : 0) + (i < b.length ? b[i] : 0); if (ret[i] >= mod) { ret[i] -= mod; } } return ret; } @Override public long[] neg(long[] a) { long[] ret = new long[a.length]; for (int i = 0; i < a.length; i++) { if (a[i] != 0) { ret[i] = mod - a[i]; } } return ret; } @Override public boolean equals(long[] a, long[] b) { int n = Math.max(a.length, b.length); for (int i = 0; i < n; i++) { long va = (i < a.length) ? a[i] : 0; long vb = (i < b.length) ? b[i] : 0; if (va != vb) { return false; } } return true; } @Override public boolean equals(Object o) { if (this == o) { return true; } if ((o == null) || (getClass() != o.getClass())) { return false; } PolynomialFpDynamic that = ((PolynomialFpDynamic) (o)); return mod == that.mod; } @Override public int hashCode() { return Objects.hash(mod); } @Override public long norm(long[] a) { return deg(a) + 1; } @Override public long[] canonicalUnit(long[] a) { int d = deg(a); if (d == (-1)) { return one(); } return new long[]{ a[d] }; } @Override public long[] sub(long[] a, long[] b) { if ((a.length == 1) && (b.length == 1)) { long diff = a[0] - b[0]; if (diff < 0) { diff += mod; } return new long[]{ diff }; } long[] ret = new long[Math.max(a.length, b.length)]; for (int i = 0; i < ret.length; ++i) { ret[i] = (i < a.length ? a[i] : 0) - (i < b.length ? b[i] : 0); if (ret[i] < 0) { ret[i] += mod; } } return ret; } public int deg(long[] a) { int n = a.length; if (n == 0) { return -1; } if (a[n - 1] != 0) { return n - 1; } for (int i = n - 2; i >= 0; i--) { if (a[i] != 0) { return i; } } return -1; } public boolean isZero(long[] f) { return deg(f) == (-1); } public long[] resize(long[] a) { return Arrays.copyOf(a, Math.max(0, deg(a)) + 1); } public long[] monic(long[] a) { int deg = deg(a); if (deg == (-1)) { return new long[0]; } if (deg == 0) { return new long[]{ 1 }; } long inv = fp.inv(a[deg]); long[] res = new long[deg + 1]; for (int i = 0; i < deg; i++) { res[i] = (a[i] * inv) % mod; } res[deg] = 1; return res; } public long[] divNaive(long[] a, long[] b) { int degA = deg(a); int degB = deg(b); if (degB == (-1)) { throw new ArithmeticException("division by zero polynomial"); } if (degA < degB) { return new long[]{ 0 }; } long[] r = Arrays.copyOf(a, degA + 1); long[] q = new long[(degA - degB) + 1]; long invB = fp.inv(b[degB]); for (int i = degA; i >= degB; i--) { if (r[i] == 0) { continue; } long c = (r[i] * invB) % mod; q[i - degB] = c; for (int j = 0; j <= degB; j++) { r[(i - degB) + j] -= (c * b[j]) % mod; if (r[(i - degB) + j] < 0) { r[(i - degB) + j] += mod; } } } return resize(q); } public long[] modNaive(long[] a, long[] b) { int degA = deg(a); int degB = deg(b); if (degB == (-1)) { throw new ArithmeticException("division by zero polynomial"); } if (degA < degB) { return resize(a); } long[] r = Arrays.copyOf(a, degA + 1); long invB = fp.inv(b[degB]); for (int i = degA; i >= degB; i--) { if (r[i] == 0) { continue; } long c = (r[i] * invB) % mod; for (int j = 0; j <= degB; j++) { r[(i - degB) + j] -= (c * b[j]) % mod; if (r[(i - degB) + j] < 0) { r[(i - degB) + j] += mod; } } } return resize(r); } public static class DivModResult { public long[] q; public long[] r; public DivModResult(long[] q, long[] r) { this.q = q; this.r = r; } } public DivModResult divmod(long[] a, long[] b) { var q = div(a, b); var r = sub(a, mul(q, b)); r = resize(r); return new DivModResult(q, r); } public long[] gcdNaive(long[] a, long[] b) { a = resize(a); b = resize(b); while (deg(b) != (-1)) { long[] r = modNaive(a, b); a = b; b = r; } return monic(a); } public long[] differentiate(long[] a) { long[] ret = new long[a.length]; for (int i = 1; i < a.length; ++i) { ret[i - 1] = (i * a[i]) % mod; } return ret; } public long[] integrate(long[] a) { long[] ret = new long[a.length]; for (int i = 0; (i + 1) < a.length; ++i) { ret[i + 1] = (MathUtils.modInv(i + 1, mod) * a[i]) % mod; } return ret; } public long[] invNaive(long[] a) { long[] g = new long[a.length]; long inv0 = MathUtils.modInv(a[0], mod); g[0] = inv0; for (int i = 1; i < a.length; i++) { long sum = 0; for (int j = 1; j <= i; j++) { if (j < a.length) { sum = (sum + (a[j] * g[i - j])) % mod; } } g[i] = (sum == 0) ? 0 : (((mod - sum) % mod) * inv0) % mod; } return g; } public long[] invFFT(long[] a) { long[] g = new long[]{ MathUtils.modInv(a[0], mod) }; for (int len = 1; len < a.length; len *= 2) { long[] fftG = Arrays.copyOf(g, len * 4); long[] fftA = new long[4 * len]; System.arraycopy(a, 0, fftA, 0, Math.min(2 * len, a.length)); prepareRoots(4 * len); fftToBitReversed(fftG); fftToBitReversed(fftA); for (int i = 0; i < fftG.length; ++i) { fftG[i] = (((fftG[i] * fftG[i]) % mod) * fftA[i]) % mod; } ifftFromBitReversed(fftG); for (int i = 0; i < len; ++i) { fftG[i] = g[i]; } for (int i = len; i < (2 * len); ++i) { if (fftG[i] != 0) { fftG[i] = mod - fftG[i]; } } g = Arrays.copyOf(fftG, Math.min(a.length, 2 * len)); } return g; } public long[] inv(long[] a) { return isNTTFriendly ? invFFT(a) : invNaive(a); } public long[] log(long[] a) { return integrate(Arrays.copyOf(mul(differentiate(a), inv(a)), a.length)); } public long[] expNaive(long[] a) { if (a[0] != 0) { throw new AssertionError(); } int n = a.length; long[] g = new long[n]; long[] df = new long[n]; for (int i = 0; i < n; i++) { df[i] = (i * a[i]) % mod; } g[0] = 1; for (int i = 1; i < n; i++) { long sum = 0; for (int j = 1; j <= i; j++) { sum = (sum + (df[j] * g[i - j])) % mod; } g[i] = (sum * fp.inv(i)) % mod; } return g; } public long[] expFFT(long[] a) { if (a[0] != 0) { throw new AssertionError(); } long[] g = new long[]{ 1 }; for (int len = 1; len < a.length; len *= 2) { long[] tmp = sub(Arrays.copyOf(a, Math.min(2 * len, a.length)), log(Arrays.copyOf(g, Math.min(2 * len, a.length)))); tmp[0] = addMod(tmp[0], 1); g = Arrays.copyOf(mul(g, tmp), Math.min(2 * len, a.length)); } return Arrays.copyOf(g, a.length); } public long[] exp(long[] a) { return isNTTFriendly ? expFFT(a) : expNaive(a); } public long[] mul(long[] a, long b) { long[] ret = new long[a.length]; for (int i = 0; i < a.length; i++) { ret[i] = (a[i] * b) % mod; } return ret; } public long[] pow(long[] a, long m) { int len = a.length; if (m == 0) { long[] ret = new long[len]; ret[0] = 1; return ret; } if (m == 1) { return a.clone(); } if (m == 2) { return squared(a); } int s = 0; while ((s < a.length) && (a[s] == 0)) { ++s; } if (s == a.length) { return a.clone(); } long[] aa = (s != 0) ? Arrays.copyOfRange(a, s, a.length) : a.clone(); long b = MathUtils.modInv(aa[0], mod); for (int i = 0; i < aa.length; i++) { aa[i] = (b * aa[i]) % mod; } aa = log(aa); for (int i = 0; i < aa.length; i++) { aa[i] = ((m % mod) * aa[i]) % mod; } aa = exp(aa); b = MathUtils.modPow(MathUtils.modInv(b, mod), m % (mod - 1), mod); for (int i = 0; i < aa.length; i++) { aa[i] = (b * aa[i]) % mod; } long[] ret = new long[len]; if (s <= ((len - 1) / m)) { for (long i = ((long) (s)) * m; (i < len) && ((i - (s * m)) < aa.length); ++i) { ret[((int) (i))] = aa[((int) (i - (s * m)))]; } } return ret; } public long[] negatedX(long[] f) { long[] g = Arrays.copyOf(f, f.length); for (int i = 1; i < g.length; i += 2) { if (g[i] != 0) { g[i] = mod - g[i]; } } return g; } public long[] divideByX(long[] f, int repeat) { return Arrays.copyOfRange(f, repeat, f.length); } public long[] evenMul(long[] f) { if (f.length == 1) { return new long[]{ (f[0] * f[0]) % mod }; } if ((!isNTTFriendly) || (((f.length + f.length) - 1) <= FFT_NAIVE_THRESHOLD)) { long[] ret = new long[(2 * f.length) - 1]; for (int i = 0; i < f.length; i += 2) { for (int j = i + 2; j < f.length; j += 2) { ret[i + j] = (ret[i + j] + ((f[i] * f[j]) * 2)) % mod; } } for (int i = 0; i < f.length; i += 2) { ret[2 * i] = (ret[2 * i] + (f[i] * f[i])) % mod; } for (int i = 1; i < f.length; i += 2) { for (int j = i + 2; j < f.length; j += 2) { ret[i + j] = ((ret[i + j] + mod) - (((f[i] * f[j]) * 2) % mod)) % mod; } } for (int i = 1; i < f.length; i += 2) { ret[2 * i] = ((ret[2 * i] + mod) - ((f[i] * f[i]) % mod)) % mod; } return ret; } int n = 1; while (n < ((f.length + f.length) - 1)) { n *= 2; } prepareRoots(n); prepareRoots(n / 2); long[] fft = Arrays.copyOf(f, n); for (int i = 0; i < fft.length; i++) { fft[i] = fp.reduce(fft[i]); } fftToBitReversed(fft); long[] fft2 = new long[n]; for (int i = 0; i < fft.length; i++) { fft2[i] = fft[(i ^ 1) % n]; } for (int i = 0; (2 * i) < fft.length; i++) { fft[i] = (fft[2 * i] * fft2[2 * i]) % mod; } fft = Arrays.copyOf(fft, n / 2); ifftFromBitReversed(fft); long[] ret = new long[(f.length + f.length) - 1]; for (int i = 0; (2 * i) < ret.length; i++) { ret[2 * i] = fft[i]; } return ret; } public long nth(long n, long[] numerator, long[] denominator) { if (numerator.length == 0) { return 0; } if (denominator[0] != 1) { throw new AssertionError(); } while (n != 0) { long[] a = Arrays.copyOf(denominator, denominator.length); for (int i = 1; i < a.length; i += 2) { if (a[i] != 0) { a[i] = mod - a[i]; } } numerator = mul(numerator, a); denominator = evenMul(denominator); long[] num2 = new long[(numerator.length + 1) / 2]; long[] den2 = new long[(denominator.length + 1) / 2]; for (int i = ((int) (n % 2)); i < numerator.length; i += 2) { num2[i / 2] = numerator[i]; } for (int i = 0; i < denominator.length; i += 2) { den2[i / 2] = denominator[i]; } numerator = num2; denominator = den2; n /= 2; } return (numerator[0] + mod) % mod; } @Override public long[] exactDiv(long[] a, long[] b) { return div(a, b); } @Override public long[] div(long[] a, long[] b) { if (b.length == 1) { if (b[0] == 0) { throw new ArithmeticException("Division by zero"); } if (b[0] == 1) { return a.clone(); } long inv = fp.inv(b[0]); long[] res = new long[a.length]; for (int i = 0; i < a.length; i++) { res[i] = (a[i] * inv) % mod; } return res; } int degA = deg(a); int degB = deg(b); if ((isNTTFriendly && (((degA - degB) + 1) > FFT_NAIVE_THRESHOLD)) && (degB >= 10)) { return divFast(a, b); } return divNaive(a, b); } public long[] divFast(long[] a, long[] b) { int degA = deg(a); int degB = deg(b); if (degA < degB) { return new long[]{ 0 }; } long[] ra = resize(a); long[] rb = resize(b); ArrayUtils.reverse(ra); ArrayUtils.reverse(rb); rb = Arrays.copyOf(rb, degA + 1); long[] q = mul(ra, inv(rb)); q = Arrays.copyOf(q, (degA - degB) + 1); ArrayUtils.reverse(q); return q; } @Override public long[] mod(long[] a, long[] b) { int degA = deg(a); int degB = deg(b); if (isNTTFriendly && (((degA - degB) + 1) > FFT_NAIVE_THRESHOLD)) { return modFast(a, b); } return modNaive(a, b); } public long[] modFast(long[] a, long[] b) { long[] q = divFast(a, b); return resize(sub(a, mul(b, q))); } public class HalfGcdResult { public long[] p00; public long[] p01; public long[] p10; public long[] p11; public HalfGcdResult(long[] p00, long[] p01, long[] p10, long[] p11) { this.p00 = p00; this.p01 = p01; this.p10 = p10; this.p11 = p11; } public long[][] apply(long[] a, long[] b) { return new long[][]{ resize(add(mul(p00, a), mul(p01, b))), resize(add(mul(p10, a), mul(p11, b))) }; } HalfGcdResult swapColumns() { return new HalfGcdResult(p01, p00, p11, p10); } } HalfGcdResult identityMatrix() { return new HalfGcdResult(new long[]{ 1 }, new long[]{ 0 }, new long[]{ 0 }, new long[]{ 1 }); } HalfGcdResult leftMulEuclideanStep(HalfGcdResult mat, long[] q) { return new HalfGcdResult(mat.p10, mat.p11, sub(mat.p00, mul(q, mat.p10)), sub(mat.p01, mul(q, mat.p11))); } HalfGcdResult multiplyMatrix(HalfGcdResult a, HalfGcdResult b) { return new HalfGcdResult(resize(add(mul(a.p00, b.p00), mul(a.p01, b.p10))), resize(add(mul(a.p00, b.p01), mul(a.p01, b.p11))), resize(add(mul(a.p10, b.p00), mul(a.p11, b.p10))), resize(add(mul(a.p10, b.p01), mul(a.p11, b.p11)))); } HalfGcdResult halfGcdNaiveOrdered(long[] a, long[] b) { int threshold = deg(a) / 2; HalfGcdResult mat = identityMatrix(); while (deg(b) > threshold) { DivModResult dm = divmod(a, b); mat = leftMulEuclideanStep(mat, dm.q); a = b; b = dm.r; } return mat; } public HalfGcdResult halfGcd(long[] a, long[] b) { a = resize(a); b = resize(b); int degA = deg(a); int degB = deg(b); if (degB == (-1)) { return identityMatrix(); } if (degA < degB) { return halfGcd(b, a).swapColumns(); } if (degB <= (degA / 2)) { return identityMatrix(); } if (degA <= 128) { return halfGcdNaiveOrdered(a, b); } int threshold = degA / 2; int shift = (degA + 1) / 2; HalfGcdResult mat = halfGcd(divideByX(a, shift), divideByX(b, shift)); long[][] cd = mat.apply(a, b); long[] c = cd[0]; long[] d = cd[1]; if (deg(d) <= threshold) { return mat; } DivModResult dm = divmod(c, d); mat = leftMulEuclideanStep(mat, dm.q); c = d; d = dm.r; if (deg(d) <= threshold) { return mat; } int secondShift = (2 * threshold) - deg(c); if (secondShift < 0) { throw new AssertionError(); } return multiplyMatrix(halfGcd(divideByX(c, secondShift), divideByX(d, secondShift)), mat); } @Override public long[] gcd(long[] a, long[] b) { int degA = deg(a); int degB = deg(b); if (degA == (-1)) { return monic(b); } if (degB == (-1)) { return monic(a); } if ((degA == 0) || (degB == 0)) { return new long[]{ 1 }; } a = resize(a); b = resize(b); if (degA < degB) { long[] t = a; a = b; b = t; } while (deg(b) != (-1)) { if ((!isNTTFriendly) || (Math.max(deg(a), deg(b)) <= 3072)) { return gcdNaive(a, b); } HalfGcdResult mat = halfGcd(a, b); long[][] cd = mat.apply(a, b); a = cd[0]; b = cd[1]; if (deg(b) == (-1)) { break; } DivModResult dm = divmod(a, b); a = b; b = dm.r; if (deg(a) < deg(b)) { long[] t = a; a = b; b = t; } } return monic(a); } @Override public EuclideanDomainStrategy.ExtGCDResult extgcd(long[] f, long[] g) { f = resize(f); g = resize(g); long[] a = f; long[] b = g; long[] x0 = new long[]{ 1 }; long[] y0 = new long[]{ 0 }; long[] x1 = new long[]{ 0 }; long[] y1 = new long[]{ 1 }; if (Math.max(deg(a), deg(b)) <= 3072) { while (deg(b) != (-1)) { DivModResult dm = divmod(a, b); long[] nx = sub(x0, mul(dm.q, x1)); long[] ny = sub(y0, mul(dm.q, y1)); a = b; b = dm.r; x0 = x1; y0 = y1; x1 = resize(nx); y1 = resize(ny); } } else { if (deg(a) < deg(b)) { long[] t = a; a = b; b = t; t = x0; x0 = x1; x1 = t; t = y0; y0 = y1; y1 = t; } while (deg(b) != (-1)) { HalfGcdResult mat = halfGcd(a, b); long[][] cd = mat.apply(a, b); long[] nx0 = resize(add(mul(mat.p00, x0), mul(mat.p01, x1))); long[] ny0 = resize(add(mul(mat.p00, y0), mul(mat.p01, y1))); long[] nx1 = resize(add(mul(mat.p10, x0), mul(mat.p11, x1))); long[] ny1 = resize(add(mul(mat.p10, y0), mul(mat.p11, y1))); a = cd[0]; b = cd[1]; x0 = nx0; y0 = ny0; x1 = nx1; y1 = ny1; if (deg(b) == (-1)) { break; } DivModResult dm = divmod(a, b); nx1 = resize(sub(x0, mul(dm.q, x1))); ny1 = resize(sub(y0, mul(dm.q, y1))); a = b; b = dm.r; x0 = x1; y0 = y1; x1 = nx1; y1 = ny1; if (deg(a) < deg(b)) { long[] t = a; a = b; b = t; t = x0; x0 = x1; x1 = t; t = y0; y0 = y1; y1 = t; } } } a = resize(a); int d = deg(a); if (d == (-1)) { return new EuclideanDomainStrategy.ExtGCDResult<>(new long[]{ 0 }, new long[]{ 0 }, new long[]{ 0 }); } long inv = MathUtils.modInv(a[d], mod); return new EuclideanDomainStrategy.ExtGCDResult<>(resize(mul(x0, inv)), resize(mul(y0, inv)), resize(mul(a, inv))); } public long[] solveLinearCongruence(long[] a, long[] b, long[] c) { a = resize(a); b = resize(b); c = resize(c); if (isZero(c)) { throw new ArithmeticException("modulo by zero polynomial"); } EuclideanDomainStrategy.ExtGCDResult res = extgcd(a, c); long[] g = res.gcd(); DivModResult dmB = divmod(b, g); if (!isZero(dmB.r)) { return null; } long[] q = dmB.q; long[] modulus = div(c, g); if (deg(modulus) == 0) { return new long[0]; } long[] particular = mul(q, res.x()); return mod(particular, modulus); } public record Term(int d, long v) {} public long[] sparseMul(long[] a, ArrayList sparseTerms, int sparseLen) { if ((a.length == 0) || (sparseLen == 0)) { return new long[0]; } long[] res = new long[(a.length + sparseLen) - 1]; if (sparseTerms.isEmpty()) { return res; } for (int i = 0; i < a.length; i++) { long v = fp.reduce(a[i]); if (v == 0) { continue; } for (Term t : sparseTerms) { res[i + t.d] = (res[i + t.d] + (v * t.v)) % mod; } } return res; } public ArrayList getTerms(long[] p, int initialCapacity) { ArrayList terms = new ArrayList<>(initialCapacity); for (int i = 0; i < p.length; i++) { long v = fp.reduce(p[i]); if (v != 0) { terms.add(new Term(i, v)); } } return terms; } private long[] consecutiveTermsOfInv(long[] F, long l, long r) { if ((F.length == 0) || (fp.reduce(F[0]) != 1)) { throw new IllegalArgumentException("F[0] must be congruent to 1"); } int d = deg(F); if (d == (-1)) { throw new ArithmeticException("Division by zero polynomial"); } if (l >= r) { return new long[0]; } if ((d == 0) || (r <= 0)) { long[] res = new long[((int) (r - l))]; if (r == 0) { res[((int) ((-l) - 1))] = 1; } return res; } long[] F_neg = negatedX(F); long[] prod = mul(F, F_neg); long[] V = new long[(prod.length + 1) / 2]; for (int i = 0; i < V.length; i++) { V[i] = prod[2 * i]; } long lPrime = Math.ceilDiv((l - 1) - d, 2); long rPrime = Math.floorDiv(r, 2); long[] C = consecutiveTermsOfInv(V, lPrime, rPrime); int bSize = ((int) (r - (l - d))); long[] B = new long[bSize]; for (int i = 0; i < C.length; i++) { long degree = 2 * ((lPrime + 1) + i); int idx = ((int) (degree - ((l - d) + 1))); if ((0 <= idx) && (idx < bSize)) { B[idx] = C[i]; } } int size = ((int) (r - l)); long[] G = negatedX(F); long[] conv = mul(G, B); if (conv.length == size) { return conv; } else { return Arrays.copyOfRange(conv, d, d + size); } } public long[] xpowMod(long k, long[] F) { if (k < 0) { throw new IllegalArgumentException("k must be non-negative"); } int d = deg(F); if (d == (-1)) { throw new ArithmeticException("Division by zero polynomial"); } if (d == 0) { return new long[0]; } if (k < d) { long[] res = new long[((int) (k)) + 1]; res[((int) (k))] = 1; return res; } long[] F_rev = F.clone(); ArrayUtils.reverse(F_rev); long fd = fp.reduce(F[d]); long inv_fd = fp.inv(fd); long[] G = new long[d + 1]; for (int i = 0; i <= d; i++) { G[i] = (F_rev[i] * inv_fd) % mod; if (G[i] < 0) { G[i] += mod; } } long[] terms = consecutiveTermsOfInv(G, k - d, k + 1); for (int i = 0; i < terms.length; i++) { terms[i] = (terms[i] * inv_fd) % mod; } long[] R_rev = mul(terms, F_rev); long[] ret = new long[d]; for (int i = 0; i < R_rev.length; i++) { if (((d - 1) - i) >= 0) { ret[i] = R_rev[(d - 1) - i]; } } return resize(ret); } public long[] geometricSumMod(long k, long[] F) { if (k < 0) { throw new IllegalArgumentException("k must be non-negative"); } int d = deg(F); if (d == (-1)) { throw new ArithmeticException("Division by zero polynomial"); } if (d == 0) { return new long[0]; } if (k < d) { long[] res = new long[((int) (k))]; Arrays.fill(res, 1); return res; } long[] P = new long[d + 2]; for (int i = 0; i <= (d + 1); i++) { long f_prev = ((i > 0) && ((i - 1) <= d)) ? F[i - 1] : 0; long f_curr = (i <= d) ? F[i] : 0; P[i] = (f_prev - f_curr) % mod; if (P[i] < 0) { P[i] += mod; } } long[] A = xpowMod(k, P); int len = Math.max(1, A.length); long[] R = new long[len]; System.arraycopy(A, 0, R, 0, A.length); R[0] = ((R[0] - 1) + mod) % mod; long[] T = new long[R.length - 1]; long cur = 0; for (int i = 0; i < T.length; i++) { cur = (cur - R[i]) % mod; if (cur < 0) { cur += mod; } T[i] = cur; } return resize(T); } public long solveCyclicCongruencePoint(long[] A, long N, long[] BA, long[] BleqR_modA, long R) { if (N <= 0) { throw new IllegalArgumentException("N must be positive"); } if ((A.length == 0) || (fp.reduce(A[0]) == 0)) { throw new IllegalArgumentException("A(0) must be non-zero"); } int k = deg(A); if (k <= 0) { throw new IllegalArgumentException("deg(A) must be positive"); } if ((R < 0) || (R >= N)) { return 0L; } long[] r = xpowMod(N, A); long[] r_minus_1 = r.clone(); if (r_minus_1.length == 0) { r_minus_1 = new long[]{ 0 }; } r_minus_1[0] = subMod(r_minus_1[0], 1); long[] Q = solveLinearCongruence(r_minus_1, neg(BA), A); if (Q == null) { throw new ArithmeticException("No solution for Q exists"); } long[] R0 = sub(BleqR_modA, Q); long invA0 = fp.inv(fp.reduce(A[0])); long[] normA = mul(A, invA0); long[] normR0 = mul(R0, invA0); return nth(R, normR0, normA); } } interface RingStrategy extends SemiRingStrategy { T neg(T a); default T sub(T a, T b) { return add(a, neg(b)); } } interface SemiRingStrategy { T zero(); T one(); T add(T a, T b); T mul(T a, T b); boolean equals(T a, T b); default T pow(T a, long n) { if (n < 0) { throw new IllegalArgumentException("Exponent must be non-negative"); } T res = one(); T base = a; while (n > 0) { if ((n & 1) == 1) { res = mul(res, base); } base = mul(base, base); n >>= 1; } return res; } default boolean isZero(T a) { return equals(zero(), a); } default boolean isOne(T a) { return equals(one(), a); } default int hashCode(T a) { return Objects.hashCode(a); } } interface UFDStrategy extends GCDDomainStrategy {} class Zn implements LongCommutativeRingStrategy { final long mod; public Zn(long mod) { this.mod = mod; } @Override public long zero() { return 0; } @Override public long one() { return 1 % mod; } public long pow(long a, long n) { if (n < 0) { throw new AssertionError(); } return MathUtils.modPow(a, n, mod); } public static long crt(long[] a, long[] m) { int N = a.length; long fac = 1; long x = 0; for (int i = 0; i < N; i++) { Fp mo = new Fp(m[i]); long c = mo.reduce((a[i] % m[i]) - (x % m[i])) * MathUtils.modInv(fac % m[i], m[i]); c = mo.reduce(c); x = x + (fac * c); fac *= m[i]; } return x; } @Override public long add(long a, long b) { long ret = (a + b) % mod; if (ret < 0) { ret += mod; } return ret; } public long sub(long a, long b) { long ret = (a - b) % mod; if (ret < 0) { ret += mod; } return ret; } @Override public long mul(long a, long b) { return (a * b) % mod; } @Override public long neg(long a) { return a == 0 ? 0 : mod - a; } @Override public boolean equals(long a, long b) { return a == b; } public long reduce(long a) { a %= mod; if (a < 0) { a += mod; } return a; } }