問題と答え

twitterで出した問題
問題:k<=100, N<=10^18, M<=1,000,000,007のとき、1^k+2^k+…+N^k mod M を1秒で求めるプログラムを作成せよ
の想定解法。


まず、1^K mod M, 2^K mod M, …, N^k mod M を計算して足し算するという方法が普通に思いつくと思いますが、当然間に合いません。
次に、Σi^k の式を導出して、それに当てはめるという方法が思いつきますが、これはかなり大変です。係数が大変なことになります(おまけ参照)。
さらに、この方法の場合、多倍長演算が必須となります。実は64bitの範囲で計算しきる方法があります。

1^k, 2^k, … , N^kを考えるのに、それぞれを「たった1つの整数」と思わずに、「(i^0, i^1, …, i^k)の片割れ」と思うことにします。
すると、(i^0, i^1, …, i^k) -> ((i+1)^0, (i+1)^1, …, (i+1)^k) の一次変換が存在します。
k=2の場合は、((i+1)^0, (i+1)^1, (i+1)^2) = (1, i+1, i^2+2i+1) のようになります。この一次変換行列のi行j列成分は(0-origin)、iCjになっていることがわかります(i ((i+1)^0, (i+1)^1, …, (i+1)^k) の一次変換を表す行列をSとします。
すると、(S + S^2 + S^3 + … + S^N) . (1, 0, 0, …, 0)を計算すると、0<=i<=kについての答えを表すベクトルが求められます。
この(S + S^2 + S^3 + … + S^N)は、バイナリ法の要領でO(k^3 log N)で求められるので、この方法で間に合うだろうということがわかります。

おまけ

k = 100 の場合の式。Mathematicaによる。ちなみに横が欠けているようです。
(1/7365930)(-20906166358025989690245262708358112533867327978807792885478238861152818386758591744131 N +
137557060998053063780068383273859695429104777440177341411371119932294150844290959249250 N^3 -
271526754925460375939402892699167005770418367491814037763270375918276117880045400685095 N^5 +
255224919564365277058990072820560271053016244391693441958072447278567409821755859362000 N^7 -
139942721633340949307633805624021906773714246359176928292504751825697599207672393283700 N^9 +
50224701866648902906954018304514209704736042738729105160129474087651854088194724289000 N^11 -
12710203553499459929294176966622970540528106302179132781279183003135609396427212642900 N^13 +
2389422493911616664836779243894724842240655767196057138221882634675914210951903687600 N^15 -
346803746499588279914206269588107219220145529701446756696165232486828108652911639950 N^17 +
40032933131773271765182465298879862489639859879507043787801643485644055249470618500 N^19 -
3762944885960623652955837495745226808453251910977880325399207835506839937546325670 N^21 +
293587173182084315222746839055242855198576266567076751379780768030615924494998000 N^23 -
19317261710275220507857305536477658038676759298069353253511132330734056425257524 N^25 +
1086346046682206649981174238119049428500340967792881450076255194466186738029000 N^27 -
52816776962762229474656148378493836433453154946253754168169571295149831998900 N^29 +
2242067502636701709164091147297234581622383956917258238119088560741823578000 N^31 -
83819391255919008895135145598205353342826434627838377116422016319019589925 N^33 +
2780720110373274462297712120317225285535422081908456994384941034405935550 N^35 -
82416238557170440526782248207578723483418676557618113800099842474551025 N^37 +
2195453902253926778163367246330521207960545317456334303939477321652000 N^39 -
52849418283107712906850708021164964350433395055064802149942091092520 N^41 +
1155266558764058249487792070417380223571229632489123955251647746000 N^43 -
23034391743048702146501114229324871023492482185484729894647161960 N^45 +
420611163965975694281664633353032579031247207551875474083564000 N^47 -
7059975841880803147184261282983483515340623242468228268420900 N^49 +
109300656691155091709412204719881150027752184764339561747000 N^51 -
1565681048361370422361573852778143937697238790603265659700 N^53 +
20811653287016494480134557104826009045952719006548466400 N^55 -
257396973528186174839230093964501393897010047439942600 N^57 +
2969498892765985532374691547974297108746616400738000 N^59 -
32030359935620917647849279378095194001090104866120 N^61 +
323734747967171553365585645640664627322536404000 N^63 -
3072244128228098234286487275872333745665970045 N^65 +
27428162976485483540852645494857822863274750 N^67 -
230780151810561633119001526604600278143225 N^69 +
1833166039102261941559847115582999882000 N^71 -
13769119983258859450691640597927209700 N^73 +
97942896363747870114638577311052360 N^75 -
660736626263222043038653275855300 N^77 +
4233177722036019126254899014000 N^79 -
25790011956737117676508391470 N^81 +
149596205898338659194234500 N^83 -
827155349897338962489990 N^85 +
4364646526295493438000 N^87 -
22004747918212359300 N^89 +
106149290488241000 N^91 -
491294051348100 N^93 +
2200645246800 N^95 -
9925590675 N^97 +
61382750 N^99 +
3682965 N^100 +
72930 N^101)