(*Minimal redundant expansions in Gaussian Integers*)

BeginPackage["MinimalRedundantGauss`"]

CNS::usage="CNS[d_0, ..., d_l] represents a digit expansion with \
respect to some basis."

ExpandCNS::usage ="ExpandCNS[beta][d_0,...,d_l] gives Sum[d_j beta^j, \
{j,0,l}]."

StandardCNS::usage="StandardCNS[alpha, beta] gives the standard \
expansion of alpha in the canonical number system defined by beta."

MinimumCNS::usage="MinimumCNS[alpha,beta] gives all minimal expansions \
of alpha in base beta with respect to the costs (4). Previous results \
are stored in order to improve running time. This may consume much memory."

ShortestCNS::usage="ShortestCNS[alpha,beta] gives a shortest expansion \
of alpha in base beta using digits of absolute value < Norm[beta]. Previous results \
are stored in order to improve running time. This may consume much memory."

FindRule::usage="FindRule[CNS[d_0,...,d_l],beta] finds the set of r_0 \
such that the conditions of Lemma 4 are fulfilled."

Begin["`Private`"]


(*Utilities for complex numbers*)

norm[z_Complex] := Re[z]^2 + Im[z]^2


modbase[z_?NumberQ, beta_?NumberQ] := 
  With[{c = Re[z], d = Im[z], a = -Re[beta]}, Mod[c + a*d, norm[beta]]]


(*Utilities*)


ExpandCNS[beta_?NumberQ][u__] := Module[{i, l = {u}},
    Sum[l[[i + 1]]*beta^i, {i, 0, Length[l] - 1}]]

cost[CNS[u___]] := Module[{l = {u}},
    Plus @@ Map[Abs, l]
    ]

(*Bounds for the Digits *)

