# Iterated polynomials modulo $p$ and divisibility properties

Looking around Math.SE, I found this interesting question about the compositeness of the product of a list of numbers (specifically primes) plus one. As the questioner points out there, this recursive application (starting after the first application) can also be represented by the polynomial $$x(x-1)+1 = x^2-x+1$$. Furthermore, a user comments that for inputs not congruent to $$0,1 \mod 7$$, this will always result in a composite number in a finite number of steps. One can easily see that $$0$$ and $$1$$ $$\mod p$$ for any prime $$p$$ will always be an exception.

An intuitive way to represent this mathematically is a graph structure that maps an input residue to its output residue and checking if there is a path from the input to zero (signalling compositeness for a certain residue class after exactly [path length] iterations). Let's call these residues "good". Let's also call residue systems for that every input ends up at zero "good systems".

Naturally, the first question to ask after looking at the comments is: Can we do "better"?.

Consulting Mathematica$${}^*$$ and looking at residue systems that maximize the relative amount of starting residues that eventually end up at zero, there is only one example within reasonable search space that hardly does any better ($$\mod 139$$ with $$100$$ good residues). I have gone beyond probably $$p > 1400$$ without having found a "better" prime $$p$$ with a ratio bigger than $$\frac{100}{139}$$ and my calculations regularly time out, which is why I am turning to this community in hopes of clever optimizations and bigger crunching power:

Is there a reasonably sized residue system for $$x^2-x+1$$ that has a larger amount of "good" residues relative to its size than $$139$$ with $$\frac{100}{139}$$?

EDIT: This question has been severly edited to split up/remove cluttering extra questions and (currently) unnecessary exposition. Sorry for that.

*Relevant Mathematica code:

PolyModRec[poly_,var_,mod_]:=Table[k \[DirectedEdge] Mod[poly /. var -> k,mod],{k,0,mod-1}]

A search through the first thousand primes reveals the rock star $$p=6469$$, for which the proportion of good residues is $$6413/6469 \approx 0.991343$$. Indeed, the only non-good residues modulo $$6469$$ are: $$\{0,1, 16, 20, 241, 305, 348, 381, 1128, 1135, 1145, 1380, 1658, 1717, 1720, 2019, 2155, 2184, 2264, 2454, 2463, 2562, 2857, 2872, 2977, 2978, 3026, 3137, 3143, 3327, 3333, 3444, 3492, 3493, 3598, 3613, 3908, 4007, 4016, 4206, 4286, 4315, 4451, 4750, 4753, 4812, 5090, 5325, 5335, 5342, 6089, 6122, 6165, 6229, 6450, 6454\}$$
Based on this limit data I would certainly conjecture that there are good-residue proportions that are arbitrarily close to $$1$$. Indeed an optimistic conjecture would be that if we looked at the set of all such proportions, ranging over all primes, then its closure was all of $$[0,1]$$.
Clearly $$0$$ and $$1$$ are never good residues themselves; it's natural to wonder where there is a prime larger than $$7$$ for which those are the only two non-good residues. On the other hand, one could also conjecture that the number of non-good residues tends to infinity (even if the proportion itself can be close to $$1$$).
• I basically found the solutions $x$ to $x^2-x+1\equiv0\pmod p$, then for each one found the solutions $w$ to $w^2-w+1\equiv x\pmod p$, then the solutions $v$ to $v^2-v+1\equiv w\pmod p$, and so on (just using Mathematica's Solve over and over) until there were no more solutions. Basically doing a backwards search through the directed graph in your OP. The first 1000 primes took about 5 minutes on my computer. There's probably even a more efficient way of doing it, especially precomputing the function for each prime modulus like you did.... Aug 27, 2022 at 15:31
• Taking the first $2400$ primes, we also find $\frac{20587}{20641} \approx 0.9973$ which has even less (to be precise $2$ less) non-good residues than $6469$! The calculation is still expensive at $O(p^2)$, since even by precomputing the values we still need to search at most $p$ entries for at most $p$ branches of solutions to previous values. Aug 29, 2022 at 0:25
• I feel like this should be closer to $O(p)$ after (hence including) the precomputation. After precomputing, sort the directed edges by their terminal vertex, so that lookups are faster—we don't have to check all $p$ edges for every new vertex (and we can remove vertices/edges from the graph once they've been accounted for, if that helps). Aug 29, 2022 at 6:32
• I have benchmarked all possible variations of your suggestion and it turns out that deleting list entries actually makes the whole process slower (by a factor of around $2.5$) but that sorting entries slightly speeds it up (e.g. for $p=6469$ it takes around $5.8s$ as opposed to $6.03s$). Because we can always have zero or two entries pointing to the same number, I don't see a nice way of figuring out the index we need to look at in $O(1)$ time.... The biggest improvement still seems to be the reverse tree search quickly discarding primes with small trees. Aug 29, 2022 at 10:03