Sum of Powers
巾乗和[2001.11.19]
整数 k>=0 に対して、1からnまでの巾乗(k乗)の和
Fk(n) = 1k+2k+3k+....+nk
をPARI/GP 2.1.1を使って求めてみよう。
Fk(x)は、xの有理数係数の(k+1)次多項式になることは、すぐにわかる。
さらに、Fk(x)は、以下の再帰式を満たす。
Fm(x) = ((x+1)m+1-1)/(m+1)-(Σi=0m-1 binomial(m+1,i)Fi(x))/(m+1)
ここで、binomial(k,i)は、２項係数kCi=k!/{i!(k-i)!}である。
これを使って、F1(x),...,F100(x)を求めるプログラムをPARI/GPで作成する。
高校生の時に、これらを手で計算しようとしたが、m=6,7ぐらいで挫折した思い出がある。
原理的には計算できることがわかっているのに、実際には、計算が面倒で実行できず、悔しかったが、今ならパソコンが以下のように計算してくれる。
(22:14) gp > read("sumofpower.gp");
(22:21) gp > f(1,x)
%25 = 1/2*x^2 + 1/2*x
(22:21) gp > f(2,x)
%26 = 1/3*x^3 + 1/2*x^2 + 1/6*x
(22:21) gp > f(3,x)
%27 = 1/4*x^4 + 1/2*x^3 + 1/4*x^2
(22:21) gp > f(4,x)
%28 = 1/5*x^5 + 1/2*x^4 + 1/3*x^3 - 1/30*x
(22:21) gp > f(5,x)
%29 = 1/6*x^6 + 1/2*x^5 + 5/12*x^4 - 1/12*x^2
(22:21) gp > f(6,x)
%30 = 1/7*x^7 + 1/2*x^6 + 1/2*x^5 - 1/6*x^3 + 1/42*x
(22:22) gp > f(8,x)
%32 = 1/9*x^9 + 1/2*x^8 + 2/3*x^7 - 7/15*x^5 + 2/9*x^3 - 1/30*x
(22:22) gp > f(9,x)
%33 = 1/10*x^10 + 1/2*x^9 + 3/4*x^8 - 7/10*x^6 + 1/2*x^4 - 3/20*x^2
(22:22) gp > f(10,x)
%34 = 1/11*x^11 + 1/2*x^10 + 5/6*x^9 - x^7 + x^5 - 1/2*x^3 + 5/66*x
(22:22) gp > f(11,x)
%35 = 1/12*x^12 + 1/2*x^11 + 11/12*x^10 - 11/8*x^8 + 11/6*x^6 - 11/8*x^4 + 5/12*x^2
(22:22) gp > f(12,x)
%36 = 1/13*x^13 + 1/2*x^12 + x^11 - 11/6*x^9 + 22/7*x^7 - 33/10*x^5 + 5/3*x^3 - 691/2730*x
(22:22) gp > f(13,x)
%37 = 1/14*x^14 + 1/2*x^13 + 13/12*x^12 - 143/60*x^10 + 143/28*x^8 - 143/20*x^6 + 65/12*x^4 - 691/420*x^2
(22:22) gp > f(14,x)
%38 = 1/15*x^15 + 1/2*x^14 + 7/6*x^13 - 91/30*x^11 + 143/18*x^9 - 143/10*x^7 + 91/6*x^5 - 691/90*x^3 + 7/6*x
(22:23) gp > f(15,x)
%39 = 1/16*x^16 + 1/2*x^15 + 5/4*x^14 - 91/24*x^12 + 143/12*x^10 - 429/16*x^8 + 455/12*x^6 - 691/24*x^4 + 35/4*x^2
(22:23) gp > f(16,x)
%40 = 1/17*x^17 + 1/2*x^16 + 4/3*x^15 - 14/3*x^13 + 52/3*x^11 - 143/3*x^9 + 260/3*x^7 - 1382/15*x^5 + 140/3*x^3 - 3617/510*x
(22:23) gp > f(17,x)
%41 = 1/18*x^18 + 1/2*x^17 + 17/12*x^16 - 17/3*x^14 + 221/9*x^12 - 2431/30*x^10 + 1105/6*x^8 - 11747/45*x^6 + 595/3*x^4 - 3617/60*x^2
(22:23) gp > f(18,x)
%42 = 1/19*x^19 + 1/2*x^18 + 3/2*x^17 - 34/5*x^15 + 34*x^13 - 663/5*x^11 + 1105/3*x^9 - 23494/35*x^7 + 714*x^5 - 3617/10*x^3 + 43867/798*x
(22:23) gp > f(19,x)
%43 = 1/20*x^20 + 1/2*x^19 + 19/12*x^18 - 323/40*x^16 + 323/7*x^14 - 4199/20*x^12 + 4199/6*x^10 - 223193/140*x^8 + 2261*x^6 - 68723/40*x^4 + 43867/84*x^2
(22:23) gp > f(20,x)
%44 = 1/21*x^21 + 1/2*x^20 + 5/3*x^19 - 19/2*x^17 + 1292/21*x^15 - 323*x^13 + 41990/33*x^11 - 223193/63*x^9 + 6460*x^7 - 68723/10*x^5 + 219335/63*x^3 - 174611/330*x
(22:23) gp > f(21,x)
%45 = 1/22*x^22 + 1/2*x^21 + 7/4*x^20 - 133/12*x^18 + 323/4*x^16 - 969/2*x^14 + 146965/66*x^12 - 223193/30*x^10 + 33915/2*x^8 - 481061/20*x^6 + 219335/12*x^4 - 1222277/220*x^2
(22:23) gp > f(22,x)
%46 = 1/23*x^23 + 1/2*x^22 + 11/6*x^21 - 77/6*x^19 + 209/2*x^17 - 3553/5*x^15 + 11305/3*x^13 - 223193/15*x^11 + 124355/3*x^9 - 755953/10*x^7 + 482537/6*x^5 - 1222277/30*x^3 + 854513/138*x
(22:23) gp > f(23,x)
%47 = 1/24*x^24 + 1/2*x^23 + 23/12*x^22 - 1771/120*x^20 + 4807/36*x^18 - 81719/80*x^16 + 37145/6*x^14 - 5133439/180*x^12 + 572033/6*x^10 - 17386919/80*x^8 + 11098351/36*x^6 - 28112371/120*x^4 + 854513/12*x^2
(22:23) gp > f(24,x)
%48 = 1/25*x^25 + 1/2*x^24 + 2*x^23 - 253/15*x^21 + 506/3*x^19 - 14421/10*x^17 + 29716/3*x^15 - 10266878/195*x^13 + 208012*x^11 - 17386919/30*x^9 + 22196702/21*x^7 - 28112371/25*x^5 + 1709026/3*x^3 - 236364091/2730*x
(22:23) gp > f(25,x)
%49 = 1/26*x^26 + 1/2*x^25 + 25/12*x^24 - 115/6*x^22 + 1265/6*x^20 - 24035/12*x^18 + 185725/12*x^16 - 25667195/273*x^14 + 1300075/3*x^12 - 17386919/12*x^10 + 277458775/84*x^8 - 28112371/6*x^6 + 21362825/6*x^4 - 1181820455/1092*x^2
(22:24) gp > f(26,x)
%50 = 1/27*x^27 + 1/2*x^26 + 13/6*x^25 - 65/3*x^23 + 16445/63*x^21 - 16445/6*x^19 + 142025/6*x^17 - 10266878/63*x^15 + 2600150/3*x^13 - 20548177/6*x^11 + 3606964075/378*x^9 - 52208689/3*x^7 + 55543345/3*x^5 - 1181820455/126*x^3 + 8553103/6*x
(22:24) gp > f(27,x)
%51 = 1/28*x^28 + 1/2*x^27 + 9/4*x^26 - 195/8*x^24 + 4485/14*x^22 - 29601/8*x^20 + 142025/4*x^18 - 15400317/56*x^16 + 1671525*x^14 - 61644531/8*x^12 + 721392815/28*x^10 - 469878201/8*x^8 + 166630035/2*x^6 - 3545461365/56*x^4 + 76977927/4*x^2
(22:24) gp > f(28,x)
%52 = 1/29*x^29 + 1/2*x^28 + 7/3*x^27 - 273/10*x^25 + 390*x^23 - 9867/2*x^21 + 52325*x^19 - 905901/2*x^17 + 3120180*x^15 - 33193209/2*x^13 + 65581165*x^11 - 365460823/2*x^9 + 333260070*x^7 - 709092273/2*x^5 + 179615163*x^3 - 23749461029/870*x
(22:24) gp > f(29,x)
%53 = 1/30*x^30 + 1/2*x^29 + 29/12*x^28 - 609/20*x^26 + 1885/4*x^24 - 26013/4*x^22 + 303485/4*x^20 - 8757043/12*x^18 + 22621305/4*x^16 - 137514723/4*x^14 + 1901853785/12*x^12 - 10598363867/20*x^10 + 4832271015/4*x^8 - 6854558639/4*x^6 + 5208839727/4*x^4 - 23749461029/60*x^2
(22:24) gp > f(30,x)
%54 = 1/31*x^31 + 1/2*x^30 + 5/2*x^29 - 203/6*x^27 + 1131/2*x^25 - 16965/2*x^23 + 216775/2*x^21 - 2304485/2*x^19 + 19959975/2*x^17 - 137514723/2*x^15 + 731482225/2*x^13 - 31795091601/22*x^11 + 8053785025/2*x^9 - 102818379585/14*x^7 + 15626519181/2*x^5 - 23749461029/6*x^3 + 8615841276005/14322*x
.... (省略) .....
(22:29) gp > f(50,x)
%1 = 1/51*x^51 + 1/2*x^50 + 25/6*x^49 - 490/3*x^47 + 75670/9*x^45 - 416185*x^43 + 56941675/3*x^41 - 92183691110/117*x^39 + 88715129650/3*x^37 - 33921877388371/34*x^35 + 541322044801475/18*x^33 - 2413278263674516/3*x^31 + 56849978179309700/3*x^29 - 45602759590076311910/117*x^27 + 6930788058683050926*x^25 - 316006038824222942020/3*x^23 + 12150871151987276191100/9*x^21 - 1465057369388478597658465/102*x^19 + 248857551895105646012275/2*x^17 - 100303446943247952675791038/117*x^15 + 13680778417903412433848650/3*x^13 - 594659492189651782474535017/33*x^11 + 451887027095744497869967025/9*x^9 - 91571456915900369733225670*x^7 + 292261044681016960006200890/3*x^5 - 196329117914923619018719464145/3978*x^3 + 495057205241079648212477525/66*x
.... (省略) .....
(22:29) gp > f(100,x)
%2 = 1/101*x^101 + 1/2*x^100 + 25/3*x^99 - 2695/2*x^97 + 298760*x^95 - 66698170*x^93 + 43232541100/3*x^91 - 2987368590010*x^89 + 592545208316600*x^87 - 1909010939318560231/17*x^85 + 60927624576260699950/3*x^83 - 10503770178403996919537/3*x^81 + 574696979476592789539800*x^79 - 89701724868851868404757210*x^77 + 13296745470530926863904296852*x^75 - 1869298239768618416234153813290*x^73 + 248870955751990847260270884407400*x^71 - 1065245686771269279784908613651828005/34*x^69 + 3723652407297582727619274890591931075*x^67 - 10844299000116828980379757772973769420469/26*x^65 + 43950288418050613210495571589828389262800*x^63 - 4348447505694585428839166185138223415249684*x^61 + 403139711179170251736670257248480111641926600*x^59 - 34944260063316672143127900206016265956506516820*x^57 + 2825393845314372316887963516463774302220183874480*x^55 - 46975128963737486419489164499794297560673231041202090/221*x^53 + 14838677083702274079364344314958348779821717660137900*x^51 - 958463607702055700662952442255558159708363131670845130*x^49 + 57102248319760803358389861613269821873306861123018474800*x^47 - 3127153223428501512572222954783017354698250212191091945572*x^45 + 156839198684220220595062954768424384099662857573873761392200*x^43 - 93273006623793637434656479802977293641893710056414115793830132/13*x^41 + 298055222117767447988694875776788702575308931452828672542296400*x^39 - 380420681562789081339436627697748498619486609696130138245054547645/34*x^37 + 377511069257143967197314136886615170865786408764196373626268649635*x^35 - 22758671683254934243234770245768111655371809025564559292966948184145/2*x^33 + 304383493005866429515905139920856508495517057169598168611307541714600*x^31 - 93215398532963113031284566771666289746833192047887884379325411297276490/13*x^29 + 147482537396120605270641214092320919218664984298368495230915199366025300*x^27 - 13112574861745373977119865065563790341936971501269597493806710307275562234/5*x^25 + 39857448167724145521712375634202721882854746999642509687138048831663608600*x^23 - 8684587426344073606489504709883050170678418405635672010429984157277666084023/17*x^21 + 16304634906294224258925539055711850026937478314146500355475130833028845746350/3*x^19 - 612067818251686839120746668057583204003010059302601007211600982133792394509295/13*x^17 + 324388433491984944852419075920450620931865462636226130063940688368734730163320*x^15 - 1725539552167813151807602972961047761861449443882732089672204732211086637590530*x^13 + 225010984573490896048358130480328881791747873028668541553377868767760647319540900/33*x^11 - 56995948223784756021697384698478769187481110881793715780290371409596995575985270/3*x^9 + 34649381621107623485288357725441359210991177541965976048926944361209977534643400*x^7 - 16293234618989521508515025064456465992824384487957638029599182473343901462949018943/442*x^5 + 18674771685049011296614057325260991542019103825338734064995339343748060441015725*x^3 - 94598037819122125295227433069493721872702841533066936133385696204311395415197247711/33330*x
(22:31) gp > quit
Good bye!
PARI/GPを使って、巾乗和の多項式F1(x),...,F10(x)を因数分解すると、
F1(x)=x(x+1)/2
F2(x)=x(x+1)(2*x+1)/6
F3(x)=x2(x+1)2/4
F4(x)=x(x+1)(2*x+1)(3*x2+3*x-1)/30
F5(x)=x2(x+1)2(2*x2+2*x-1)/12
F6(x)=x(x+1)(2*x+1)(3*x4+6*x3-3*x+1)/42
F7(x)=x2(x+1)2(3*x4+6*x3-x2-4*x+2)/24
F8(x)=x(x+1)(2*x+1)(5*x6+15*x5+5*x4-15*x3-x2+9*x-3)/90
F9(x)=x2(x+1)2(x2+x-1)(2*x4+4*x3-x2-3*x+3)/20
F10(x)=x(x+1)(2*x+1)(x2+x-1)(3*x6+9*x5+2*x4-11*x3+3*x2+10*x-5)/66
となる。
参考文献
- [1]Joseph. H. Silverman, "A Friendly Introduction to Number Theory Second Edition", Prentice Hall Inc., p317-327, 2001, ISBN0-13-030954-0.
| Last Update: 2005.06.12 |
| H.Nakao |