The Function FP

Function to find all nonempty SSRS periods within a closed area Q⊂D3. 
The function may return more periods. These additional ones are nonempty too, but the polyhedron do not intersect with Q.
Usage:
FP[Q]
Q....description of the set Q by using x,y and z 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:
FP[x-z≤-1∧x-y+z≤1/2∧-x+y≤1/2∧x>-1/2]
Here Q equals the area {(x,y,z) |x - z ≤ -1, x - y + z ≤ (1/2), -x + y ≤ (1/2), x > -(1/2)}.
FP[{x,y,z}=={1/2,1/2,-1/2}]  
Here Q equals only the point ((1/2),(1/2),-(1/2)). When the input is just one point, the usage of any real nuber is allowed but works more slowly.

<< DiscreteMath`Combinatorica`; 
$RecursionLimit = Infinity;
Off[General::spell1];

\(\(SetOfWit = 
        Function[{r, bound}, 
          Block[\ {S, S1, S2, GenerateList, Tau, 
              TauR}, 

 (*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})\) &; 

 (*\(\(\(\ \)\(w\  -> \ {w, Min \((rw)\), Max \((rw)\)\)\()\)\)\(}\)\), \ 
              r\  ∈ \ R*) 
GenerateList = 
              Function[
                w, 
{w, \(Minimize[{x, y, z} . w, 
                      r, {x, y, z}]\)[\([1]\)], \(Maximize[{x, y, z} . w, 
                      r, {x, y, 
                        z}]\)[\([1]\)]}]; 

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

 (*Initialization*) 
S = 
              Map[GenerateList, 
                IdentityMatrix[
                    3]\ \  ⋃ \ \(-IdentityMatrix[
                      3]\)]; 
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, \({#[\([1]\)], \(-Floor[#[\([2]\)] + 
                          1/2]\), \(-Floor[#[\([3]\)] + 1/2]\)} &\) /@ 
                S, {}]]];\)
\)
\(\(BuildGraph = 
        Function[S, 
          Block[{mx, my, mz, Vx, AL, AdjazenzList, 
              gr}, 

 (*Generates\ the\ Graph\ for\ a\ set\ \((using\ the\ Combinatorica - 
                    package)\)\ *) 
 (*A\ Mathematica\ graph\ \((Combinatorica - package)\)\ works\ with\ two - 
                dimensional\ \(\(points\)\(.\)\)*) 
 (*Therefore\ the\ output\ of\ this\ function\ is\ a\ pair\ {g, 
                  V}\ where\ V\ is\ a\ set\ of\ three\ dimensional\ points\ and\ g\ is\ a\ graph\ with\ vertices\ of\ the\ form\ {1, 
                  k}\ such\ that\ V[\([k]\)]\ is\ the\ corresponding\ \(\(point\)\(.\)\)*) \ 

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

  

l = 
              Length[Seq]; 
UGS = 
              And @@ Table[
                  1/2 > x*Seq[[i, 1]] + y*Seq[[i, 2]] + 
                      z*Seq[[i, 3]] + 
                      Seq[[i + 1, 3]] ≥ -1/2, {i, 1, 
                    l - 1}]; 
FindInstance[
                UGS, {x, y, 
                  z}] ≠ {}]];

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

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

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


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

 (*chooses\ one\ representative\ out\ of\ the\ set\ of\ all\ \((equivalent)\)\ rotations\ of\ a\ period*) 

l = 
              Length[cy]; 
Li = 
              Table[RotateLeft[cy, i], {i, 0, 
                  l - 1}]; 
\(Sort[
                Li]\)[\([1]\)]]];\)
\)
FindAllCycles = 
        Function[{g, Vxx}, 
          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[Vxx, 
                              Last /@ 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[Vxx, 
                      Last /@ Part[Vertices[g], 
                          Append[SC, i]]]], 
If [
                    MemberQ[SC, i], 
If [
                      SC[[1]] == 
                        i, {NormCyc[
                          Last /@ Part[Vxx, 
                              Last /@ 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[
              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, 
                      z}] ≠ {}, 
in = 
                  Append[in, (List @@ 
                          inq)[[i]]]]; 
i++]; 
Head[inq] @@ 
              in
]];

\(\(FP = 
        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, zm, dx, dy, dz, r1, rz, cc, Vxx, res, 
                    bound}, 

 (*minimal\ x, y\ and\ z\ values\ of\ R*) \ 
xm = \(Minimize[x, 
                        r, {x, y, 
                          z}]\)[\([1]\)]; 
ym = \(Minimize[
                        y, r, {x, y, 
                          z}]\)[\([1]\)]; 
zm = \(Minimize[
                        z, r, {x, y, 
                          z}]\)[\([1]\)]; 

 (*expanse\ in\ x, \ 
                    y\ and\ z\ direction*) \ 
dx = \(Maximize[x, r, {x, y, z}]\)[\([1]\)] - 
                      xm; 
dy = \(Maximize[y, 
                          r, {x, y, z}]\)[\([1]\)] - 
                      ym; 
dz = \(Maximize[z, 
                          r, {x, y, z}]\)[\([1]\)] - 
                      zm; 

 (*\(-- -> \) >>> \ 
                      calculation\ of\ the\ bound\ for\ the\ set\ of\ \(\(\(witnesses\ !!\)!!\)!!\)*) 
If[
                    Max[dx, dy, dz] > 0, 
bound = 
                      IntegerPart[
                        20/Max[dx, dy, dz]], 
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, Vxx} = 
                      BuildGraph[
                        S]; 
cc = \(InduceSubgraph[
                            g, #] &\) /@ \((Select[
                            StronglyConnectedComponents[g], 
                            Length[#] ≥ 
                                2 &])\); 
res = 
                      Union @@ \((\(FindAllCycles[#, Vxx] &\) /@ 
                            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\  ∧ \ 
                        dx ≥ dz, 
rz = 
                        Rationalize[xm + dx/2, 
                          dx/100]; 
r1 = {x ≤ 
                            rz, x ≥ rz}, 
If[
                        dy ≥ dz, 
rz = 
                          Rationalize[ym + dy/2, 
                            dy/100]; 
r1 = {y ≤ 
                              rz, y ≥ 
                              rz}, 
rz = 
                          Rationalize[zm + dz/2, 
                            dz/100]; 
r1 = {z ≤ 
                              rz, z ≥ 
                              rz}]]; 
Union @@ \((\(Appr[$ModuleNumber + 1, RedInq[r\  ∧ \ #]] &\) /@ 
                          r1)\)]]]; 

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