(* 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]];