Problem:
 Max Score: $1400$ Points   Problem Statement E869120 defined a sequence $a$ like this:   $a_1=a_2=1$, $a_{k+2}=a_{k+1}+a_k \ (k \ge 1)$   He also defined sequences $d_1, d_2, d_3, \dots , d_n$, as the following recurrence relation :   $d_{1, j} = a_j$  $d_{i, j} = \sum_{k = 1}^j d_{i - 1, k} \ (i \ge 2)$   You are given integers $n$ and $m$. Please calculate the value of $d_{n, m}$. Since the answer can be large number, print the answer modulo $998,244,353$. Can you solve this problem???        Input The input is given from standard input in the following format.   $n \quad m$      Output   Print $d_{n, m}$ modulo $998,244,353$.     Constraints   $1 \le n \le 200,000$  $1 \le m \le 10^{18}$     Subtasks Subtask 1 [ $100$ points ]   The testcase in this subtask satisfies $1 \le n, m \le 3,000$.  Subtask 2 [ $170$ points ]   The testcase in this subtask satisfies $1 \le m \le 200,000$.  Subtask 3 [ $230$ points ]   The testcase in this subtask satisfies $1 \le n \le 3$.  Subtask 4 [ $420$ points ]   The testcase in this subtask satisfies $1 \le n \le 1000$.  Subtask 5 [ $480$ points ]   There are no additional constraints.         Sample Input 1  4 7      Sample Output 1  176 The sequence is following:     $d_{i, 1}$  $d_{i, 2}$  $d_{i, 3}$  $d_{i, 4}$  $d_{i, 5}$  $d_{i, 6}$  $d_{i, 7}$    $d_1$  1  1  2  3  5  8  13    $d_2$  1  2  4  7  12  20  33    $d_3$  1  3  7  14  26  46  79    $d_4$  1  4  11  25  51  97  176    As a result, the answer is $d_{4, 7} = 176$.       Sample Input 2  12 20      Sample Output 2  174174144  Interesting number. It only contains $1$, $4$ and $7$.      Sample Input 3  16 30      Sample Output 3  102292850 You can calculate that $d_{16, 30} = 1,193,004,294,685$. Please don't forget printing the answer mod $998,244,353$.    Writer：square1001

Solution (Python):
import sys
input = sys.stdin.readline

import numpy as np

MOD = 998244353  
N,M = map(int,input().split())

def cumprod_mod(arr):
    L = len(arr); Lsq = int((L**.5)+1)
    arr = np.resize(arr,Lsq**2).reshape(Lsq,Lsq)
    for n in range(1,Lsq):
        arr[:,n] *= arr[:,n-1]; arr[:,n] %= MOD
    for n in range(1,Lsq):
        arr[n] *= arr[n-1,-1]; arr[n] %= MOD
    return arr.ravel()[:L]


Nsq = int((N**.5)+1)
a = np.zeros((Nsq,Nsq+1),dtype=np.int64)
b = np.zeros((Nsq,Nsq+1),dtype=np.int64)
a[0,0] = 1
for n in range(1,Nsq+1):
    a[0,n] = -b[0,n-1] % MOD; b[0,n] = (a[0,n-1]+3*b[0,n-1]) % MOD
u,v = a[0,Nsq],b[0,Nsq]
for n in range(1,Nsq):
    a[n] = (a[n-1]*u-b[n-1]*v) % MOD
    b[n] = (a[n-1]*v+b[n-1]*u+3*b[n-1]*v) % MOD
a = a[:,:-1].ravel()
b = b[:,:-1].ravel()
G = (-a-2*b) % MOD
G = G[:N+1]

A,B = -G[N-1],-G[N]+3*G[N-1]
A,B = A+B,-B
A %= MOD; B %= MOD

