Math Challenge



#23

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


#24

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?


#25

I have updated my original post to include some of th formulas for a few summations.

Namir

#26

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


#27

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'

http://www.mathpages.com/home/kmath279.htm


#28

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


#29

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.


#30

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

#31

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


#32

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


#33

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

#34

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.


#35

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 (<grin>).


Namir

Edited: 29 Oct 2006, 5:59 p.m.


#36

Quote:
I guess we can call it the Barbosa-Shammas Summation Series Approximation... (<grin>).

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.

#37

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.


#38

Cool tables!!!

You additional terms make the approximation give exact values for p = 2, and reduces the error for higher p values.

Namir


#39

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.


#40

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<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 re-executing the command XEQ SSA and without assigning SSA to a key.

Namir

#41

"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.


#42

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.


#43

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

#44

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,867 12-13-2013, 02:23 AM
Last Post: Didier Lachieze
  OT: a math competition site Pier Aiello 0 512 09-16-2013, 06:03 AM
Last Post: Pier Aiello
  Simple Math Question Namir 2 707 08-09-2013, 06:13 PM
Last Post: Eddie W. Shore
  Cool math clock Bruce Bergman 28 3,421 04-10-2013, 03:13 AM
Last Post: Siegfried (Austria)
  Math Challenge I could not solve Meindert Kuipers 22 2,801 01-05-2013, 04:43 PM
Last Post: Thomas Klemm
  Math Question Namir 0 451 11-06-2012, 07:43 AM
Last Post: Namir
  Survey for Special Math Problem Namir 7 1,300 06-03-2012, 09:46 PM
Last Post: Namir
  math question Don Shepherd 22 2,535 04-25-2012, 11:38 PM
Last Post: Don Shepherd
  Math Help!! Namir 10 1,556 04-16-2012, 11:32 PM
Last Post: Allen
  go41cx overlay Math module Alexander Oestert 12 1,848 02-04-2012, 05:09 AM
Last Post: Alexander Oestert

Forum Jump: