r"""
Polynomial system generator for the 'Data Encryption Standard' (DES).

The Data Encryption Standard (DES) is a cipher selected as an official
Federal Information Processing Standard (FIPS) for the United States
in 1976, and which has subsequently enjoyed widespread use
internationally. DES came under intense academic scrutiny, and
motivated the modern understanding of block ciphers and their
cryptanalysis.

The specification of DES can be found at

   \url{http://www.itl.nist.gov/fipspubs/fip46-2.htm}

One set of equations for the S-Boxes S1,...,S8 are based on a bitslice
DES C implementation by Matthew Kwan available at:

   \url{http://www.darkside.com.au/bitslice/nonstd.c}

These low gate count S-Box representations are described and discussed
in:

    Matthew Kwan; Reducing the Gate Count of Bitslice DES; Cryptology
    ePrint Archive, Report 2000/051; 2000;
    \url{http://eprint.iacr.org/}

If the \code{DES} constructor sees the parameter \code{sbox_eq='opns'}
then exactly these equations are used. If \code{sbox_eq='opns_gb'}
which is the default then a 'degrevlex' Gr\"obner basis is
pre-computed for each S-Box in the boolean ring

  $$GF(2)[y_0,...,y_4,x_0,...,x_5,t_0,...,t_{52}]/I$$

where $I$ is the 'field ideal' spanned by polynomials of the form
$x_i^2 + x_i$ (c.f. \code{sage.rings.ideal.FieldIdeal}). This
pre-computation is executed once for each Sage session when the first
DES polynomial system is constructored with \code{sbox_eq='opns_gb'}.

The constructor also supports \code{sbox_eq='cubic'} which means that
fully cubic equations with no intermediate variables are used to
represent each S-Box. For Groebner basis computation these equations
are probably the best option.

More equation systems for the S-Boxes S1,...,S8 can be found in

  Nicolas T. Courtois and Gregory V. Bard; Algebraic Cryptanalysis of
  the Data Encryption Standard; Cryptology ePrint Archive, Report
  2006/402; 2006; \url{http://eprint.iacr.org/}

Note that selecting an optimal algebraic representation for the DES
S-Boxes is an open research problem.

Once an algebraic representation for each S-Box S1,...,S8 is found
constructing a full DES is straight forward, because the cipher does
not contain any other non-linear components.

EXAMPLE:

    We encrypt the plaintext $M$ under key $K$, this example is from:

       \url{http://www.aci.net/kalliste/des.htm}

    sage: execfile('des.py')
    sage: des = DES()
    sage: M = '0000000100100011010001010110011110001001101010111100110111101111'
    sage: K = '0001001100110100010101110111100110011011101111001101111111110001'
    sage: M = [int(b) for b in M]
    sage: K = [int(b) for b in K]
    sage: ''.join(map(str, des(M,K)))
    '1000010111101000000100110101010000001111000010101011010000000101'

AUTHOR: 
    Martin Albrecht (2008-09) initial implementation
"""
import re
import sage.crypto.mq as mq

from sage.crypto.mq.mpolynomialsystemgenerator import MPolynomialSystemGenerator
from sage.misc.misc import union
from sage.rings.integer_ring import ZZ
from sage.structure.sage_object import SageObject

