power sum evaluations

**danbaron** · 30-08-2011, 08:43

This program calculates a power sum twice.

By a power sum, I mean, for given values of x and n,

power sum = x^0 + x^1 + x^2 + x^3 + .. + x^n.

As it is, in the program I have set x=2, and n=100.

The program calculates the given power sum twice, only the first time it calls the function, "slowevaluation()", and the second time it calls the function, "fastevaluation()".

(Actually, both of the functions are very slow, but one is much faster than the other.)

Both evaluations give the same sum.

The question is, why is the fast function so much faster than the slow function (or vice versa)?

' code ----------------------------------------------------------------------------------------------
Uses "console"
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%neg = -1
%pos = +1

%intsize = 500

Type genint
digit(%intsize) As Byte
sign As Integer
error As Byte
End Type

'--------------------------------------------------------------

Function TBMain()
Local s As String
Local t1, t2, tt As Integer
Local x, ps As genint
Local n As DWord

'Set value of x (string s).
s = "2"

'Set number of terms, n.
n = 20

genintsetvalue(x,s)

PrintL "power sum = x^0 + x^1 + x^2 + .. + x^n"
PrintL
PrintL "x = ", Format$(s)
PrintL "n = ", n
PrintL
PrintL "-------------------------------------------------------------"
PrintL

PrintL "Calculating power sum by slow method."
t1 = Timer
slowevaluation(x,n, ps)
s = genintgetstringrep(ps)
t2 = Timer
tt = t2 - t1
PrintL
PrintL "sum = ", s
PrintL
PrintL "elapsed time = ", tt, " seconds"
PrintL
PrintL "-------------------------------------------------------------"
PrintL

PrintL "Calculating power sum by fast method."
t1 = Timer
fastevaluation(x,n, ps)
s = genintgetstringrep(ps)
t2 = Timer
tt = t2 - t1
PrintL
PrintL "sum = ", s
PrintL
PrintL "elapsed time = ", tt, " seconds"
PrintL
PrintL "-------------------------------------------------------------"
PrintL

PrintL "Done. Press a key to quit."

WaitKey
End Function

'--------------------------------------------------------------

Function slowevaluation(ByRef x As genint, n As DWord, psum As genint)
Local xp As genint
Local s As String
Local i As Integer

s = "0"
genintsetvalue(psum,s)
If n=0 Then Return
For i = 0 To n
genintpower(x,i,xp)
genintadd(psum,xp,psum)
Next
End Function

'--------------------------------------------------------------

Function fastevaluation(ByRef x As genint, n As DWord, psum As genint)
Local temp, xp As genint
Local s As String
Local i As Integer

s = "0"
genintsetvalue(psum,s)
If n=0 Then Return
s = "1"
genintsetvalue(psum,s)
genintsetvalue(xp,s)
For i = 1 To n
genintmult(xp,x,temp)
genintsetequal(xp,temp)
genintadd(psum,xp,psum)
Next
End Function

'--------------------------------------------------------------

Function genintpower(ByRef x As genint, p As DWord, ByRef xp As genint)
Local product, temp As genint
Local s As String
Local i As Integer

s = "1"
genintsetvalue(xp,s)
If p=0 Then Return
For i = 1 To p
genintmult(xp,x,temp)
genintsetequal(xp,temp)
Next
End Function

'--------------------------------------------------------------

Function genintfactorial(n As Integer) As String
Local num1, num2, num3 As genint
Local s As String
Local i As Integer

s = "1"
genintsetvalue(num1,s)
For i = n To 2 Step -1
PrintL i
s = Format$(i)
genintsetvalue(num2,s)
genintmult(num1,num2,num3)
genintsetequal(num1,num3)
Next
Return genintgetstringrep(num3)
End Function

'--------------------------------------------------------------

Function genintsetequal(gn2 As genint, gn1 As genint)
'Set gn2 = gn1.
Local i As DWord
For i = 1 To %intsize
gn2.digit(i) = gn1.digit(i)
Next
End Function

'--------------------------------------------------------------

Function genintsetvalue(gn As genint, vs As String)
Local i, l As DWord
genintzero(gn)
l = Len(vs)
For i = 1 To l
gn.digit(l - i + 1) = Val(Mid$(vs, i, 1))
Next
End Function

