# What is the next “Tribonacci-like” pseudoprime?

Given the three roots of $x^3=x^2+x+1$. Then we get the tribonacci-like sequence,

$B_n = x_1^n+x_2^n+x_3^n = 3, 1, 3, 7, 11, 21, 39, 71, 131,\dots$

where $B_n = B_{n-1}+B_{n-2}+B_{n-3}$, and the index starts with n = 0. This is A001644 of the OEIS.

I noticed that, like for the Lucas numbers, if n is prime then n divides $B_n-1$. (For example $B_7=71$ minus 1 is div by 7.) But it also has pseudoprimes, one of which is n = 182.

Question: What is the next pseudoprime?

(If I did my Mathematica session correctly, then this is the ONLY one for n < 5000, versus the Lucas numbers which has 5 pseudoprimes within that range.)

-
I don't have Mathematica to work with and doubt I could put something together quickly enough, but how have you been doing your calculations? It shouldn't be hard to extend this up into the $10^6$ range with a slightly-smarter algorithm: instead of computing the sequence once and thus dealing with exponentially-large terms, instead for each number $n$ run the recurrence mod $n$ - and to avoid taking $O(n)$ time to compute the recurrence (and thus having an overall $O(n^2)$ algorithm) you can find the $n$th term in $O(\log n)$ time via the usual matrix exponentiation methods... –  Steven Stadnicki May 22 '12 at 16:58
(For the details of the matrix method, if you haven't seen them before, have a look at my answer to math.stackexchange.com/questions/41667/… ) –  Steven Stadnicki May 22 '12 at 17:06
I can confirm your result that there is only one pseudoprime for $n<5000$, and additionally say that there are no pseudoprimes for $n<16000$, which is as far as my code has got! It's been running for about ten minutes so far. I'm going to leave it running overnight. It's a pretty quick and dirty implementation in Haskell - it computes $B_n$ in $O(n)$ time and has a fairly slow prime number checker. I put the code here if you want to use it. –  Chris Taylor May 22 '12 at 17:28
@Steven I did it in a most brute force manner. $B_n \approx x_1^n$ (the real root) for large n, since the contribution of the complex conjugates $x_2^n, x_3^n$ rapidly diminishes. But the calculation gets longer when n gets into the thousands. –  Tito Piezas III May 22 '12 at 17:32
Are you computing $B_n$ using symbolic algebra or floating point arithmetic? –  Chris Taylor May 22 '12 at 17:35

I got together a quick C++ program using the matrix method mentioned in the comments (and described in detail in Fibonacci, tribonacci and other similar sequences) to search for more pseudoprimes; the next several after $182$ are $25201 (=11\cdot 29\cdot 79)$, $54289 (=233^2!)$ and $63618 (=2\cdot 3\cdot 23\cdot 461)$. There are another three between $10^5$ and $10^6$, and then five more between $10^6$ and $10^7$ (for a total of twelve $\leq 10^7$); it seems likely that the count is growing logarithmically or polylogarithmically, but the numbers involved are just too small to make a reasonable conjecture on the total number.
ADDED: while I wait for my OEIS application to go through, here's the full list less than $10^9$. This took about 3 hours for my less-than-optimal C++ code to compute, so it wouldn't be too hard to go to $10^{10}$ if someone really wanted:
$182, 25201, 54289, 63618, 194390, 750890, 804055, 1889041, 2487941, 3542533, 3761251, 6829689, 12032021, 12649337, 18002881, 22586257, 28250321, 68355001, 72374401, 74458790, 79351441, 100595461, 116406374, 123872111, 191239529, 221265526, 225853633, 248947777, 338458807, 358313761, 379732501, 381427201, 509551201, 517567051, 813015901, 859481921$
(I'm currently waiting for my OEIS registration to be approved; as soon as it is, I'll submit the full sequence for $n\leq 10^7$ there) –  Steven Stadnicki May 22 '12 at 19:02
For example, given the roots of $x^4=x^3+x^2+x+1$ and the sequence $F_n = x_1^n+x_2^n+x_3^n+x_4^n$, then the first 6 pseudoprimes n such that n divides $F_n-1$ are $n = 5^2, 7^2, 5^3, 13^2, 7^3, 23^2$. Interesting that they are all powers. –  Tito Piezas III May 22 '12 at 19:43
(Note that if my hypothesis about polylogarithmic growth is correct, then heuristically there should only be a finite number of perfect powers among the pseudoprimes; the probability that a number $n$ is a perfect power is $\Theta(1/n)$ and the probability that it's a pseudoprime is $\Theta(\log^k(n)/n)$, so the probability that it's both is $\Theta(\log^k(n)/n^2)$ and summing this over all $n$ yields a finite result.) –  Steven Stadnicki May 22 '12 at 22:15