|
%I #19 Jul 27 2012 12:32:12
%S 1,1,7,1,23,1,47,1,79,1,113,89,23,73,31,1,283,89,113,139,173,67,47,1,
%T 619,131,109,83,113,211,191,1,1087,1,1223,1,1367,1,1511,367,83,1,1847,
%U 1,2017,317,571,241,199,1,2593,367,211,271,223,1,3229,1117,239,1
%N Numerators of a series with denominators n! and sum 1.
%C a(n) is the greatest integer x such that gcd(x,n!) = 1 and x/n! < 1 - sum_{j=2}^{n-1} a(j)/j!.
%C The infinite series sum_{n=2}^infinity a(n)/n! = 1
%C For each n, either a(n) = 1 or a(n) >= n+1.
%H Robert Israel and T. D. Noe, <a href="/A214686/b214686.txt">Table of n, a(n) for n = 2..10000</a> (first 1000 terms from Robert Israel)
%e 1 - a(2)/2! = 1/2 = 3/3!, gcd(2,3!)>1 so a(3) = 1.
%e 1 - a(2)/2! - a(3)/3! = 8/4! so a(4) = 7.
%p N:= 100; a[2]:= 1; R[2]:= 1/2;
%p for j from 3 to N do
%p T:= R[j-1] *j!;
%p for x from T-1 by -1 while igcd(x,j!) > 1 do end do;
%p a[j]:= x;
%p R[j]:= R[j-1] - x/j!
%p end do:
%p seq(a[j],j=2..N);
%t s = 1; Table[d = n! s; q = If[d <= n, 1, If[d <= n^2, Prime[PrimePi[d]], Print["d > n^2"]; Abort[]]]; s = s - q/n!; q, {n, 2, 100}] (* _T. D. Noe_, Jul 27 2012 *)
%o (Sage)
%o def A214686_list(n) :
%o a = [1]; R = 1/2
%o for j in (3..n+1) :
%o J = factorial(j)
%o T = R * J
%o for x in range(T-1, -1, -1) :
%o if gcd(x, J) == 1 : break
%o a.append(x)
%o R = R - x / J
%o return a
%o A214686_list(51) # _Peter Luschny_, Jul 27 2012
%K nonn
%O 2,3
%A _Robert Israel_, Jul 25 2012
|
