OFFSET
1,3
COMMENTS
If n is a perfect number then a(n) = 0. But if a(n) = 0, n needs not be perfect, e.g., a(24) = 0, but 24 is not a perfect number. See A064510.
LINKS
Nathaniel Johnston, Table of n, a(n) for n = 1..10000
FORMULA
a(1) = 0, a(2^n) = 1.
a(p) = p-1, a(p^n) = (p^(n+1) - 2*p^n + 1)/(p-1), if p is a prime.
a(n) = n - A117552(n). - Ridouane Oudra, Jan 25 2024
EXAMPLE
a(14) = 4: 14-1 = 13, 13-2 = 11, 11-7 = 4.
a(6) = 0: 6-1 = 5, 5-2 = 3, 3-3 = 0. 6 is a perfect number.
a(35) = 22: 35-1 = 34, 34-5 = 29, 29-7 = 22.
MAPLE
A109883:=proc(n)local d, j, k, m:if(n=1)then return 0:fi:j:=1:m:=n:d:=divisors(n); k:=nops(d):for j from 1 to k do m:=m-d[j]:if(m<d[j+1])then return m:fi:od:end: # Nathaniel Johnston, Apr 15 2011
MATHEMATICA
subtract = If[ #1 < #2, Throw[ #1], #1 - #2]&;
a[n_] := Catch @ Fold[subtract, n, Divisors @ n]
Table[ a[n], {n, 80}] (* Bobby R. Treat (DrBob(AT)bigfoot.com), Jul 14 2005 *)
PROG
(PARI) a(n) = {my(r = n); fordiv(n, d, if (r < d, return (r)); r -= d; ); 0; } \\ Michel Marcus, Dec 28 2018
(Python)
from sympy import divisors
def A109883(n):
if n == 1: return 0
s = n
for d in divisors(n)[:-1]:
if s < d: break
s -= d
return s
print([A109883(n) for n in range(1, 80)]) # Michael S. Branicky, Mar 31 2024
CROSSREFS
KEYWORD
easy,nonn
AUTHOR
Amarnath Murthy, Jul 11 2005
EXTENSIONS
More terms from Jason Earls and Robert G. Wilson v, Jul 12 2005
STATUS
approved