def fibonacci(n):
    
    if n == 0:
        return 1,0
    a,b = fibonacci(n//2)
    a,b = a*a+b*b,2*a*b+b*b
    a,b = a%MOD,b%MOD
    return (b,a+b) if n&1 else (a,b)

x = (np.arange(N,dtype=np.int64)+M)%MOD
x[0] = 1
num = cumprod_mod(x)

x = np.arange(N,dtype=np.int64)
x[0] = 1
fact = cumprod_mod(x)

x = np.arange(N,0,-1,dtype=np.int64)
x[0] = pow(int(fact[-1]),MOD-2,MOD)
fact_inv = cumprod_mod(x)[::-1]

if N == 1:
    x = 0
else:
    x = ((num * fact_inv % MOD)[:N-1] * G[N-2::-1] % MOD).sum() % MOD
y,z = fibonacci(M+1)
answer = (z*A + y*B + x)%MOD
print(answer)

Solution (Python):
def f(n):
	e=[1,0,0,1];z=[1,1,1,0]
	while n:e=[e,g(e,z)][n%2];z=g(z,z);n/=2
	return e[1]
M=998244353;g=lambda a,b:[(a[x]*b[y]+a[z]*b[w])%M for x,y,z,w in[[0,0,1,2],[0,1,1,3],[2,0,3,2],[2,1,3,3]]];n,m=map(int,raw_input().split());I=[1,1];r=f(m+2*n-2);c=1
for i in range(2,n):I+=[(M-M/i)*I[M%i]%M]
for i in range(n-1):r=(r+c*(M-f(2*n-2-2*i)))%M;c=c*(m+i)*I[i+1]%M
print r

Solution (Python):
def f(n):
	e=[1,0,0,1];z=[1,1,1,0]
	while n:e=[e,g(e,z)][n%2];z=g(z,z);n//=2
	return e[1]
M=998244353;g=lambda a,b:[(a[0]*b[0]+a[1]*b[2])%M,(a[0]*b[1]+a[1]*b[3])%M,(a[2]*b[0]+a[3]*b[2])%M,(a[2]*b[1]+a[3]*b[3])%M];n,m=map(int,input().split());I=[1,1];r=f(m+2*n-2);c=1
for i in range(2,n):I+=[(M-M//i)*I[M%i]%M]
for i in range(n-1):r=(r+c*(M-f(2*n-2-2*i)))%M;c=c*(m+i)*I[i+1]%M
print(r)

Solution (Python):
def f(n):
	i=e=[1,0,0,1];z=[1,1,1,0]
	while n:e=[e,g(e,z)][n%2];z=g(z,z);n//=2
	return e[1]
M=998244353;g=lambda a,b:[(a[0]*b[0]+a[1]*b[2])%M,(a[0]*b[1]+a[1]*b[3])%M,(a[2]*b[0]+a[3]*b[2])%M,(a[2]*b[1]+a[3]*b[3])%M];n,m=map(int,input().split());I=[1,1];r=f(m+2*n-2);c=1
for i in range(2,n):I+=[(M-M//i)*I[M%i]%M]
for i in range(n-1):r=(r+c*(M-f(2*n-2-2*i)))%M;c=c*(m+i)*I[i+1]%M
print(r)

Solution (Python):
def f(n):
	e=[1,0,0,1];z=[1,1,1,0]
	while n:
		if n%2>0:e=g(e,z)
		z=g(z,z);n//=2
	return e[1]
M=998244353;g=lambda a,b:[(a[0]*b[0]+a[1]*b[2])%M,(a[0]*b[1]+a[1]*b[3])%M,(a[2]*b[0]+a[3]*b[2])%M,(a[2]*b[1]+a[3]*b[3])%M];n,m=map(int,input().split());I=[1,1];r=f(m+2*n-2);c=1
for i in range(2,n):I+=[(M-M//i)*I[M%i]%M]
for i in range(n-1):r=(r+c*(M-f(2*n-2-2*i)))%M;c=c*(m+i)*I[i+1]%M
print(r)

Solution (Python):
M=998244353
g=lambda a,b:[(a[0]*b[0]+a[1]*b[2])%M,(a[0]*b[1]+a[1]*b[3])%M,(a[2]*b[0]+a[3]*b[2])%M,(a[2]*b[1]+a[3]*b[3])%M]
def f(n):
	e=[1,0,0,1];z=[1,1,1,0]
	while n:
		if n%2>0:e=g(e,z)
		z=g(z,z);n//=2
	return e[1]
n,m=map(int,input().split());I=[1,1];r=f(m+2*n-2);c=1
for i in range(2,n):I+=[(M-M//i)*I[M%i]%M]
for i in range(n-1):r=(r+c*(M-f(2*n-2-2*i)))%M;c=c*(m+i)*I[i+1]%M
print(r)