class DES(MPolynomialSystemGenerator):
    def __init__(self, Nr=16, **kwds):
        """
        Return a new polynomial system generator for the 'Data
        Encryption Standard' (DES). See \code{sage.crypto.mq.des} for
        a brief description of the S-Box representation.

        INPUT:
            Nr -- number of rounds $1 <= Nr <= 16$ (default: 16)
            order -- a string to specify the term ordering of the
                     variables (default: 'degrevlex')
            postfix -- a string which is appended after the variable
                       name except for key variables (default: '')
            polybori -- use PolyBoRi as implementation for the
                        polynomials (probably deprecated eventually,
                        since PolyBoRi will the standard, default: True)
            sbox_eq -- string; algebraic representation of the S-Boxes
                       (default: 'opns_gb')

        S-BOX REPRESENTATIONS:
            cubic -- fully cubic equations
            opns -- sparse quadratic equations with new intermediate
                    variables
            opns_gb -- same as above but the Groebner basis for each
                       S-Box is precomputed
                       
        EXAMPLE:
            sage: execfile('des.py')
            sage: des = DES()
            sage: des
            DES(Nr=16)

        REFERENCES:
              FIPS Publication 46-1: Data Encryption standard, NIST,
              Washington, D.C., Jan. 22, 1988; originally issued by
              the National Bureau of Standards.
              \url{http://www.itl.nist.gov/fipspubs/fip46-2.htm}
              
              Nicolas T. Courtois and Gregory V. Bard; Algebraic
              Cryptanalysis of the Data Encryption Standard;
              Cryptology ePrint Archive, Report 2006/402; 2006;
              \url{http://eprint.iacr.org/}
        """
        self._params = {}
        self._params["Nr"] = Nr
        self._params["postfix"] = kwds.get("postfix","")
        self._params["polybori"] = kwds.get("polybori", True)
        self._params["order"] = kwds.get("order","degrevlex")
        self._params["sbox_eq"] = kwds.get("sbox_eq", "opns_gb")

        for name,value in self._params.iteritems():
            setattr(self, "_" + name, value)

        self.Bs = 64
        self._base = GF(2)
        self._S = [DESSBox(i) for i in range(1,9)]
        self._tlength = 376 # the number of additional intermediate
                            # variables t in each round

    def new_generator(self, **kwds):
        """
        Return a new polynomial system generator for DES matching this
        generator except for the keywords provided.

        INPUT:
            **kwds -- see constructor for allowed keywords

        EXAMPLE:
            sage: execfile('des.py')
            sage: des = DES(Nr=4, polybori=True)
            sage: type(des.ring())
            <type 'sage.rings.polynomial.pbori.BooleanPolynomialRing'>

            sage: des = des.new_generator(polybori=False)
            sage: type(des.ring())
            <type 'sage.rings.polynomial.multi_polynomial_libsingular.MPolynomialRing_libsingular'>
        """
        params = {}
        
        for name, value in self._params.iteritems():
            params[name] = kwds.get(name,value)

        return self.__class__(**params)

    def __getattr__(self, attr):
        """
        EXAMPLE:
            sage: execfile('des.py')
            sage: des = DES(Nr=2)
            sage: des.k[0]
            k00

            sage: des = DES(Nr=1)
            sage: des.R.ngens()
            464

            sage: des.sbox_eq
            'opns_gb'
        """
        if attr in self._params.keys():
            return getattr(self, "_" + attr)
        if attr == "R":
            self.R = self.ring()
            return self.R
        if attr == "k":
            return self.vars("k",0)
        else:
            raise AttributeError("'%r' object has no attribute '%s'"%(type(self),attr))
            
    def _repr_(self):
        """
        EXAMPLE:
            sage: execfile('des.py')
            sage: des = DES(Nr=2)
            sage: des # indirect doctest
            DES(Nr=2)
        """
        return "DES(Nr=%d)"%self.Nr

    def __cmp__(self, other):
        """
        EXAMPLE:
            sage: execfile('des.py')
            sage: des1 = DES(Nr=2)
            sage: des2 = DES(Nr=5)
            sage: des1 == des2
            False
        """
        return cmp(self._params, other._params)

    def __reduce__(self):
        """
        TESTS:
            sage: execfile('des.py')
            sage: des = DES()
            sage: des == loads(dumps(des)) # not tested since it is broken!
            True
        """
        return unpickle_DES, (self._params,)

    def random_element(self, length=64):
        """
        Return random elements suitable as keys, plaintexts etc.

        INPUT:
            length -- number of bits (default: 64)

        EXAMPLE:
            sage: execfile('des.py')
            sage: des = DES()
            sage: des.random_element() # e.g. a plaintext
            [0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0]
            sage: des.random_element(56) # e.g. a key
            [1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0]
        """
        re = self._base.random_element
        return [re() for _ in xrange(length)]

    def base_ring(self):
        """
        Return GF(2)

        EXAMPLE:
            sage: execfile('des.py')
            sage: des = DES()
            sage: des.base_ring()
            Finite Field of size 2
        """
        return self._base

    def IP(self, X):
        """
        Perform the initial permutation $IP$ as specified in the DES.

        INPUT:
            X -- a list/vector of length 64 for IP to act on

        EXAMPLE:
            sage: execfile('des.py')
            sage: des = DES()
            sage: e = range(64)
            sage: des.IP(e)
            [57, 49, 41, 33, 25, 17, 9, 1, 59, 51, 43, 35, 27, 19, 11,
             3, 61, 53, 45, 37, 29, 21, 13, 5, 63, 55, 47, 39, 31, 23,
             15, 7, 56, 48, 40, 32, 24, 16, 8, 0, 58, 50, 42, 34, 26,
             18, 10, 2, 60, 52, 44, 36, 28, 20, 12, 4, 62, 54, 46, 38,
             30, 22, 14, 6]
        """
        IP = [57, 49, 41, 33, 25, 17, 9, 1, 59, 51, 43, 35, 27, 19,
              11, 3, 61, 53, 45, 37, 29, 21, 13, 5, 63, 55, 47, 39,
              31, 23, 15, 7, 56, 48, 40, 32, 24, 16, 8, 0, 58, 50, 42,
              34, 26, 18, 10, 2, 60, 52, 44, 36, 28, 20, 12, 4, 62,
              54, 46, 38, 30, 22, 14, 6]

        return [ X[ i ] for i in IP ]

    def E(self, X):
        """
        Perform the expansion permutation $E$ from a list/vector of
        length 32 to a list of length 48 as specificed in the DES.

        INPUT:
            X -- a list/vector of length 32 for E to act on

        EXAMPLE:
            sage: execfile('des.py')
            sage: des = DES()
            sage: e = range(32)
            sage: des.E(e)
            [31, 0, 1, 2, 3, 4, 3, 4, 5, 6, 7, 8, 7, 8, 9, 10, 11, 12, 11, 12, 13, 14, 15, 16, 15, 16, 17, 18, 19, 20, 19, 20, 21, 22, 23, 24, 23, 24, 25, 26, 27, 28, 27, 28, 29, 30, 31, 0]
        """

        E = [31, 0, 1, 2, 3, 4, 3, 4, 5, 6, 7, 8, 7, 8, 9, 10, 11, 12,
             11, 12, 13, 14, 15, 16, 15,16, 17, 18, 19, 20, 19, 20,
             21, 22, 23, 24, 23, 24, 25, 26, 27, 28, 27, 28, 29, 30,
             31, 0]
        
        return [ X[ i ] for i in E ]

    def add(self, X, Y, Z=None):
        """
        Component-wise addition of the elements in X and Y (and Z if
        provided).

        INPUT:
            X -- a list/vector
            Y -- a list/vector
            Z -- an optional list/vector
        """
        if Z is not None:
            return [ X[i]  + Y[i] + Z[i] for i in xrange(len(X)) ]
        else:
            return [ X[i]  + Y[i] for i in xrange(len(X)) ]

    def S(self, X):
        """
        Perfom the S-Box operation on X.

        INPUT:
            X -- a list/vector of length 48
        """
        Y = sum([self._S[i/6][X[i:i+6]] for i in xrange(0,48,6)],[])
        return Y

    def P(self, X):
        """
        Perform the compression permutation $P$ as specified in the
        DES.

        INPUT:
            X -- list/vector of length 32
            

        EXAMPLE:
            sage: execfile('des.py')
            sage: des = DES()
            sage: e = range(32)
            sage: des.P(e)
            [15, 6, 19, 20, 28, 11, 27, 16, 0, 14, 22, 25, 4, 17, 30, 9, 1, 7, 23, 13, 31, 26, 2, 8, 18, 12, 29, 5, 21, 10, 3, 24]
        """
        P = [15, 6, 19, 20, 28, 11, 27, 16, 0, 14, 22, 25, 4, 17, 30,
             9, 1, 7, 23, 13, 31, 26, 2, 8, 18, 12, 29, 5, 21, 10, 3,
             24]

        return [ X[ i ] for i in P ]

    def f(self, R, K):
        r"""
        Perform the round function $f$ of the DES specified as

        $$f(R_{i-1},K_i) = P( S( E(R_{i-1}) \oplus K) ).$$

        INPUT:
            R -- a vector/list of length 32
            K -- a vector/list (round key) of length 48
        """
        return self.P( self.S( self.add( self.E(R), K ) ) )

    def IPm1(self, X):
        """
        Perform the final permutation $IP^{-1}$ as specified in the
        DES.

        INPUT:
            X -- a list/vector of length 64 for $IP^{-1}$ to act on

        EXAMPLE:
            sage: execfile('des.py')
            sage: des = DES()
            sage: des.IPm1(range(64))
            [39, 7, 47, 15, 55, 23, 63, 31, 38, 6, 46, 14, 54, 22, 62,
            30, 37, 5, 45, 13, 53, 21, 61, 29, 36, 4, 44, 12, 52, 20,
            60, 28, 35, 3, 43, 11, 51, 19, 59, 27, 34, 2, 42, 10, 50,
            18, 58, 26, 33, 1, 41, 9, 49, 17, 57, 25, 32, 0, 40, 8,
            48, 16, 56, 24]
        """
        IPm1 = [39, 7, 47, 15, 55, 23, 63, 31, 38, 6, 46, 14, 54, 22,
                62, 30, 37, 5, 45, 13, 53, 21, 61, 29, 36, 4, 44, 12,
                52, 20, 60, 28, 35, 3, 43, 11, 51, 19, 59, 27, 34, 2,
                42, 10, 50, 18, 58, 26, 33, 1, 41, 9, 49, 17, 57, 25,
                32, 0, 40, 8, 48, 16, 56, 24]

        return [ X[ i ] for i in IPm1 ]

    def PC1(self, X):
        """
        Perform the $PC-1$ permutation as specified in the DES.

        INPUT:
            X -- a list/vector of length 64 to act on

        EXAMPLE:
            sage: execfile('des.py')
            sage: des = DES()
            sage: des.PC1(range(64))
            [56, 48, 40, 32, 24, 16, 8, 0, 57, 49, 41, 33, 25, 17, 9,
            1, 58, 50, 42, 34, 26, 18, 10, 2, 59, 51, 43, 35, 62, 54,
            46, 38, 30, 22, 14, 6, 61, 53, 45, 37, 29, 21, 13, 5, 60,
            52, 44, 36, 28, 20, 12, 4, 27, 19, 11, 3]
        """
        PC_1 = [56, 48, 40, 32, 24, 16, 8, 0, 57, 49, 41, 33, 25, 17, 9,
                1, 58, 50, 42, 34, 26, 18, 10, 2, 59, 51, 43, 35, 62, 54,
                46, 38, 30, 22, 14, 6, 61, 53, 45, 37, 29, 21, 13, 5, 60,
                52, 44, 36, 28, 20, 12, 4, 27, 19, 11, 3]
        return [ X[ i ] for i in PC_1 ]

    def PC2(self, X):
        """
        Perform the $PC-2$ permutation as specified in the DES.

        INPUT:
            X -- a list/vector of length 56 to act on

        EXAMPLE:
            sage: execfile('des.py')
            sage: des = DES()
            sage: des.PC2(range(64))
            [13, 16, 10, 23, 0, 4, 2, 27, 14, 5, 20, 9, 22, 18, 11, 3,
            25, 7, 15, 6, 26, 19, 12, 1, 40, 51, 30, 36, 46, 54, 29,
            39, 50, 44, 32, 47, 43, 48, 38, 55, 33, 52, 45, 41, 49,
            35, 28, 31]
        """
        PC_2 = [13, 16, 10, 23, 0, 4, 2, 27, 14, 5, 20, 9, 22, 18, 11, 3,
                25, 7, 15, 6, 26, 19, 12, 1, 40, 51, 30, 36, 46, 54, 29,
                39, 50, 44, 32, 47, 43, 48, 38, 55, 33, 52, 45, 41, 49,
                35, 28, 31]
        return [ X[ i ] for i in PC_2 ]


    def Ks(self, K):
        """
        Return list of subkeys $K_1, ..., K_{Nr}$ derrived from $K$
        according to the key schedule specification of the DES.

        If $K$ has length 64 it is assumed that it contains parity
        bits which are ignored. If $K$ has length 56 it is assumed it
        only contains real key bits.

        INPUT:
            K -- a list/vector of either length 56 or 64

        EXAMPLE:
            sage: execfile('des.py')
            sage: des = DES()
            sage: K = '0001001100110100010101110111100110011011101111001101111111110001'
            sage: K = [int(b) for b in K]
            sage: for k in des.Ks(K):
            ...    print "".join(map(str,k))
            000110110000001011101111111111000111000001110010
            011110011010111011011001110110111100100111100101
            010101011111110010001010010000101100111110011001
            011100101010110111010110110110110011010100011101
            011111001110110000000111111010110101001110101000
            011000111010010100111110010100000111101100101111
            111011001000010010110111111101100001100010111100
            111101111000101000111010110000010011101111111011
            111000001101101111101011111011011110011110000001
            101100011111001101000111101110100100011001001111
            001000010101111111010011110111101101001110000110
            011101010111000111110101100101000110011111101001
            100101111100010111010001111110101011101001000001
            010111110100001110110111111100101110011100111010
            101111111001000110001101001111010011111100001010
            110010110011110110001011000011100001011111110101
        """
        if len(K) == 56:
            #reintroduce parity bits
            K = K[ 0: 7]+[0]+\
                K[ 7:14]+[0]+\
                K[14:21]+[0]+\
                K[21:28]+[0]+\
                K[28:35]+[0]+\
                K[35:42]+[0]+\
                K[42:49]+[0]+\
                K[49:56]+[0]

        K = self.PC1(K)
        KL, KR = K[0:28], K[28:56]

        Ki = []

        for i in range(self.Nr):
            ki = self.L(i, KL, KR)
            Ki.append(ki)

        return Ki

    def L(self, i, k0, k1=None):
        """
	Return a subkey for round i from the initial key k0 or
	[k0,k1], if k1 is provided.

        INPUT:
            i -- round counter
            k0 -- list of key bits
            k1 -- additional optional list of key bits
	"""
        if k1 is None and len(k0) == 56:
            k0, k1 = k0[0:28], k0[28:]

        k0bar, k1bar = k0, k1
        for j in range(i+1):
            k0bar, k1bar = self.LS(k0bar,j+1), self.LS(k1bar,j+1)
        return self.PC2(k0bar + k1bar)

    def v(self, i):
        """
        Return the cyclic left-shift constant for $i$

        INPUT:
            i -- integer
        """
        if i in (1,2,9,16):
            return 1
        else:
            return 2
        
    def LS(self, X, i):
        """
        Perform cyclic left-shift on $X$ according to $i$.

        INPUT:
            X -- a list/vector of length 28
            i -- an integer $0<i<=16$.
        """
        i = self.v(i)
        return X[i:]+X[:i]

    def __call__(self, P, K):
        """
        Encrypt the plaintext $P$ under the key $K$.
        
        INPUT:
            P -- a list/vector of length 64
            K -- a list/vector of length either 56 or 64

        EXAMPLE:
            sage: execfile('des.py')
            sage: des = DES()
            sage: M = '0000000100100011010001010110011110001001101010111100110111101111'
            sage: K = '0001001100110100010101110111100110011011101111001101111111110001'
            sage: M = [int(b) for b in M]
            sage: K = [int(b) for b in K]
            sage: ''.join(map(str, des(M,K)))
            '1000010111101000000100110101010000001111000010101011010000000101'
        """
        K = [self._base(e) for e in K]
        P = [self._base(e) for e in P]

        Ki = self.Ks(K)
        x = self.IP(P)
        L,R = x[:32],x[32:]
        for i in range(self.Nr):
            L,R = R, self.add(L, self.f(R, Ki[i]) )
            #print to_str(L+R)
        return self.IPm1(R+L)

    def varformatstr(self, name):
        r"""
        Return format string for variable named \code{name}.

        INPUT:
            name -- variable name
        """
        l = str(max([len(str(self.Nr)), 2]))

        if name.startswith("k"):
            return name + "%0" + l + "d"
        if name.startswith("P") or name.startswith("C"): 
            name += self._postfix
            return name + "%0" + l + "d"
        else:
            name += self._postfix
            return name + "%0" + l + "d" + "%0" + l + "d"

    def varstrs(self, name, r):
        r"""
        Return variable string for variable named \code{name} and
        round \code{r}.

        INPUT:
            name -- variable name
            r -- integer $0 < i <= 16$
        """
        s = self.varformatstr(name)
        if name == "k":
            return [s%(i) for i in xrange(56) ]
        if name.startswith("t"):
            return [s%(r,i) for i in xrange(self._tlength) ]
        if name.startswith("P") or name.startswith("C"):
            return [s%(i) for i in xrange(64)]
        else:
            return [s%(r,i) for i in xrange(32) ]

    def vars(self, name, i, *args):
        r"""
        Return variables for variable named \code{name} and round
        \code{i}.

        INPUT:
            name -- variable name
            i -- integer $0 < i <= 16$
        """
        gd = self.R.gens_dict()
        return [gd[e] for e in self.varstrs(name, i)]

    def ring(self, order=None, n=1):
        """
        Return a new ring (and cache it) for this equation system
        generator.

        INPUT:
            order -- term ordering (default: None)
            n -- number of plaintext-ciphertext pairs
        """
        if order is None:
            order = self._order
        if order == "block":
            order = self.block_order(n)

        var_names = []
        var_names +=[self.varformatstr("k")%(j) for j in xrange(56)]
        if n == 1:
            for nr in range(1,self.Nr+1):
                if nr < self.Nr-1:
                    var_names += [self.varformatstr("z")%(nr,j) for j in xrange(32)]
                var_names += [self.varformatstr("y")%(nr,j) for j in xrange(32)]
                if self.sbox_eq.startswith("opns"):
                    var_names += [self.varformatstr("t")%(nr,j) for j in xrange(self._tlength)]
        elif n > 1:
            for nr in range(1,self.Nr+1):
                for i in xrange(n):
                    if nr < self.Nr-1:
                        var_names += [self.varformatstr("z%d"%i)%(nr,j) for j in xrange(32)]
                    var_names += [self.varformatstr("y%d"%i)%(nr,j) for j in xrange(32)]
                    if self.sbox_eq.startswith("opns"):
                        var_names += [self.varformatstr("t%d"%i)%(nr,j) for j in xrange(self._tlength)]
        
        if self._polybori:
            R = BooleanPolynomialRing(len(var_names), var_names, order=order )
        else:
            R = PolynomialRing(self._base, len(var_names), var_names, order=order)

        if n == 1:
            self.R = R

        return R

    def t(self, i):
        if self.sbox_eq.startswith("opns"):
            T = self.vars("t",  i )
        else:
            T = []
        return T

    def _load_sbox_eq(self):
        """
        """
        if self._sbox_eq == "cubic":
            if "s1cubic" not in globals():
                print "Pre-computing fully cubic S-Box equations."
                load_sXcubic()

        elif self._sbox_eq == "opns_gb":
            if "s1opns_gb" not in globals():
                print "Pre-computing Groebner bases for S-Boxes, this might take a moment."
                load_sXopns_gb()

    def polynomial_system(self, P=None, K=None, C=None):
        """
        Return a DES polynomial system and the plaintext $P$ and the
        key $K$.

        INPUT:
            P -- plaintext (default: random)
            K -- key (default: random)

        OUTPUT:
            MPolynomialSystem, dictionary of solutions
        """

        symbolic = False
        self._load_sbox_eq()
        
        if not P:
            P = self.random_element()
        if not K and not C:
            K = self.random_element(56)
            Kval, Kvar = list(K), list(self.k)
        else:
            symbolic = True
            Kvar = list(self.k)
        if not C:
            C = self(P, Kval)

        rounds = []
        K = iter(self.Ks(Kvar))
        P = self.IP(P)

        L,R = P[:32], P[32:]

        for i in xrange(1,self.Nr+1):
            X = self.add( self.E(R), K.next() )
            T = self.t(i)
            Y = self.vars("y", i)

            if i == self.Nr-1:
                Z = self.IP(C)[32:]
            elif i == self.Nr:
                Z = self.IP(C)[:32]
            else:
                Z = self.vars("z", i)

            # Y = S(E(R_{i-1}) + K_i)
            sbox = self.sbox_polynomials( X, Y, T)

            # R_i = P(Y) + L_{i-1}
            lin = self.lin_polynomials(Z, Y, L)

            # field polynomials
            fld = self.field_polynomials(T)
            fld += self.field_polynomials(Y)
            if i < self.Nr - 1:
                fld += self.field_polynomials(Z)

            L,R  = R,Z

            rounds.append(mq.MPolynomialRoundSystem(self.R, sbox+lin+fld))

        L = [L[i] + self.IP(C)[32:][i] for i in range(len(L))]
        R = [R[i] + self.IP(C)[:32][i] for i in range(len(R))]

        rounds.append(mq.MPolynomialRoundSystem(self.R, L+R))

        if not symbolic:
            return mq.MPolynomialSystem(self.R, rounds), dict(zip(Kvar,Kval))
        else:
            return mq.MPolynomialSystem(self.R, rounds), {}

    def field_polynomials(self, X):
        r"""
        Return $x_i^2 + x_i$ for each $x_i$ in X.
        
        INPUT:
            X -- a list

        NOTE: This method returns an empty list if
        \code{self.polybori} is \code{True} because in that case field
        equations are redundant.

        INPUT:
            X -- a list 
        """
        if self._polybori:
            return []
        else:
            return [x**2 + x for x in X]

    def lin_polynomials(self, Z, Y, L):
        """
        Return linear polynomials for the equation $Z = P(Y) + L$.

        INPUT:
            Z -- list of length 32
            Y -- list of length 32
            L -- list of length 32

        EXAMPLE:
            sage: execfile('des.py')
            sage: Z = [' z%02d '%i for i in range(32)]
            sage: Y = [' y%02d '%i for i in range(32)]
            sage: L = [' l%02d '%i for i in range(32)]
            sage: des = DES()
            sage: des.lin_polynomials(Z, Y, L)
            [' z00  y15  l00 ',
             ' z01  y06  l01 ',
             ' z02  y19  l02 ',
             ' z03  y20  l03 ',
             ' z04  y28  l04 ',
             ' z05  y11  l05 ',
             ' z06  y27  l06 ',
             ' z07  y16  l07 ',
             ' z08  y00  l08 ',
             ' z09  y14  l09 ',
             ' z10  y22  l10 ',
             ' z11  y25  l11 ',
             ' z12  y04  l12 ',
             ' z13  y17  l13 ',
             ' z14  y30  l14 ',
             ' z15  y09  l15 ',
             ' z16  y01  l16 ',
             ' z17  y07  l17 ',
             ' z18  y23  l18 ',
             ' z19  y13  l19 ',
             ' z20  y31  l20 ',
             ' z21  y26  l21 ',
             ' z22  y02  l22 ',
             ' z23  y08  l23 ',
             ' z24  y18  l24 ',
             ' z25  y12  l25 ',
             ' z26  y29  l26 ',
             ' z27  y05  l27 ',
             ' z28  y21  l28 ',
             ' z29  y10  l29 ',
             ' z30  y03  l30 ',
             ' z31  y24  l31 ']
        """
        return self.add(Z, self.P(Y), L)

    def sbox_polynomials(self, X, Y, T=None):
        """
        Return S-Box polynomials for all S-Boxes in 48 input variables
        $X_i$, 32 output variables $Y_i$ (and $376$ intermediate variables
        $T_i$).

        INPUT:
            X -- list of length 48
            Y -- list of length 32
            T -- optional list of length 376

        EXAMPLE:
            sage: execfile('des.py')
            sage: des = DES(Nr=4)
            sage: K = des.k
            sage: Ki = des.Ks(K)
            sage: X = des.E(des.vars('z',1)) + Ki[1]
            sage: Y = des.vars('y',2)
            sage: T = des.vars('t',2)
            sage: sb = des.sbox_polynomials(X, Y, T)
            Pre-computing Groebner bases for S-Boxes, this might take a moment.

            sage: sb[:2]
            [z0104*z0103 + t0233 + t0235 + t0238 + t0243 + t0244 + t0248 + t0249 + t0250,
             z0104*z0102 + t0229 + t0239 + z0102 + z0103]

            sage: len(sb)
            3222
        """
        self._load_sbox_eq()

        if self._sbox_eq == "opns":
            ret  = s1opns( X[ 0: 6], Y[ 0: 4], T[  0: 52])
            ret += s2opns( X[ 6:12], Y[ 4: 8], T[ 52: 98])
            ret += s3opns( X[12:18], Y[ 8:12], T[ 98:147])
            ret += s4opns( X[18:24], Y[12:16], T[147:182])
            ret += s5opns( X[24:30], Y[16:20], T[182:234])
            ret += s6opns( X[30:36], Y[20:24], T[234:283])
            ret += s7opns( X[36:42], Y[24:28], T[283:330])
            ret += s8opns( X[42:48], Y[28:32], T[330:376])

        elif self._sbox_eq == "cubic":
            ret  = s1cubic( X[ 0: 6], Y[ 0: 4])
            ret += s2cubic( X[ 6:12], Y[ 4: 8])
            ret += s3cubic( X[12:18], Y[ 8:12])
            ret += s4cubic( X[18:24], Y[12:16])
            ret += s5cubic( X[24:30], Y[16:20])
            ret += s6cubic( X[30:36], Y[20:24])
            ret += s7cubic( X[36:42], Y[24:28])
            ret += s8cubic( X[42:48], Y[28:32])

        elif self._sbox_eq == "opns_gb":
            ret  = s1opns_gb( X[ 0: 6], Y[ 0: 4], T[  0: 52])
            ret += s2opns_gb( X[ 6:12], Y[ 4: 8], T[ 52: 98])
            ret += s3opns_gb( X[12:18], Y[ 8:12], T[ 98:147])
            ret += s4opns_gb( X[18:24], Y[12:16], T[147:182])
            ret += s5opns_gb( X[24:30], Y[16:20], T[182:234])
            ret += s6opns_gb( X[30:36], Y[20:24], T[234:283])
            ret += s7opns_gb( X[36:42], Y[24:28], T[283:330])
            ret += s8opns_gb( X[42:48], Y[28:32], T[330:376])
        else:
            raise TypeError, "sbox_eq='%s' parameter unknown."%self.sbox_eq
            

        return ret

