MILP Model for Division Property Analysis of Clyde-128

milp-clyde-dp.sage

from sage.crypto.boolean_function import BooleanFunction
from sage.crypto.sbox import SBox


def algebraic_normal_form(self):
    """
    Computes the algebraic normal forms (ANFs) of every coordinate.
    """
    n = self.input_size()
    return [self.component_function(i).algebraic_normal_form()
            for i in [1<<j for j in range(n)]]


SBox.algebraic_normal_form = algebraic_normal_form


def division_trail(self, k):
    """
    Computes the output division property for the starting input dp k.

    INPUT:

        - ``k``, the input division property
    """

    def gt(a, b):
        """
        check whether a >= b
        """
        from operator import ge
        return all(map(lambda x: ge(*x), zip(a, b)))

    n = self.input_size()

    S = set()
    for e in range(2^n):
        kbar = ZZ(e).digits(base=2, padto=n)
        if gt(kbar, k):
            S.add(tuple(kbar))

    ys = self.algebraic_normal_form()[::-1]
    P = ys[0].ring()
    x = P.gens()[::-1]

    F = set()
    for kbar in S:
        F.add(P(prod([x[i] for i in range(n) if kbar[i] == 1])))

    Kbar = set()
    for e in range(2^n):
        u = ZZ(e).digits(base=2, padto=n)
        puy = prod([ys[i] for i in range(n) if u[i] == 1])
        puyMon = P(puy).monomials()
        contains = False
        for mon in F:
            if mon in puyMon:
                contains = True
                break

        if contains:
            Kbar.add(tuple(u))

    K = []
    for kbar in Kbar:
        greater = False
        for kbar2 in Kbar:
            if(kbar != kbar2 and gt(kbar, kbar2)):
                greater = True
                break
        if not greater:
            K.append(kbar)

    return sorted(K)


SBox.division_trail = division_trail


def division_trail_table(self):
    """
    Return a dict containing all possible division propagation of the SBOX,
    where y is a list containing the ANF of each output bits
    """
    n = self.input_size()

    D = dict()
    for c in range(2^n):
        k = tuple(ZZ(c).digits(base=2, padto=n))
        D[k] = self.division_trail(k)

    return D

SBox.division_trail_table = division_trail_table

def sbox_inequalities(sbox, analysis="differential", algorithm="greedy", big_endian=False):
    """
    Computes inequalities for modeling the given S-box.

    INPUT:

    - ``sbox`` - the S-box to model
    - ``analysis`` - string, choosing between 'differential' and 'linear' cryptanalysis
        or 'division_property'
        (default: ``differential``)
    - ``algorithm`` - string, choosing the algorithm for computing the S-box model,
        one of ['none', 'greedy', 'milp'] (default: ``greedy``)
    - ``big_endian`` - representation of transitions vectors (default: little endian)
    """
    ch = convex_hull(sbox, analysis, big_endian)

    if algorithm is "greedy":
        return cutting_off_greedy(ch)
    elif algorithm is "milp":
        return cutting_off_milp(ch)
    elif algorithm is "none":
        return list(ch.inequalities())
    else:
        raise ValueError("algorithm (%s) has to be one of ['greedy', 'milp']" % (algorithm,))


SBox.milp_inequalities = sbox_inequalities


def convex_hull(sbox, analysis="differential", big_endian=False):
    """
    Computes the convex hull of the differential or linear behaviour of the given S-box.

    INPUT:

    - ``sbox`` - the S-box for which the convex hull should be computed
    - ``analysis`` - string choosing between differential and linear behaviour
      (default: ``differential``)
    - ``big_endian`` - representation of transitions vectors (default: little endian)
    """
    from sage.geometry.polyhedron.constructor import Polyhedron

    if analysis is "differential":
        valid_transformations_matrix = sbox.difference_distribution_table()
    elif analysis is "linear":
        valid_transformations_matrix = sbox.linear_approximation_table()
    elif analysis is "division_property":
        valid_transformations = sbox.division_trail_table()
    else:
        raise TypeError("analysis (%s) has to be one of ['differential', 'linear']" % (analysis,))

    if analysis is "division_property":
        points = [tuple(x) + tuple(y)
                  for x, ys in valid_transformations.iteritems() for y in ys]
    else:
        n, m = sbox.input_size(), sbox.output_size()

        if big_endian:
            def to_bits(x):
                return ZZ(x).digits(base=2, padto=sbox.n)
        else:
            def to_bits(x):
                return ZZ(x).digits(base=2, padto=sbox.n)[::-1]

        points = [to_bits(i) + to_bits(o)
                for i in range(1 << n)
                for o in range(1 << m)
                if valid_transformations_matrix[i][o] != 0]

    return Polyhedron(vertices=points)


