(* Super Power Mod by Ilan Vardi, from "Computational Recreations in Mathematica," Addison-Wesley Publishing Co., Redwood City, CA, 1991, pages 226-229. *) SuperPowerMod[a_, k_, n_] := 0 /; Mod[a, n] == 0; SuperPowerMod[a_, k_, n_] := Mod[a, n] /; k == 1; SuperPowerMod[a_, k_, n_] := PowerMod[a, a, n] /; k == 2; SuperPowerMod[a_, k_, n_] := Mod[a, 2] /; n == 2; (* It was incorrectly stated as Mod[2, a] on page 227. *) SuperPowerMod[a_, k_, n_] := SuperPowerMod[a, k, n] = Block[{fi = FactorInteger[n], m = Mod[a, n], i, rprimem, gcdn}, rprimem = m; gcdlist = Block[{i = 0}, While[ Mod[ rprimem, #[[1]]] == 0, rprimem /= #[[1]]; i++]; {#[[1]], #[[2]], i}] & /@ Select[fi, Mod[m, #[[1]]] == 0 &]; If[ gcdlist == {}, PowerMod[m, SuperPowerMod[a, k - 1, EulerPhi[n]], n], gcdn = Times @@ (#[[1]]^#[[2]] & /@ gcdlist); If[ LogStar[a, Max[ #[[2]] - #[[3]] & /@ gcdlist]] > k - 1, PowerMod[m, SuperPower[a, k-1], n], If[ gcdn == n, 0, Mod[ gcdn PowerMod[ gcdn, -1, n/gcdn] PowerMod[ m/rprimem, SuperPowerMod[a, k - 1, EulerPhi[n/gcdn]], n/gcdn] PowerMod[ rprimem, SuperPowerMod[a, k - 1, EulerPhi[n]], n], n]]]]]; LogStar[a_, n_] := 0 /; n < a; LogStar[a_, n_] := LogStar[a, Log[a, N[n]]] + 1; LogStarStar[a_, n_] := 0 /; n < a; LogStarStar[a_, n_] := LogStarStar[a, LogStar[a, n]] + 1; SuperSuperPowerMod[a_, k_, n_] := Mod[a, n] /; k == 1; SuperSuperPowerMod[a_, k_, n_] := SuperPowerMod[a, a, n] /; k == 2; SuperSuperPowerMod[a_, k_, n_] := SuperPowerMod[a, SuperSuperPower[a, k - 1], n] /; LogStarStar[Log[2, n]] > k; SuperSuperPowerMod[a_, k_, n_] := SuperPowerMod[a, n, n] SuperPower[a_, 1] := a; SuperPower[a_, k_] := a^SuperPower[a, k - 1]; SuperSuperPower[a_, 1] := 1; SuperSuperPower[a_, k_] := SuperPower[a, SuperSuperPower[a, k - 1]];