▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
Write a program that adds the sum of integers raised to a given power (integer or noninteger).
The program should take as input the number of integers and the power used to raise each integer.
The summation starts with the integer 1 and ends with the specified integer.
The challenge is to write a fast routine to add integers raised to integer or noninteger powers.
Your program MAY calculate the summation as a close approximation to the actual sum.
The following website discusses the summation of integers taken to various powers:
http://planetmath.org/encyclopedia/Summing.html
If S(n,p) is the sum of the first n integers, each integer rased to power p, then:
S(n,1) = 1 + 2 + 3 + 4 + ... + n = n(n + 1)/2 which is Gauss's formula
S(n,2) = 1^2 + 2^2 + 3^2 + 4^2 + ... + n^2 = n(n + 1)(2n + 1)/6
S(n,3) = 1^3 + 2^3 + 3^3 + 4^3 + ... + n^3 = n^2 (n + 1)^2 /4
S(n,4) = 1^4 + 2^4 + 3^4 + 4^4 + ... + n^4 = n(n + 1)(2n + 1)(3n^2 + 3n  1)/30
S(n,5) = 1^5 + 2^5 + 3^5 + 4^5 + ... + n^5 = n^2 (n + 1)^2 (2n^2 + 2n  1)/12
and so on ... the link metion above shows the summations for S(n,1) through S(n,10). I was unable to locate formulas for noninteger summations.
The input for the routine you write is:
1. The value for n
2. The value for the powes p.
You routine should output the value for summation of integers S(n,p).
You can code in RPN, RPL, or BASIC. Use your favorite machine.
May the best code win!
Namir
Edited: 28 Oct 2006, 4:06 p.m. after one or more responses were posted
▼
Posts: 562
Threads: 27
Joined: Oct 2006
greetings, can you give an example of the inputs and outputs?? Looks interesting, but why not just [n!] [y^x]? Is there a specific platform? 15C, 42S, 48G?
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
I have updated my original post to include some of th formulas for a few summations.
Namir
Posts: 562
Threads: 27
Joined: Oct 2006
Here is a slow one I wrote on the 48G for testing the fast one:
<< { } 1 ROT
FOR I I + NEXT
SWAP ^ SUMLIST
>>
Checksum #17Eh
Size (bytes): 47
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
For testing purpose, this is extremely fast on the HP50G in exact mode:
<<
> n p 'SIGMA(k=1,n,k^p)'
>>
Less than 3.5 seconds for n = 1000000000 and p = 10. Awkwardly this is also the time to compute the sum of 10th powers of the first 1000 numbers:
91409924241424243424241924242500
This is the same result obtained by Jacques Bernoulli himself around 1690, in about seven minutes and thirty seconds, using '(n*(n+1)*(2*n+1)*(n^2+n1)*(3*n^6+9*n^5+2*n^411*n^3+3*n^2+10*n5))/66'
http://www.mathpages.com/home/kmath279.htm
▼
Posts: 1,041
Threads: 15
Joined: Jan 2005
Umm, actually, I think that
91409924241424243424241924242500
is the result for n=1000 and p=10. For n=1000000000 (10 ^{9}) and p=10, I get:
90909091409090909924242424242424241424242424242424243424242424242424241924242424242424242500000000
With either set of inputs, your program runs in less than 3.5 seconds on my 49g+.
I suspect that, barring some relatively littleknown mathematical "trick", the obvious builtin tools will be hard to beat on an RPL model.
Regards, James
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
Hello James,
Quote:
Umm, actually, I think that
91409924241424243424241924242500
is the result for n=1000 and p=10.
That's what I stated. I had both results, but I chose to present the one calculated by Bernoulli.
Quote:
your program runs in less than 3.5 seconds on my 49g+.
This was not meant to be the program Namir is looking for. I am not working on this but I am curious to know what he's got under his sleeves :)
Regards,
Gerson.
▼
Posts: 1,041
Threads: 15
Joined: Jan 2005
Hello Gerson,
Quote:
That's what I stated. I had both results, but I chose to present the one calculated by Bernoulli.
Okay, now I see. Sorry, I misunderstood your post.
Quote:
This was not meant to be the program Namir is looking for.
But it is pretty fast, even with "real numbers" instead of "exact
integers" for the integers and fractional powers instead of
integer powers. For example, with n=1000. and p=9.9, it
returns (in approximate mode) 4.6231487245E31 in about 11.5
seconds on my 49g+.
Quote:
I am not working on this but I am curious to know what he's got under his sleeves :)
I've considered methods such as loops instead of the builtin
summation function, but nothing that would be quite as fast occurs
to me. I'm curious as to what he's got up his sleeve too;
presumably something other than a fairly obvious "bruteforce"
method.
Regards, James
Posts: 2,247
Threads: 200
Joined: Jun 2005
Good work guys!! It’s nice to see some cool RPL code. Here is my solution.
I started looking at summations of integers raised to realvalue powers. I did a basic study using Excel and
linearized regression. I noticed that the loglog of S(n,p) vs n fits gave good straight lines. That is,
ln[S(n,p)] = a + b ln(n) had a high value for the coefficient of determination. This model worked for integer
and real values of p. I realized that there I can apply some approximations to the various series.
Armed with the results of the power fit, I focused on the approximation approach. Using the analytical equations
for the series S(n,p) it is easy to see that the term with the highest power of n is the most relevant. It is
these terms that appeared in the power curve fir that I performed. My first step is to deduce the first approximation
for S(n,p), which is:
S(n,p) = n^(p+1) / (p+1)
This approximation works well as the values for n and p increase.
To improve the above approximation, I decided to look at including the second most relevant term.
For a simple series, where p = 1, S(n,1) is:
n^2 / + n / 2
The second term is n / 2.
For the case of p = 2, S(n, 2) is:
n^3/3 + n^2/2 + n/6
The second term is n^2 / 2.
For the case of p = 3, S(n, 3) is:
n^4/4 + n^3/2 + n^2/4
The second term is n^3 / 2.
For the case of p = 4, S(n, 4) is:
n^5/5 + n^4/2 + n^3/3 – n/30
The second term is n^4 / 2.
For the case of p = 5, S(n, 5) is:
n^6/6 + n^5/2 + 5n^4/12 – n^2/12
T second term is n^5 / 2.
For the case of p = 6, S(n, 6) is:
n^7/7 + n^6/2 + n^5/2 – n^3/6 + n/42
The second term is n^6 / 2.
By induction, I write the general formula for the second most relevant term for S(n.p) as:
n^p / 2
Thus, a refined approximation for S(n,p) is:
S(n,p) = n^p/(p+1) + n^p / 2
Which I rewrite as:
S(n,p) = n^p (n/(p+1) + 1/2)
The summation approximation formula works for p = 1,2,3,4, and on. For p = 1, the above equation gives an exact
result since it matches the Gauss formula. The approximation also does well for positive real values of the power p.
The summation approximation formula gives good approximations at a fraction of the time needed in using loops
(especially for positive noninteger values of p). Moreover the calculation time is independent of the value of n.
The duration of loops used to calculate the summations is directly proportional to the value of n.
And finally, here is a sample results that compare the exact summations withe the approximated values. The error is taken
per billion instead of a percent.
N S(n,1.5) Approx S(n,1.5) Error (ppb)