'--------------------------------------------------------------

Function genintzero(gn As genint)
Local i As DWord
For i = 1 To %intsize
gn.digit(i) = 0
Next
End Function

'--------------------------------------------------------------

Function genintsetsign(gn As genint, s As Integer)
gn.sign = s
End Function

'--------------------------------------------------------------

Function genintgetstringrep(gn As genint) As String
Local i, j As DWord
Local vs As String

vs = ""

If gn.sign > -1 Then
vs = vs & "+"
Else
vs = vs & "-"
EndIf

i = 0

Do
If %intsize - i = 0 Then Return "0"
If gn.digit(%intsize - i) > 0 Then Exit Do
i += 1
Loop
 
For j = i + 1 To %intsize
vs = vs & Chr$(48 + gn.digit(%intsize - j + 1))
Next

Return vs
End Function

'--------------------------------------------------------------

Function genintgetsize(gn As genint) As DWord
Local i As DWord

i = 0

Do
If %intsize - i = 0 Then Return 0
If gn.digit(%intsize - i) > 0 Then Exit Do
i += 1
Loop

Return %intsize - i
End Function

'--------------------------------------------------------------

Function genintseterror(gn As genint)
gn.error = TRUE
End Function

'--------------------------------------------------------------

Function genintshiftleft(gn As genint, nplaces As DWord)
Local sz, i, j As DWord

sz = genintgetsize(gn)

For i = 1 To sz
j = sz - i + 1
gn.digit(j + nplaces) = gn.digit(j)
Next

For i = 1 To nplaces
gn.digit(i) = 0
Next

End Function

'--------------------------------------------------------------

Function genintadd(a As genint, b As genint, c As genint)
Local i As DWord
Local s, v, k As Byte
Local asize, bsize, maxsize
Local temp As genint

asize = genintgetsize(a)
bsize = genintgetsize(b)
maxsize = asize
If bsize > maxsize Then maxsize = bsize

genintzero(temp)

k = 0
For i = 1 To maxsize + 1
s = a.digit(i) + b.digit(i) + k
v = Mod(s, 10)
k = s \ 10
temp.digit(i) = v
Next

For i = 1 To %intsize
c.digit(i) = temp.digit(i)
Next

End Function

'--------------------------------------------------------------

Function genintmult(a As genint, b As genint, c As genint)
Local i, j As DWord
Local p, v, k As Byte
Local asize, bsize, tsize As DWord
Local temp As genint

genintzero(c)
genintzero(temp)

asize = genintgetsize(a)
bsize = genintgetsize(b)
tsize = asize + bsize
If tsize > %intsize Then
genintseterror(c)
Return
EndIf

For i = 1 To asize
k = 0
For j = 1 To bsize + 1
p = a.digit(i) * b.digit(j) + k
v = Mod(p, 10)
k = p \ 10
temp.digit(j) = v
Next

genintshiftleft(temp, i - 1)
genintadd(c, temp, c)
genintzero(temp)
Next

End Function

' output --------------------------------------------------------------------------------------------

power sum = x^0 + x^1 + x^2 + .. + x^n

x =  2
n =  100

-------------------------------------------------------------

Calculating power sum by slow method.

sum =  +2535301200456458802993406410751

elapsed time =  506  seconds

-------------------------------------------------------------

Calculating power sum by fast method.

sum =  +2535301200456458802993406410751

elapsed time =  16  seconds

-------------------------------------------------------------

Done. Press a key to quit.

**danbaron** · 31-08-2011, 08:13

I changed the program a little bit, and ran it again, this time with x=2 and n=300.

Theoretically, the ratio between the slow time and the fast time, when only considering the multiplications, should be
(n+1)/2.

For n=300, that would be 150.5

Here, I got, 2931/30 = 97.7.

But, there are a lot of additions too, and the number of additions are the same for both methods, so, I think that
should tend to make the ratio smaller.

'changed code -----------------------------------------------------------------------------------------------------------------
Uses "console"
'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%neg = -1
%pos = +1

%intsize = 100

Type genint
digit(%intsize) As Byte
sign As Integer
error As Byte
End Type

'--------------------------------------------------------------

Function TBMain()
Local s As String
Local t1, t2, tt As Integer
Local x, ps As genint
Local n As DWord