def cutting_off_greedy(poly):
    """
    Computes a set of inequalities that is cutting-off equivalent to the
    H-representation of the given convex hull.

    INPUT:

    - ``poly`` - the polyhedron representing the convex hull
    """
    from sage.modules.free_module import VectorSpace
    from sage.modules.free_module_element import vector
    from sage.rings.finite_rings.finite_field_constructor import GF
    from sage.modules.free_module_element import vector

    chosen_ineqs = []

    poly_points = poly.integral_points()
    remaining_ineqs = list(poly.inequalities())
    impossible = [vector(poly.base_ring(), v)
                  for v in VectorSpace(GF(2), poly.ambient_dim())
                  if v not in poly_points]

    while impossible != []:

        if len(remaining_ineqs) == 0:
            raise ValueError("no more inequalities to choose, but still "\
                             "%d impossible points left" % len(impossible))

        # find inequality in remaining_ineqs that cuts off the most
        # impossible points and add this to the chosen_ineqs
        ineqs = []
        for i in remaining_ineqs:
            cnt = sum(map(lambda x: not(i.contains(x)), impossible))
            ineqs.append((cnt, i))
        chosen_ineqs.append(sorted(ineqs, reverse=True)[0][1])

        # remove ineq from remaining_ineqs
        remaining_ineqs.remove(chosen_ineqs[-1])

        # remove all cut off impossible points
        impossible = [v
                      for v in impossible
                      if chosen_ineqs[-1].contains(v)
                      ]
    return chosen_ineqs


def cutting_off_milp(poly, number_of_ineqs=None, **kwargs):
    """
    Computes a set of inequalities that is cutting-off equivalent to the
    H-representation of the given convex hull by solving a MILP.

    The representation can either be computed from the minimal number of
    necessary inequalities, or by a given number of inequalities. This
    second variant might be faster, because the MILP solver that later
    uses this representation can do some optimizations itself.

    INPUT:

    - ``poly`` - the polyhedron representing the convex hull
    - ``number_of_ineqs`` - integer; either `None` (default) or the number
        of inequalities that should be used for representing the S-box.

    REFERENCES:

    - [SasTod17]_ "New Algorithm for Modeling S-box in MILP Based
        Differential and Division Trail Search"
    """

    from sage.matrix.constructor import matrix
    from sage.modules.free_module import VectorSpace
    from sage.modules.free_module_element import vector
    from sage.numerical.mip import MixedIntegerLinearProgram
    from sage.rings.finite_rings.finite_field_constructor import GF

    ineqs = list(poly.inequalities())
    poly_points = poly.integral_points()
    impossible = [vector(poly.base_ring(), v)
                  for v in VectorSpace(GF(2), poly.ambient_dim())
                  if v not in poly_points]

    # precompute which inequality removes which impossible point
    precomputation = matrix(
        [[int(not(ineq.contains(p)))
          for p in impossible]
         for ineq in ineqs]
    )

    milp = MixedIntegerLinearProgram(maximization=False, **kwargs)
    var_ineqs = milp.new_variable(binary=True, name="ineqs")

    # either use the minimal number of inequalities for the representation
    if number_of_ineqs is None:
        milp.set_objective(sum([var_ineqs[i] for i in range(len(ineqs))]))
    # or the given number
    else:
        milp.add_constraint(sum(
            [var_ineqs[i]
             for i in range(len(ineqs))]
        ) == number_of_ineqs)

    nrows, ncols = precomputation.dimensions()
    for c in range(ncols):
        lhs = sum([var_ineqs[r]
                   for r in range(nrows)
                   if precomputation[r][c] == 1])
        milp.add_constraint(lhs >= 1)

    milp.solve()

    remaining_ineqs = [
        ineq
        for ineq, (var, val) in zip(ineqs, milp.get_values(var_ineqs).iteritems())
        if val == 1
    ]

    return remaining_ineqs


