cp's OEIS Frontend

A337748 T(m,n) is the least k such that the partial sum of the series Sum_{j=0..k} 1/(m*j+1) is > n, read by ascending antidiagonals.

%I A337748 #19 Oct 19 2020 11:44:48
%S A337748 1,1,3,1,7,10,1,17,56,30,1,43,353,418,82,1,111,2374,7100,3091,226,1,
%T A337748 289,16497,129644,142624,22845,615,1,762,116759,2448436,7078315,
%U A337748 2864677,168803,1673,1,2021,835695,47104189,363380155,386462913,57538579,1247297,4549,1,5389,6025478,916451548,19003186159,53930396770,21100160132,1155693253,9216353,12366
%N A337748 T(m,n) is the least k such that the partial sum of the series Sum_{j=0..k} 1/(m*j+1) is > n, read by ascending antidiagonals.
%C A337748 H(m,k) = Sum_{j=0..k} 1/(m*j+1) may be thought of as a generalized harmonic sequence: H(1,k) = H(k+1), where H(k) = Sum_{j=1..k} 1/j is the harmonic and H(2,k) the odd harmonic sequence. As a consequence, T(1,n) = A002387(n)-1 and T(2,n) = A092315(n). For m > 2, the sequence T(m,n) is not currently in the OEIS. Therefore, I give a detailed derivation of the formulas below, see link "Generalized harmonic series". The parameter c(m) in the formula for T(m,n): c(1) = gamma (Euler's constant), c(2) = gamma+1-log(2). For m > 2, c(m) has to be determined individually. The algorithm is described in the formula section though it is not only one formula.
%C A337748 In the Maxima code, two tests are implemented, see Formula. The first test will probably be passed without exception, the second one limits the terms. For greater T(m,n), the precision must be increased.
%H A337748 Gerhard Kirchner, <a href="/A337748/b337748.txt">Table of n, a(n) for n = 1..496, 31 antidiagonals</a>
%H A337748 Gerhard Kirchner, <a href="/A337748/a337748.pdf">Generalized harmonic sequence</a>
%F A337748 T(m,n) = floor(k1+1) with k1 = exp(m*n-c(m))-(m+2)/(2*m).
%F A337748 For avoiding bad terms, k1 has to pass two tests:
%F A337748 Test 1: floor(k1) = floor(k2) with k2 = k1-1/(24*k1).
%F A337748 Test 2: k1 < d*10^p, where p is the precision of c(m) (number of digits) and d the distance between k1 and the next integer floor(k1) or floor(k1)+1.
%F A337748 The parameter c(m):
%F A337748 c(m) = gamma+m*(1-beta(m)) with (P1) beta(m) = lim_{k->oo} S(m,k) and S(m,k) = Sum_{j=1..k}(1/(m*j)-1/(m*j+1)).
%F A337748 For faster convergence, the infinite sum must be split up: (P2) beta(m) = S(m,k-1) + F(m,k).
%F A337748 F(m,k) can be written as:
%F A337748    (P3) F(m,k) = Sum_{r=1..s-1} b(r)/k^r + R(s) with R(s) ~ b(s)/k^s.
%F A337748 Example: For the precision p = 100 digits, k = 100 is large enough to find s with b(s)/k^s < 10^(-p) where s >= 73.
%F A337748 The coefficients b(r) are defined by the recurrence b(1) = 1/m^2 and, for r>1: (P4) b(r) = ((-1/m)^(r+1)-Sum_{j=1..r-1} (-1)^j*b(r-j)*binomial(r,j+1))/r.
%F A337748               1            m-1           m^2-3*m+2            m^2-2*m+1
%F A337748 (P5) b(1) = -----, b(2) = -----, b(3) = -----------,  b(4) = -----------,
%F A337748              m^2          2*m^3            6*m^4                4*m^5
%e A337748 T(m,n) begins:
%e A337748 m=1: 1   3    10      30        82 ...
%e A337748 m=2: 1   7    56     418      3091 ...
%e A337748 m=3: 1  17   353    7100    142624 ...
%e A337748 m=4: 1  43  2374  129644   7078315 ...
%e A337748 m=5: 1 111 16497 2448436 363380155 ...
%e A337748      .............................
%e A337748 Evaluation of c(3):
%e A337748 1. S(3,99) = 0.1472791427382, see formulas (P1) and (P2).
%e A337748 2. F(3,100) = 10^(-2)/9 + 10^(-4)/27 + 10^(-6)/243 - 10^(-8)/243 ..., see (P5) = 0.001114818889, see formulas (P3) and (P4).
%e A337748 The next summand is 10^(-10)/729 ~ 10^(-13), the precision is approximately  p = 13 digits.
%e A337748 3. beta(3) = 0.1483939616270.
%e A337748 4. c(3) = 3.1320337800204.
%e A337748 Evaluation of T(3,5):
%e A337748 5. Summing up directly with k = 142624: H(3,k-1) = 4.9999999, H(3,k) = 5.0000022. => T(3,5) = 142624.
%e A337748 6. Using the formula with k1 = exp(3*5-c(3))-5/6 = 142623.0446139:
%e A337748    T(3,5) = floor(k1+1) = 142624.
%e A337748 7,1. floor(k1) = floor(k2) is satisfied with k2 = k1-1/(24*k1) = 142623.0446136.
%e A337748 7,2. With d = 0.0446 and p = 13, k1 < d*10^(-13) is satisfied, too.
%o A337748 (Maxima)
%o A337748 block(p:100, km: 100, eps:10^(-p), lim:10^(p), cc:[],
%o A337748 /*"harm-bfile.txt" with terms < lim is saved*/
%o A337748 gam: bfloat(%gamma),
%o A337748 fpprec:p, load(newton1),
%o A337748 fl: openw("harm-bfile.txt"),
%o A337748 c(m):= block(x:0, delta: 1, r:0,
%o A337748 if m=1 then return(gam) else
%o A337748   (for k from 1 thru km-1 do x: x+bfloat(1/(m*k*(m*k+1))),
%o A337748    while delta=0 or delta >eps do
%o A337748       (r: r+1, su:0,
%o A337748       for j from 1 thru r-1 do su: su+a[r-j]*(-1)^j*binomial(r,j+1),
%o A337748       ar: bfloat(-((-1)^r/m^(r+1)+su)/r),  ad: ar/km^r, a: append(a,[ar]),
%o A337748       x: x+bfloat(ad), delta: abs(ad)), return(gam+m*(1-x)))),
%o A337748 mn:0, ok: true,  nr:0,
%o A337748 while ok do (mn: mn+1, a:[], cc: append(cc, [c(mn)]), m:mn+1,
%o A337748    while m>1 and ok do  (m: m-1, n: mn+1-m,
%o A337748    k1: exp(m*n-cc[m])-(m+2)/(2*m), k2: k1-1/(24*k1),
%o A337748    ka: floor(k1+1), kb: floor(k2+1),
%o A337748    if ka>1 then
%o A337748      (d: kb-k2, if d>1-d then d: 1-d, if ka#kb or kb> d*lim then  ok: false),
%o A337748   if ok then (nr: nr+1,  printf(fl, concat(nr, " ", ka)), newline(fl) ) ) ),
%o A337748 close(fl));
%Y A337748 Cf. A002387, A092315.
%K A337748 nonn,tabl
%O A337748 1,3
%A A337748 _Gerhard Kirchner_, Sep 18 2020