100 40501.22452 40500 30234.0
250 397263.082 397261.1311 4910.9
500 2241660.917 2241658.147 1235.5
1000 12664925.96 12664922.03 310.1
5000 707283566.7 707283557.9 12.5
10000 4000500012 4000500000 3.1
N S(n,2) Approx S(n,2) Error (ppb)

100 338350 338333.3333 49258.7
250 5239625 5239583.333 7952.2
500 41791750 41791666.67 1994.0
1000 333833500 333833333.3 499.3
5000 41679167500 41679166667 20.0
10000 3.33383E+11 3.33383E+11 5.0
N S(n,2.5) Approx S(n,2.5) Error (ppb)

100 2907351.199 2907142.857 71660.3
250 71081484.32 71080660.8 11585.6
500 801393120.5 801390791.2 2906.5
1000 9050897005 9050890417 727.9
5000 2.52627E+12 2.52627E+12 29.2
10000 2.85764E+13 2.85764E+13 7.3
N S(n,3) Approx S(n,3) Error (ppb)

100 25502500 25500000 98029.6
250 984390625 984375000 15872.8
500 15687562500 15687500000 3984.0
1000 2.505E+11 2.505E+11 998.0
5000 1.56313E+14 1.56313E+14 40.0
10000 2.5005E+15 2.5005E+15 10.0
N S(n,3.5) Approx S(n,3.5) Error (ppb)