def milp_spook_sbox_constraints(milp, sbox, xi, yi):
    sbox_ineqs = sbox.milp_inequalities(analysis="division_property", algorithm="greedy")

    permuted_bits = matrix(ZZ, 4, 32, range(128)).columns()
    in_outs = [([xi[i] for i in sbox_indices], [yi[i] for i in sbox_indices])
               for sbox_indices in permuted_bits]

    for ineq in sbox_ineqs:
        for inputs, outputs in in_outs:
            milp.add_constraint(sum([inputs[i] * ineq[i+1]
                                     for i in range(len(inputs))] +
                                    [outputs[i] * ineq[i+1+len(inputs)]
                                     for i in range(len(outputs))]
                                   ) + ineq[0] >= 0)

def rotate(x, n):
    return x[n:] + x[:n]


def copy2(milp, x, y0, y1):
    for i in range(len(x)):
        milp.add_constraint(x[i] - y0[i] - y1[i] == 0)


def copy3(milp, x, y0, y1, y2):
    for i in range(len(x)):
        milp.add_constraint(x[i] - y0[i] - y1[i] - y2[i] == 0)


def xor2(milp, x0, x1, y):
    for i in range(len(x0)):
        milp.add_constraint(x0[i] + x1[i] - y[i] == 0)


def milp_spook_llayer_constraints(milp, xi, yi, ai, bi, rnd=0):
    s = milp.new_variable(binary=True, name="tmp_s")
    t = milp.new_variable(binary=True, name="tmp_t")
    u = milp.new_variable(binary=True, name="tmp_u")
    v = milp.new_variable(binary=True, name="tmp_v")

    s0 = [s[rnd, 0, i] for i in range(32)]
    s1 = [s[rnd, 1, i] for i in range(32)]
    s2 = [s[rnd, 2, i] for i in range(32)]
    s3 = [s[rnd, 3, i] for i in range(32)]
    s4 = [s[rnd, 4, i] for i in range(32)]
    s5 = [s[rnd, 5, i] for i in range(32)]
    s6 = [s[rnd, 6, i] for i in range(32)]
    s7 = [s[rnd, 7, i] for i in range(32)]

    copy3(milp, xi, s0, s1, s2)
    xor2(milp, s1, rotate(s2, 12), s3)
    copy2(milp, s3, s4, s5)
    xor2(milp, s4, rotate(s5, 3), s6)
    xor2(milp, s6, rotate(s0, 17), s7)

    t0 = [t[rnd, 0, i] for i in range(32)]
    t1 = [t[rnd, 1, i] for i in range(32)]
    t2 = [t[rnd, 2, i] for i in range(32)]
    t3 = [t[rnd, 3, i] for i in range(32)]
    t4 = [t[rnd, 4, i] for i in range(32)]
    t5 = [t[rnd, 5, i] for i in range(32)]
    t6 = [t[rnd, 6, i] for i in range(32)]
    t7 = [t[rnd, 7, i] for i in range(32)]

    copy3(milp, yi, t0, t1, t2)
    xor2(milp, t1, rotate(t2, 12), t3)
    copy2(milp, t3, t4, t5)
    xor2(milp, t4, rotate(t5, 3), t6)
    xor2(milp, t6, rotate(t0, 17), t7)

    s8 = [s[rnd, 8, i] for i in range(32)]
    s9 = [s[rnd, 9, i] for i in range(32)]
    u0 = [u[rnd, 0, i] for i in range(32)]
    u1 = [u[rnd, 1, i] for i in range(32)]
    u2 = [u[rnd, 2, i] for i in range(32)]
    u3 = [u[rnd, 3, i] for i in range(32)]
    u4 = [u[rnd, 4, i] for i in range(32)]

    copy3(milp, s7, s8, u0, u1)
    xor2(milp, rotate(u0, 31), u1, u2)
    copy2(milp, u2, u3, u4)
    xor2(milp, s8, rotate(u3, 15), s9)

    t8 = [t[rnd, 8, i] for i in range(32)]
    t9 = [t[rnd, 9, i] for i in range(32)]
    v0 = [v[rnd, 0, i] for i in range(32)]
    v1 = [v[rnd, 1, i] for i in range(32)]
    v2 = [v[rnd, 2, i] for i in range(32)]
    v3 = [v[rnd, 3, i] for i in range(32)]
    v4 = [v[rnd, 4, i] for i in range(32)]

    copy3(milp, t7, t8, v0, v1)
    xor2(milp, rotate(v0, 31), v1, v2)
    copy2(milp, v2, v3, v4)
    xor2(milp, t8, rotate(v3, 15), t9)

    xor2(milp, s9, rotate(v4, 26), ai)
    xor2(milp, t9, rotate(u4, 25), bi)


