Pillai's arithmetical function is defined by $P(n)=\sum_{k=1}^n\gcd(k,n)$.
Is it true that $1+P(n)\equiv 0\pmod{n}$ iff $n$ is prime?
It was stated without proof in 2014 on the OEIS A018804 and has since been reverted as a conjecture. A question posted last week on MSE didn't get a full proof: the case of squarefree composites with more than 3 factors is completely open.
Do you have a reference or a proof for that case?
Copying my answer from MSE
Surprisingly, this is not true! A counterexample is given by $n=3\cdot37\cdot43\cdot42307\cdot116341=23492890653051$ which has $$\frac{P(n)+1}{n}=26$$
An answer to the linked question shows that $P(n)$ is equal to $$\sum_{d|n}\phi(n/d)d$$ and that $P(n)$ is multiplicative. And as other answers say, it's enough to consider $n$ square-free. So let's say $n=p_1p_2\cdots p_k$ with $k\ge 2$. Then $$P(n)=\prod_{i=1}^k(2p_i-1)$$
Consider $$\frac{1+P(n)}{n}$$ which we are trying to show is never an integer. We can bound $\frac{1+P(n)}{n}$ as being strictly between $\prod_{i=1}^k\frac{2p_i-1}{p_i}$ and $2^k$ when $k\ge 2$.
Write $p_1<\cdots < p_k$. Suppose we fix the $m$ smallest prime factors, and then want to bound the $(m+1)$-th smallest one. Call $r=p_{m+1}\cdots p_k$. We have $$p_1\cdots p_mr|[(2p_1-1)\cdots(2p_m-1)P(r)+1]$$ Since $P(r)\le 2^{k-m}r-1$, we have $$[(2p_1-1)\cdots(2p_m-1)P(r)+1]\le (2p_1-1)\cdots(2p_m-1)[2^{k-m}r-1]+1$$
If $c\cdot p_1\cdots p_mr\le (2p_1-1)\cdots(2p_m-1)[2^{k-m}r-1]+1< (2p_1-1)\cdots (2p_m-1)2^{k-m}r$, then $$c< \frac{(2p_1-1)\cdots(2p_m-1)}{p_1\cdots p_m}2^{k-m}$$
Note that this cannot be an integer unless $p_1=2$ and $m=1$. If it were, then $p_1\cdots p_m$ would divide $P(n)$ so couldn't divide $P(n)+1$. More generally, if $p_i|(2p_j-1)$ for any $i,j$, it is already impossible for $n|[P(n)+1]$ to be the case.
But then suppose $p_{m+1}$ and subsequent primes are too large (e.g. all at least $N$). $\frac{2x-1}{x}$ is increasing in $x$ for $x>0$, so we have $$\frac{(2p_1-1)\cdots(2p_m-1)P(r)}{p_1\cdots p_mr}\ge \frac{(2p_1-1)\cdots(2p_m-1)}{p_1\cdots p_m}(2-\epsilon)^{k-m}$$
Since it has to lie in between these bounds, we only have to check a finite amount of $p_{m+1}$.
For example, with $m=1$, $p_1=3$, and $k=4$, we have that $c\le 40/3$, so $c\le 13$. If $\frac{2N-1}{N}>\sqrt[3]{39/5}$, which happens when $N\ge 60$, there is no solution. So we would only need to check $p_2\le 59$.
If $m=0$ (i.e. no primes are given) and $k\ge 2$, we have the strict inequality of $c<2^k$. Thus, if $\left(\frac{2N-1}{N}\right)^k>2^{k}-1$, there is no solution. For $k=4$, this tells us $p_1\le 31$.
In general, for a given $k$, there are only a finite amount of cases that we need to check. Using a Python program, I checked it for $k=4$, and it holds.
But it doesn't hold for $k=5$, and we have a counterexample! With the primes $\{3,37,43,42307,116341\}$, we have $n=23492890653051$, $P(n)=610815156979325$, and $\frac{P(n)+1}{n}=26$. I don't know if this is the smallest counterexample though.
