OFFSET
1,1
COMMENTS
Motivated by Erdős Problem 307, which asks for finite sets of primes P and Q with (Sum_{p in P} 1/p)*(Sum_{q in Q} 1/q) = 1. If A' = B and B' = A is a 2-cycle of the arithmetic derivative (with A, B squarefree, this gives such a solution), then N = A*B satisfies N' + 2*N = (A+B)^2 and N' - 2*N = (A-B)^2, so both N' - 2*N and N' + 2*N are perfect squares. This sequence imposes only the plus condition. The minus condition N' - 2*N = square requires N' >= 2*N, i.e., Sum_{p|N} 1/p >= 2 for squarefree N, which forces N >= 2*3*5*...*277 > 10^112; so no term of the intersection can be found by direct search.
Theorem: for every odd prime l dividing a term k, the cofactor k/l is a nonzero quadratic residue mod l. Proof: for squarefree k = l*m, k' = m + l*m' == m (mod l), so s^2 = k' + 2*k == m (mod l), and m !== 0 (mod l) since k is squarefree.
Semiprime terms k = p*q satisfy (2*p+1)*(2*q+1) = 2*s^2 + 1 where s^2 = k' + 2*k, and necessarily p == q == 1 (mod 4). Proof: (2*s)^2 == -2 (mod r) for every prime r dividing 2*p+1, so r == 1 or 3 (mod 8), hence 2*p+1 == 1 or 3 (mod 8) and p == 0 or 1 (mod 4); p = 2 is impossible since s^2 = 5*q + 2 == 2 (mod 5) has no solution (and q = 5 fails directly).
Conjecturally infinite: for terms of the form 5*q one needs 11*q + 5 = s^2, i.e., q = (s^2 - 5)/11 with s == +-4 (mod 11); s even is equivalent to q == 1 (mod 4). Bunyakovsky's conjecture applied to this quadratic gives infinitely many prime values q, hence infinitely many terms 5*q.
The counting function appears to be ~ 0.25*sqrt(X) (square-density heuristic; the observed ratio is 0.2549 at X = 2*10^5 and 0.2510 at X = 10^6), suggesting a(n) ~ 16*n^2.
LINKS
Thomas F. Bloom, Problem 307, Erdős Problems.
FORMULA
For semiprime terms k = p*q with s^2 = k' + 2*k: (2*p+1)*(2*q+1) = 2*s^2 + 1.
EXAMPLE
130 = 2*5*13 is a term: 130' = 101 and 101 + 2*130 = 361 = 19^2.
145 = 5*29 is a term: 145' = 34 and 34 + 2*145 = 324 = 18^2.
MAPLE
with(numtheory): readlib(issqr): P:=proc(k); if issqrfree(k) and not isprime(k)
then if issqr(k*add(op(2, p)/op(1, p), p=ifactors(k)[2])+2*k) then k; fi; fi; end:
seq(P(i), i=2..10^5); # Paolo P. Lava, Jun 10 2026
PROG
(Python)
from math import isqrt
from sympy import factorint
def ok(n):
f = factorint(n)
if len(f) < 2 or max(f.values()) > 1: return False
d = sum(n//p for p in f)
return isqrt(d + 2*n)**2 == d + 2*n
print([n for n in range(2, 18000) if ok(n)])
(PARI)
ad(n) = my(f=factor(n)); n*sum(i=1, #f~, f[i, 2]/f[i, 1]);
isok(k) = issquarefree(k) && omega(k) > 1 && issquare(ad(k) + 2*k);
CROSSREFS
KEYWORD
nonn
AUTHOR
Vico Bonfioli, Jun 09 2026
STATUS
approved