The Function Br

Function to find all nonempty SRS periods within a closed area Q⊂D2.
The function may return more periods. These additional ones are nonempty too, but the polyhedron do not intersect with Q.
Usage:
Br[Q]
Q....description of the set Q by using x and y in logical expressions.
For numberical expressions only use integers and rationals. Do not use strict inequalities.
The function returns a list of  periods that yield nonempty polyhedra. The list may contain periods, whose cutouts do not intersect with Q.
Examples:
Br[x≥5/6∧x≤99/100∧y≤x-1∧y≥-x/2]
Here Q equals the area {(x,y) | (5/6) ≤ x ≤ (99/100), -(x/2) ≤ y ≤ x - 1}.
Br[{x,y}=={9/10,17/10}]  
Here Q equals only the point ((9/10),(17/10)). When the input is just one point, the usage of any real nuber is allowed but works more slowly.
Additionally the list of returned periods can include more periods with polyhedra not containing the point: this example returns 2 periods while the same input using 0.9 an 1.7 instead gives 3 periods (an additional one whose polygon is "near" the point (0.9,1.7)).

<< DiscreteMath`Combinatorica`; 
$RecursionLimit = Infinity;
\(\(\(Off[General::spell1];\)\(
\)
    \)\)
\(\(SetOfWit = 
        Function[{R, bound}, 
          Block[\ {S, S1, S2, GenerateList, Tau, TauR, 
              ma}, 

 (*Generates\ the\ Set\ of\ Witnesses\ for\ a\ Set\ R\ given\ as\ System\ of\ \((In)\) 
                equalities*) 

 (*LeftShift\  + \ #2\ on\ the\ last\ Position*) 
Tau = \(({Apply[\((##2)\) &, #1], #2})\) &; 

 (*\ z\  -> \ {z, \(-Floor\) \((Min \((rz)\))\), \(-Floor\) \((Max \((rz)\))\), Floor \((Min \((\(-rz\))\))\)}, \ 
              r\ \in\ R*) 
GenerateList = 
              Function[
                z, 
ma = \(Maximize[{x, y} . z, 
                      R, {x, y}]\)[\([1]\)]; 
{z, \(-Floor[\(Minimize[{x, y} . z, R, {x, y}]\)[\([1]\)]]\), \(-Floor[ma]\), \(-Ceiling[
                      ma]\)}]; 

 (*Generates\ Tau_ {r} \((z)\)\ and\  - 
                  Tau_ {r} \((\(-z\))\)\ for\ all\ r\  ∈ \ 
                R*) 
TauR = 
              Function[z, 
                Table[GenerateList[Tau[z[\([1]\)], i]], {i, z[\([4]\)], 
                    z[\([2]\)]}]]; 

 (*Initialization*) 
S = 
              Map[GenerateList, 
                IdentityMatrix[
                    2]\ \  ⋃ \ \(-IdentityMatrix[
                      2]\)]; 
S1 = {}; 

 (*S_ {i + 1} = 
                TauR \((S_ {i}\S_ {i - 1})\)\  ⋃ \ 
                  S_ {i - 
                        1}\ up\ to\ stabilization\ or\ overflow*) 
While[\((S\  ≠ \ 
                    S1)\)\  ∧ \((Length[S]\  ≤ 
                    bound)\), 
S2 = 
                Complement[S, S1]; 
S1 = 
                S; 
S = 
                S ⋃ 
                  Flatten[TauR /@ S2, 
                    1]]; 

If\ [
              Length[S]\  ≤ \ bound, 
              S, {}]]];\)
\)
\(\(BuildGraph = 
        Function[S, 
          Block[{mx, my, Vx, AL, AdjazenzList, 
              gr}, 

 (*Generates\ the\ Graph\ for\ a\ set\ \((using\ the\ Combinatorica - 
                    package)\)\ *) 

 (*The\ List\ of\ Edges\ is\ Created\ without\ using\ the\ function\ "Position"\ \((very\ slow)\)*) 
AdjazenzList = 
              Function[
                z, 
Table[{{\((z[\([1, 1]\)] + 
                              mx)\)*\((2*my + 1)\) + z[\([1, 2]\)] + my + 
                        1, \((z[\([1, 2]\)] + mx)\)*\((2*my + 1)\) + my + 1 + 
                        i}}, {i, z[\([3]\)], 
                    z[\([2]\)]}]]; \ 
mx = 
              Max @@ \((Abs /@ \((First /@ \((First /@ 
                            S)\))\))\); 
my = 
              Max @@ \((Abs /@ \((Last /@ \((First /@ 
                            S)\))\))\); 
Vx = 
              Flatten[Table[{{j, i}}, {j, \(-mx\), mx}, {i, \(-my\), my}], 
                1]; 
AL = 
              Join @@ \((AdjazenzList /@ S)\); 
gr = 
              Graph[AL, Vx, 
                EdgeDirection -> 
                  True]; 
InduceSubgraph[gr, 
              Union[First /@ Edges[gr]]]]];\)
\)
CheckSeq = 
        Function[Seq, 
          Block[{UGS, 
              l}, 

  

l = 
              Length[Seq]; 
UGS = (x < 
                    1) ∧ (And @@ 
                    Table[1 > 
                        x*Seq[[i, 1]] + y*Seq[[i, 2]] + 
                          Seq[[i + 1, 2]] ≥ 0, {i, 1, 
                        l - 1}]); 
FindInstance[
                UGS, {x, y}] ≠ {}]];

CheckCyc = 
        Function[Cy, 
          Block[{UGS, 
              l}, 

 (*Returns TRUE if the given List of Integers induces a nonempty cutout*) 

l = 
              Length[Cy]; 
UGS = (x < 1) ∧ (And @@ 
                    Table[1 > 
                        x*Cy[[i + 1]] + y*Cy[[Mod[i + 1, l] + 1]] + 
                          Cy[[Mod[i + 2, l] + 1]] ≥ 0, {i, 
                        0, 
                        l - 1}]); 
Cy == {} ∨ (FindInstance[UGS, {x, 
                      y}] ≠ {})]];

NormCyc = 
        Function[cy, 
          Block[{l, 
              Li}, 

  

l = 
              Length[cy]; 
Li = 
              Table[RotateLeft[cy, i], {i, 0, 
                  l - 1}]; 
Sort[
                Li][[1]]]];

FindAllCycles = 
        Function[g, 
          Block[{IL, FC, 
              FAC}, 

 (*Returns the list of all cycles within the graph g that give nonemty cutouts*) 
 (*The cycles are rotated such that the leftmost element is the lexicographical lowest*) 

 (*1.  way :  
                  Recursive calculation of all cycles .  
                      Test afterwards*) 
FC = 
              Function[{$ModuleNumber, i, SC}, 
                Block[{erg}, 
If[
                    MemberQ[SC, i], 
NormCyc /@ 
                      Select[{First /@ 
                            Part[Vertices[g], 
                              Take[SC, 
                                Append[Flatten[Position[SC, i]], -1]]]}, 
                        CheckCyc], 
Union @@ (FC[ $ModuleNumber + 1, #, Append[SC, 
                                i]] & /@   IL[[i]])]]]; 

 (*2.  way :  
                  Recursive calculation of all cycles .  
                      A path is not further persecuted,  
              if it cannot be part of a valid cycle*) 
FAC = 
              Function[{$ModuleNumber, i, SC}, 
If[
                  CheckSeq[
                    Part[Vertices[g], 
                      Append[SC, i]]], 
If [
                    MemberQ[SC, i], 
If [
                      SC[[1]] == 
                        i, {NormCyc[
                          First /@ 
                            Part[Vertices[g], 
                              SC]]}, {}], 
Union @@ (FAC[ $ModuleNumber + 1, #, 
                              Append[SC, 
                                i]] & /@   IL[[i]])], 
{}]]; 

IL = 
              ToAdjacencyLists[
                g]; 
 (*#  cycles  ≤ 
                  2^10  ->  1.  way,  
              otherwise  ->  2.  way*) 
If[
              cc == {} ∨ 
                Times @@ (Length /@ IL) ≤ 
                  2^10 , 
FC[$ModuleNumber + 1, 
                1, {}], 
Union @@ (FAC[$ModuleNumber + 
                          1, #, {}] & /@ 
                    Table[i, {i, 1, 
                        V[g]}])]
]];

RedInq = 
        Function[inq, 
          Block[{in, i, 
              l}, 

 (*Removes redundant Inequalities of a system (given as List or linked with And)*) 
 (*In this context it gives clearer results as the built in command Reduce*) 

in = 
              List @@ inq; 
i = 1; 
l = 
              Length[in]; 
While[
              i ≤ l, 
in = 
                Complement[
                  in, {(List @@ inq)[[i]]}]; 
If[
                FindInstance[(And @@ in) ∧ 
                      Not[(List @@ inq)[[i]]], {x, 
                      y}] ≠ {}, 
in = 
                  Append[in, (List @@ 
                          inq)[[i]]]]; 
i++]; 
Head[inq] @@ 
              in
]];

\(\(Br = 
        Function[R, 
          Block[{Appr, 
              per}, 

 (*Calculated\ all\ cycles\ that\ give\ nonemty\ cutouts\ within\ a\ closed\ set\ R*) 

 (*Main\ Routine, \ 
              Recursive*) 
Appr = 
              Function[{$ModuleNumber, r}, 
                Block[{S, g, xm, ym, dx, dy, r1, rz, cc, bound, 
                    res}, 

 (*minimal\ x\ and\ y\ values\ of\ R*) \ 
xm = \(Minimize[x, 
                        r, {x, 
                          y}]\)[\([1]\)]; 
ym = \(Minimize[
                        y, r, {x, 
                          y}]\)[\([1]\)]; 
 (*expanse\ in\ x\ and\ y\ direction*) \ 
dx = \(Maximize[x, 
                          r, {x, y}]\)[\([1]\)] - 
                      xm; 
dy = \(Maximize[y, 
                          r, {x, y}]\)[\([1]\)] - 
                      ym; 
 (*\(-- -> \) >>> \ 
                      calculation\ of\ the\ bound\ for\ the\ set\ of\ \(\(\(witnesses\ !!\)!!\)!!\)*) 
If[
                    Max[dx, dy] > 0, 
bound = 
                      IntegerPart[
                        100/Max[dx, dy]], 
bound = 
                      Infinity]; 
Print["Maximal size of the set of witnesses: ", 
                    bound, " elements"]; 

 (*Calculation\ of\ the\ set\ of\ witnesses*) 
S = SetOfWit[r, bound]; 
If[
                    S ≠ {}, 

 (*go\ on\ if\ the\ set\ of\ witnesses\ is\ small\ enough ... *) 
g = 
                      BuildGraph[
                        S]; 
cc = \(InduceSubgraph[
                            g, #] &\) /@ \((Select[
                            StronglyConnectedComponents[g], 
                            Length[#] ≥ 
                                2 &])\); 
res = 
                      Union @@ \((FindAllCycles /@ 
                            cc)\); 
per = 
                      per + 1/
                          2^\(($ModuleNumber - 1)\); 
Print[
                      per*100\  // 
                        N, " % of the area investigated"]; 
res, 

 (*\(...\) \(othwerwise\)\(\ \)\(split\)\(\ \)\(the\)\(\ \)\(area\)\(\ \)\(and\)\(\ \)\(apply\)\(\ \)\(Appr\)\(\ \)\(on\)\(\ \)\(each\)\(\ \)\(of\)\(\ \)\(them\)*) 
If[
                      dx > dy, 
rz = 
                        Rationalize[xm + dx/2, 
                          dx/100]; 
r1 = {x ≤ 
                            rz, x ≥ 
                            rz}, 
rz = 
                        Rationalize[ym + dy/2, 
                          dy/100]; 
r1 = {y ≤ 
                            rz, y ≥ 
                            rz}]; 
Union @@ \((\(Appr[$ModuleNumber + 1, RedInq[r\  ∧ \ #]] &\) /@ 
                          r1)\)]]]; 

per = 0; 
x =. ; 
            y =. ; 
If[
              FindInstance[R, {x, y}] ≠ \ {}, 
              Appr[1, 
                R], "Expression equals the empty set"]]];\)\)