結果
| 問題 |
No.980 Fibonacci Convolution Hard
|
| ユーザー |
 uwi
|
| 提出日時 | 2020-01-31 22:30:19 |
| 言語 | Java (openjdk 23) |
| 結果 |
AC
|
| 実行時間 | 1,411 ms / 2,000 ms |
| コード長 | 14,948 bytes |
| コンパイル時間 | 4,505 ms |
| コンパイル使用メモリ | 90,196 KB |
| 実行使用メモリ | 63,564 KB |
| 最終ジャッジ日時 | 2024-12-24 16:13:22 |
| 合計ジャッジ時間 | 30,693 ms |
|
ジャッジサーバーID （参考情報） |
judge4 / judge2 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 17 |
ソースコード
package contest200131;
import java.io.ByteArrayInputStream;
import java.io.IOException;
import java.io.InputStream;
import java.io.PrintWriter;
import java.util.Arrays;
import java.util.InputMismatchException;
public class F {
InputStream is;
PrintWriter out;
String INPUT = "";
// 1 6 29 126
// 1 8 50 280
// 1 10 77 530
// 0 1 p p^2+1 p^3+2p p^4+3p^2+1 p^5+4p^3+3p
// 2p
// 3p^2+2
// 4p^3+6p
// 2p^4+6p^2+2+2p^4+4p^2+p^4+2p^2+1
// 5p^4+12p^2+3
// C(n-m,m)*(n+1-m)*p^2~
// (n-m)!/(n-2m)!/m!*(n+1-m)*p^4+..
void solve()
{
int P = ni();
long[] a = new long[200];
// long[] a = new long[10];
a[1] = 1;
int mod = 1000000007;
for(int i = 2;i < a.length;i++){
a[i] = (a[i-1] * P + a[i-2]) % mod;
}
long[] b = convolute(a, a, 3, mod);
long[][] f = Arrays.copyOf(guessLeaned(mod, Arrays.copyOf(b, 20)), 3);
long[] ret = Arrays.copyOf(b, 2000001);
for(int i = f.length;i <= 2000000;i++){
long[] q = Arrays.copyOfRange(ret, i-(f.length-1), i);
ret[i] = f(f, q, i, mod);
}
for(int Q = ni();Q > 0;Q--){
out.println(ret[ni()-2]);
}
}
public static long[][] guessLeaned(int mod, long... a)
{
int n = a.length;
// #formula >= #variable
// n-r+2 >= r(r+1)/2
for(int r = n;r >= 1;r--){
if(n-r+2 < r*(r+1)/2)continue;
int[][] M = new int[n-r+2][r*(r+1)/2];
for(int i = 0;i < n-r+1;i++){
int p = 0;
for(int j = 0;j < r;j++){
long prod = 1;
for(int k = 0;k <= r-j-1;k++){
M[i][p++] = (int)(prod*a[i+j]%mod);
prod = prod * i % mod;
}
}
}
M[n-r+1][0] = 1;
int[] v = new int[n-r+2];
v[n-r+1] = 1;
Result res = gaussElimination(M, v, mod);
if(res.exists){
long[][] ret = new long[r][];
int p = 0;
for(int i = 0;i < r;i++){
ret[i] = new long[r-i];
for(int j = 0;j < r-i;j++){
ret[i][j] = res.sol[p++];
}
}
return ret;
}
}
return null;
}
public static Result gaussElimination(int[][] M, int[] v, int mod)
{
int n = M.length, m = M[0].length;
int[] head = new int[n];
// if not needed, comment out.
for(int[] row : M){
for(int i = 0;i < row.length;i++){
row[i] %= mod;
if(row[i] < 0)row[i] += mod;
}
}
// Forward Elimination
int row = 0;
for(int col = 0;col < m;col++){
// select pivot
boolean pivotFound = false;
out:
for(int prow = row;prow < n;prow++){
if(M[prow][col] != 0){
// pivot found
if(prow != row){
// swap rows
for(int k = 0;k < m;k++){
int u = M[prow][k]; M[prow][k] = M[row][k]; M[row][k] = u;
}
int dum = v[prow]; v[prow] = v[row]; v[row] = dum;
}
pivotFound = true;
break out;
}
}
if(!pivotFound)continue;
head[row] = col;
// diag to 1
long imul = invl(M[row][col], mod);
for(int k = 0;k < m;k++){
M[row][k] = (int)(M[row][k] * imul % mod);
}
v[row] = (int)(v[row] * imul % mod);
for(int j = row+1;j < n;j++){
if(M[j][col] != 0){
long mul = mod-M[j][col];
for(int k = col;k < m;k++){
M[j][k] = (int)((M[j][k] + M[row][k] * mul) % mod);
}
v[j] = (int)((v[j] + v[row] * mul) % mod);
}
}
row++;
}
Result ret = new Result();
ret.mat = M;
for(int i = row;i < n;i++){
if(v[i] != 0){
ret.rank = row;
ret.exists = false;
return ret;
}
}
for(int i = row-1;i >= 0;i--){
for(int j = i-1;j >= 0;j--){
if(M[j][head[i]] != 0){
long mul = mod-M[j][head[i]];
for(int k = head[i];k < m;k++){
M[j][k] = (int)((M[j][k] + M[i][k] * mul) % mod);
}
v[j] = (int)((v[j] + v[i] * mul) % mod);
}
}
}
int[] retv = new int[m];
for(int i = 0;i < row;i++){
retv[head[i]] = v[i];
}
ret.sol = retv;
ret.rank = row;
ret.exists = true;
return ret;
}
public static class Result
{
public int[][] mat;
public int[] sol;
public int rank;
public boolean exists;
}
/**
*
* @param ged guessしたやつ
* @param prevs 求めたい項の直前len(ged)-1項。末尾が求めたい項に近いほう
* @param x 求めたい項の番号
* @param mod
* @return
*/
public static long f(long[][] ged, long[] prevs, long x, int mod)
{
int n = ged.length;
assert prevs.length == n-1;
x -= n-1;
long s = 0;
long tar = 0;
for(int i = 0;i < n;i++){
long co = 0;
for(int j = ged[i].length-1;j >= 0;j--){
co = (co * x + ged[i][j]) % mod;
}
if(i < n-1){
s += co * prevs[i];
s %= mod;
}else{
tar = co;
}
}
long ret = -invl(tar, mod) * s % mod;
if(ret < 0)ret += mod;
return ret;
}
public static long f(long x, long mod, long[] f)
{
long ret = 0;
for(int i = f.length-1;i >= 0;i--)ret = (ret * x + f[i]) % mod;
return ret;
}
public static long[] guess(long mod, long... y)
{
int n = y.length;
long[] dp = new long[n+1];
dp[0] = 1;
// (x-x0)(x-x1)...(x-x{n-1})
for(int i = 0;i < n;i++){
for(int j = i;j >= 0;j--){
dp[j+1] += dp[j];
if(dp[j+1] >= mod)dp[j+1] -= mod;
dp[j] = dp[j]*-i%mod;
if(dp[j] < 0)dp[j] += mod;
}
}
long[] f = new long[n+1];
f[0] = 1;
for(int i = 1;i <= n;i++)f[i] = f[i-1] * i % mod;
long[] ret = new long[n];
for(int i = 0;i < n;i++){
long den = f[i]*f[n-1-i]%mod; // (-1)^(n-1-i)*i!*(n-1-i)!
if(((i^n-1)&1) == 1){
den = mod - den;
}
long iden = invl(den, mod) * y[i] % mod;
long minus = 0;
for(int j = n-1;j >= 0;j--){
minus = (dp[j+1] + minus * i) % mod;
ret[j] = (ret[j] + minus*iden)%mod;
}
}
return ret;
}
// public static final int[] NTTPrimes = {1053818881, 1051721729, 1045430273, 1012924417, 1007681537, 1004535809, 998244353, 985661441, 976224257, 975175681};
// public static final int[] NTTPrimitiveRoots = {7, 6, 3, 5, 3, 3, 3, 3, 3, 17};
public static final int[] NTTPrimes = {1012924417, 1004535809, 998244353, 985661441, 975175681, 962592769, 950009857, 943718401, 935329793, 924844033};
public static final int[] NTTPrimitiveRoots = {5, 3, 3, 3, 17, 7, 7, 7, 3, 5};
public static long[] convoluteSimply(long[] a, long[] b, int P, int g)
{
int m = Math.max(2, Integer.highestOneBit(Math.max(a.length, b.length)-1)<<2);
long[] fa = nttmb(a, m, false, P, g);
long[] fb = a == b ? fa : nttmb(b, m, false, P, g);
for(int i = 0;i < m;i++){
fa[i] = fa[i]*fb[i]%P;
}
return nttmb(fa, m, true, P, g);
}
public static long[] convolute(long[] a, long[] b)
{
int USE = 2;
int m = Math.max(2, Integer.highestOneBit(Math.max(a.length, b.length)-1)<<2);
long[][] fs = new long[USE][];
for(int k = 0;k < USE;k++){
int P = NTTPrimes[k], g = NTTPrimitiveRoots[k];
long[] fa = nttmb(a, m, false, P, g);
long[] fb = a == b ? fa : nttmb(b, m, false, P, g);
for(int i = 0;i < m;i++){
fa[i] = fa[i]*fb[i]%P;
}
fs[k] = nttmb(fa, m, true, P, g);
}
int[] mods = Arrays.copyOf(NTTPrimes, USE);
long[] gammas = garnerPrepare(mods);
int[] buf = new int[USE];
for(int i = 0;i < fs[0].length;i++){
for(int j = 0;j < USE;j++)buf[j] = (int)fs[j][i];
long[] res = garnerBatch(buf, mods, gammas);
long ret = 0;
for(int j = res.length-1;j >= 0;j--)ret = ret * mods[j] + res[j];
fs[0][i] = ret;
}
return fs[0];
}
public static long[] convolute(long[] a, long[] b, int USE, int mod)
{
int m = Math.max(2, Integer.highestOneBit(Math.max(a.length, b.length)-1)<<2);
long[][] fs = new long[USE][];
for(int k = 0;k < USE;k++){
int P = NTTPrimes[k], g = NTTPrimitiveRoots[k];
long[] fa = nttmb(a, m, false, P, g);
long[] fb = a == b ? fa : nttmb(b, m, false, P, g);
for(int i = 0;i < m;i++){
fa[i] = fa[i]*fb[i]%P;
}
fs[k] = nttmb(fa, m, true, P, g);
}
int[] mods = Arrays.copyOf(NTTPrimes, USE);
long[] gammas = garnerPrepare(mods);
int[] buf = new int[USE];
for(int i = 0;i < fs[0].length;i++){
for(int j = 0;j < USE;j++)buf[j] = (int)fs[j][i];
long[] res = garnerBatch(buf, mods, gammas);
long ret = 0;
for(int j = res.length-1;j >= 0;j--)ret = (ret * mods[j] + res[j]) % mod;
fs[0][i] = ret;
}
return fs[0];
}
// static int[] wws = new int[270000]; // outer faster
// Modifed Montgomery + Barrett
private static long[] nttmb(long[] src, int n, boolean inverse, int P, int g)
{
long[] dst = Arrays.copyOf(src, n);
int h = Integer.numberOfTrailingZeros(n);
long K = Integer.highestOneBit(P)<<1;
int H = Long.numberOfTrailingZeros(K)*2;
long M = K*K/P;
int[] wws = new int[1<<h-1];
long dw = inverse ? pow(g, P-1-(P-1)/n, P) : pow(g, (P-1)/n, P);
long w = (1L<<32)%P;
for(int k = 0;k < 1<<h-1;k++){
wws[k] = (int)w;
w = modh(w*dw, M, H, P);
}
long J = invl(P, 1L<<32);
for(int i = 0;i < h;i++){
for(int j = 0;j < 1<<i;j++){
for(int k = 0, s = j<<h-i, t = s|1<<h-i-1;k < 1<<h-i-1;k++,s++,t++){
long u = (dst[s] - dst[t] + 2*P)*wws[k];
dst[s] += dst[t];
if(dst[s] >= 2*P)dst[s] -= 2*P;
// long Q = (u&(1L<<32)-1)*J&(1L<<32)-1;
long Q = (u<<32)*J>>>32;
dst[t] = (u>>>32)-(Q*P>>>32)+P;
}
}
if(i < h-1){
for(int k = 0;k < 1<<h-i-2;k++)wws[k] = wws[k*2];
}
}
for(int i = 0;i < n;i++){
if(dst[i] >= P)dst[i] -= P;
}
for(int i = 0;i < n;i++){
int rev = Integer.reverse(i)>>>-h;
if(i < rev){
long d = dst[i]; dst[i] = dst[rev]; dst[rev] = d;
}
}
if(inverse){
long in = invl(n, P);
for(int i = 0;i < n;i++)dst[i] = modh(dst[i]*in, M, H, P);
}
return dst;
}
// Modified Shoup + Barrett
private static long[] nttsb(long[] src, int n, boolean inverse, int P, int g)
{
long[] dst = Arrays.copyOf(src, n);
int h = Integer.numberOfTrailingZeros(n);
long K = Integer.highestOneBit(P)<<1;
int H = Long.numberOfTrailingZeros(K)*2;
long M = K*K/P;
long dw = inverse ? pow(g, P-1-(P-1)/n, P) : pow(g, (P-1)/n, P);
long[] wws = new long[1<<h-1];
long[] ws = new long[1<<h-1];
long w = 1;
for(int k = 0;k < 1<<h-1;k++){
wws[k] = (w<<32)/P;
ws[k] = w;
w = modh(w*dw, M, H, P);
}
for(int i = 0;i < h;i++){
for(int j = 0;j < 1<<i;j++){
for(int k = 0, s = j<<h-i, t = s|1<<h-i-1;k < 1<<h-i-1;k++,s++,t++){
long ndsts = dst[s] + dst[t];
if(ndsts >= 2*P)ndsts -= 2*P;
long T = dst[s] - dst[t] + 2*P;
long Q = wws[k]*T>>>32;
dst[s] = ndsts;
dst[t] = ws[k]*T-Q*P&(1L<<32)-1;
}
}
// dw = dw * dw % P;
if(i < h-1){
for(int k = 0;k < 1<<h-i-2;k++){
wws[k] = wws[k*2];
ws[k] = ws[k*2];
}
}
}
for(int i = 0;i < n;i++){
if(dst[i] >= P)dst[i] -= P;
}
for(int i = 0;i < n;i++){
int rev = Integer.reverse(i)>>>-h;
if(i < rev){
long d = dst[i]; dst[i] = dst[rev]; dst[rev] = d;
}
}
if(inverse){
long in = invl(n, P);
for(int i = 0;i < n;i++){
dst[i] = modh(dst[i] * in, M, H, P);
}
}
return dst;
}
static final long mask = (1L<<31)-1;
public static long modh(long a, long M, int h, int mod)
{
long r = a-((M*(a&mask)>>>31)+M*(a>>>31)>>>h-31)*mod;
return r < mod ? r : r-mod;
}
private static long[] garnerPrepare(int[] m)
{
int n = m.length;
assert n == m.length;
if(n == 0)return new long[0];
long[] gamma = new long[n];
for(int k = 1;k < n;k++){
long prod = 1;
for(int i = 0;i < k;i++){
prod = prod * m[i] % m[k];
}
gamma[k] = invl(prod, m[k]);
}
return gamma;
}
private static long[] garnerBatch(int[] u, int[] m, long[] gamma)
{
int n = u.length;
assert n == m.length;
long[] v = new long[n];
v[0] = u[0];
for(int k = 1;k < n;k++){
long temp = v[k-1];
for(int j = k-2;j >= 0;j--){
temp = (temp * m[j] + v[j]) % m[k];
}
v[k] = (u[k] - temp) * gamma[k] % m[k];
if(v[k] < 0)v[k] += m[k];
}
return v;
}
private static long pow(long a, long n, long mod) {
// a %= mod;
long ret = 1;
int x = 63 - Long.numberOfLeadingZeros(n);
for (; x >= 0; x--) {
ret = ret * ret % mod;
if (n << 63 - x < 0)
ret = ret * a % mod;
}
return ret;
}
private static long invl(long a, long mod) {
long b = mod;
long p = 1, q = 0;
while (b > 0) {
long c = a / b;
long d;
d = a;
a = b;
b = d % b;
d = p;
p = q;
q = d - c * q;
}
return p < 0 ? p + mod : p;
}
void run() throws Exception
{
is = INPUT.isEmpty() ? System.in : new ByteArrayInputStream(INPUT.getBytes());
out = new PrintWriter(System.out);
long s = System.currentTimeMillis();
solve();
out.flush();
if(!INPUT.isEmpty())tr(System.currentTimeMillis()-s+"ms");
// Thread t = new Thread(null, null, "~", Runtime.getRuntime().maxMemory()){
// @Override
// public void run() {
// long s = System.currentTimeMillis();
// solve();
// out.flush();
// if(!INPUT.isEmpty())tr(System.currentTimeMillis()-s+"ms");
// }
// };
// t.start();
// t.join();
}
public static void main(String[] args) throws Exception { new F().run(); }
private byte[] inbuf = new byte[1024];
public int lenbuf = 0, ptrbuf = 0;
private int readByte()
{
if(lenbuf == -1)throw new InputMismatchException();
if(ptrbuf >= lenbuf){
ptrbuf = 0;
try { lenbuf = is.read(inbuf); } catch (IOException e) { throw new InputMismatchException(); }
if(lenbuf <= 0)return -1;
}
return inbuf[ptrbuf++];
}
private boolean isSpaceChar(int c) { return !(c >= 33 && c <= 126); }
private int skip() { int b; while((b = readByte()) != -1 && isSpaceChar(b)); return b; }
private double nd() { return Double.parseDouble(ns()); }
private char nc() { return (char)skip(); }
private String ns()
{
int b = skip();
StringBuilder sb = new StringBuilder();
while(!(isSpaceChar(b))){ // when nextLine, (isSpaceChar(b) && b != ' ')
sb.appendCodePoint(b);
b = readByte();
}
return sb.toString();
}
private char[] ns(int n)
{
char[] buf = new char[n];
int b = skip(), p = 0;
while(p < n && !(isSpaceChar(b))){
buf[p++] = (char)b;
b = readByte();
}
return n == p ? buf : Arrays.copyOf(buf, p);
}
private int[] na(int n)
{
int[] a = new int[n];
for(int i = 0;i < n;i++)a[i] = ni();
return a;
}
private long[] nal(int n)
{
long[] a = new long[n];
for(int i = 0;i < n;i++)a[i] = nl();
return a;
}
private char[][] nm(int n, int m) {
char[][] map = new char[n][];
for(int i = 0;i < n;i++)map[i] = ns(m);
return map;
}
private int[][] nmi(int n, int m) {
int[][] map = new int[n][];
for(int i = 0;i < n;i++)map[i] = na(m);
return map;
}
private int ni() { return (int)nl(); }
private long nl()
{
long num = 0;
int b;
boolean minus = false;
while((b = readByte()) != -1 && !((b >= '0' && b <= '9') || b == '-'));
if(b == '-'){
minus = true;
b = readByte();
}
while(true){
if(b >= '0' && b <= '9'){
num = num * 10 + (b - '0');
}else{
return minus ? -num : num;
}
b = readByte();
}
}
private static void tr(Object... o) { System.out.println(Arrays.deepToString(o)); }
}