def unpickle_DES(kwds):
    """
    TESTS:
        sage: execfile('des.py')
        sage: des = DES()
        sage: des == loads(dumps(des)) #  # not tested since it is broken!
        True
    """
    return DES(**kwds)


class DESDC(DES):
    def __init__(self, r=1, *args, **kwds):
        r"""
        Return a polynomial system generator for a \code{DESDC} equation
        system.

        \code{DESDC} is a cipher that is related to DES as
        follows. Given two plaintexts $P'$ and $P''$ related by the
        difference $\delta_0$ assume that after $r$ rounds the
        difference $\delta_r$ holds. Also denote $C'$ and $C''$ the
        ciphertexts under the key $K$ for $P'$ and $P''$
        respectively.\code{DESDC} is the cipher that decrypts $N_r -
        r$ rounds from $C'$ to the assumed difference $\delta_r$
        relates the state variables via that difference and encrypts
        to $C2$ again. If $\delta_r$ holds with probability $p$ then
        the probability that $C2 = C''$ is $p$ too.
        
        """
        DES.__init__(self, *args, **kwds)
        self._params["r"] = r
        D = kwds.get("D", DESCharacteristic(self))
        if isinstance(D,str):
            D = from_str(D)
        elif isinstance(D, DifferentialCharacteristicIterator):
            it = iter(D)
            d = it.next()
            for _ in range(r):
                d = it.next()
            D = d
        self._params["D"] = D
        self._D = D
        self._r = r

    def _repr_(self):
        return "DESDC(Nr=%d,r=%d)"%(self.Nr,self.r)

    def __call__(self, C, K):
        Nr = self.Nr
        r = self.r
        K = [self._base(e) for e in K]
        C = [self._base(e) for e in C]
        D = self.D

        Ki = self.Ks(K)
        C = self.IP(C)
        L,R = C[:32], C[32:]

        # decrypt last Nr-r rounds
        for i in range(Nr-r):
            L,R = R, self.add(L, self.f(R, Ki[Nr-1-i]) )
            #print to_str(L+R)
            
        # relate via difference
        L,R = self.add(R, D[32:]), self.add(L, D[:32])

        # encrypt last Nr-r rounds
        for i in range(Nr-r):
            L,R = R, self.add(L, self.f(R, Ki[r+i]) )
            #print to_str(L+R)
                         
        return self.IPm1(R+L)

    def polynomial_system(self, C, Cbar):
        """
        """
        Nr, r = self.Nr, self.r
        length = Nr-r

        D = self.D

        self._load_sbox_eq()

        rounds = []
        Ki = self.Ks(self.k)

        C = self.IP(C)
        L,R = C[:32], C[32:]

        # going backward from C to difference
        for i in xrange(length):
            X = self.add( self.E(R), Ki[Nr-1-i] )
            T = self.t(i)
            Y = self.vars("y", i )
            Z = self.vars("z", i )
            
            sbox = self.sbox_polynomials(X, Y, T)
            lin = self.lin_polynomials(Z, Y, L)

            fld  = self.field_polynomials(T)
            fld += self.field_polynomials(Y)
            fld += self.field_polynomials(Z)

            rounds.append(mq.MPolynomialRoundSystem(self.R, [self.R(f) for f in sbox+lin+fld]))

            L, R = R, Z

        # relate via difference
        L,R = self.add(R, D[32:]), self.add(L, D[:32])

        # going forward from difference to Cbar
        for i in xrange(length):
            X = self.add( self.E(R), Ki[r+i] )
            T = self.t( length + i )
            Y = self.vars("y", length + i )

            if r+i+1 == self.Nr-1:
                Z = self.IP(Cbar)[32:]
            elif r+i+1 == self.Nr:
                Z = self.IP(Cbar)[:32]
            else:
                Z = self.vars("z", length + i )

            sbox = self.sbox_polynomials(X, Y, T)
            lin = self.lin_polynomials(Z, Y, L)

            fld  = self.field_polynomials(T)
            fld += self.field_polynomials(Y)
            fld += self.field_polynomials(Z)

            rounds.append(mq.MPolynomialRoundSystem(self.R, [self.R(f) for f in sbox+lin+fld]))

            L, R = R, Z


        # additional linear information about the variables
        L,R = C[:32], C[32:]
        # decrypt last Nr-r rounds
        for i in range(Nr-r):
            L,R = R, self.add(L, self.P(self.vars("y", i )))
            
        # relate via difference
        L,R = self.add(R, D[32:]), self.add(L, D[:32])

        # encrypt last Nr-r rounds
        for i in range(Nr-r):
            L,R = R, self.add(L, self.P(self.vars("y", length + i)))

        lin = self.add(self.IP(Cbar), R+L)
        rounds.append(mq.MPolynomialRoundSystem(self.R, lin ))

        return mq.MPolynomialSystem(self.R, rounds)

    def ring(self, order=None):
        """
        Return a new ring (and cache it) for this equation system
        generator.

        INPUT:
            order -- term ordering (default: None)
        """
        if order is None:
            order = self._order
        if order == "block":
            order = self.block_order()

        Nr = self.Nr
        r = self.r
        length = Nr - r

        var_names = []
        for nr in reversed(range(2*length)):
            if nr < 2*length - 2 or nr==0:
                var_names += [self.varformatstr("z")%(nr,j) for j in xrange(32)]
                
            var_names += [self.varformatstr("y")%(nr,j) for j in xrange(32)]
            if self.sbox_eq.startswith("opns"):
                var_names += [self.varformatstr("t")%(nr,j) for j in xrange(self._tlength)]
        
        var_names +=[self.varformatstr("k")%(j) for j in xrange(56)]

        if self._polybori:
            R = BooleanPolynomialRing(len(var_names), var_names, order=order )
        else:
            R = MPolynomialRing(self._base, len(var_names), var_names, order=order)

        self.R = R
        return R