100 227251388.7 227222222.2 128344.6
250 13848978155 13848689927 20812.2
500 3.11964E+11 3.11963E+11 5226.5
1000 7.0431E+12 7.0431E+12 1309.6
5000 9.82535E+15 9.82535E+15 52.5
10000 2.22272E+17 2.22272E+17 13.1
N S(n,6) Approx S(n,6) Error (ppb)

100 1.47907E+13 1.47857E+13 338038.7
250 8.84187E+15 8.84138E+15 55223.5
500 1.1239E+18 1.12388E+18 13902.5
1000 1.43358E+20 1.43357E+20 3487.8
5000 1.11685E+25 1.11685E+25 139.9
10000 1.42907E+27 1.42907E+27 35.0
Edited: 29 Oct 2006, 5:33 p.m. after one or more responses were posted
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
Quote:
I noticed that the loglog of S(n,p) vs n fits gave good straight lines.
Excellent insight!
The HP34, my slowest calculator, gives the results below in about two seconds.
1000, 10 => 9.140909E31
1000000000, 10 => 9.090909E97
These match the HP50G results at this precision. Although exact results are always beautiful, for real world applications nice approximations might be preferable. Really excellent work!
Gerson.

01 LBL A
02 x<>y
03 ENTER
04 ENTER
05 R^
06 y^x
07 x<>y
08 LST x
09 1
10 +
11 /
12 .5
14 +
15 *
16 RTN
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
Gerson,
Thank you for your kind complements and for your clever listing that does not use any memory register.
The approximation route gives the user the ability to quickly estimate the sum of integers raised to integer or real powers, if such approximation works for his.her application.
Namir
Posts: 2,761
Threads: 100
Joined: Jul 2005
Hello Namir,
Quote:
S(n,p) = n^p (n/(p+1) + 1/2)
After a little work, I made a slight modification in your approximation formula. It's not exact for p=1 anymore but it appears to be somewhat more accurate for all other powers:
S(n,p) = n^p (n/(p+1) + (6n + p)/(12n))
So far I have tested it only with integer powers up to 5, and n up to 10000.
Regards,
Gerson.
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
Gerson,
Good work!!!! The new terms makes the approximation for p >= 2 even better.
Many Thanks!
I guess we can call it the BarbosaShammas Summation Series Approximation of BSSSA for short (<grin>).
Namir
Edited: 29 Oct 2006, 5:59 p.m.
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
Quote:
I guess we can call it the BarbosaShammas Summation Series Approximation... (<grin>).
Better call it Shammas^{10}Barbosa Summation Series Approximation since you did the hard work. All I did was noticing the error for p=0 was always exactly 1/2. So I guessed a correction factor could be introduced next to that constant to minimize the error. I computed the correction factors for p=2, 3, 4 and 5 and found values very close to (1 + 1/(k*n)), k = 3, 2, 1.5 and 1.2, respectively, that is, k = 6/p. The error for p=0 still remains, but the approximation is good enough for practical applications (if any :).
Gerson.
Edited: 29 Oct 2006, 7:23 p.m.
Posts: 2,761
Threads: 100
Joined: Jul 2005
Some tables:

n p ~ =
100 0 100.50 100.00
250 0 250.50 250.00
500 0 500.50 500.00
1,000 0 1,000.50 1,000.00
5,000 0 5,000.50 5,000.00
10,000 0 10,000.50 10,000.00

