login
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.
1

%I #19 Oct 19 2020 11:44:48

%S 1,1,3,1,7,10,1,17,56,30,1,43,353,418,82,1,111,2374,7100,3091,226,1,

%T 289,16497,129644,142624,22845,615,1,762,116759,2448436,7078315,

%U 2864677,168803,1673,1,2021,835695,47104189,363380155,386462913,57538579,1247297,4549,1,5389,6025478,916451548,19003186159,53930396770,21100160132,1155693253,9216353,12366

%N 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 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 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 Gerhard Kirchner, <a href="/A337748/b337748.txt">Table of n, a(n) for n = 1..496, 31 antidiagonals</a>

%H Gerhard Kirchner, <a href="/A337748/a337748.pdf">Generalized harmonic sequence</a>

%F T(m,n) = floor(k1+1) with k1 = exp(m*n-c(m))-(m+2)/(2*m).

%F For avoiding bad terms, k1 has to pass two tests:

%F Test 1: floor(k1) = floor(k2) with k2 = k1-1/(24*k1).

%F 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 The parameter c(m):

%F 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 For faster convergence, the infinite sum must be split up: (P2) beta(m) = S(m,k-1) + F(m,k).

%F F(m,k) can be written as:

%F (P3) F(m,k) = Sum_{r=1..s-1} b(r)/k^r + R(s) with R(s) ~ b(s)/k^s.

%F 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 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 1 m-1 m^2-3*m+2 m^2-2*m+1

%F (P5) b(1) = -----, b(2) = -----, b(3) = -----------, b(4) = -----------,

%F m^2 2*m^3 6*m^4 4*m^5

%e T(m,n) begins:

%e m=1: 1 3 10 30 82 ...

%e m=2: 1 7 56 418 3091 ...

%e m=3: 1 17 353 7100 142624 ...

%e m=4: 1 43 2374 129644 7078315 ...

%e m=5: 1 111 16497 2448436 363380155 ...

%e .............................

%e Evaluation of c(3):

%e 1. S(3,99) = 0.1472791427382, see formulas (P1) and (P2).

%e 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 The next summand is 10^(-10)/729 ~ 10^(-13), the precision is approximately p = 13 digits.

%e 3. beta(3) = 0.1483939616270.

%e 4. c(3) = 3.1320337800204.

%e Evaluation of T(3,5):

%e 5. Summing up directly with k = 142624: H(3,k-1) = 4.9999999, H(3,k) = 5.0000022. => T(3,5) = 142624.

%e 6. Using the formula with k1 = exp(3*5-c(3))-5/6 = 142623.0446139:

%e T(3,5) = floor(k1+1) = 142624.

%e 7,1. floor(k1) = floor(k2) is satisfied with k2 = k1-1/(24*k1) = 142623.0446136.

%e 7,2. With d = 0.0446 and p = 13, k1 < d*10^(-13) is satisfied, too.

%o (Maxima)

%o block(p:100, km: 100, eps:10^(-p), lim:10^(p), cc:[],

%o /*"harm-bfile.txt" with terms < lim is saved*/

%o gam: bfloat(%gamma),

%o fpprec:p, load(newton1),

%o fl: openw("harm-bfile.txt"),

%o c(m):= block(x:0, delta: 1, r:0,

%o if m=1 then return(gam) else

%o (for k from 1 thru km-1 do x: x+bfloat(1/(m*k*(m*k+1))),

%o while delta=0 or delta >eps do

%o (r: r+1, su:0,

%o for j from 1 thru r-1 do su: su+a[r-j]*(-1)^j*binomial(r,j+1),

%o ar: bfloat(-((-1)^r/m^(r+1)+su)/r), ad: ar/km^r, a: append(a,[ar]),

%o x: x+bfloat(ad), delta: abs(ad)), return(gam+m*(1-x)))),

%o mn:0, ok: true, nr:0,

%o while ok do (mn: mn+1, a:[], cc: append(cc, [c(mn)]), m:mn+1,

%o while m>1 and ok do (m: m-1, n: mn+1-m,

%o k1: exp(m*n-cc[m])-(m+2)/(2*m), k2: k1-1/(24*k1),

%o ka: floor(k1+1), kb: floor(k2+1),

%o if ka>1 then

%o (d: kb-k2, if d>1-d then d: 1-d, if ka#kb or kb> d*lim then ok: false),

%o if ok then (nr: nr+1, printf(fl, concat(nr, " ", ka)), newline(fl) ) ) ),

%o close(fl));

%Y Cf. A002387, A092315.

%K nonn,tabl

%O 1,3

%A _Gerhard Kirchner_, Sep 18 2020