　特性方程式は
$$t^2-xt-y=0$$
である．$x,y$ を固定したとき，これの解を $\alpha,\beta$ とすると，以下が成り立つ：
$$a_n=\begin{cases}
\dfrac{\beta^n-\alpha^n}{\beta-\alpha} & (\alpha\neq \beta) \\\\
n\alpha^{n-1} & (\alpha=\beta)
\end{cases}$$
　以下，特性方程式を $\mathbb{F}_p$ 係数として考え，以下の三つのパターンに場合分けして考える．ただし，$\mathbb{F}_p$ は位数 $p$ の有限体である．

+ (1) $\alpha=\beta$であるような場合 $(p-1)$通り
+ (2) $\alpha\neq \beta$かつ解が$\mathbb{F}_p$上に存在する場合 $\dfrac{(p-1)(p-2)}{2}$通り
+ (3) 解 $\alpha,\beta$ が  $\mathbb{F}\_p$ 上に存在しない場合．このとき解は$\mathbb{F}_{p^2}\setminus \mathbb{F}_p$の元である．$\dfrac{p(p-1)}{2}$通り


(1)について

このとき,一般項は
$$a_n=n\alpha^{n-1}$$
である.そして, $k$を周期としたとき
$$n\alpha^{n-1}\equiv (n+k)\alpha^{n+k-1}$$
十分大きな整数$m$について,$n=mp$としたときに等号が成立する.

$$\begin{aligned}
mp\alpha^{mp-1}&\equiv (mp+k)\alpha^{mp+k-1}\\\\
0&\equiv k\alpha^{mp+k-1}\\\\
\end{aligned}$$

よって,$k$が$p$の倍数であることが必要である.さらに,$n=mp+1$の場合について考察すると
$$\begin{aligned}
(mp+1)\alpha^{mp}&\equiv (mp+1+k)\alpha^{mp+k}\\\\
\alpha^{mp}&\equiv \alpha^{mp+k}\\\\
\end{aligned}$$
よって,$\alpha^k=1$であることが必要である.

よって,周期となるためには,$p$と$\alpha$の位数で両方割り切れることが必要である.逆にこの場合$a_{n+k}\equiv a_n \forall n\in\mathbb{N}$なので周期としての条件を満たしている.よって,$f(x,y)$は$p$と${\rm ord}(\alpha)$の最小公倍数である.

ここで,$\alpha$の位数は$p-1$の約数となる.これは$p$とは互いに素なので$f(x,y)=p\cdot {\rm ord}(\alpha)$である.さらに,$\mathbb{F}\_p$の乗法群は$\mathbb{Z}\/(p-1)\mathbb{Z}$に同型であることに気をつけると,
位数が$d$であるような$\alpha\in \mathbb{F}\_p^{\times}$の個数はトーシェント関数を使って$\varphi(d)$と書ける.
よって,(1)の範囲での総和は
$$\sum_{d|(p-1)}dp\varphi(d)$$
である.



(2)について

このとき,一般項は

$$a_n=\dfrac{\beta^n-\alpha^n}{\beta-\alpha}$$

である.

さらに,
$$(\alpha,\beta);1\leq \alpha\lt \beta\leq p-1$$
のペアがそれぞれ1つずつ登場するので,これらについての総和を求める.

周期が$k$だとした場合,ある整数$m$が存在して,
$$a_{m(p-1)}\equiv a_{m(p-1)+k}\pmod{p}$$
であるため,
$$\beta^{m(p-1)}-\alpha^{m(p-1)}=\beta^{m(p-1)+k}-\alpha^{m(p-1)+k}$$
である.ここで,フェルマーの小定理より,
$$0=\beta^{k}-\alpha^{k}$$

さらに,
$$a_{m(p-1)+1}\equiv a_{m(p-1)+k+1}\pmod{p}$$
であるため,
$$\beta^{m(p-1)+1}-\alpha^{m(p-1)+1}=\beta^{m(p-1)+k+1}-\alpha^{m(p-1)+k+1}$$
であるが,
$$\beta-\alpha=\beta^{k+1}-\alpha^{k+1}$$
が必要だが,
$$\beta^{k+1}-\alpha^{k+1}=\beta^{k+1}-\alpha \cdot \beta^{k}=\beta^k (\beta-\alpha)$$
である.よって,$k$が周期であるためには
$$\alpha^k\equiv \beta^k\equiv 1\pmod{p}$$
が必要である.逆にこのときに周期となる.
よって,
$$f(x,y)={\rm lcm}({\rm ord}(\alpha),{\rm ord}(\beta))$$

