Cut the Cards - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: Not HP Calculators (/forum-7.html) +--- Forum: Not quite HP Calculators - but related (/forum-8.html) +--- Thread: Cut the Cards (/thread-15399.html) Cut the Cards - David Hayden - 07-30-2020 08:00 PM PPC Journal V5N4 Page 13 poses a question about a deck of cards. Choose a card at random from a deck of 52 and put the card back. What is the average number of choices required to see all 52 cards? I simulated the experiment in C++ and got the right answer (236) but I can't figure how to do this analytically. Is it possible? If you simplify the problem down to a deck of 2 cards then it's easier: you choose the first card, and then it's just a matter of the average number of choices to get the second: 1*1/2 + 2*1/4 + 3*1/8 ..., but can that be extended to more cards? Dave RE: Cut the Cards - Albert Chan - 07-30-2020 08:58 PM (07-30-2020 08:00 PM)David Hayden Wrote:  Choose a card at random from a deck of 52 and put the card back. What is the average number of choices required to see all 52 cards? First card picked must never been seen, thus only 1 try needed to get 1st card. Probability of getting the next "new" card = 51/52, or averaged 52/51 tries to get 2nd card ... We can get averaged picks either way: XCas> sum(1/t, t=1 .. 52) * 52.0 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 235.978285436 XCas> int((t^52-1)/(t-1), t=0 .. 1) * 52.0 ﻿ ﻿ ﻿ → 235.978285436 → averaged 236 picks to see all 52 cards RE: Cut the Cards - Jim Horn - 07-30-2020 09:49 PM So, it approaches x*(ln(x)+0.57721), where the latter constant is the Euler-Mascheroni constant. This isn't as accurate as doing it analytically but is good for fast calculations. For instance, if you picked a random second of a year every second, how many years would it take? About 17.84 total. RE: Cut the Cards - Albert Chan - 07-30-2020 10:21 PM Here is simulations, using Lua Note: math.random(n) return 1 .. n, each with same probability. Code: function run(tests,cards)   -- averaged picks     local n = 0             -- bad tries     for i=1,tests do         for t=1,cards do             while math.random(cards) < t do n=n+1 end         end     end     return n/tests + cards end lua> run(1e6, 52) 235.970842 RE: Cut the Cards - John Keith - 07-31-2020 12:24 AM (07-30-2020 09:49 PM)Jim Horn Wrote:  So, it approaches x*(ln(x)+0.57721), where the latter constant is the Euler-Mascheroni constant. Or in other words, n*H(n), where H(n) is the nth harmonic number. (Albert's first formula in post #2) (07-30-2020 09:49 PM)Jim Horn Wrote:  This isn't as accurate as doing it analytically but is good for fast calculations. Actually it is very accurate as the first image in the Wikipedia link shows: LN(n)+gamma is the asymptotic limit of H(n) as n approaches infinity. RE: Cut the Cards - Albert Chan - 08-21-2020 11:00 PM (07-30-2020 08:58 PM)Albert Chan Wrote:  First card picked must never been seen, thus only 1 try needed to get 1st card. Probability of getting the next "new" card = 51/52, or averaged 52/51 tries to get 2nd card ... I might jumped too quickly to say averaged 1/p tries to get the next card. Let probability of getting the card = p ; of not getting it = q = 1-p Expected tries = p + 2qp + 3q²p + 4q³p + ... = (1-q) (1 + 2q + 3q² + 4q³ + ...) = (1 + 2q + 3q² + 4q³ + ...) - (q + 2q² + 3q³ + ...) = 1 + q + q² + q³ + ... = 1/(1-q) = 1/p For n > 0, we have: ﻿ ﻿ ﻿ ﻿ $$\Large {1 \over (1-x)^n}\normalsize = \displaystyle \sum_{k=0}^∞ \binom{-n}{k}(-x)^k = \sum_{k=0}^∞ \binom{k+n-1}{n-1} x^k$$ Examples: XCas> series(1/(1-x)^2, x, polynom) → 1+2*x + 3*x^2 + 4*x^3 + 5*x^4 + 6*x^5 XCas> series(1/(1-x)^3, x, polynom) → 1+3*x + 6*x^2+10*x^3+15*x^4+21*x^5 ﻿ ﻿ // triangular numbers coefficients XCas> series(1/(1-x)^4, x, polynom) → 1+4*x+10*x^2+20*x^3+35*x^4+56*x^5 ﻿ ﻿ // tetrahedral numbers coefficients RE: Cut the Cards - Gerson W. Barbosa - 08-24-2020 01:57 PM (07-31-2020 12:24 AM)John Keith Wrote:   (07-30-2020 09:49 PM)Jim Horn Wrote:  So, it approaches x*(ln(x)+0.57721), where the latter constant is the Euler-Mascheroni constant. Or in other words, n*H(n), where H(n) is the nth harmonic number. (Albert's first formula in post #2) The H(n) function can be easily defined on the HP 50g: « 1. + Psi 1. Psi - » 'H' STO For instance, 52 ENTER ENTER H × → 235.978285436 RE: Cut the Cards - pinkman - 08-24-2020 09:49 PM Thanks for all those formulas and explanations. Here is an implementation of the brute force algorithm for the experimental analysis: Code: // Forward declaration CRDS(N); // Main func: // N is the number of cards // L is the number of tests // Returns average number of peeks EXPORT CARDS(N,L) BEGIN  LOCAL I, S := 0;  FOR I FROM 1 TO L DO   S := S + CRDS(N);   END;  RETURN S/L; END; // Peek a card until all cards //   are peeked at least once EXPORT CRDS(N) BEGIN  LOCAL LL, C := 0;  LL:=MAKELIST(0,A,1,N);  WHILE ΠLIST(LL)=0 DO   LL(RANDINT(1,N)):=1;   C := C+1;  END;   RETURN C; END; Usage: CARDS(number-of-cards, number-of-analysis) Example: CARDS(52, 1000) Returns: 235.951 RE: Cut the Cards - Albert Chan - 08-25-2020 06:14 PM (07-30-2020 09:49 PM)Jim Horn Wrote:  So, it approaches x*(ln(x)+0.57721), where the latter constant is the Euler-Mascheroni constant. Replacing ln(x) by ln(x+0.5) gives much better estimate for ψ(x+1). ψ(x+1) = slope of ln(gamma(t)) at t = x+1. When x is large: $$\psi(x+1) ≈ {\ln\Gamma(x+2)\;-\;\ln\Gamma(x)\over 2}=\ln\sqrt{(x+1)(x)} ≈ \ln(x+0.5)﻿$$ >>> x = 52 >>> H = log(x+0.5) + 0.57722 >>> print format(x*H, 'g') 235.978 see https://en.wikipedia.org/wiki/Digamma_function#Computation_and_approximation RE: Cut the Cards - Gerson W. Barbosa - 08-25-2020 11:41 PM FWIW, « 1. Psi + EXP 6. OVER INV 12. + / + 1. - » 'H­-¹' STO Example: 1100 H → 7.58073560027 H­-¹ → 1100 H-­¹ should be accurate to 11 or 12 digits for arguments equal or greater than 7.5 RE: Cut the Cards - Albert Chan - 08-26-2020 03:06 AM (08-25-2020 11:41 PM)Gerson W. Barbosa Wrote:  FWIW, « 1. Psi + EXP 6. OVER INV 12. + / + 1. - » 'H­-¹' STO Example: 1100 H → 7.58073560027 H­-¹ → 1100 H-­¹ should be accurate to 11 or 12 digits for arguments equal or greater than 7.5 Here is another version, using asymptotic estimate: $$\quad\exp(\psi(x+0.5)) ≈ x + \large{1\over24x}$$ Let k = exp(ψ(x+0.5)), n = x - 0.5 → k = exp(H(n) + ψ(1)) ≈ x + 1/(24x) If x is big, x ≈ k, we have: ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ $$n ≈ k - \large{1\over 24k} - {1\over2}$$ HP50g inverse harmonic function: ﻿ ﻿ ﻿ ﻿ ﻿ « 1. Psi + EXP DUP 24. * INV - .5 - » Both formulas perform well for large H, but this new version is also good for smaller H >>> from mpmath import exp, euler, harmonic >>> f1 = lambda k: k + 6/(1/k+12) - 1 >>> f2 = lambda k: k - 1/(24*k) - .5 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ # new version >>> for n in range(50,100,10): ... ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ h = harmonic(n) ... ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ k = exp(h - euler) ... ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿e1, e2 = n - f1(k), n - f2(k) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ # under-estimated errors ... ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿﻿print '%d\t%.10f\t%.10f\t%g' % (n, e1, e2, abs(e1/e2)) ... 50 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿-0.0000013228 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿0.0000000364 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿36.3546 60 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿-0.0000009261 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿0.0000000212 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿43.7606 70 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿-0.0000006844 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿0.0000000134 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿51.1669 80 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿-0.0000005263 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿0.0000000090 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿58.5736 90 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿-0.0000004172 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿0.0000000063 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿65.9804 RE: Cut the Cards - Gerson W. Barbosa - 08-26-2020 08:23 AM (08-26-2020 03:06 AM)Albert Chan 0.247866177136� 54 Wrote:  HP50g inverse harmonic function: ﻿ ﻿ ﻿ ﻿ ﻿ « 1. Psi + EXP DUP 24. * INV - .5 - » Both formulas perform well for large H, but this new version is also good for smaller H For H less than 2.0701076 the first formula is slightly better: « DUPDUP H H­¹ - SWAP DUP H H­¹2 - / ABS » 50 → 36.4435261708 10 → 6.75851700917 2.0701076 → 1. 1 → 0.290013346857 0.5131 → 2.30966647121E-5 0.5 → 7.2984301462E-3 0.2 → 0.162258684815 0 → 0.247866177136 RE: Cut the Cards - Albert Chan - 08-26-2020 02:13 PM I added a tiny correction to compensate errors for smaller H. HP50g H-1 ﻿: ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ « 1. Psi + EXP DUPDUP SQ 24. * 2.7 + / - .5 - » Code: from mpmath import exp, euler, harmonic, psi f1 = lambda k: k + 6/(1/k+12) - 1 f2 = lambda k: k - 1/(24*k) - .5 f3 = lambda k: k - k/(24*k*k+2.7) - .5    # new formula def test(lst):     print 'n\t err(f1)\t err(f2)\t err(f3)'     for n in lst:         h = harmonic(n)         k = exp(h - euler)         print n,'\t% .10f\t% .10f\t% .10f' % (n-f1(k), n-f2(k), n-f3(k)) >>> test(range(10) + range(10,60,10)) n ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ err(f1) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ err(f2) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ err(f3) 0 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.0031607566 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.0127518672 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0067666262 1 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0003177730 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.0010957186 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0001620997 2 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0002588484 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.0002719596 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000171164 3 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0001712995 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.0001037256 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000035307 4 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0001178786 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.0000497922 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000010529 5 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000852046 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.0000275594 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000003957 6 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000641850 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.0000167992 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000001742 7 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000499820 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.0000109785 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000000859 8 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000399758 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.0000075616 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000000463 9 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000326776 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.0000054263 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000000266 10 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000271983 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.0000040243 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000000162 20 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000076840 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.0000005432 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000000006 30 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000035570 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.0000001651 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000000001 40 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000020419 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.0000000705 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000000000 50 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000013228 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 0.0000000364 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -0.0000000000 If needed, we could apply Newton's method to get better estimate. >>> h = harmonic(10) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ # h = 2.9289682539682538 >>> x = f3(exp(h-euler)) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ # x = 10.000000016205583 >>> x - (harmonic(x) - h) / psi(1, x+1) mpf('10.0') RE: Cut the Cards - Gerson W. Barbosa - 08-26-2020 06:13 PM D’oh, my latest improvement ended up being equivalent to your f2 program, same number of bytes as yours (56.5). That’s an old attempt at exact results, but I wasn’t pleased with it. Perhaps I can make it shorter some day: https://hpmuseum.org/forum/thread-6157-post-55046.html RE: Cut the Cards - Gerson W. Barbosa - 08-27-2020 10:07 PM How about « 1. Psi + EXP DUP 1. + 4. * INV OVER 1. EXP / + INV OVER 24. * + INV - .5 - » ? Again, best for smaller H. A previous attempt ( « 1. Psi + EXP DUP SQ 24. * 1. EXP + INV NEG 1. + * .5 - » ) gave the same results as yours when I replaced 1 EXP ( e ) with 2.7. RE: Cut the Cards - Albert Chan - 08-28-2020 09:26 PM (08-26-2020 02:13 PM)Albert Chan Wrote:  I added a tiny correction to compensate errors for smaller H. HP50g H-1 ﻿: ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ « 1. Psi + EXP DUPDUP SQ 24. * 2.7 + / - .5 - » (08-27-2020 10:07 PM)Gerson W. Barbosa Wrote:  How about « 1. Psi + EXP DUP 1. + 4. * INV OVER 1. EXP / + INV OVER 24. * + INV - .5 - » ? Just to be clear, both corrections are not based on actual continued fraction of e^ψ(x+1/2) My correction (and probably Gerson's version) are based on curve fitting, for small H. I picked 2.7 instead of e to get better accuracy at bigger H, it is simply a trade-off. To build e^ψ(x+1/2) continued fraction with XCas: XCas> k := exp(Psi(x+0.5)) XCas> limit(k/x, x=inf) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 1 XCas> limit(1/(k -x)/x, x=inf) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 24 XCas> limit(1/(1/(k -x) -24x)/x , x=inf) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 10/37 XCas> limit(1/(1/(1/(k -x) -24x) -10/37*x)/x , x=inf) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 0 ??? The last one is likely a XCas bug, Mathematica get the coefficient: 689976/74381 $$\qquad\qquad\exp( \psi(x+1/2)) = x +\frac{1}{24 \cdot x +\Large\frac{1}{\frac{10}{37} \cdot x +\frac{1}{\frac{689976}{74381} \cdot x \; +\; \cdots}}}$$ To confirm, lets try H-1(H(20)): XCas> f(x) := 1/(24*x + 1/(10/37*x + 1/(689976/74381*x))) XCas> h := float(harmonic(20)) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 3.59773965714 XCas> k := exp(h - euler_gamma)﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 20.5020317757 XCas> fsolve(x+f(x) = k, x=k) - 0.5 ﻿ ﻿ ﻿ → 20.0 If f(x) is small, x ≈ k. We can also iterate for x, like this: XCas> x0 := k XCas> x0 := k - f(x0)﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 20.5000002012 XCas> x0 := k - f(x0)﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 20.5 XCas> n := x0 - 0.5 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 20.0 RE: Cut the Cards - Gerson W. Barbosa - 08-28-2020 11:39 PM For the inverse digamma function I’ve been using a continued fraction with increasing number of terms. I am sure about only the first half of them, though. It’s difficult to be sure about the other half when one is limited to only 12 digits. $$\rm{e}^{{x}}+\frac{1}{2+\frac{1}{6\rm{e}^{{x}}-\frac{1}{2+\frac{1}{\rm{e}^{{x}-1}-\frac{1}{2+\frac{1}{\rm{e}^{{x}}+\frac{1}{2}}}}}}}$$ RE: Cut the Cards - Albert Chan - 08-29-2020 04:02 PM (08-27-2020 10:07 PM)Gerson W. Barbosa Wrote:  How about « 1. Psi + EXP DUP 1. + 4. * INV OVER 1. EXP / + INV OVER 24. * + INV - .5 - » ? FYI, above H-1 is equvialent to Ψ-1 (previous post) Since H(n) = Ψ(n+1) - Ψ(1), we flip the sign of -0.5 to +0.5 XCas> invPsi1(k) := k - 1/(24k + 1/(k/e + 1/(4k+4))) + 1/2 $$k \rightarrow k-\frac{1}{24\cdot k+ \large\frac{1}{\frac{k}{e}+\frac{1}{4\cdot k+4}}}+\frac{1}{2}$$ XCas> invPsi2(k) := k + 1/(2 + 1/(6k - 1/(2 + 1/(k/e - 1/(2 + 1/(k+1/2)))))) $$k \rightarrow k+\frac{1}{2+\frac{1}{6\cdot k-\frac{1}{2+\frac{1}{\frac{k}{e}-\frac{1}{2+\frac{1}{k+\frac{1}{2}}}}}}}$$ XCas> simplify(invPsi1(k) - invPsi2(k)) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 0 Let's try a few x, to see if it round-trip. XCas> invPsi(x) := invPsi1(exp^float(x)) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // Note: k = e^x XCas> map(range(1,6), x->[x, Psi(invPsi(x))]) $$\qquad\left(\begin{array}{cc} 1 & 0.999999666188 \\ 2 & 2.00000000607 \\ 3 & 3.00000000018 \\ 4 & 4.0 \\ 5 & 5.0 \end{array}\right)$$ RE: Cut the Cards - Albert Chan - 06-23-2021 12:08 AM An integral, from a very old (2006) thread: So your HP can INTEGRATE ... I revived the old thread recently, for an integral that HP71B failed, see HP71B IBOUND fooled Code:             / Inf             |     FRAC(x)      I7 =   |     -------  .dx             |       x²            / 1 I put the solution here because integral is closely related to harmonic number I7 = ∫((x-1)/x^2, x=1..2) + ∫((x-2)/x^2, x=2..3) + ∫((x-3)/x^2, x=3..4) + ... ﻿ ﻿ ﻿ ﻿ =﻿ (ln(2) - ln(1) - 1/2)﻿ ﻿ ﻿ ﻿ + (ln(3) - ln(2) - 1/3) ﻿ ﻿ + (ln(4) - ln(3) - 1/4) ﻿ ﻿ ﻿ + ... $$I_7 - 1 = \displaystyle \lim_{n \to ∞} (\,\ln(n) - H_n\,) = - \gamma$$ $$I_7 = 1 - \gamma ≈ 0.422784$$