Math Challenge « Next Oldest | Next Newest »

 ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-28-2006, 11:07 AM 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 ▼ Allen Senior Member Posts: 562 Threads: 27 Joined: Oct 2006 10-28-2006, 02:29 PM 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? ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-28-2006, 03:46 PM I have updated my original post to include some of th formulas for a few summations. Namir Allen Senior Member Posts: 562 Threads: 27 Joined: Oct 2006 10-28-2006, 08:25 PM 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 ``` ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-28-2006, 09:06 PM 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'` ▼ James M. Prange (Michigan) Posting Freak Posts: 1,041 Threads: 15 Joined: Jan 2005 10-28-2006, 10:52 PM 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 ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-28-2006, 11:13 PM 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. ▼ James M. Prange (Michigan) Posting Freak Posts: 1,041 Threads: 15 Joined: Jan 2005 10-28-2006, 11:51 PM 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 Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-29-2006, 10:19 AM ```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 ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-29-2006, 11:35 AM 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 ``` ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-29-2006, 01:20 PM 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 Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-29-2006, 05:08 PM 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. ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-29-2006, 05:36 PM 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. ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-29-2006, 06:25 PM 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. Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-29-2006, 06:12 PM 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. ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-29-2006, 10:13 PM Cool tables!!! You additional terms make the approximation give exact values for p = 2, and reduces the error for higher p values. Namir ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-29-2006, 10:46 PM 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. ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-30-2006, 07:42 AM 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= 2. Namir Edited: 30 Oct 2006, 7:45 a.m. ▼ GE Member Posts: 230 Threads: 11 Joined: Jan 1970 11-15-2006, 04:34 AM 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 ▼ GE Member Posts: 230 Threads: 11 Joined: Jan 1970 11-15-2006, 04:36 AM Oooops, I meant C(p,p+1) = 1/(p+1)

 Possibly Related Threads... Thread Author Replies Views Last Post Need help understanding math.... cyrille de Brébisson 9 1,721 12-13-2013, 02:23 AM Last Post: Didier Lachieze OT: a math competition site Pier Aiello 0 487 09-16-2013, 06:03 AM Last Post: Pier Aiello Simple Math Question Namir 2 662 08-09-2013, 06:13 PM Last Post: Eddie W. Shore Cool math clock Bruce Bergman 28 3,133 04-10-2013, 03:13 AM Last Post: Siegfried (Austria) Math Challenge I could not solve Meindert Kuipers 22 2,633 01-05-2013, 04:43 PM Last Post: Thomas Klemm Math Question Namir 0 417 11-06-2012, 07:43 AM Last Post: Namir Survey for Special Math Problem Namir 7 1,181 06-03-2012, 09:46 PM Last Post: Namir math question Don Shepherd 22 2,399 04-25-2012, 11:38 PM Last Post: Don Shepherd Math Help!! Namir 10 1,443 04-16-2012, 11:32 PM Last Post: Allen go41cx overlay Math module Alexander Oestert 12 1,723 02-04-2012, 05:09 AM Last Post: Alexander Oestert

Forum Jump: