Elektrotehniški vestnik 80(1-2): 34-38, 2013 Pregledni znanstveni članek
Algoritmi za vsoto potenc naravnih števil
Jurij Mihelic
Univerza v Ljubljani, Fakulteta za računalništvo in informatiko, Tržaška 25, 1000 Ljubljana, Slovenija E-pošta: jurij.mihelic@fri.uni-lj.si
Povzetek. Zaporedja imajo ključni pomen v matematiki, teoretičnem računalništvu, analizi algoritmov, teoriji izračunljivosti, računski zahtevnosti in številnih drugih natančnih znanostih. V članku se posvetimo zaporedjem k-tih potenč prvih n naravnih števil. Za takšna zaporedja za poljubna k in n izpeljemo več različnih formul za izračun delne vsote, tj. vsote prvih n členov. Na podlagi izpeljanih formul razvijemo pet algoritmov, ki za podana n in K izračunajo vse delne vsote za k med 0 in K. Prva dva izmed predstavljenih algoritmov temeljita neposredno na definičiji delne vsote, preostali trije pa na rekurenčni enačbi, ki vsoto izračuna iz vsot nizjih redov. Za vse algoritme zapišemo implementačijo in podamo asimptotično časovno zahtevnost.
Ključne besede: zaporedje, vrsta, delna vsota, rekurenčna enačba, algoritem
Algorithms for sum of powers of integers
Sequences are of key importance in mathematics, theoretical computer science, analysis of algorithms, computability theory, computational complexity and many other exact sciences. In this paper we focus on the sequences of k-th powers of the first n natural numbers. For any k and n we derive several formulas to calculate the partial sum of the sequence, i.e., the sum of the first n terms. Based on these formulas, we develop five algorithms which for any n and K compute all partial sums for any k between 0 and K. The first two algorithms are based on the definition of the partial sum and the remaining three on the recurrence equation which calculates the sum from the lower order sums. For all the algorithms we give implementations and analyze asymptotic time complexity.
nega n G N izračunamo po dobro znani formuli [3]:
En . n(n + 1)
4 =
i= 1
2
(1)
Pri zaporednih n G N delne vsote ^n=i i tvorijo novo zapredje 1, 3,6,10,15,..., tako imenovanih trikotniških števil. Ta pomenijo število objektov, ki jih lahko razporedimo v obliko enakostraničnega trikotnika (glej Sliko 1).
cro


1 Uvod
Zapredja spadajo med najbolj raziskovane objekte v matematiki. Prav tako imajo velik pomen v (teoretičšnem) računalništvu, še zlasti pri analizi algoritmov, teoriji izračšunljivosti in zahtevnosti.
Z zaporedjem lahko opišemo časovno ali prostorsko zahtevnost algoritma, pri čemer n-ti člen zaporedja podaja zahtevnost algoritma pri vhodnih podatkih velikosti n [1]. Zaporedje lahko opišemo z rodovno funkcijo, področje, ki se s tem poglobljeno ukvarja, se imenuje analitična kombinatorika [2].
V nadaljevanju se bomo ukvarjali s posplosšeno obliko enega najbolj znanih zaporedij, tj. z zaporedjem potenč reda k prvih n naravnih števil. Ce je k čelo število, so tudi vsi čšleni taksšnega zaporedja čela sštevila. V čšlanku bomo namenili pozornost algoritmom za izračšun vsote taksšnega zaporedja za poljubna k in n.
Najbolj znano in tudi enostavno je verjetno zaporedje naravnih števil 1,2, 3,4, 5 ..., katerega vsoto do poljub-
Prejet 4. april, 2013 Odobren 15. april, 2013
Slika 1: Trikotniška števila 1, 3, 6 in 10
Formula (1) koristi pri analizi algoritmov. Na primer: preštejmo, kolikokrat se izvede tretja vrstiča v naslednjem algoritmu.
Obe zanki for preprosto spremenimo v vsoti in ju seštejemo:
nin
E E i = E i =
i=1 j= 1 i= 1
2
Dobro znano je tudi zaporedje kvadratov naravnih števil oz. kvadratnih števil, tj. 1,4, 9,16,25,..., ki podajajo število objektov, razmeščenih v obliko kvadrata (glej sliko 2). Vsoto prvih n kvadratnih števil izračunamo po formuli:
E
i=1
,2 _ n(n + 1)(2n + 1) 4 = 6 '
O O o o o o
Slika 2: Kvadratna števila 1,4, 9 in 16
Nadalje, za zaporedne n G N vsote po enačbi (2) tvorijo novo zaporedje 1, 5,14, 30, 55,... tako imenovanih kvadratnih piramidnih števil, ki podajajo število objektov, razmeščenih v obliko štiristrane piramide s kvadratno osnovno ploskvijo (glej sliko 3). Prav tako n-to kvadratno piramidno število npr. podaja število vseh mogočih kvadratov v šahovnici velikosti n x n.
M A»
Slika 3: Kvadratna piramidna števila 1, 5,14 in 30
V nadaljevanju članka se bomo ukvarjali s posplošeno vrsto prvih n potenc reda k naravnih števil. Ta je natančneje definirana s formulo:
Sfc(n) = £ ik = 1k + 2k + 3k + ... + '
(3)
i=1
Vsote različnih redov so med seboj večkrat povezane z enačbami, npr. S3(n) = (S1(n))2 (vsoti kubov zato pravimo tudi kvadrat trikotniškega števila).
V analizi algoritmov pogosto srečamo tudi harmonično vrsto:
H (n) = S-1(n)
n 1
£ I.
i=1
formulo (3). Izračun po formuli je treba izvesti (K +1)-krat; za vsak k od 0 do K. Implementačija samega algoritma je preprosta. Potrebujemo dve zanki: prvo prek parametra k in drugo za sam izračšun vsote.
1	sums = [n, 0, 0, ...,
2	for k = 1 to K do
3	sum = 1
4	for i = 2 to n do
•h
5	sum += i
6	sums [k] = sum
0]
In ne nazadnje neskončšna hiperharmoničšna vrsta, ki je v matematični analizi znana kot Riemanova zeta funkcija:
Z (s) = s_s(^) = £ rs.
i=i
V literaturi [4] zasledimo tudi posplosšitev vrste (3), kjer namesto zaporedja prvih n sštevil nastopa aritmetičšno zaporedje (b + ia).
V nadaljevanju čšlanka bomo predstavili večš algoritmov, ki za podana n in K izračunajo vse vsote Sk(n) za k med 0 in K. Najprej predstavimo algoritma, ki temeljita na definičiji Sk (n) (enačba (3)). Nato sledita algoritma, zasnovana na rekurenčnih enačbah, ki Sk(n) podajata kot vsoto predhodno izračunanih Si(n), kjer je i < k. In kot zadnjega opišemo še algoritem, ki ravno tako temelji na rekurenčni enačbi, le da za izračun potrebuje vsak drugi člen. Na konču opišemo še nekaj mogočših optimizačij predstavljenih algoritmov.
V vrstiči 1 ustvarimo tabelo sums dolzine K +1, pri čemer prvi element tabele neposredno nastavimo na S0(n) = n (seštevek n enič). Vrstiča 5 vsebuje izračun ik, ki ga po potrebi lahko implementiramo tudi z dodatno zanko. Po izvedbi algoritma k-ti element tabele sums vsebuje vrednost Sk(n).
Asimptotičšna čšasovna zahtevnost algoritma je O(K2n), prostorska pa O(K). Pri tem smo upoštevali, da operačija potenčiranja zahteva O (k) mnozenj. Opozorimo še, da gre pravzaprav za eksponentno časovno zahtevnost, ker za zapis vhoda, tj. števil K in n, potrebujemo le log(Kn) bitov.
3 Optimizacija potenciranja
Operačijo potenčiranja v algoritmu iz predhodnega razdelka je mogoče pohitriti. Vsote računamo postopoma po naraščajočih k: potenče reda k pa seveda ni treba vsakičš račšunati znova, ampak je mogočše k-to potenčo izračunati iz (k - 1)-ve. Velja namreč:
xk = x • xk—1
Nov algoritem potrebuje dodaten prostor, tj. tabelo velikosti n - 1 za hranjenje potenč reda k za vse 2, 3, 4, . . . , n. Implementačija je podobna algoritmu iz predhodnega razdelka.
1	sums = [ n , 0 , 0 , . . . ,	0]
2	pows = [ 2 , 3 , 4 , . . . , n ]	
3	for k = 1 to K do	
4	sum = 1	
5	for i = 2 to n do	
6	sum += pows[i -	2]
7	pows[i — 2] *=	i
8	sums [ k] = sum	
V drugi vrstiči se nahaja iničializačija tabele pows (dolzine n - 1), ki hrani potenče reda k. Dodana je še vrstiča 7, ki iz predhodne potenče reda k - 1 izračuna potenčo reda k.
Asimptotičšna čšasovna zahtevnost algoritma je O(Kn). S preprosto optimizačijo smo za red velikosti zmanjšali časovno zahtevnost. Prostorska zahtevnost algoritma je O(K + n).
2 IZRAČUN po definiciji
Najprej si oglejmo algoritem, ki ga dobimo, če vsote Sk (n) računamo neposredno po definičiji, podani s
4 Alternirajoča vrsta
V tem razdelku razvijemo algoritem za zadani problem s pomočjo enačbe, ki med seboj povezuje vsote Si(n),
k
kjer je 0 < i < k. Enačba temelji na alternirajoči vrsti členov Si(n). Najprej zapišimo enačbo, nato jo bomo izpeljali:
k
T
¿=0
k +1
( —1) ¿Si(n) =
- „k + 1
(4)
V nadaljevanju bomo večkrat potrebovali binomski izrek, zato ga zapisimo:
(x + y)n = T
0
x y
pri čemer je binomski koeficient definiran kot:
n!
(n — i)! i!
(5)
(6)
n-1
Sk+i(n) =
k+1 + n
k+1
k + — 1)k+1-iSi(n) + nk+1.
Sk (n)
1
k + 1
nk+1 +
k-1
T
¿=0
k + 1
(—1)
k+1-
¿Si(n)
(7)
Preden dokončno zapišemo psevdokodo algoritma, naštejmo nekaj ugotovitev. Potenčo nk+1 lahko izračunavamo sproti ob samem računanju vsote. Pri računanju vsote je treba upoštevati vse predhodno izračunane Si(n), kjer je 0 < i < k. Vrednosti Si(n) alternirajoče enkrat prištevamo, drugič odštevamo; to je odvisno od ( —1)k+1-i. Začetni člen pri i = 0 je lahko pozitiven ali negativen, kar je odvisno od k. Zadnji člen pri i = k — 1 pa je ( — 1)2 = 1, torej vedno pozitiven. Cš e vsoto začšnemo račšunati od k — 1 proti 0, lahko vedno začnemo s prištevanjem. Poleg tega je tako tudi lazše račšunati binomske koefičiente, kot bomo videli v nadaljevanju. Pri seštevanju je treba posamezne
člene Si(n), kjer 0 < i < k še mnoziti z binomskim koefičientom (k+^.
Oglejmo si še izračun binomskih koefičientov. Izračun po definičiji (6) se pogosto izjalovi zaradi faktoriele oz. prevelikih sštevil. Poleg tega je počšasen, ker izračšun n! vsebuje n—2 mnozenj. Spomnimo se, daje binomski ko-efičient mogoče izračunati tudi po naslednji rekurenčni enačbi [3]:
k + 1
i
+
k
i1
Izpeljavo rekurenčne formule za Sk (n) začnemo iz definičije za vsoto Sk+1(n). Najprej izrazimo zadnji, tj. n-ti člen, nato premaknemo sumačijski indeks za ena in upoštevamo, da je člen pri j = 1 enak 0. Nadaljujemo z uporabo binomskega izreka, nato obrnemo vrstni red vsot.
Takšen izračun je za naš algoritem pravzaprav zelo priročen. Rekurenčna enačba nam da dobro znani Pa-scalov trikotnik (glej tabelo 1). Pri računanju Sk(n) potrebujemo prvih k koefičientov vrstiče (k +1) trikotnika. Vedno potrebujemo le eno vrstičo, ki jo z uporabo rekurenčšne enačšbe izračšunamo iz predhodne. Celotne vrstiče ni treba naračšunati vnaprej, ampak lahko ustrezni koefičient izračšunamo natanko takrat, ko ga potrebujemo.
56
j=1 n
= T(j — 1)k+1 + nk+1 j=1
n k+1 /, . , \
= TT (k + —1)k+1-i + nk+1 j = 1 ¿=0 ^ i ' k+1
= T+1
¿=0
V zadnjem izrazu je čšlen v vsoti pri i = k + 1 enak Sk+1(n), kar se odšteje s členom na levi strani enačbe. Preostanek je lepše zapisan v obliki enačbe (4).
Lotimo se implementačije algoritma po izpeljani enačbi. Najprej izrazimo člen Sk (n), ki ga zelimo izračšunati:
k\i	0	1	2	3	4
0	1				
1	1	1			
2	1	2	1		
3	1	3	3	1	
4	1	4	6	4	1
5	1	5	10	10	5
6	1	6	15	20	15
Pasčalov trikotnik za				0<	k
1
6
vsebuje števila (k) za 0 < i < k
Na voljo imamo dovolj informačij, da lahko zapišemo implementačijo. Zopet potrebujemo dve zanki: prvo prek parametra k in drugo za sam izračšun vsote po formuli (7).
1	binoms = [1 , 0, 0, ... , 0]
2	sums = [n, 0, 0, ... , 0]
3	nk1 = n
4	for k = 1 to K do
5	sum = 0 , neg = 1
6	for i = k - 1 to 0 by -1 do
7	if i > 0 then
8	binoms[i] += binoms[i —1]
9	sum += neg *binoms [ i ] * sums [ i ]
10	neg = —neg
11	nk1 *= n
12	sums[k] = (sum+nk1) / (k + 1)
13	binoms [ k] = k + 1
Vrstiče od 1 do 3 vsebujejo iničializačijo spremenljivk. Tabela binoms (dolzine K +1) pomeni aktivno vrstičo Pasčalovega trikotnika. Spremenljivka nk1 vsebuje potenčo nk+1, ki jo, kot je ze bilo rečeno, računamo sproti v prvi zanki. Obe tabeli, tako binoms kot sums, imata svoj prvi element ze nastavljen na pravo vrednost.
Zunanja zanka v vrstiči 4 je namenjena izračšunu Sk (n) za vse k od 1 do K .V akumulatorju sum se
k
n
1
računa skupna vsota, spremenljivka neg je indikator prištevanja oz. odštevanja. Notranja zanka z začetkom v vrstici 6 je implementacija zgoraj zapisane rekurenčne formule. Najprej v vrstici 8 izračunamo nov binomski koeficient, nato ga pomnozimo (vrstica 9) z ustrezno ze predhodno izračunano vsoto, rezultat pa prištejemo oz. odštejemo v akumulator. Po zanki dokončamo izračun vsote in rezultat shranimo, nato sledi le še priprava na morebitni naslednji korak zanke.
Asimptotična časovna zahtevnost algoritma je O(K2), prostorska pa O(k). Algoritem je torej asimptotično hitrejsši od algoritma, opisanega v predhodnem razdelku.
5 Nealternirajoca vrsta
Algoritem po alternirajoči vrsti zahteva nekaj dodatnega dela in pazljivosti pri seštevanju členov. Ze Blaise Pascal naj bi poznal rekurenčno enačbo, v kateri se členi v vrsti le prištevajo:
+ ')s,(n) = (n + 1)'+' - 1.
^A + 1
,=0
(8)
Veljavnost enačbe je mogoče pokazati s pomočjo indukcije, kar pa je ne pomaga izpeljati. To lahko storimo na več načinov: z rodovnimi funkcijami [7], z uporabo kombinatorike [8] ali z uporabo algebraičnih metod [9], [10]. Izpeljavo povzemamo po [10]. Začnimo z naslednjo enačbo:
E
j=i
((j + 1)'+' - j'+') =(n +1)'+' - 1
Da enačba res drzi, vidimo, če zapišemo nekaj členov vrste:
(2'+'-1'+') + (3'+1 -2'+') + - • •+((n+1)'+1 -n'+')
Večina členov se odšteje, ostaneta le: (-1)'+' in (n + 1)'+'.
Lotimo se leve strani identitete. Najprej člen (j + 1)'+' razvijemo po binomskem izreku, nato se člen j'+' razvoja odšteje z -j'+'. V preostanku obrnemo vrstni red seštevanja vrst, upoštevamo definicijo za S,(n) in enačšba (8) je izpeljana:
£ (£ k+1 j' - j*+')-
j = ' i=0
£ £ (k+1)j,
j = 'i=0 '
E(k + 1 )S,(n).
i=0
Lotimo se še implementacije algoritma. Najprej enačbo (8) preoblikujmo, da bo S' (n) na levi strani,
vse drugo pa na desni: 1

S' (n) =
k + 1
(n +1)'+' - 1	(k +1 )S,(n)
(9)
Dobljena enačšba je podobna enačšbi (7), le da vsota ne vsebuje člena (-1)'+'-i, ki povzroči alternacijo predznaka členov, kar malce poenostavi algoritem. Pri implementaciji se zgledujemo po algoritmu iz predhodnega razdelka.
1	binoms = [ 1 , 0, 0, ...	, 0]
2	sums = [n, 0, 0, ...,	0]
3	n1k1 = n + 1	
4	for k = 1 to K do	
5	sum = 0	
6	for i = k - 1 to 0	by —1 do
7	if i > 0 then	
8	binoms [ i ] +=	binoms [i — 1]
9	sum += binoms [ i ]	* sums [ i ]
10	n1k1 *= n + 1	
11	sums[k] = (n1k1 —1—sum) / (k + 1)	
12	binoms [ k] = k + 1	
Asimptotičšna čšasovna zahtevnost algoritma je enaka kot pri algoritmu iz predhodnega razdelka, tj. O(K2).
6 Prepolovitev dela
V tem razdelku najprej izpeljemo rekurenčno enačbo, ki za izračun S'(n) ne potrebuje vseh vrednosti S,(n), kjer i < k, ampak le vsako drugo. Nato zapišemo še implementačijo algoritma po izpeljani enačšbi.
Vzemimo enačbi (4) in (8) in ju seštejmo, tako dobimo:
ik + M S,(n)((-1)'-' + 1) = (»+1)'+'+»'+1 —1.
+1
0
V vsoto vpeljimo nov indeks j = k - i in upoštevajmo simetričnost binomskega koeficienta (') = (nnJ:
k + 1N
(k + S'-j(n) ((-1)j+1) = (n+1)'+'+n'+'-1.
'
jSv j + 1,
Opazimo, da je vsak drugi čšlen vrste enak ničš oz. natančšneje, da velja:
0 j je lih 2 j je sod.
Vpeljimo nov indeks i, pri čemer j = 2i, in dobimo:
L'/2J
(-1)j + 1
2
k + 1 2i + 1
S'-2,(n) = (n +1)'+' + n'+' - 1.
Dobljeno enačbo preoblikujmo, da bo primerna za uporabo v algoritmu:
1 ((n +1)'+' + n'+' - 1
S' (n) =
k+1
L'/2J
2
L'/2J	)
- E 2i ++1 S'-2i(»>). (10)
i=' x 7
Zapišimo se implementacijo. Na vsakem koraku je treba izračunati celotno vrstico Pascalovega trikotnika, Čeprav gre vsota v enacbi (10) le do [k/2j. Zato v algoritmu nastopata dve notranji zanki (vrstici 9 in 12). Prva skrbi za izracun binomskih koeficientov, druga pa dodaja se izracun vsote po enacbi (10).
1	binoms = [1 , 1 , 0, 0, ... , 0]	
2	sums = [n, 0, 0, ... , 0]	
3	nlkl = n + 1	
4	nkl = n	
5	for k = 1 to K do	
6	sum = 0	
7	binoms [ k + 1 ] = 1	
8	i = k	
9	while i > k/2 do	
10	binoms [i] += binoms [ i —	1]
11	i —= 1	
12	while i >= 1 do	
13	binoms [ i ] += binoms [ i —	1]
14	sum += binoms [2* i + 1] *	
15	sums [k — 2* i ]	
16	i —= 1	
17	n1 k1 *= n + 1	
18	nk1 *= n	
19	tmp = (n1k1 + nk1 — 1) / 2	
20	sums [k] = (tmp — sum) / (k	+ 1)
Asimptoticna casovna zahtevnost je enaka O (K2), vendar je število mnozenj prepolovljeno. Opozorimo, da gre pravzaprav za mnozenja velikih števil - casovna zahtevnost takega mnozsenja pa ni konstantna, ampak je odvisna od velikosti števila.
7 Optimizacije
Omenimo še nekaj mogocih splošnih optimizacij, ki jih lahko uporabimo v zgoraj opisanih algoritmih. Vsote za majhne k lahko vnajprej izracunamo po dobro znanih formulah. S tem lahko prihranimo nekaj iteracij zunanje zanke v vseh omenjenih algoritmih.
Vrstice v Pascalovem trikotniku so simetricne po pravilu (?) = (n„J. Z upoštevanjem simetricnosti lahko izracšunamo le polovico vrstice Pascalovega trikotnika. Druge polovice nam niti ni treba hraniti, v algoritmu pa moramo zaradi tega dodati preverjanje parametra i binomskega koeficienta.
Simetricnost lahko izrabimo tudi tako, da zmanjšamo število mnozenj velikih števil. Pri izracunu vsote namrec lahko upoštevamo obe Si in Sk+1-i hkrati. Kot primer, predelajmo enacšbo (9) tako, da zmanjsšamo sštevilo mnozenj za polovico:
= kri (<"+ 1)k+1 -1
- So(n) - (k +1)Si(n) -
Sk (n)

k/2
E
i=2
k + 1
(Si(n) + Sk+i_i(n))
Zapisana enacšba velja pri sodih vrednostih k, pri lihih k ravnamo podobno, le srednji cšlen moramo posebej uposštevati.
8 Sklep
Zaporedja in njihove vsote so pomembni na sštevilnih znanstvenih področjih. V clanku smo se ukvarjali z vsoto zaporedja potenc reda k prvih n naravnih sštevil. Opisali smo vec algoritmov za izracun te vsote. Vse algoritme smo matematicšno utemeljili. Poleg tega smo zapisali tudi njihovo implementacijo, na podlagi katere smo podali tudi asimptoticšno cšasovno zahtevnost.
Literatura
[1]
Robert Sedgewick, Philippe Flajolet, An Introduction to the Analysis of Algorithms, 2nd ed. Addison-Wesley, 2013.
[2]	Philippe Flajolet, Robert Sedgewick, Analytic Combinatorics, Cambridge University Press, 2009.
[3]	Ronald Graham, Donald Knuth, and Oren Patashnik, Concrete Mathematics: A Foundation for Computer Science, Addi-son-Wesley, 1994.
[4]	N. Gauthier, "Sum of the m-th powers of n successive terms of an arithmetic sequence: bm + (a + b)m + (2a + b)m + ... + ((n — l)a + b)m", International Journal of Mathematical Education in Science and Technology, vol. 37, no. 3, pp.
[5]	Michael Z. Spivey, "The Euler-Maclaurin Formula and Sums of Powers", Mathematics Magazine, vol. 79, no. 1, pp. 61-65, 2006.
[6]	David M. Bloom, "An Old Algorithm for the Sums of Integer Powers", Mathematics Magazine, vol. 66, no. 5, pp. 304-305, 1993.
[7]	J. Riordan, CombinatoriaI ldentities, Wiley, 1968. J. L. Paul, "On the Sum of the fc-th Powers of the First n Integers", The American Mathematical Monthly, vol. 78, no. 3, pp. 271-272, 1971.
Clive Kelly, "An Algorithm for Sums of Integer Powers", Mathematics Magazine, vol. 57, no. 5, pp. 296-297, 1984.
[10] Dumitru Acu, "Some Algorithms for the Sums of Integer Powers", Mathematics Magazine, vol. 61, no. 3, pp. 189-191, 1988.
Jurij Mihelic je leta 2006 doktoriral s podrocja prilagodljivosti v optimizacijskih problemih na Fakulteti za racunalništvo in informatiko Univerze v Ljubljani. Je asistent in raziskovalec na omenjeni fakulteti. Njegovo podrocje raziskovanja vkljucuje nacrtovanje in analizo algoritmov ter podatkovnih struktur, kombinatoricšno optimizacijo, teorijo grafov.
[8]
[9]