'Set value of x (string s).
s = "2"

'Set number of terms, n.
n = 300

'--------------------------------------------------------------

Function slowevaluation(ByRef x As genint, n As DWord, psum As genint)
Local xp As genint
Local s As String
Local i As Integer

s = "0"
genintsetvalue(psum,s)
If n=0 Then Return
For i = 0 To n
genintpower(x,i,xp)
genintadd(psum,xp,psum)
Next
End Function

'--------------------------------------------------------------

Function fastevaluation(ByRef x As genint, n As DWord, psum As genint)
Local xp As genint
Local s As String
Local i As Integer

s = "0"
genintsetvalue(psum,s)
If n=0 Then Return
s = "1"
genintsetvalue(psum,s)
genintsetvalue(xp,s)
For i = 1 To n
genintmult(xp,x,xp)
genintadd(psum,xp,psum)
Next
End Function

'--------------------------------------------------------------

Function genintpower(ByRef x As genint, p As DWord, ByRef xp As genint)
Local s As String
Local i As Integer

s = "1"
genintsetvalue(xp,s)
If p=0 Then Return
For i = 1 To p
genintmult(xp,x,xp)
Next
End Function

'--------------------------------------------------------------

Function genintfactorial(n As Integer) As String
Local num1, num2 As genint
Local s As String
Local i As Integer

s = "1"
genintsetvalue(num1,s)
For i = n To 2 Step -1
PrintL i
s = Format$(i)
genintsetvalue(num2,s)
genintmult(num1,num2,num1)
Next
Return genintgetstringrep(num1)
End Function

'--------------------------------------------------------------

Function genintmult(a As genint, b As genint, c As genint)
Local i, j As DWord
Local p, v, k As Byte
Local asize, bsize, tsize As DWord
Local temp, temp1 As genint

genintzero(temp)
genintzero(temp1)

asize = genintgetsize(a)
bsize = genintgetsize(b)
tsize = asize + bsize
If tsize > %intsize Then
genintseterror(c)
Return
EndIf

For i = 1 To asize
k = 0
For j = 1 To bsize + 1
p = a.digit(i) * b.digit(j) + k
v = Mod(p, 10)
k = p \ 10
temp.digit(j) = v
Next

genintshiftleft(temp, i - 1)
genintadd(temp1, temp, temp1)
genintzero(temp)
Next
genintsetequal(c,temp1)
End Function

'output -----------------------------------------------------------------------------------------------------------------------

power sum = x^0 + x^1 + x^2 + .. + x^n

x =  2
n =  300

-------------------------------------------------------------

sum = +4074071952668972172536891376818756322102936787331872501272280898708762599526673412366794751

elapsed time = 2931 seconds

-------------------------------------------------------------

Calculating power sum by fast method.

sum = +4074071952668972172536891376818756322102936787331872501272280898708762599526673412366794751

elapsed time = 30 seconds

-------------------------------------------------------------

Done. Press a key to quit.

**REDEBOLT** · 31-08-2011, 10:32

sum = x^0 + x^1 + x^2 + x^3 + .. + x^n

Then x*sum = x^1 + x^2 + x^3 + .. + x^(n+1)

Subtract second equation from the first

sum - x*sum = x^0 - x^(n+1)

This is called "telescoping cancellation"

So, sum*(x-1) = x^(n+1) - x^0

           x^(n+1) - x^0
And sum = ---------------
               (x-1)

Let set x=2, and n=100.
 
             2^(101) - 1
Then sum = --------------
                  1

So sum = 2,535,301,200,456,458,802,993,406,410,751

And for n = 300, sum = 4.0740719526689721725368913768188e+90

The results are from the windows calculator.

Regards,
Bob

**peter** · 31-08-2011, 12:51

is this mathematically forum?
where is the fun at all ? Games, OpenGl and drolly programs.

**henksq49** · 31-08-2011, 17:45

I appreciate the nice programming aspects but cannot help by noting that taking the limit of n to infinity is even nicer:

sum= 1 + x^1 + x^2 + x^3 + ... = 1 + x*( 1 + x^1 + x^2 + ...)

= 1 + x*sum i.e. sum = 1/(1-x)