である.
$1\leq \alpha,\beta\leq p-1$の範囲内での${\rm lcm}({\rm ord}(\alpha),{\rm ord}(\beta))$の総和は
$$\sum_{d_1|(p-1)}\sum_{d_2|(p-1)}\varphi(d_1)\varphi(d_2){\rm lcm}(d_1,d_2)$$
である.さらに,$1\leq \alpha=\beta\leq p-1$の範囲内での${\rm lcm}({\rm ord}(\alpha),{\rm ord}(\beta))$の総和は
$$\sum_{d|(p-1)}d\varphi(d)$$
である.

よって,(2)の場合の答えの総和は

$$\dfrac{1}{2}\sum_{d_1|(p-1)}\sum_{d_2|(p-1)}\varphi(d_1)\varphi(d_2){\rm lcm}(d_1,d_2)-\dfrac{1}{2}\sum_{d|(p-1)}d\varphi(d)$$
である.

(3)について

このときも(2)と同様に一般項は

$$a_n=\dfrac{\beta^n-\alpha^n}{\beta-\alpha}$$

である.ただし,$\alpha,\beta$は$\mathbb{F}_{p^2}\setminus \mathbb{F}_p$の範囲内を動く.ただし,$\alpha$と$\beta$を入れ替えているものは同一視しているため,該当の$(x,y)$の組み合わせは$\dfrac{p(p-1)}{2}$通りである.

(2)と同様にすると,
$$f(x,y)={\rm lcm}({\rm ord}(\alpha),{\rm ord}(\beta))$$
だが,このときのordとは,$\mathbb{F}_{p^2}^{\times}$における位数である.

さらにこのとき,$\alpha,\beta$は$\mathbb{F}\_{p^2}$の元の中で共役元の関係になっている.よって位数は等しい.
つまり${\rm ord}(\alpha)={\rm ord}(\beta)$ということである.なので,結局周期の総和は
$$\dfrac{1}{2}\sum_{\alpha\in\mathbb{F}\_{p^2}\setminus\mathbb{F}\_p}{\rm ord}(\alpha)$$
ここで,実は$\mathbb{F}_{p^2}^{\times}$は群$\mathbb{Z}\/(p^2-1)\mathbb{Z}$に同型であることに注意.
さらにこのとき,$\mathbb{Z}\/(p^2-1)\mathbb{Z}$の元の中で$(p+1)$の倍数であるようなもの$(p-1)$個が$\mathbb{F}_p$の元に対応している.
よって,(3)についての総和は

$$\dfrac{1}{2}\sum_{d|(p^2-1)}d\varphi(d)-\dfrac{1}{2}\sum_{d|(p-1)}d\varphi(d)$$

である.

結局合計すると,最終的な答えは

$$(p-1)\sum_{d|(p-1)}d\varphi(d) + \dfrac{1}{2}\sum_{d|(p^2-1)}d\varphi(d)+\dfrac{1}{2}\sum_{d_1|(p-1)}\sum_{d_2|(p-1)}\varphi(d_1)\varphi(d_2){\rm lcm}(d_1,d_2)$$

となる.この値を$p-1$や$p^2-1$の約数全てについて値を求めるのは人力では大変である.

便宜上以下のように関数$F(n),G(n)$を定める.

$$F(n):=\sum_{d|n} d\varphi(d); G(n):=\sum_{d_1|n}\sum_{d_2|n}\varphi(d_1)\varphi(d_2){\rm lcm}(d_1,d_2)$$

このとき,求めるべき答えは

$$(p-1)F(p-1) + \dfrac{1}{2}F(p^2-1)+\dfrac{1}{2}G(p-1)$$

つまり,$p=2017$であるため,

$$2016\cdot F(2016) + \dfrac{1}{2}F(2016\cdot 2018)+\dfrac{1}{2}G(2016)$$

実は$F(n),G(n)$は乗法的関数である.つまり,$n,m$が互いに素な正の整数とした場合,$F(nm)=F(n)F(m),G(nm)=G(n)G(m)$となるという性質がある.

この性質を使うことで,求めるべき答えは

$$2016\cdot F(2^5)F(3^2)F(7) + \dfrac{1}{2}F(2^6)F(3^2)F(7)F(1009)+\dfrac{1}{2}G(2^5)G(3^2)G(7)$$

であるとわかる.それぞれについて答えを求める.

