sum of large sequence numbers

[ericyuan]

Hi,

I have a puzzle that needs to be solved.

sum(k/(t+1) + k*2/(t+2) + k*3/(t+3) + k*4/(t+4) ... k*n/(t+n))
the n can be very large number, such as 10million
k and t are constants.

Is there any better approach to calculate the sum?
I can solve it using a loop in a computer, but it's not efficient.

Thank you so much for reading.

Cheers,
Eric

 

[ericyuan]

sorry, it should be
sum(k/(t+1) + k/(t+2) + k/(t+3) + k/(t+4) ... k/(t+n))

 

[Dr.Peterson]

sorry, it should be
sum(k/(t+1) + k/(t+2) + k/(t+3) + k/(t+4) ... k/(t+n))

I think the word "sum" there is irrelevant, and you just mean the sum

k/(t+1) + k/(t+2) + k/(t+3) + k/(t+4) + ... + k/(t+n) = k(1/(t+1) + 1/(t+2) + 1/(t+3) + 1/(t+4) + ... + 1/(t+n))​

That is the difference of two harmonic numbers: https://en.wikipedia.org/wiki/Harmonic_number

There is no simple formula, but there are ways to approximate these numbers for large n.

 

[ericyuan]

Thank you so much, Dr Peterson for your prompt response.

I forgot one constant
1/(t+1*k) + 1/(t+2*k) + 1/(t+3*k) + 1/(t+4*k) ... +1/(t+n*k))

for example
k = 50
t = 1000
n = 5000
1/ (1000+50) + 1/(1000+50+50) + 1(1000+50+50+50) .... + 1/(1000 + 5000 x 50)

What's the most efficient approach to calculate instead of a loop?
It's used to calculate an investment return for a customer, so it needs to be precise.
Any suggestion would be much appreciated.

 

[JeffM]

Thank you so much, Dr Peterson for your prompt response.

I forgot one constant
1/(t+1*k) + 1/(t+2*k) + 1/(t+3*k) + 1/(t+4*k) ... +1/(t+n*k))

for example
k = 50
t = 1000
n = 5000
1/ (1000+50) + 1/(1000+50+50) + 1(1000+50+50+50) .... + 1/(1000 + 5000 x 50)

What's the most efficient approach to calculate instead of a loop?
It's used to calculate an investment return for a customer, so it needs to be precise.
Any suggestion would be much appreciated.

First, if you are using a computer, as is implied by the mention of a loop, programming a loop is scarcely arduous and is exactly what computers are good at.

Second, you have not explained the nature of the underlying investment that generated this series. Perhaps there is an alternative formula that involves far fewer calculations. We cannot answer that question without more information than you have given. You have provided an answer to an unknown question. Not much that we can do without the question itself.

Third, I greatly doubt that you need an answer that is exact beyond $0.005 if the investment is in dollars. No investor today expects to be paid in bits.

 

[Cubist]

I forgot one constant
1/(t+1*k) + 1/(t+2*k) + 1/(t+3*k) + 1/(t+4*k) ... +1/(t+n*k))

The following website suggests a solution:- WolframAlpha(click). I've never heard of the "digamma" function before. It's related to the harmonic numbers, which makes sense given @Dr.Peterson 's reply.

C++ Boost libraries have an implementation of the digamma function. If you're not using boost then the following pseudo code works...

define digamma(double z) {
  # this returns an approximation, but hopefully it should be accurate enough for your needs
  # see https://en.wikipedia.org/wiki/Digamma_function#Asymptotic_expansion
  if (z<20) {
    return digamma(z+1) - 1.0/z; # recursion is unlikely to occur if the numbers used are similar to your example
  } else {
    double x=1.0/(z*z);
    return ln(z) - 1.0/(2*z) + x*(-1.0/12 + x*(1.0/120 + x*(-1.0/252 + x*(1.0/240 + x*(-1.0/132 + x*(691.0/32760 +x*(-1.0/12 + x*(3617.0/8160 + x*(-43867.0/14364 + x*174611.0/6600)))))))));
  }
}

k=50
t=1000
n=5000

# Calculate with a loop
s=0
for (j=n;j>=1;--j) s+=1/(t+j*k); # Using reverse order will add the smallest numbers first, perhaps giving a more accurate result (less underflow)
print s;

# Calculate without a loop
print (digamma(n+t/k+1)-digamma(t/k+1))/k

# Results...
0.1100152163746085104465404192569
0.1100152163746085104465404192123
# First difference here       ^

It's used to calculate an investment return for a customer, so it needs to be precise.

I completely agree with @JeffM that you're unlikely to require a lot of precision. But I'll leave it to you to judge if the above code is good enough (obviously you need to test it with a range of inputs greater than you'd expect in practice).