now take x = 2 and find 1 + 2 + 4 + 8 + 16 + ... = -1 !!!

Juggling with infinite sums can give funny results.

Kind regards, henks

**danbaron** · 31-08-2011, 20:15

A geometric series,

a + ar + ar² + ar³ + .. (to infinity)

only converges to the value,

sum = a/(1-r),

if abs(r) < 1, right?

--------------------------------------------------------------------------------

For a conditionally convergent series, you can make the sum anything you want, depending on the order in which you add the terms.

The series,

1 - 1/2 + 1/3 - 1/4 + .., is conditionally convergent.

For this series, by rearranging the terms, you can make the sum any number, you can even make it diverge (become infinite).

The reason is that, the positive terms alone form a divergent series, and so do the negative terms alone.

--------------------------------------------------------------------------------

Here is what I wanted to show in the program.

The slow method calculates the sum like this:

1 + x + x*x + x*x*x + x*x*x*x + .. + x*x*x*x*x* (n times for the last term)

You can see that in the program, during each iteration, it calls the power function, which calculates x^i by multiplying x*x*x.., i times.

The number of multiplications of x is, n*(n+1)/2, or, (n^2+n)/2.

The fast method calculates the sum like this:

xp = 1
sum = 0
for i = 0 to n
xp = xp*i
sum = sum + xp
next

So, in order to find the sum, it only calculates x*x*x*x.., up to n, one time.

The number of multiplications of x, is n.

The time ratio of the slow method to the fast method, when only taking into account the x multiplications is,

(n^2+n)/2/n, = (n+1)/2.

--------------------------------------------------------------------------------

For x=2, and n,

you can calculate the sum instantaneously.

It is,

x^(n+1) - 1, (what Bob said).

**jack** · 31-08-2011, 23:50

Dan, try this inpower function, of course you need to change it to use your math functions.

Function odd(ByVal n As Long) As Long
  Function = n And 1
End Function
Function even(ByVal n As Long) As Long
  Function = 1-(n And 1)
End Function
Function intpower(ByRef x As Ext, p As DWord, byref xp as ext)
Local xp, xt As Ext
Local rm As Integer
rm = p
xp = 1
xt = x
While rm>0
  While even(rm)
    rm = rm\2
    xt = xt*xt
  Wend
  rm -= 1
  xp=xp*xt
Wend
End Function

**danbaron** · 01-09-2011, 07:09

Concerning fun and games, post them.

-------------------------------------------------------------------------------------------------------------------------------

I see what Bob did with the series.

That is really nice, and works for any value of x.

-------------------------------------------------------------------------------------------------------------------------------

Wow, what a difference your function makes, Johan!

Very smart, continually squaring x.

For x=2 and n=300, the time went from 2931 seconds, to 46 seconds!!

(I changed your "even" function, because, I didn't understand it.)

'changed code -----------------------------------------------------------------------------------------------------------------

Function slowevaluation(ByRef x As genint, n As DWord, psum As genint)
Local xp As genint
Local s As String
Local i As Integer

s = "0"
genintsetvalue(psum,s)
If n=0 Then Return
For i = 0 To n
PrintL i
genintpowerjack(x,i,xp)
genintadd(psum,xp,psum)
Next
End Function

'--------------------------------------------------------------

Function even(ByVal n As Long) As Long
If Mod(n,2)= 0 Then Return TRUE
Return FALSE
End Function

'--------------------------------------------------------------

Function genintpowerjack(ByRef x As genint, p As DWord, ByRef xp As genint)
Local xt As genint
Local rm As Integer
Local s As String
rm = p
s = "1"
genintsetvalue(xp,s)
genintsetequal(xt,x)
While rm>0
While even(rm)
rm = rm\2
genintmult(xt,xt,xt)
Wend
rm -= 1
genintmult(xp,xt,xp)
Wend
End Function

'output -----------------------------------------------------------------------------------------------------------------------

sum =  +4074071952668972172536891376818756322102936787331872501272280898708762599526673412366794751

elapsed time =  46  seconds

-------------------------------------------------------------

Calculating power sum by fast method.

sum =  +4074071952668972172536891376818756322102936787331872501272280898708762599526673412366794751

elapsed time =  32  seconds

-------------------------------------------------------------

Done. Press a key to quit.