$$\begin{aligned}
F(7)&=&1\cdot \varphi(1)+7\cdot \varphi(7)&=&1\cdot 1+7\cdot 6&=&43\\\\
F(1009)&=&1\cdot \varphi(1)+1009\cdot \varphi(1009)&=&1\cdot 1+1009\cdot 1008&=&1017073\\\\
F(3^2)&=&1\cdot \varphi(1)+3\cdot \varphi(3)+9\cdot \varphi(9)&=&1+6+54&=&61\\\\
F(2^5)&=&1\cdot 1+2\cdot 1+4\cdot 2+\cdots+32\cdot 16&=&1+2(1+4+\cdots +4^{4})&=&683\\\\
F(2^6)&=&1\cdot 1+2\cdot 1+4\cdot 2+\cdots+64\cdot 32&=&1+2(1+4+\cdots +4^{5})&=&2731\\\\
G(7)&=&1\cdot 1+(7^2-1)\cdot 7&&&=&337\\\\
G(3^2)&=&1\cdot 1+(3^2-1)\cdot 3+(3^4-3^2)\cdot 9&&&=&673\\\\
G(2^5)&=&1\cdot 1+ (2^2-1)\cdot 2+\cdots+ (2^{10}-2^8)\cdot 2^5&=& 1+6\cdot(1+8+\cdots+8^4)&=&28087\\\\
\end{aligned}$$

よって答えは

$$\begin{aligned}
&2016\cdot 683\cdot 61\cdot 43+\dfrac{1}{2}\cdot (2731\cdot 61\cdot 43\cdot 1017073+28087\cdot 673\cdot 337)\\\\
&=3611682144+\dfrac{1}{2}\cdot(7285713950149+6370159687)\\\\
&=\bm{3649653737062}\\\\
\end{aligned}$$
公式解説
https://onlinemathcontest.com/contests/nfhai2023/editorial/9862
ここの解説では, $f(x,y)$ がwell-definedであることの証明を行います.

$f(x,y)$の定義は,「～～～を満たす最小の正の整数$k$」と書かれています. ここで,「～～～を満たす正の整数$k$」が1つでも存在することが言えれば,well-definedであることが言えます. 

つまり,

$$X=\lbrace k\in\mathbb{Z}_{\gt 0} ; \exists N,\forall n\geq N, s.t. a\_{n+k}\equiv a\_{n}\pmod{p}\rbrace$$

という集合に対して,

$$f(x,y):=\min X$$

としているわけですが,

$X$が空集合でなかったら$X$は(正の整数からなる集合なので)最小値の存在が保証されるわけです.

逆に言えば,$X$が空集合だったらまずいので,$X$が空集合でないことを証明すればよいことになります.

というわけで証明

集合$A$を以下で定めます.

$$A=\lbrace (a,b);0\leq a,b\leq p-1\rbrace$$

このとき,$|A|=p^2\lt \infty$となります.

$\lbrace (a_n\mod{p} ,a_{n+1}\mod{p})\rbrace_{n=0,1,\ldots}$

という列を考えると,各項は全て$A$の元となります.

よって,

$$(a_0,a_1),(a_1,a_2),\ldots,(a_{p^2},a_{p^2+1})$$

の$p^2+1$個の組を考えると,鳩の巣論法によって,一致するペアが存在することが言えます.

つまり,ある$0\leq i\lt j\leq p^2$という整数$i,j$が存在して,

$$(a_i,a_{i+1}) \equiv (a_j,a_{j+1}) \Leftrightarrow a_i\equiv a_j \land a_{i+1}\equiv a_{j+1}$$

となります.このとき,数列$\lbrace a\rbrace$は三項間漸化式で定義されるので,2項さえ一致していればそれ以降はずっと同じです.

よって,

$$a_{i+2}\equiv a_{j+2} \land a_{i+3}\equiv a_{j+3}\land a_{i+4}\equiv a_{j+4}\land \cdots $$

となります.

このとき,$N=i$とすれば,$n\geq N$なる任意の整数$n$について

$$a_{n+(j-i)}\equiv a_{n}\pmod{p}$$

が成立します.

よって,$(j-i)\in X$となるため,$X$が空でないことが言えました.

----

余談の余談

周期の定義をこのようにしていますが,実は$0\leq x\lt p,1\leq y\lt p$の条件ならば$N=0$としても問題ないです.つまり,周期の定義を

「任意の非負整数$n$に対して,$a_{n+k}\equiv a_n$が恒等的に成り立つような最小の正の整数$k$」

としても問題ないわけです. つまり,この問題で問われている状況の範囲内では,正直どっちでもいいのですが,この定義にすると$y=0$のときに問題が発生します.

$y=0,x\neq 0$ならば,$a_n$は

$$0,1,x,x^2,x^3,\ldots $$

となります.$a_n=0$となるような$n$は1個しか無いので周期が定義できません.

特に,$x=0$かつ$y=0$だと

$$0,1,0,0,0,\ldots$$

となると,こっちもうまく周期が定義できなくなってしまいます.

