The Museum of HP Calculators

HP Forum Archive 16

 Math ChallengeMessage #1 Posted by Namir on 28 Oct 2006, 11:07 a.m. Write a program that adds the sum of integers raised to a given power (integer or non-integer). 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 non-integer 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 non-integer 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

 Re: Math ChallengeMessage #2 Posted by Allen on 28 Oct 2006, 2:29 p.m.,in response to message #1 by Namir 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?

 Re: Math ChallengeMessage #3 Posted by Namir on 28 Oct 2006, 3:46 p.m.,in response to message #2 by Allen I have updated my original post to include some of th formulas for a few summations. Namir

 Re: Math ChallengeMessage #4 Posted by Allen on 28 Oct 2006, 8:25 p.m.,in response to message #1 by Namir 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 ```

 Re: Math ChallengeMessage #5 Posted by Gerson W. Barbosa on 28 Oct 2006, 9:06 p.m.,in response to message #4 by Allen For testing purpose, this is extremely fast on the HP-50G 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+n-1)*(3*n^6+9*n^5+2*n^4-11*n^3+3*n^2+10*n-5))/66'`

 Re: Math ChallengeMessage #6 Posted by James M. Prange (Michigan) on 28 Oct 2006, 10:52 p.m.,in response to message #5 by Gerson W. Barbosa Umm, actually, I think that ```91409924241424243424241924242500 ``` is the result for n=1000 and p=10. For n=1000000000 (109) 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 little-known mathematical "trick", the obvious built-in tools will be hard to beat on an RPL model. Regards,James

 Re: Math ChallengeMessage #7 Posted by Gerson W. Barbosa on 28 Oct 2006, 11:13 p.m.,in response to message #6 by James M. Prange (Michigan) 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.

 Re: Math ChallengeMessage #8 Posted by James M. Prange (Michigan) on 28 Oct 2006, 11:51 p.m.,in response to message #7 by Gerson W. Barbosa 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 built-in 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 "brute-force" method. Regards,James

 Re: Math ChallengeMessage #9 Posted by Namir on 29 Oct 2006, 10:19 a.m.,in response to message #1 by Namir ```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 real-value powers. I did a basic study using Excel and linearized regression. I noticed that the log-log 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 non-integer 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

 Re: Math ChallengeMessage #10 Posted by Gerson W. Barbosa on 29 Oct 2006, 11:35 a.m.,in response to message #9 by Namir Quote: I noticed that the log-log of S(n,p) vs n fits gave good straight lines. Excellent insight! The HP-34, my slowest calculator, gives the results below in about two seconds. ``` 1000, 10 => 9.140909E31 1000000000, 10 => 9.090909E97 ``` These match the HP-50G 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 ```

 Re: Math ChallengeMessage #11 Posted by Namir on 29 Oct 2006, 1:20 p.m.,in response to message #10 by Gerson W. Barbosa 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

 Re: Math ChallengeMessage #12 Posted by Gerson W. Barbosa on 29 Oct 2006, 5:08 p.m.,in response to message #9 by Namir 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.

 Re: Math ChallengeMessage #13 Posted by Namir on 29 Oct 2006, 5:36 p.m.,in response to message #12 by Gerson W. Barbosa Gerson, Good work!!!! The new terms makes the approximation for p >= 2 even better. Many Thanks! I guess we can call it the Barbosa-Shammas Summation Series Approximation of BSSSA for short (). Namir Edited: 29 Oct 2006, 5:59 p.m.

 Re: Math ChallengeMessage #14 Posted by Gerson W. Barbosa on 29 Oct 2006, 6:25 p.m.,in response to message #13 by Namir Quote: I guess we can call it the Barbosa-Shammas Summation Series Approximation... (). Better call it Shammas10-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.

 Re: Math ChallengeMessage #15 Posted by Gerson W. Barbosa on 29 Oct 2006, 6:12 p.m.,in response to message #12 by Gerson W. Barbosa 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.

 Re: Math ChallengeMessage #16 Posted by Namir on 29 Oct 2006, 10:13 p.m.,in response to message #15 by Gerson W. Barbosa Cool tables!!! You additional terms make the approximation give exact values for p = 2, and reduces the error for higher p values. Namir

 Re: Math ChallengeMessage #17 Posted by Gerson W. Barbosa on 29 Oct 2006, 10:46 p.m.,in response to message #16 by Namir 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 HP-12C (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.

 Re: Math ChallengeMessage #18 Posted by Namir on 30 Oct 2006, 7:42 a.m.,in response to message #17 by Gerson W. Barbosa Here is an HP-41C 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

 Re: Math ChallengeMessage #19 Posted by GE on 30 Oct 2006, 4:48 a.m.,in response to message #9 by Namir "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.

 Re: Math ChallengeMessage #20 Posted by Namir on 30 Oct 2006, 7:36 a.m.,in response to message #19 by GE 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.

 Re: Math ChallengeMessage #21 Posted by GE on 15 Nov 2006, 4:34 a.m.,in response to message #20 by Namir 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+1-i)*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

 Re: Math ChallengeMessage #22 Posted by GE on 15 Nov 2006, 4:36 a.m.,in response to message #21 by GE Oooops, I meant C(p,p+1) = 1/(p+1)

Go back to the main exhibit hall