class DESSBox(mq.SBox):
    def __init__(self, i=1):
        """
        Return one of the eight DES S-Boxes depending on the parameter
        $i$.

        INPUT:
            i -- integer $1<=i<=8$

        EXAMPLE:
            sage: execfile('des.py')
            sage: DESSBox(1)
            S1
        """
        self._F = GF(2)
        self.m = 6
        self.n = 4
        self._big_endian = True

        S = [ [[14,  4, 13,  1,  2, 15, 11,  8,  3, 10,  6, 12,  5,  9,  0,  7],
               [ 0, 15,  7,  4, 14,  2, 13,  1, 10,  6, 12, 11,  9,  5,  3,  8],
               [ 4,  1, 14,  8, 13,  6,  2, 11, 15, 12,  9,  7,  3, 10,  5,  0],
               [15, 12,  8,  2,  4,  9,  1,  7,  5, 11,  3, 14, 10,  0,  6, 13]],
              
              [[15,  1,  8, 14,  6, 11,  3,  4,  9,  7,  2, 13, 12,  0,  5, 10],
               [ 3, 13,  4,  7, 15,  2,  8, 14, 12,  0,  1, 10,  6,  9, 11,  5],
               [ 0, 14,  7, 11, 10,  4, 13,  1,  5,  8, 12,  6,  9,  3,  2, 15],
               [13,  8, 10,  1,  3, 15,  4,  2, 11,  6,  7, 12,  0,  5, 14,  9]],

              [[10,  0,  9, 14,  6,  3, 15,  5,  1, 13, 12,  7, 11,  4,  2,  8],
               [13,  7,  0,  9,  3,  4,  6, 10,  2,  8,  5, 14, 12, 11, 15,  1],
               [13,  6,  4,  9,  8, 15,  3,  0, 11,  1,  2, 12,  5, 10, 14,  7],
               [ 1, 10, 13,  0,  6,  9,  8,  7,  4, 15, 14,  3, 11,  5,  2, 12]],
              
              [[ 7, 13, 14,  3,  0,  6,  9, 10,  1,  2,  8,  5, 11, 12,  4, 15],
               [13,  8, 11,  5,  6, 15,  0,  3,  4,  7,  2, 12,  1, 10, 14,  9],
               [10,  6,  9,  0, 12, 11,  7, 13, 15,  1,  3, 14,  5,  2,  8,  4],
               [ 3, 15,  0,  6, 10,  1, 13,  8,  9,  4,  5, 11, 12,  7,  2, 14]],

              [[ 2, 12,  4,  1,  7, 10, 11,  6,  8,  5,  3, 15, 13,  0, 14,  9],
               [14, 11,  2, 12,  4,  7, 13,  1,  5,  0, 15, 10,  3,  9,  8,  6],
               [ 4,  2,  1, 11, 10, 13,  7,  8, 15,  9, 12,  5,  6,  3,  0, 14],
               [11,  8, 12,  7,  1, 14,  2, 13,  6, 15,  0,  9, 10,  4,  5,  3]],

              [[12,  1, 10, 15,  9,  2,  6,  8,  0, 13,  3,  4, 14,  7,  5, 11],
               [10, 15,  4,  2,  7, 12,  9,  5,  6,  1, 13, 14,  0, 11,  3,  8],
               [ 9, 14, 15,  5,  2,  8, 12,  3,  7,  0,  4, 10,  1, 13, 11,  6],
               [ 4,  3,  2, 12,  9,  5, 15, 10, 11, 14,  1,  7,  6,  0,  8, 13]],

              [[ 4, 11,  2, 14, 15,  0,  8, 13,  3, 12,  9,  7,  5, 10,  6,  1],
               [13,  0, 11,  7,  4,  9,  1, 10, 14,  3,  5, 12,  2, 15,  8,  6],
               [ 1,  4, 11, 13, 12,  3,  7, 14, 10, 15,  6,  8,  0,  5,  9,  2],
               [ 6, 11, 13,  8,  1,  4, 10,  7,  9,  5,  0, 15, 14,  2,  3, 12]],

              [[13,  2,  8,  4,  6, 15, 11,  1, 10,  9,  3, 14,  5,  0, 12,  7],
               [ 1, 15, 13,  8, 10,  3,  7,  4, 12,  5,  6, 11,  0, 14,  9,  2],
               [ 7, 11,  4,  1,  9, 12, 14,  2,  0,  6, 10, 13, 15,  3,  5,  8],
               [ 2,  1, 14,  7,  4, 10,  8, 13, 15, 12,  9,  0,  3,  5,  6, 11]] ]

        self._S = [[self.to_bits(x,4) for x in r] for r in S[i-1]]

    def __call__(self, X):
        r"""
        Apply subsitution to X.

        If X is a list it is interpreted as sequence of bits depending
        on the bit order of this S-Box.

        INPUT:
            X -- either an integer or a tuple of GF(2) elements of
                 length \code{self.m}
        """
        ret = lambda x: x
        if isinstance(X, Integer) or isinstance(X, int):
            X = self.to_bits(X,6)
            ret = lambda x: ZZ(map(int, list(reversed(x))),2)
        if isinstance(X, (list,tuple)) and len(X) == 6:
            b1,b2,b3,b4,b5,b6 = map(int, X)
            return ret( self._S[ZZ([b6,b1],2)][ZZ([b5,b4,b3,b2],2)] )
        else:
            raise TypeError

    def _repr_(self):
        """
        Return name of S-Box ('S1' to 'S8')
        """
        if self._S[0][0] == [1,1,1,0]:
            return "S1"
        elif self._S[0][0] == [1,1,1,1]:
            return "S2"
        elif self._S[0][0] == [1,0,1,0]:
            return "S3"
        elif self._S[0][0] == [0,1,1,1]:
            return "S4"
        elif self._S[0][0] == [0,0,1,0]:
            return "S5"
        elif self._S[0][0] == [1,1,0,0]:
            return "S6"
        elif self._S[0][0] == [0,1,0,0]:
            return "S7"
        elif self._S[0][0] == [1,1,0,1]:
            return "S8"
        else:
            return "Unknown S-Box"