ちなみに,問題で採用している,「ある正整数$N$が存在し～」だと,$y=0$の場合でも問題なく周期を定義することができます.

実はこの問題の原案はシグマを取る範囲が$1\leq y\lt p$ではなく$0\leq y\lt p$でした.さらに,周期の定義も「ある正整数$N$が存在し～」ではないものでした.もしこのままだと$y=0$の場合に$f(x,y)$を定義できなくなるため,問題不備になってしまいます. 

有志コン作問のメンバーに上記の点を指摘されたとき,「とりあえず周期の定義を変えれば問題不備はなくなる」と思ったため現在のような周期の定義に変えたのですが,その上で特性方程式の解に0がある場合の例外処理があまり好きじゃないと思ったため,$y=0$の場合を省くことにしたわけです.

結果,より広い場合でも適用できるような定義だけが残ったという裏話でした.
余談　f(x,y)がwell-definedであることの証明
https://onlinemathcontest.com/contests/nfhai2023/editorial/9862/343
（公式解説の補足）

周期を特性方程式の解$\alpha,\beta$の位数で表す部分は，有限体上の行列を考えると楽である．

以下，$x$を$\mathbb{F_p}$の元とし，$y$を$\mathbb{F_p^\times}$の元とする．
$\mathbb{F_p}$上の行列$A_n$および$B$を
$$
A_n=\begin{pmatrix}a_{n+2}&a_{n+1}\\\ a_{n+1}&a_n\end{pmatrix},\quad B=\begin{pmatrix}x&y\\\ 1&0\end{pmatrix}
$$
により定めると，$a_n$の初期値及び漸化式は次のように言い換えられる：
$$
A_0=\begin{pmatrix}x&1\\\ 1&0\end{pmatrix},\quad A_{n+1}=BA_n.
$$
$A_0,B$は正則行列なので$A_n$は常に正則行列である．
数列$a_n$のmod $p$での周期は，任意の$n$に対して$B^kA_n=A_n$が成り立つような最小の正整数$k$に等しい．
$A_n$は正則なのでこれは$B$の$\mathrm{GL_2}(\mathbb{F_p})$における位数に等しいことがわかる．

$B$の固有方程式$t^2-xt-y=0$（これは数列$a_n$の特性方程式でもある）の解を$\alpha,\beta$とする．$y\in \mathbb{F_p^\times}$より$\alpha,\beta\in \mathbb{F_{p^2}^\times}$である．

- $\alpha,\beta\in \mathbb{F_p^\times}$の場合：$\alpha\neq \beta$ならば$B$は
$$
\begin{pmatrix}\alpha&0\\\ 0&\beta\end{pmatrix}
$$
と対角化されるため，その位数は$\mathrm{lcm}(\mathrm{ord}(\alpha),\mathrm{ord}(\beta))$に等しい．$\alpha=\beta$の場合，$B$はスカラー行列ではないことに注意すると，$B$のJordan標準形は
$$
\begin{pmatrix}\alpha&1\\\ 0&\alpha\end{pmatrix}
$$
となり，その位数は$p\cdot \mathrm{ord}(\alpha)$である．

- $\alpha,\beta\in \mathbb{F_{p^2}^\times} \setminus\mathbb{F_p^\times}$の場合，$\alpha$と$\beta$は$\mathbb{F_p}$上共役である．$\alpha\neq \beta$なので$B$は以下のように対角化される：
$$
\begin{pmatrix}\alpha&0\\\ 0&\beta\end{pmatrix}.
$$
よって$B$の位数は$\mathrm{lcm}(\mathrm{ord}(\alpha),\mathrm{ord}(\beta))=\mathrm{ord}(\alpha)$である．
行列を使って周期を求める方法
https://onlinemathcontest.com/contests/nfhai2023/editorial/9862/363
　$p=2017$ とします（これは素数です）．\
　整数 $x,y$ について，数列 $\\{a_n\\}\_{n=0,1,\ldots}$ を以下の漸化式で定義します：
$$a_0=0, \quad a_1=1, \quad a_{n+2}=xa_{n+1}+ya_n \quad (n=0,1,\ldots)$$
このとき，$\mathrm{mod}~p$ における $a_n$ の周期を $f(x,y)$ として，総和
$$\sum_{x=0}^{p-1}\sum_{y=1}^{p-1}f(x,y)$$
を求めてください.

<details><summary>$f(x,y)$ の定義<\/summary>

　以下を満たす最小の正の整数 $k$ をさします：

- ある正整数 $N$ が存在して，$n\geq N$ なる任意の整数 $n$ について $a_{n+k}\equiv a_n\pmod{p}$ が成り立つ．

なお，このような $k$ は $x,y$ によらず存在します．
<\/details>