def milp_model_spook(initial_dp, rnds=1):
    sbox = SBox([0, 8, 1, 15, 2, 10, 7, 9, 4, 13, 5, 6, 14, 3, 11, 12])

    from itertools import product

    # initialise MILP object
    milp = MixedIntegerLinearProgram(maximization=False, solver="CPLEX")

    # sbox layer inputs
    xs = milp.new_variable(binary=True, name="x", indices=product(range(rnds), range(128)))
    # sbox layer outputs / linear layer inputs
    ys = milp.new_variable(binary=True, name="y", indices=product(range(rnds), range(128)))
    # linear layer outputs
    zs = milp.new_variable(binary=True, name="z", indices=product(range(rnds), range(128)))

    # model for each round the sbox layer and linear layer transitions
    for r in range(rnds):
        xi = [xs[(r, i)] for i in range(128)]
        yi = [ys[(r, i)] for i in range(128)]
        milp_spook_sbox_constraints(milp, sbox, xi, yi)

        yi = [[ys[(r, 32*j+i)] for i in range(32)] for j in range(4)]
        zi = [[zs[(r, 32*j+i)] for i in range(32)] for j in range(4)]
        milp_spook_llayer_constraints(milp, yi[0], yi[1], zi[0], zi[1], rnd=(r, 0))
        milp_spook_llayer_constraints(milp, yi[2], yi[3], zi[2], zi[3], rnd=(r, 1))

        # link each rounds output with next rounds input
        if r < rnds-1:
            for i in range(128):
                milp.add_constraint(zs[(r, i)] == xs[(r+1, i)])

    # Set input variables to initial division property
    from sage.crypto.sbox import integer_types
    if type(initial_dp) in integer_types + (Integer, ):
        initial_dp = ZZ(initial_dp).digits(base=2, padto=128)
    for i in range(128):
        milp.add_constraint(xs[(0, i)] == initial_dp[i])

    # Objective function is to minimize the weight of the output division property
    milp.set_objective(sum(zs[(rnds-1, i)] for i in range(128)))

    return milp, xs, ys, zs


def check_dp(rnds=1, milp_model=milp_model_spook, block_size=128):
    from sage.numerical.mip import MIPSolverException

    for i in range(block_size):
        k = ((1<<block_size) - 1)^^(1<<i)
        milp, xs, ys, zs = milp_model(initial_dp=k, rnds=rnds)

        cnt = 0
        found_unit_vector = True
        while found_unit_vector:
            try:
                obj = int(milp.solve())
                inp = [int(x)
                       for x in milp.get_values([xs[(   0  , j)]
                                                 for j in range(block_size)])]
                out = [int(x)
                       for x in milp.get_values([zs[(rnds-1, j)]
                                                 for j in range(block_size)])]

                inpstr = "".join(map(lambda x:"%d" % x, inp))
                outstr = "".join(map(lambda x:"%d" % x, out))

                cnt += 1
                print("%3d/%3d: %3d %s -> %s" % (i, cnt, obj, inpstr, outstr))

                if obj > 1:
                    print("found a distinguisher:")
                    print("%3d: %3d %s -> %s" % (i, obj, inpstr, outstr))
                    return inp, out
                else:
                    idx = out.index(1)
                    milp.add_constraint(zs[(rnds-1, idx)] == 0)

            except MIPSolverException as e:
                print("i = %d: no feasible solution" % i)
                found_unit_vector = False

    return None, None


if __name__ == "__main__":
    import sys
    if len(sys.argv) < 2:
        print("Usage:\n%s rounds" % (sys.argv[0]))
        sys.exit(1)

    rounds = int(sys.argv[1])
    check_dp(rounds, milp_model_spook, block_size=128)