def s1opns(x, y, t):
    r"""
    Return sparse quadratic equations for the S-Box S1 in the input
    variables $x_0,...,x_5$, the output variables $y_0,...,y_3$ and 52
    intermediate variables $t_i$.

    These equations are based on

       http://www.darkside.com.au/bitslice/nonstd.c

    INPUT:
        x -- a list of length 6
        y -- a list of length 4
        t -- a list of length 52
    """
    t = list(t)

    # we add these silly entries, to have uniform notation with
    # nonstd.c
    
    t.insert(21,None)
    t.insert(34,None)
    t.insert(45,None)
    t.insert(55,None)

    return [ t[0] + x[2] * (x[4] + 1), t[1] + t[0] + x[3], t[2] + x[2]
        * (x[3] + 1), t[3] + t[2] * x[4] + t[2] + x[4], t[4] + x[5] *
        t[3], t[5] + t[1] + t[4], t[6] + x[3] * (x[4] + 1), t[7] +
        x[2] + x[3], t[8] + x[5] * (t[7] + 1), t[9] + t[6] + t[8],
        t[10] + x[1] * t[9] + x[1] + t[9], t[11] + t[5] + t[10], t[12]
        + x[4] + t[4], t[13] + t[12] * t[7], t[14] + x[4] * (x[3] +
        1), t[15] + t[2] + t[13], t[16] + x[5] * t[15] + x[5] + t[15],
        t[17] + t[14] + t[16], t[18] + x[1] * t[17] + x[1] + t[17],
        t[19] + t[13] + t[18], t[20] + x[0] * t[19], y[1] + t[11] +
        (t[20] + 1), t[22] + t[0] * t[4] + t[0] + t[4], t[23] + t[22]
        + t[7], t[24] + t[17] * (t[1] + 1), t[25] + x[1] * (t[24] +
        1), t[26] + t[23] + t[25], t[27] + t[5] * t[6] + t[5] + t[6],
        t[28] + t[27] + t[24], t[29] + t[8] + t[23], t[30] + t[17] *
        (t[29] + 1), t[31] + x[1] * t[30], t[32] + t[28] + t[31],
        t[33] + x[0] * t[32], y[3] + t[26] + t[33], t[35] + x[2] *
        t[27], t[36] + t[17] * (t[35] + 1), t[37] + x[1] * t[2] + x[1]
        + t[2], t[38] + t[36] + t[37], t[39] + x[2] * t[30] + x[2] +
        t[30], t[40] + t[23] * (t[36] + 1), t[41] + t[40] * t[2] +
        t[40] + t[2], t[42] + t[41] * (x[1] + 1), t[43] + t[39] +
        t[42], t[44] + x[0] * (t[43] + 1), y[0] + t[38] + (t[44] + 1),
        t[46] + t[32] * (t[8] + 1), t[47] + t[46] + t[38], t[48] +
        t[3] + t[35], t[49] + t[48] * (t[4] + 1), t[50] + t[41] *
        t[17] + t[41] + t[17], t[51] + t[50] + x[4], t[52] + x[1] *
        (t[51] + 1), t[53] + t[49] + t[52], t[54] + x[0] * t[53] +
        x[0] + t[53], y[2] + t[47] + (t[54] + 1), ]


