Let $\mathbb{Z}^+$ denote the natural numbers. Find the number of functions $f:\mathbb{Z}^+\to \{0,1,\cdots, 16\}$ so that $f(x+17)=f(x), f(x^2)\equiv f(x)^2+15\mod 17$ for all $x\ge 1$.
By plugging in $x=0,$ we get $f(0) = 16$ or $2$ (by the quadratic formula, the roots are $(1\pm 3)\cdot 9$). Also, the relation $f(x+17)=f(x)$ implies that it suffices to only consider the values $\{f(0),\cdots, f(16)\}$; these values uniquely determine f. I think it could be useful to come up with some sort of inductive formula by finding some sort of "inductive parameter". We have $f(1) \equiv f(1)^2 + 15\mod 17$, so $f(1)=16$ or $2$. I'm not sure about how to determine the values of f on the other elements of $\{0,1,\cdots, 16\}$; this seems to involve some sort of nonlinear system. For instance, $f(4)\equiv f(2)^2 +15, f(9)\equiv f(3)^2 + 15,f(16)\equiv f(4)^2 + 15, f(1)\equiv f(16)^2 + 15,$ etc. At least from $f(1)\equiv f(16)^2+15$, we can deduce that $f(16)^2\equiv 1,4$ so $f(16)\equiv \pm 1, \pm 2$. How can I solve this problem?
As you said, $f(1)$ is either $2$ or $16$.
First assume $f(1)=2$. Then, because of $2+f(1) \equiv f^2(1)\equiv f^2(16)$, we conclude that $f(16)$ is either $2$ or $-2$, and moreover $f(4), f(13)\in \{ 2,0,-2\}.$ Note that:
$$\\ 2+f(16) \equiv f^2(4)\equiv f^2(13)\\2+f(13) \equiv f^2(8)\equiv f^2(9)\\2+f(8) \equiv f^2(5)\equiv f^2(12)\\2+f(9) \equiv f^2(3)\equiv f^2(14),$$ and, similarly; $$\\ 2+f(16) \equiv f^2(4)\equiv f^2(13)\\2+f(4) \equiv f^2(2)\equiv f^2(15)\\2+f(2) \equiv f^2(6)\equiv f^2(11)\\2+f(15) \equiv f^2(7)\equiv f^2(10).$$ Three cases happen:
Case $1$: $f(13)=2.$
In this case, for choosing $f(8), f(9), f(5), f(12), f(3), f(14)$, there are $25$ options. Because:
if $(f(8), f(9))=(2,2)$, there are $16$ options for choosing $f(5), f(12), f(3), f(14)$;
if $(f(8), f(9))=(2,-2)$, there are $4$ options for choosing $f(5), f(12), f(3), f(14)$;
if $(f(8), f(9))=(-2,2)$, there are $4$ options for choosing $f(5), f(12), f(3), f(14)$;
if $(f(8), f(9))=(-2,-2)$, there are $1$ option for choosing $f(5), f(12), f(3), f(14)$.
Similarly, if $f(4)=2$, we have $25$ options for choosing $f(2), f(15), f(6), f(11), f(7), f(10).$
Case $2$: $f(13)=-2$.
In this case we have $16$ options for choosing $f(8), f(9), f(5), f(12), f(3), f(14)$.
Similarly, if $f(4)=-2$, we have $16$ options for choosing $f(2), f(15), f(6), f(11), f(7), f(10).$
Case $3$: $f(13)=0$.
In this case we have $64$ options for choosing $f(8), f(9), f(5), f(12), f(3), f(14)$.
Similarly, if $f(4)=0$, we have $64$ options for choosing $f(2), f(15), f(6), f(11), f(7), f(10).$
While still assuming $f(1)=2$;
we have $25\times 25$ options if $f(16)=2$ and $(f(4), f(13))=(2,2);$
we have $25\times 16$ options if $f(16)=2$ and $(f(4), f(13))=(2,-2);$
we have $25\times 16$ options if $f(16)=2$ and $(f(4), f(13))=(-2,2);$
we have $16\times 16$ options if $f(16)=2$ and $(f(4), f(13))=(-2,-2);$
we have $64\times 64$ options if $f(16)=-2$.
On the other hand assuming $f(1)=16$, it is very easy to observe that $f(16)=16$ (otherwise $2+f(16) \equiv f^2(4)\equiv f^2(13)$ has no solution.) Similarly, $f(4)=f(13)=16$, and consequently $f(8)=f(9)=f(2)=f(15)=16$, and for choosing $f(5), f(12), f(3), f(4), f(6), f(11), f(7), f(10)$ we have $2^8$ options.
In total for choosing $f(1), ..., f(16)$, we have: $(25\times 25)+ 2(25 \times 16)+(16\times 16)+(64\times 64) +2^8$ options.
Finally, note that there exist $2$ options for $f(17)$ as well.