cp's OEIS Frontend

A058972 For a rational number p/q let f(p/q) = sum of aliquot divisors of p+q divided by number of divisors of p+q; sequence gives numbers k such that, starting at k/1 and iterating f, an integer is eventually reached.

%I A058972 #30 Feb 09 2022 09:04:28
%S A058972 3,9,15,24,25,29,33,35,50,51,55,57,59,63,73,79,80,81,85,87,89,90,95,
%T A058972 99,105,119,120,121,128,131,139,143,145,169,177,179,181,183,193,195,
%U A058972 201,203,204,211,215,217,218,219,221,225,227,233,247,248,255,273,275,288
%N A058972 For a rational number p/q let f(p/q) = sum of aliquot divisors of p+q divided by number of divisors of p+q; sequence gives numbers k such that, starting at k/1 and iterating f, an integer is eventually reached.
%H A058972 Reinhard Zumkeller, <a href="/A058972/b058972.txt">Table of n, a(n) for n = 1..10000</a>
%H A058972 P. Schogt, <a href="http://arxiv.org/abs/1211.6583">The Wild Number Problem: math or fiction?</a>, arXiv preprint arXiv:1211.6583 [math.HO], 2012. - From _N. J. A. Sloane_, Jan 03 2013
%e A058972 f(9/1) = 8/4 = 2, an integer, so 9 is in the sequence;
%e A058972 f(10/1) = 1/2 and f(1/2)=1/2, so 10 is not in the sequence.
%t A058972 f[r_] := If[init == False && IntegerQ[r], r, init = False; p = Numerator[r]; q = Denominator[r]; d = Most[Divisors[p+q]]; Total[d]/(Length[d]+1)]; ok[n_] := IntegerQ[ init = True; FixedPoint[f, n/1]]; ok[1] = False; A058972 = Select[ Range[300], ok] (* _Jean-François Alcover_, Dec 21 2011 *)
%o A058972 (Haskell)
%o A058972 import Data.Ratio ((%), numerator, denominator)
%o A058972 a058972 n = a058972_list !! (n-1)
%o A058972 a058972_list = map numerator $ filter ((f [])) [1..] where
%o A058972    f ys q = denominator y == 1 || not (y `elem` ys) && f (y : ys) y
%o A058972             where y = a001065 q' % a000005 q'
%o A058972                   q' = numerator q + denominator q
%o A058972 -- _Reinhard Zumkeller_, Jun 15 2013
%o A058972 (PARI) f2(p,q) = (sigma(p+q)-p-q)/numdiv(p+q);
%o A058972 f1(r) = f2(numerator(r), denominator(r));
%o A058972 loop(list) = {my(v=Vecrev(list)); for (i=2, #v, if (v[i] == v[1], return(1)););}
%o A058972 ff(n) = {my(ok=0, m=f2(n,1), list=List()); while(denominator(m) != 1, m = f1(m); listput(list, m); if (loop(list), return (0));); return(m);}
%o A058972 isok(m) = ff(m) > 0; \\ _Michel Marcus_, Feb 09 2022
%Y A058972 Cf. A058971, A058973.
%Y A058972 Cf. A001065, A000005.
%K A058972 nonn,easy,nice
%O A058972 1,1
%A A058972 _N. J. A. Sloane_, Jan 14 2001
%E A058972 Corrected and extended by _Matthew Conroy_, Apr 18 2001