def s2opns(x, y, t):
    r"""
    Return sparse quadratic equations for the S-Box S1 in the input
    variables $x_0,...,x_5$, the output variables $y_0,...,y_3$ and 46
    intermediate variables $t_i$.

    These equations are based on

       http://www.darkside.com.au/bitslice/nonstd.c

    INPUT:
        x -- a list of length 6
        y -- a list of length 4
        t -- a list of length 46
    """
    t = list(t)

    # we add these silly entries, to have uniform notation with
    # nonstd.c

    t.insert(13,None)
    t.insert(30,None)
    t.insert(39,None)
    t.insert(49,None)

    return [ t[0] + x[0] + x[5], t[1] + t[0] + x[4], t[2] + x[5] *
        x[4], t[3] + x[0] * (t[2] + 1), t[4] + x[1] * (t[3] + 1), t[5]
        + t[1] + t[4], t[6] + t[2] * t[4] + t[2] + t[4], t[7] + t[6] *
        (t[0] + 1), t[8] + x[2] * t[7] + x[2] + t[7], t[9] + t[5] +
        t[8], t[10] + x[4] * (t[3] + 1), t[11] + t[10] * x[1] + t[10]
        + x[1], t[12] + x[3] * t[11], y[0] + t[9] + (t[12] + 1), t[14]
        + t[3] + y[0], t[15] + t[14] * (x[1] + 1), t[16] + t[1] +
        t[15], t[17] + x[5] * (t[3] + 1), t[18] + t[5] + t[10], t[19]
        + x[1] * t[18], t[20] + t[17] + t[19], t[21] + x[2] * t[20],
        t[22] + t[16] + t[21], t[23] + x[4] + x[1], t[24] + t[23] *
        (t[7] + 1), t[25] + t[5] * x[0] + t[5] + x[0], t[26] + t[25] +
        x[1], t[27] + x[2] * (t[26] + 1), t[28] + t[24] + t[27], t[29]
        + x[3] * t[28] + x[3] + t[28], y[2] + t[22] + t[29], t[31] +
        t[17] * t[24] + t[17] + t[24], t[32] + t[31] + t[9], t[33] +
        t[26] * t[19] + t[26] + t[19], t[34] + x[2] * t[33], t[35] +
        t[32] + t[34], t[36] + t[23] * t[33], t[37] + t[11] * (t[36] +
        1), t[38] + x[3] * t[37] + x[3] + t[37], y[3] + t[35] + (t[38]
        + 1), t[40] + x[1] + t[1], t[41] + t[40] * (t[32] + 1), t[42]
        + t[41] + t[28], t[43] + x[2] * (t[42] + 1), t[44] + t[40] +
        t[43], t[45] + t[2] * t[19] + t[2] + t[19], t[46] + x[2] *
        t[2], t[47] + t[45] + t[46], t[48] + x[3] * (t[47] + 1), y[1]
        + t[44] + (t[48] + 1), ]


def s3opns(x, y, t):
    r"""
    Return sparse quadratic equations for the S-Box S3 in the input
    variables $x_0,...,x_5$, the output variables $y_0,...,y_3$ and 49
    intermediate variables $t_i$.

    These equations are based on

       http://www.darkside.com.au/bitslice/nonstd.c

    INPUT:
        x -- a list of length 6
        y -- a list of length 4
        t -- a list of length 49
    """
    t = list(t)

    # we add these silly entries, to have uniform notation with
    # nonstd.c

    t.insert(14,None)
    t.insert(29,None)
    t.insert(41,None)
    t.insert(52,None)
    
    return [ t[0] + x[1] + x[2], t[1] + t[0] + x[5], t[2] + x[1] *
        t[1], t[3] + x[4] * t[2] + x[4] + t[2], t[4] + t[1] + t[3],
        t[5] + x[2] + t[2], t[6] + t[5] * (x[4] + 1), t[7] + x[0] *
        t[6] + x[0] + t[6], t[8] + t[4] + t[7], t[9] + x[5] * (t[2] +
        1), t[10] + t[9] + x[4], t[11] + x[0] * t[10], t[12] + x[4] +
        t[11], t[13] + x[3] * t[12] + x[3] + t[12], y[3] + t[8] +
        t[13], t[15] + x[2] * x[5], t[16] + t[15] * t[2] + t[15] +
        t[2], t[17] + t[16] + x[4], t[18] + t[1] * (t[6] + 1), t[19] +
        t[18] + t[15], t[20] + x[0] * t[19] + x[0] + t[19], t[21] +
        t[17] + t[20], t[22] + x[1] * t[6] + x[1] + t[6], t[23] +
        t[22] + t[3], t[24] + t[10] * t[18] + t[10] + t[18], t[25] +
        t[24] + t[16], t[26] + x[0] * t[25] + x[0] + t[25], t[27] +
        t[23] + t[26], t[28] + x[3] * (t[27] + 1), y[2] + t[21] +
        (t[28] + 1), t[30] + x[2] * x[4], t[31] + t[30] + t[1], t[32]
        + t[6] * (x[2] + 1), t[33] + x[0] * t[32] + x[0] + t[32],
        t[34] + t[31] + t[33], t[35] + t[9] * t[25] + t[9] + t[25],
        t[36] + x[5] + t[16], t[37] + t[36] * (t[4] + 1), t[38] + x[0]
        * t[37], t[39] + t[35] + t[38], t[40] + x[3] * t[39], y[1] +
        t[34] + t[40], t[42] + x[1] * t[18] + x[1] + t[18], t[43] +
        t[42] + t[17], t[44] + x[5] * y[3], t[45] + t[44] + t[5],
        t[46] + t[45] * (x[0] + 1), t[47] + t[43] + t[46], t[48] +
        y[1] * (t[22] + 1), t[49] + x[0] * t[48] + x[0] + t[48], t[50]
        + t[46] + t[49], t[51] + x[3] * t[50], y[0] + t[47] + (t[51] +
        1), ]

def s4opns(x, y, t):
    r"""
    Return sparse quadratic equations for the S-Box S4 in the input
    variables $x_0,...,x_5$, the output variables $y_0,...,y_3$ and 35
    intermediate variables $t_i$.

    These equations are based on

       http://www.darkside.com.au/bitslice/nonstd.c

    INPUT:
        x -- a list of length 6
        y -- a list of length 4
        t -- a list of length 35
    """
    t = list(t)

    # we add these silly entries, to have uniform notation with
    # nonstd.c

    t.insert(21, None)
    t.insert(23, None)
    t.insert(36, None)
    t.insert(38, None)

    return [ t[0] + x[0] * x[2] + x[0] + x[2], t[1] + x[4] * t[0],
             t[2] + x[0] + t[1], t[3] + x[1] * x[2] + x[1] + x[2],
             t[4] + t[2] + t[3], t[5] + x[2] * (x[0] + 1), t[6] + t[5]
             * t[2] + t[5] + t[2], t[7] + x[1] * t[6], t[8] + x[4] +
             t[7], t[9] + x[3] * t[8], t[10] + t[4] + t[9], t[11] +
             x[2] + t[1], t[12] + x[1] * (t[11] + 1), t[13] + t[6] +
             t[12], t[14] + t[11] * t[2] + t[11] + t[2], t[15] + x[2]
             + x[4], t[16] + t[15] * (x[1] + 1), t[17] + t[14] +
             t[16], t[18] + x[3] * t[17] + x[3] + t[17], t[19] + t[13]
             + t[18], t[20] + x[5] * t[19] + x[5] + t[19], y[0] +
             t[10] + t[20], t[22] + x[5] * t[19], y[1] + t[22] +
             (t[10] + 1), t[24] + x[1] * t[8], t[25] + t[24] + t[14],
             t[26] + x[2] + t[7], t[27] + t[26] + t[16], t[28] + x[3]
             * (t[27] + 1), t[29] + t[25] + t[28], t[30] + t[10] +
             t[29], t[31] + x[1] * (t[30] + 1), t[32] + y[0] + t[31],
             t[33] + t[30] * (x[3] + 1), t[34] + t[32] + t[33], t[35]
             + x[5] * t[34] + x[5] + t[34], y[2] + t[29] + (t[35] +
             1), t[37] + t[22] + t[34], y[3] + t[37] + y[2], ]