possibleDigits[r_Integer, beta_?NumberQ] := 
  Module[{n = -Re[beta], digitbound, i},
    digitbound = Which[n >= 3, n^2, n == 2, 8, n == 1, 2];
    Select[Table[i, {i, -digitbound, +digitbound}], 
      Mod[#, n^2 + 1] == Mod[r, n^2 + 1] &]
    ]

(*Symbolic Manipulation of CNS*)

CNS /: Dot[u_CNS, v_CNS] := Join[u, v]
CNS /: Power[u_CNS, 0] := CNS[]
CNS /: Power[u_CNS, 1] := u
CNS /: Power[u_CNS, n_Integer /; n > 1] := 
  u^Floor[n/2].u^Floor[n/2].u^Mod[n, 2]

(*Expansion in the standard Canonical Number System*)


StandardCNS[0, beta_Complex] = CNS[];
StandardCNS[z_?NumberQ, beta_?NumberQ] := Module[{r, z1},
    r = modbase[z, beta];
    z1 = (z - r)/beta;
    Prepend[StandardCNS[z1, beta], r]
    ]


(* Minimal Expansion *)


MinimumCNS[0, ___] = {CNS[]};
MinimumCNS[z_?NumberQ, beta_?NumberQ, first_?IntegerQ, maxweight_?IntegerQ] :=
   Module[{expansions},
    expansions = MinimumCNS[(z - first)/beta, beta, maxweight - Abs[first]];
    expansions = Map[Prepend[#, first] &, expansions];
    expansions
    ]

MinimumCNS[z_?NumberQ, beta_?NumberQ, maxweight_?IntegerQ] := 
    Module[{r, rlist, expansions, mincost, lookup},
      If[minimumCNSTableBeta ===  beta,
        lookup = minimumCNSTable[z],
        (
          Clear[minimumCNSTable];
          minimumCNSTableBeta = beta;
          )
        ];
      If[Head[lookup] === List,
        If[cost[lookup[[1]]] <= maxweight, lookup, {}],
        ( r = modbase[z, beta];
          rlist = possibleDigits[r, beta];
          rlist = Select[rlist, Abs[#] <= maxweight &];
          
          expansions = 
            Flatten[Map[MinimumCNS[z, beta, #, maxweight] &, rlist], 1];
          mincost = Min[Map[cost, expansions]];
          expansions = Select[expansions, cost[#] <= mincost &];
          If[Length[expansions] >= 1,
            minimumCNSTable[z] = expansions
            ];
          expansions
          )
        ]
      ];

MinimumCNS[z_?NumberQ, beta_?NumberQ] := Module[{standardcost, result},
    standardcost = cost[StandardCNS[z, beta]];
    result = MinimumCNS[z, beta, standardcost];
    If[{z} == Union[result /. CNS :> ExpandCNS[beta]] && 
        Length[Union[Map[cost, result]]] == 1, result, False]
    ]



(*Shortest Expansion: Minimize Length of Expansion, digits = - \
 Norm[beta]+1, ... Norm[beta]-1*)

ShortestCNS[0, ___] = {CNS[]};

ShortestCNS[z_?NumberQ /; z != 0, beta_?NumberQ, 0] = {};

ShortestCNS[z_?NumberQ, beta_?NumberQ,  first_?IntegerQ, 
    maxlength_?IntegerQ] := Module[{expansions},
    expansions = ShortestCNS[(z - first)/beta, beta, maxlength - 1];
    expansions = Map[Prepend[#, first] &, expansions];
    expansions
    ]

ShortestCNS[z_?NumberQ, beta_?NumberQ, maxlength_?IntegerQ] := 
  Module[{r, rlist, expansions, minlength, lookup},
    If[shortestCNSTableBeta == beta,
      lookup = shortestCNSTable[z],
      (
        Clear[shortestCNSTable];
        shortestCNSTableBeta = beta;
        )
      ];
    If[Head[lookup] === List,
      If[Length[lookup[[1]]] <= maxlength, lookup, {}],
      (
        r = modbase[z, beta];
        rlist = If[r==0,{r},{r,r-norm[beta]}];
        expansions = Flatten[Map[ShortestCNS[z, beta, #, maxlength] &, rlist], 1];
        minlength = Min[Map[Length, expansions]];
        expansions = Select[expansions, Length[#] <= minlength &];
        expansions = If[Length[expansions] >= 1, {expansions[[1]]}, {}];
        If[Length[expansions] >= 1,
          shortestCNSTable[z] = expansions
          ];
        expansions
        )
      ]
    ]

ShortestCNS[z_?NumberQ, beta_?NumberQ] := Module[{standardlength, result},
    standardlength = Length[StandardCNS[z, beta]];
    result = ShortestCNS[z, beta, standardlength];
    If[Length[result] == 1 && (result /. CNS :> ExpandCNS[beta]) == {z} && 
        Length[result] <= standardlength, result[[1]], False
       ]
]



listCirclePoints[bound_?NumberQ] := Module[{a, b, result},
    result = 
      Flatten[Table[
          Union[{a + b*I, a - b*I, -a + b*I, -a - b*I}], {a, 0, 
            Floor[Sqrt[bound]]}, {b, 0, Floor[Sqrt[bound - a^2]]}]]
    ]

listRepresenters[u_CNS, beta_?NumberQ] := 
  Module[{z, bound, r, k, possibleOffsets, candidates},
    r = Length[u];
    z = u /. CNS -> ExpandCNS[beta];
    bound = 
      Floor[(2*(norm[beta] - 1)/Abs[beta]^r Sum[
                Abs[beta^k], {k, 0, r - 1}])^2];
    possibleOffsets = listCirclePoints[bound];
    candidates = Map[z + beta^r*# &, possibleOffsets];
    Select[candidates, Length[ShortestCNS[#, beta]] <= r &]
    ]



FindRule[u_CNS, beta_?NumberQ] := Module[{repr, expansions, initials},
    repr = listRepresenters[u, beta];
    expansions = Map[MinimumCNS[#, beta] &, repr] /. CNS[] -> CNS[0];
    initials = Map[Union, expansions /. CNS[a_, ___] -> a];
    Intersection @@ initials
    ]

0

End[]
EndPackage[]