A new generalisation of Fermat's little theorem? I would like to share a personal result, which I interpret as a "generalisation of Fermat's little theorem". I know that there exist already several generalisations (e.g. Euler's) but I would like to know whether the result I found is something original or not. Let me explain...
Some years ago, my researches on binary representation of numbers lead me to this sequence of naturals: 1, 2, 1, 2, 3, 6, 9, 18, 56, ... . After some researches, I found that this was already indexed as A001037 sequence in OEIS; this is known as number of "binary Lyndon words of length $n$" or "$n$-bead necklaces with beads of 2 colors". Such sequence can easily be generalised to other bases than 2. After naming $\lambda_a(n)$ the number of $n$-bead necklaces with beads of $a$ colors with $a \ge 2$, I found that
$$ a^n = \sum_{d|n}  d \lambda_a(d) $$
For example,
$$ 2^6 = 1\lambda_2(1) + 2\lambda_2(2) + 3\lambda_2(3) + 6\lambda_2(6) = 2 + 2 + 6 + 54 = 64 $$
This is NOT an original result since you can find this relationship on the afore-mentioned OEIS page, at least. 
QUESTION 1: Can someone provides references for a proof of this equation (papers, books, ...)?
Now, after defining
$$ \sigma_a(n) \triangleq \sum_{\substack{ d|n \\ 1<d<n}} d \lambda_a(d) $$
the previous equation can be rewritten
$$ a^n = a + n\lambda_a(n) + \sigma_a(n) $$
Then, using modular arithmetic, we obtain a generalisation of Fermat's little theorem:
$$ a^n \equiv a + \sigma_a(n) \pmod {n} $$
Note that this is indeed a generalisation since it applies on any natural $n$ and it boils down to 
$$ a^p \equiv a \pmod {p} $$
for any prime $p$, that is the Fermat's little theorem.
This relationship looks appealing to me, at least for studying Fermat pseudoprimes: these are the composites $n$ having the "rare" property that $\sigma_a(n) \equiv 0 \pmod {n}$.
QUESTION 2: Is this generalisation of Fermat's little theorem a known result? If yes, could you please provide references?
Thanks,
PS: Please, be indulgent: this is my very first post on MO!
 A: Here's a massive generalization, with references.
Consider a sequence $\{a_n \}_{n\ge 1}$ of integers. We'll say it satisfies the necklace congruences if $a_{n} \equiv a_{n/p} \mod n$ whenever $p \mid n$ ($p$ prime).
Here are some equivalent formulations:


*

*$a_{p^{k+1} m} \equiv a_{p^k m} \mod {p^{k+1}}$ for all $p,m,k$

*$\sum_{d \mid n} a_d \mu(n/d) \equiv 0 \mod n$ for all $n$

*The Mobius function can be replaced by any arithmetic function $f$ satisfying $\sum_{d \mid n} f(d) =0 \mod n, f(1) = \pm 1$.


So any such sequence gives rise to an integer sequence $\lambda_{\tilde{a}}(n) := \frac{1}{n}\sum_{d \mid n} a_d \mu(n/d)$ which satisfies $a_n = \sum_{d \mid n} d \lambda_{\tilde{a}}(n)$ by Mobius inversion. One can then define $\sigma_{\tilde{a}}(n):=\sum_{d \mid n, 1<d<n} d\lambda_{\tilde{a}}(d)$ and obtain the following generalization of Fermat's little theorem:
$a_n \equiv a_1 + \sigma_{\tilde{a}}(n) \mod n$
Now, let's go back to necklace congruences. A less obvious equivalence is the following:
Theorem: $\{a_n\}_{n \ge 1}$ satisfies the necklace congruences iff $\zeta_{a}(x):=\exp(\sum_{n\ge1} \frac{a_n}{n}x^n)$ has integer coefficients.
Proof: Write, $\zeta_a$ formally as $\prod_{n\ge 1} (1-x^n)^{-\frac{b_n}{n}}$. The $b_n$'s are uniquely determined, and in fact (by taking logarithmic derivative) it can be seen that they are $b_n = \sum_{d \mid n} a_d \mu(n/d)$. Now one direction is immediate and the other requires an inductive argument. $\blacksquare$
(Reference: Exercise 5.2 (and its solution) in Richard Stanley's book "Enumerative Combinatorics, vol. 2")
As Sergei remarked, for any integer square matrix $A$, the sequence $a_n := \text{Tr}(A^n)$ satisfies the necklace congruences. This is immediate from the last theorem, as the corresponding "zeta" function $\zeta_a$ is just:
$\exp( \sum_{n \ge 1} \frac{\text{Tr}((Ax)^n)}{n}) = \exp (\text{Tr} (- \ln (I -Ax)))=\det(I-Ax)^{-1} \in \mathbb{Z}[x]$
In particular, by taking a 1 on 1 matrix, we recover your observation. 
Another generalization is $a_n = [x^n]f^n(x)$, where $f \in \mathbb{Z}[[x]]$. By taking the linear polynomial $f(x)=1+ax$ we recover your example, by taking $f(x)=(1+x)^m$ we get $a_n = \binom{nm}{n}$.
The really interesting feature is that this notion has a local version. Specifically, theorems of Dieudonne-Dwork and Hazewinkel give the following result:
Theorem: Given $\{ a_n \}_{n\ge 1} \subseteq \mathbb{Q}_{p}$ ($p$-adic numbers), the following are equivalent:


*

*$\zeta_{a}(x) := \exp(\sum_{n\ge 1} \frac{a_n }{n}x^n)\in\mathbb{Z}_{p}[x]$

*$\exp( \sum_{n \ge 1} p \frac{a_n - a_{n/p}}{n} x^n) \in 1+px\mathbb{Z}_{p}[x]$ (I defnite $a_{n/p}=0$ if $p\nmid n$)

*$\sum_{n \ge 1}  \frac{a_n - a_{n/p}}{n} x^n) \in \mathbb{Z}_{p}[x]$

*$p \mid n \implies a_n \equiv a_{n/p} \mod {n\mathbb{Z}_{p}}$


Applying this simultaneously for all $p$, we recover the previous theorem.
A reference for this theorem is "A Course in p-adic Analysis" by Alain M. Robert (the theorems are stated there in a much more general context).
A: In this paper (in Russian) : http://www.mathnet.ru/php/archive.phtml?wshow=paper&jrnid=mp&paperid=238&option_lang=eng
there are discussion of this result, including theorem of Gauss, generalization to matrices by V.I.Arnold and so on.