def s5opns(x, y, t):
    r"""
    Return sparse quadratic equations for the S-Box S5 in the input variables
    $x_0,...,x_5$, the output variables $y_0,...,y_3$ and 52
    intermediate variables $t_i$.

    These equations are based on

       http://www.darkside.com.au/bitslice/nonstd.c

    INPUT:
        x -- a list of length 6
        y -- a list of length 4
        t -- a list of length 52
    """
    t = list(t)

    # we add these silly entries, to have uniform notation with
    # nonstd.c

    t.insert(20,None)
    t.insert(29,None)
    t.insert(41,None)
    t.insert(55,None)

    return [ t[0] + x[2] * (x[3] + 1), t[1] + t[0] + x[0], t[2] + x[0]
        * (x[2] + 1), t[3] + x[5] * t[2] + x[5] + t[2], t[4] + t[1] +
        t[3], t[5] + x[3] + x[0], t[6] + t[5] * t[0] + t[5] + t[0],
        t[7] + t[6] * (x[5] + 1), t[8] + x[2] + t[7], t[9] + x[4] *
        t[8] + x[4] + t[8], t[10] + t[4] + t[9], t[11] + x[2] * t[6],
        t[12] + t[11] + x[3], t[13] + t[12] * (t[2] + 1), t[14] + x[3]
        + t[2], t[15] + x[5] * t[14] + x[5] + t[14], t[16] + t[13] +
        t[15], t[17] + x[4] * t[16] + x[4] + t[16], t[18] + t[12] +
        t[17], t[19] + t[18] * (x[1] + 1), y[3] + t[10] + t[19], t[21]
        + x[3] * t[3], t[22] + t[21] + t[16], t[23] + x[0] + t[8],
        t[24] + t[1] * t[23], t[25] + x[4] * (t[24] + 1), t[26] +
        t[22] + t[25], t[27] + x[3] * t[23] + x[3] + t[23], t[28] +
        t[27] * (x[1] + 1), y[1] + t[26] + t[28], t[30] + t[16] *
        t[4], t[31] + t[6] * (t[30] + 1), t[32] + t[7] * (x[3] + 1),
        t[33] + t[32] + x[2], t[34] + x[4] * t[33], t[35] + t[31] +
        t[34], t[36] + t[12] * t[15] + t[12] + t[15], t[37] + t[8] +
        t[30], t[38] + x[4] * t[37] + x[4] + t[37], t[39] + t[36] +
        t[38], t[40] + x[1] * t[39] + x[1] + t[39], y[2] + t[35] +
        (t[40] + 1), t[42] + t[18] * (t[31] + 1), t[43] + t[42] +
        t[23], t[44] + t[26] * t[42] + t[26] + t[42], t[45] + t[44] +
        t[5], t[46] + x[4] * (t[45] + 1), t[47] + t[43] + t[46], t[48]
        + t[5] * t[37], t[49] + t[48] + t[33], t[50] + y[3] + t[37],
        t[51] + t[27] * (t[50] + 1), t[52] + x[4] * t[51], t[53] +
        t[49] + t[52], t[54] + x[1] * t[53] + x[1] + t[53], y[0] +
        t[47] + t[54], ]

def s6opns(x, y, t):
    r"""
    Return sparse quadratic equations for the S-Box S6 in the input
    variables $x_0,...,x_5$, the output variables $y_0,...,y_3$ and 49
    intermediate variables $t_i$.

    These equations are based on

       http://www.darkside.com.au/bitslice/nonstd.c

    INPUT:
        x -- a list of length 6
        y -- a list of length 4
        t -- a list of length 49
    """
    t = list(t)
    t.insert(18 ,None)
    t.insert(30 ,None)
    t.insert(39 ,None)
    t.insert(52 ,None)

    # we add these silly entries, to have uniform notation with
    # nonstd.c

    return [ t[0] + x[4] + x[0], t[1] + t[0] + x[5], t[2] + x[0] *
        x[5], t[3] + t[2] * (x[4] + 1), t[4] + x[3] * (t[3] + 1), t[5]
        + t[1] + t[4], t[6] + x[5] + t[2], t[7] + t[3] * t[6] + t[3] +
        t[6], t[8] + t[7] * (x[3] + 1), t[9] + t[6] + t[8], t[10] +
        x[1] * t[9], t[11] + t[5] + t[10], t[12] + x[5] * t[5] + x[5]
        + t[5], t[13] + t[12] * (x[4] + 1), t[14] + t[3] * t[9] + t[3]
        + t[9], t[15] + x[1] * (t[14] + 1), t[16] + t[13] + t[15],
        t[17] + t[16] * (x[2] + 1), y[0] + t[11] + (t[17] + 1), t[19]
        + y[0] * (t[0] + 1), t[20] + t[19] + t[14], t[21] + x[5] *
        (t[20] + 1), t[22] + t[21] + t[5], t[23] + x[1] * (t[22] + 1),
        t[24] + t[20] + t[23], t[25] + x[4] * x[5] + x[4] + x[5],
        t[26] + t[25] * (t[0] + 1), t[27] + x[1] * (t[23] + 1), t[28]
        + t[26] + t[27], t[29] + x[2] * (t[28] + 1), y[3] + t[24] +
        (t[29] + 1), t[31] + t[2] + t[5], t[32] + t[31] * (t[9] + 1),
        t[33] + x[5] + t[24], t[34] + x[4] * (t[33] + 1), t[35] + x[1]
        * (t[34] + 1), t[36] + t[32] + t[35], t[37] + t[20] * (x[4] +
        1), t[38] + x[2] * t[37] + x[2] + t[37], y[2] + t[36] + (t[38]
        + 1), t[40] + t[34] * t[1] + t[34] + t[1], t[41] + x[4] *
        t[6], t[42] + x[3] * (t[41] + 1), t[43] + x[1] * t[42] + x[1]
        + t[42], t[44] + t[40] + t[43], t[45] + t[22] * t[34] + t[22]
        + t[34], t[46] + t[45] + t[4], t[47] + t[25] * t[32], t[48] +
        t[47] + t[1], t[49] + x[1] * t[48], t[50] + t[46] + t[49],
        t[51] + x[2] * (t[50] + 1), y[1] + t[44] + (t[51] + 1), ]

def s7opns(x, y, t):
    r"""
    Return sparse quadratic equations for the S-Box S7 in the input
    variables $x_0,...,x_5$, the output variables $y_0,...,y_3$ and 47
    intermediate variables $t_i$.

    These equations are based on

       http://www.darkside.com.au/bitslice/nonstd.c

    INPUT:
        x -- a list of length 6
        y -- a list of length 4
        t -- a list of length 47
    """
    t = list(t)

    # we add these silly entries, to have uniform notation with
    # nonstd.c

    t.insert(18 ,None)
    t.insert(31 ,None)
    t.insert(40 ,None)
    t.insert(50 ,None)

    return [ t[0] + x[1] * x[3], t[1] + t[0] + x[4], t[2] + x[3] *
        t[1], t[3] + t[2] + x[1], t[4] + x[2] * (t[3] + 1), t[5] +
        t[1] + t[4], t[6] + x[2] + t[4], t[7] + x[5] * (t[6] + 1),
        t[8] + t[5] + t[7], t[9] + x[1] * x[3] + x[1] + x[3], t[10] +
        t[9] * x[4] + t[9] + x[4], t[11] + x[4] * (x[1] + 1), t[12] +
        x[2] * t[11] + x[2] + t[11], t[13] + t[10] + t[12], t[14] +
        t[2] + t[5], t[15] + x[5] * t[14] + x[5] + t[14], t[16] +
        t[13] + t[15], t[17] + x[0] * t[16], y[0] + t[8] + t[17],
        t[19] + x[3] * (x[2] + 1), t[20] + x[1] * (t[19] + 1), t[21] +
        x[5] * t[20], t[22] + t[8] + t[21], t[23] + x[3] + t[3], t[24]
        + x[2] * t[2] + x[2] + t[2], t[25] + t[23] + t[24], t[26] +
        x[2] + t[2], t[27] + t[26] * x[1], t[28] + x[5] * (t[27] + 1),
        t[29] + t[25] + t[28], t[30] + x[0] * t[29] + x[0] + t[29],
        y[1] + t[22] + (t[30] + 1), t[32] + t[6] + t[29], t[33] + x[1]
        * t[23] + x[1] + t[23], t[34] + t[33] + y[0], t[35] + t[34] *
        (x[5] + 1), t[36] + t[32] + t[35], t[37] + t[25] * (x[2] + 1),
        t[38] + t[37] * t[29] + t[37] + t[29], t[39] + t[38] * (x[0] +
        1), y[2] + t[36] + t[39], t[41] + x[4] * t[19] + x[4] + t[19],
        t[42] + t[41] + t[32], t[43] + x[1] + t[14], t[44] + t[23] *
        (t[43] + 1), t[45] + x[5] * t[44], t[46] + t[42] + t[45],
        t[47] + x[2] * t[21], t[48] + t[47] + t[45], t[49] + x[0] *
        t[48] + x[0] + t[48], y[3] + t[46] + t[49], ]