n p ~ =
100 1 5,050.08 5,050.00
250 1 31,375.08 31,375.00
500 1 125,250.08 125,250.00
1,000 1 500,500.08 500,500.00
5,000 1 12,502,500.08 12,502,500.00
10,000 1 50,005,000.08 50,005,000.00

n p ~ =
100 2 338,350.00 338,350.00
250 2 5,239,625.00 5,239,625.00
500 2 41,791,750.00 41,791,750.00
1,000 2 333,833,500.00 333,833,500.00
5,000 2 41,679,167,500.00 41,679,167,500.00
10,000 2 333,383,335,000.00 333,383,335,000.00

n p ~ =
100 2.5 2,907,351.19 2,907,351.20
250 2.5 71,081,484.31 71,081,484.32
500 2.5 801,393,120.46 801,393,120.47
1,000 2.5 9,050,897,005.43 9,050,897,005.42
5,000 2.5 2.52626531851E+12 2.52626531852E+12
10,000 2.5 2.85764287798E+13 2.85764287815E+13

n p ~ =
100 6 1.47907142857E+13 1.47907141190E+13
250 6 8.84186662946E+15 8.84186662689E+15
500 6 1.12389955357E+18 1.12389955354E+18
1,000 6 1.43357642857E+20 1.43357642858E+20
5,000 6 1.11685283482E+25 1.11685283482E+25
10,000 6 1.42907147857E+27 1.42907147851E+27

And
S(1000,10):
Exact: 91409924241424243424241924242500.00
Approximate: 91409924242424242424242424242424.24
Edited: 29 Oct 2006, 7:13 p.m.
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
Cool tables!!!
You additional terms make the approximation give exact values for p = 2, and reduces the error for higher p values.
Namir
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
However it's become harder to write an RPN program using only the stack :) (I prefer this method because programs using only the stack are faster).
Using the equivalent formula below it's been possible to save one step on the HP12C (another step may be saved using the 12* or 12/ instructions). By the way, my not so efficient program runs in 2.5 seconds on the 12C.
'n^p*(n/(p+1)+p/(12*n)+1/2)'
Gerson.
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
Here is an HP41C listing. No rocket science here:
01 LBL "SSA
02 LBL A
03 "N?
04 PROMPT
05 STO 00
06 "P?
07 PROMPT
08 STO 01
09 2
10 X<>Y
11 X<Y?
12 0
13 RCL 00
14 6
15 *
16 *
17 12
18 /
19 RCL 00
20 /
21 RCL 00
22 RCL 01
23 1
24 +
25 /
26 +
27 RCL 00
28 RCL 01
29 Y^X
30 *
31 "S=
32 ARCL X
33 PROMPT
34 RTN
I use LBL A because it makes it easier to rerun the program without reexecuting the command XEQ SSA and without assigning SSA to a key.
Namir
Posts: 230
Threads: 11
Joined: Jan 1970
"By induction" makes my teeth grin, but actually you are right the term is 1/2, which is (IMHO) beautiful.
Can anyone provide the next term ? I think I recall Bernoulli numbers are used.
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
I did see a gneral equation in a math handbook that specifies the first few terms. The third term did involve the Bernoulli function. The second term was n^p/2.
Gerson was able to find the third term and combine it with the second term. The result is still a compact approximation for p >= 2.
Namir
Edited: 30 Oct 2006, 7:45 a.m.
▼
Posts: 230
Threads: 11
Joined: Jan 1970
I just dug out an old note with the following result :
Sum(i=1 to n,i^p) = sum(i=1 to p+1,C(i,p)*S(p+1i)*n^i)
the second term is the wanted polynom in powers of n.
In this formula, C() is the standard number of combinations with the added twist that C(p,p+1)=1
and S() are constant (!!!!) values in the following table :
j S(J)
0 1
1 1/2
2 1/12
3 0
4 1/120
5 0
6 1/252
7 0
8 1/240
9 0
10 1/132
11 0
12 105/4978
13 0
14 1/12
I'm not sure this is funny
▼
Posts: 230
Threads: 11
Joined: Jan 1970
Oooops, I meant C(p,p+1) = 1/(p+1)