def s8opns(x, y, t):
    r"""
    Return sparse quadratic equations for the S-Box S8 in the input
    variables $x_0,...,x_5$, the output variables $y_0,...,y_3$ and 46
    intermediate variables $t_i$.

    These equations are based on

       http://www.darkside.com.au/bitslice/nonstd.c

    INPUT:
        x -- a list of length 6
        y -- a list of length 4
        t -- a list of length 46
    """
    t = list(t)

    # we add these silly entries, to have uniform notation with
    # nonstd.c

    t.insert(19 ,None)
    t.insert(31 ,None)
    t.insert(41 ,None)
    t.insert(49 ,None)

    return [ t[0] + x[2] + x[0], t[1] + x[0] * (x[2] + 1), t[2] + t[1]
        + x[3], t[3] + x[4] * t[2] + x[4] + t[2], t[4] + t[0] + t[3],
        t[5] + t[4] * (x[0] + 1), t[6] + t[5] + x[2], t[7] + t[6] *
        (x[4] + 1), t[8] + x[3] + t[7], t[9] + x[1] * (t[8] + 1),
        t[10] + t[4] + t[9], t[11] + t[5] * x[3] + t[5] + x[3], t[12]
        + t[11] + t[0], t[13] + t[12] + x[4], t[14] + t[2] * (t[13] +
        1), t[15] + t[14] + t[6], t[16] + x[1] * (t[15] + 1), t[17] +
        t[13] + t[16], t[18] + x[5] * t[17] + x[5] + t[17], y[0] +
        t[10] + (t[18] + 1), t[20] + t[4] * x[4] + t[4] + x[4], t[21]
        + t[20] + t[2], t[22] + t[10] * (x[3] + 1), t[23] + x[1] *
        (t[22] + 1), t[24] + t[21] + t[23], t[25] + x[0] * t[20],
        t[26] + x[4] * t[1], t[27] + t[26] + t[22], t[28] + x[1] *
        t[27], t[29] + t[25] + t[28], t[30] + t[29] * (x[5] + 1), y[2]
        + t[24] + t[30], t[32] + x[2] * (t[15] + 1), t[33] + t[8] *
        t[32] + t[8] + t[32], t[34] + x[1] * t[5] + x[1] + t[5], t[35]
        + t[33] + t[34], t[36] + t[1] * (t[13] + 1), t[37] + t[21] *
        y[2] + t[21] + y[2], t[38] + x[1] * (t[37] + 1), t[39] + t[36]
        + t[38], t[40] + x[5] * t[39] + x[5] + t[39], y[1] + t[35] +
        (t[40] + 1), t[42] + t[0] * (x[4] + 1), t[43] + t[42] * x[3] +
        t[42] + x[3], t[44] + x[2] + x[4], t[45] + t[44] + t[36],
        t[46] + t[45] * (x[1] + 1), t[47] + t[43] + t[46], t[48] +
        x[5] * t[47], y[3] + t[10] + (t[48] + 1), ]


def check_sbox_consistency(sbox_eq="opns"):
    """
    Check consistency of S-Box equations.

    EXAMPLE:
        sage: execfile('des.py')
        sage: check_sbox_consistency('cubic')
        S1 passed
        S2 passed
        S3 passed
        S4 passed
        S5 passed
        S6 passed
        S7 passed
        S8 passed

        sage: check_sbox_consistency('opns') # long time 
        S1 passed
        S2 passed
        S3 passed
        S4 passed
        S5 passed
        S6 passed
        S7 passed
        S8 passed
    """
    load_sXopns_gb()
    load_sXcubic()

    des = DES(Nr=1)
    t = des.vars("t",1)

    m = 6
    sopns = [eval("s%d%s"%(i+1,sbox_eq)) for i in range(8)]

    for j in range(8):
        S = DESSBox(j+1)
        for i in range(1<<m):
            i = S.to_bits(i,m)
            if 'cubic':
                F = sopns[j](i , S(i))
            else:
                F = sopns[j](i , S(i) ,t)
            F = [f for f in F if f!=0]
            if F == []:
                continue
            if Ideal(F).groebner_basis() == [1]:
                raise TypeError, "S%d failed"%(j+1)
        print "S%d passed"%(j+1)


def load_sXopns_gb():
    """
    """
    def _precomp_sbox_des_gb(f):
        """
        Given a function f which returns a list of polynomials defining an
        S-Box S1,...,S8 this function returns a list of strings
        representing a degrevlex Groebner basis for this S-Box.

        INPUT:
            f -- a function f(x,y,t) returning a list of polynomials
        """
        Y = ["y%d"%i for i in xrange(4) ]
        T = ["t%d"%i for i in xrange(56)]
        X = ["x%d"%i for i in xrange(6) ]
        var_names = Y + T + X
        P = BooleanPolynomialRing(len(var_names),var_names, order="degrevlex")
        Y = [P("y%d"%i) for i in xrange(4) ]
        T = [P("t%d"%i) for i in xrange(56)]
        X = [P("x%d"%i) for i in xrange(6) ]
        b = f(X,Y,T)
        gb = P.ideal(b).groebner_basis(red_tail=True)
        l = []
        for f in gb:
            l.append(re.sub("([a-z])([0-9]+)","\\1[\\2]",str(f)))
        return l

    fn = tmp_filename()
    fh = open(fn, "w")
    S = []
    S.append( _precomp_sbox_des_gb(s1opns) )
    S.append( _precomp_sbox_des_gb(s2opns) ) 
    S.append( _precomp_sbox_des_gb(s3opns) )
    S.append( _precomp_sbox_des_gb(s4opns) )
    S.append( _precomp_sbox_des_gb(s5opns) )
    S.append( _precomp_sbox_des_gb(s6opns) )
    S.append( _precomp_sbox_des_gb(s7opns) )
    S.append( _precomp_sbox_des_gb(s8opns) )

    for i in range(len(S)):
        fh.write("def s%dopns_gb(x, y, t):\n"%(i+1))
        fh.write("    return [\n")
        for f in S[i]:
            fh.write("        "+f+",\n")
        fh.write("           ]\n")
        fh.write("\n")
    fh.close()
    execfile(fn,globals(),globals())

def load_sXcubic():
    """
    """
    fn = tmp_filename()
    fh = open(fn, "w")

    for i in range(8):
        F = DESSBox(i+1).polynomials(degree=3)
        fh.write("def s%dcubic(x, y):\n"%(i+1))
        fh.write("    x0,x1,x2,x3,x4,x5 = x\n")
        fh.write("    y0,y1,y2,y3 = y\n")
        fh.write("    return [\n")
        for f in F:
            fh.write("        "+str(f)+",\n")
        fh.write("           ]\n")
        fh.write("\n")
    fh.close()
    execfile(fn,globals(),globals())

        
class DifferentialCharacteristicIterator(SageObject):
    def __init__(self, cipher ):
        """
        Abstract class for difference characteristic iterator.

        INPUT:
            cipher -- cipher instance
        """
        self.cipher = cipher
        self.characteristic = []
        self.__it = 0
    def __iter__(self):
        return self

    def next(self):
        if self.__it < len(self.characteristic):
            it = self.__it
            self.__it += 1
        else:
            raise StopIteration
            #it = 0
            #self.__it = 1

        return self.characteristic[it]

    def __getitem__(self, i):
        return self.characteristic[i]

def from_str(D):
    """
    """
    d = []
    for s in D:
        if not re.match("[a-fA-Z0-9]",s):
            continue
        bits = ZZ(s,16).digits(2,padto=4)
        bits = map(GF(2), bits[::-1])
        d += bits

    return d

def to_str(D):
    d = ''
    for i in range(0,len(D),8):
        bits = ZZ(list(reversed(map(int, D[i:i+8]))),2)
        d += "%02x "%bits
    return d

class DESCharacteristic(DifferentialCharacteristicIterator):
    def __init__(self, des):
        """
        """
        DifferentialCharacteristicIterator.__init__(self, des)
        c = [from_str("19 60 00 00 00 00 00 00"),
             from_str("00 00 00 00 19 60 00 00"),
             from_str("19 60 00 00 00 00 00 00"),
             from_str("00 00 00 00 19 60 00 00")]
        self.characteristic = c
     
    def _repr_(self):
        return "DES Characteristic Iterator"

#################################


class LinearStructures:
    def __init__(self, Nr, des):
        self.Nr = Nr
        self.des = des(Nr=Nr)
        self.R = 1
        self.M = Nr/2 
        self.T = Nr
        self.V = [range(4*i,4*(i+1)) for i in range(8)]

    def f(self, s):
        _F = [[1,2,3,4,5,7],
	      [0,2,3,4,6,7],
	      [1,3,4,5,6,7],
	      [0,2,4,5,6,7],
	      [0,1,2,3,5,6],
	      [0,1,2,4,6,7],
	      [0,1,2,3,5,7],
	      [0,1,3,4,5,6]]
        return reduce(union, [_F[e] for e in s], set([]))

    def w(self, i,t):
        v = [0 for _ in range(56)]
        v[t] = 1
        v = self.des.L(i-1, v)
        if max(v) == 0:
            return set()
        else:
            assert(v.count(1) == 1)
            return set([v.index(1)//6])

    def v(self, i, t):
        if i == (self.R-1):
            return set()
        if i == self.R:
            return set()
        return reduce(union, [self.v(i-2,t), self.f(self.v(i-1,t)), self.w(i-1,t)],set())

    def vbar(self, i, t):
        if i == (self.T+1):
            return set()
        if i == self.T:
            return set()
        return reduce(union, [self.vbar(i+2,t), self.f(self.vbar(i+1,t)), self.w(i+1,t)],set())

    def x(self, i, t):
        return set(reduce(union, [self.V[j] for j in self.v(i,t)], set()))

    def xbar(self, i, t):
        return set(reduce(union, [self.V[j] for j in self.vbar(i,t)], set()))

    def q(self, t):
        if self.des.__class__ == DES:
            A = map(lambda z: z+32, reduce(union, [self.x(self.M+1,t), self.xbar(self.M+1,t)], set()))
            return reduce(union, [self.x(self.M,t), self.xbar(self.M,t), A], set([]))
        elif self.des.__class__ == DESDC:
            A = map(lambda z: z+32, self.xbar(self.M+1,t))
            return reduce(union, [self.xbar(self.M,t), A], set([]))

    def __call__(self):
        return  [[t for t in range(56) if j not in self.q(t)] for j in range(64)]