Discrete antiderivatives for functions over Fpn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathop {{\mathbb {F}}}_p^n$$\end{document}

In the design of cryptographic functions, the properties of their discrete derivatives have to be carefully considered, as many cryptographic attacks exploit these properties. One can therefore attempt to first construct derivatives with the desired properties and then recover the function itself. Recently Suder developed an algorithm for reconstructing a function (also called antiderivative) over the finite field F2n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathop {{\mathbb {F}}}_{2^n}$$\end{document} given its discrete derivatives in up to n linearly independent directions. Pasalic et al. also presented an algorithm for determining a function over Fpn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathop {{\mathbb {F}}}_{p^n}$$\end{document} given one of its derivatives. Both algorithms involve solving a pn×pn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^n \times p^n$$\end{document} system of linear equations; the functions are represented as univariate polynomials over Fpn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathop {{\mathbb {F}}}_{p^n}$$\end{document}. We show that this apparently high computational complexity is not intrinsic to the problem, but rather a consequence of the representation used. We describe a simpler algorithm, with quasilinear complexity, provided we work with a different representation of the functions. Namely they are polynomials in n variables over Fp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathop {{\mathbb {F}}}_{p}$$\end{document} in algebraic normal form (for p>2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p>2$$\end{document}, additionally, we need to use the falling factorial polynomial basis) and the directions of the derivatives are the canonical basis of Fpn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathop {{\mathbb {F}}}_{p}^n$$\end{document}. Algorithms for other representations (the directions of the derivatives not being the canonical basis vectors or the univariate polynomials over Fpn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathop {{\mathbb {F}}}_{p^n}$$\end{document} mentioned above) can be obtained by combining our algorithm with converting between representations. However, the complexity of these conversions is, in the worst case, exponential. As an application, we develop a method for constructing new quadratic PN (Perfect Nonlinear) functions. We use an approach similar to the one of Suder, who used antiderivatives to give an alternative formulation of the methods of Weng et al. and Yu et al. for searching for new quadratic APN (Almost Perfect Nonlinear) functions.


Introduction
A function f with an n-bit input and n-bit output can be represented in various ways.In cryptographic applications, it is often represented as a univariate polynomial of degree up to 2 n − 1 over F2 n , the finite field with 2 n elements.More generally, functions over Fp n can be represented as univariate polynomials of degree up to p n − 1 over Fp n .Another representation, which is the one we will mostly use here, is as a function f : F n p → F m p (which includes the more general possibility of m different from n), f = ( f 1 , . . ., f m ) with each component function f i : F n p → Fp being a polynomial in n variables over Fp, of degree at most p − 1 in each variable; the polynomials are usually represented in ANF (Algebraic Normal Form), i.e. as a sum of monomials, each monomial being a constant multiplied by the product of powers of some of the variables.
Many cryptographic attacks (e. g. differential attacks, higher order differential attacks) exploit properties of the discrete derivatives of the functions.The discrete derivative (also called simply derivative when no danger of confusion exists) of a function f in the direction a is defined as the function D a f (x) = f (x + a) − f (x).
Recently Suder [4] considered a natural problem: reconstruct a function f over F2 n from its discrete derivatives in n (or fewer) linearly independent directions.This operation can be called antiderivative, or integration, in analogy to its continuous counterpart.The proposed algorithm involves solving a 2 n × 2 n system of linear equations; it is stated [4,Sect. 3.1] that the matrices are sparse and easy to compute with, but no complexity analysis is presented.Pasalic et al. [2] consider the problem of determining a function over Fp n given one of its derivatives; similarly, their algorithm involves solving a p n × p n system of linear equations.Earlier, Xiong et al. [6] determined some necessary conditions for a function to be the discrete derivative of another function.In all these papers the functions are represented as univariate polynomials over Fp n .
While the representation as univariate polynomials over Fp n is very useful for some purposes, it turns out that for computing the antiderivative the second representation above is much more convenient; when p > 2, additionally, we have to modify the algebraic normal form by using falling factorials instead of powers of variables, i.e. using x d = x(x − 1) • • • (x − (d − 1)) instead of x d ; it is well known that such a representation is more convenient for discrete differentiation.When the directions of the derivatives are vectors in the canonical basis and the input polynomials are represented as the list of their monomials, we present an algorithm which is quasilinear in the size of the input polynomials (see Sect. 3).Hence the apparently exponential computational complexity of the previous algorithms is not intrinsic to the problem, but rather a consequence of the representation used.Our algorithm is similar to the one used in calculus for reconstructing a function from its gradient/partial derivatives (see for example the textbook [3,Section 15.8]).This is possible (despite obvious differences between the classical derivative from calculus and the discrete derivative) because for multivariate polynomial functions over Fp represented as described above, the discrete derivative behaves similarly to the formal derivative.
Algorithms for other representations can be obtained by combining our algorithm with converting between representations.When the functions are represented as multivariate polynomials but the directions of the derivatives are arbitrary, our algorithm can be combined with linear changes of coordinates (using n × n matrices) which transform those directions into the canonical basis directions; changes of coordinates can, in the worst case, cause the number of monomials in the algebraic normal form representation to increase exponentially.
If the functions are given as univariate polynomials over Fp n , as considered in the work of Suder [4] and of Pasalic et al. [2], we can obtain an alternative to their antiderivative algorithm; namely we transform the function into its multivariate polynomial representation, apply our quasilinear antiderivative algorithm, and then transform the result back to the univariate representation (see Sect. 4).Unfortunately, the conversion between representations has exponential worst case complexity.
As applications of the antiderivative construction, Pasalic et al. [2] prove results regarding the degree of planar functions, whereas Suder [4] prove results regarding the construction of APN (Almost Perfect Nonlinear) functions.Both these classes of functions are defined by the properties of their derivatives.Differential attacks on cryptographic functions exploit the situation where there is a direction a such that the derivative D a f (x) takes a particular value significantly more often than other values.To withstand such attacks, we can require that D a f is bijective (only possible when p ≥ 3) or 2-to-1 (for p = 2).In the first situation the function f is called PN(Perfect Nonlinear), or planar, and in the second situation it is called APN (Almost Perfect Nonlinear).
Weng et al. [5] and of Yu et al. [7] independently presented a method of constructing quadratic APN functions by constructing n × n matrices over F2 n with certain properties.Using this characterisation and a computer search they obtain numerous new quadratic APN functions.Suder [4] applies his method to obtain an alternative approach to the methods in [5,7].He also suggested that antiderivatives could be used more generally, for other degrees and for other properties of the derivatives.Following these lines, we revisit and expand such applications in Sect. 5.
We give an analogue of [7, Theorem 1] for characterising quadratic APN functions in multivariate ANF representation.More interestingly, we present a similar result concerning quadratic PN functions.This leads to a method of constructing potentially new PN functions either from scratch (for example, we computed exhaustively all PN functions for n = 3 and p ≤ 7 up to extended affine equivalence) or by modifying known PN functions in a manner similar to Yu et al. [7].A more extensive investigation of this direction will be a topic of further research.

Preliminaries
We recall definitions and known results needed for the rest of the paper.The finite field with p n elements will be denoted Fp n , where p is a prime.When K is any field, we denote the canonical basis of the vector space K n by e 1 , . . ., e n , i.e. the vector e i has a 1 in position i and zeroes everywhere else; we denote by 0 the all-zero vector.In general we will denote vectors by boldface font.
Boolean functions in n variables can be viewed as functions f : F n 2 → F2.More generally, we can consider functions f : F n p → Fp.It is well known that any such function can be uniquely represented in its ANF (Algebraic Normal Form) i.e. as a polynomial function described by a polynomial in Fp[x1, . . ., x n ] of degree at most p − 1 in each variable, f = t c t t where c t ∈ Fp and t ranges over all the power products x i 1 1 x i 2 2 • • • x i n n with i 1 , . . ., i n ∈ {0, 1, . . ., p − 1}.We denote by deg( f ) the degree of this polynomial (also called algebraic degree of f ) and by deg x i ( f ) the degree in the variable x i .By convention, the degree of the zero polynomial is −1.
Recall that the falling factorial is defined as ), for d ≥ 0 an integer (with x 0 = 1).Note that if we view this expression as a polynomial in x, then deg(x d ) = deg(x d ) = d.It is well known that polynomials (over any field) can be represented in a basis consisting of falling factorials.In our case, multivariate polynomials of degree at most p − 1 in each variable can be represented as f = t c t t where c t ∈ Fp and t ranges over all the terms x 1 i 1 x 2 i 2 • • • x n i n with i 1 , . . ., i n ∈ {0, 1, . . ., p − 1}.We will call this representation the falling factorial ANF.Obviously, for p = 2 this is the same as the ANF.
More generally, we consider functions f : ).When we evaluate the polynomial, each variable takes values in Fp so the coefficients will only be multiplied by a scalar.By abuse of terminology we will call this the ANF of f , or the falling factorial ANF of f , respectively.This representation also allows us to define the degree of f as the highest degree of a term t for which c t = 0 (it will also equal the highest degree of the components f 1 , . . ., f m ).
We can write this compactly as a polynomial Alternatively, functions with n bits input and n bits output are often viewed as f : F2 n → F2 n represented as univariate polynomials over F2 n of degree up to 2 n − 1.In Sect. 4 we will consider this representation and the conversion between the representations.
We recall the definition of (discrete) derivative/ differentiation: Definition 1 Let K be a field and f : The differentiation operator in the direction of a associates to each function f its discrete derivative D a f : The inverse of the (discrete) differentiation has been considered in [4,6]: antiderivative (or integral) of g 1 , . . ., g k .The set of such functions f will be denoted Note that an antiderivative might not exist, and when it does exist it is not unique, e.g. the addition of a constant preserves the property.Other antiderivatives exist when k < n; this will be discussed later.
We recall a few useful well-known properties of the discrete derivative.They are straightforward to prove; for point (v) see [1].
If M is an invertible n × n matrix over Fp and g is obtained from f by the linear invertible change of coordinates described by M, i.e. g(x For any c ∈ Fp, applying inductively Proposition 1(iii) for b = ia, we obtain Furthermore, we can express the derivative in an arbitrary direction in terms of the derivatives in the canonical directions: In particular for a = (a 1 , . . ., a n ) = n i=1 a i e i we have For p = 2 the two equations above become: and a s e s .

Antiderivatives for functions over F n p
We consider the problem of retrieving a function f : F n p → F m p when we know its derivatives g 1 , . . ., g k in k ≤ n linearly independent directions a 1 , . . ., a k , i.e. we want to compute the antiderivative.Note that it suffices to consider directions which are linearly independent; if we were given, additionally, the derivative g in a direction a which is a linear combination of a 1 , . . ., a k , that would not provide any additional information about f , as g could be computed from g 1 , . . ., g k using Proposition 2 (1).
In order for a solution to exist, the g i have to satisfy certain compatibility conditions (see [4,6])).The first set of conditions is which is necessary in a field of characteristic p as g i = D a i f (see Proposition 1(viii)).The second set of conditions is as both sides should be equal to D a i D a j f = D a j D a i f , see Proposition 1(vii).The fact that these conditions are not only necessary, but also sufficient, will follow from the construction of f below.

Antiderivatives in the directions of the canonical basis
First we consider the case when the directions are vectors in the canonical basis a 1 = e 1 , . . ., a k = e k .In this case the compatibility condition (3) means that deg x i (g i ) ≤ p − 2. (see Proposition 1(vi)).Algorithm 1 constructs the antiderivative when p = 2 and Algorithm 2 for arbitrary p.These algorithms are similar to the method used in calculus for reconstructing a function from its gradient (see for example the textbook [3]).
Proof For ease of notation, and without loss of generality, let us assume 1 = 1, . . ., k = k.Note that for p = 2 Algorithm 2 becomes Algorithm 1, so it suffices to prove the correctness of Algorithm 2.
For convenience we will use the following notation and we extend A i by linearity to any polynomial expressed in falling factorial ANF.Note that for any polynomial g we have D e i (A i (g)) = g and for any j = i we have D e j (A i (g)) = A i (D e j (g)).To prove both these facts it suffices to verify for monomials (which is left as an exercise, using Proposition 1(ii)), since both operators A i and D e j are linear.
Let us denote by f (0) the initial value of f and by f (i) the value of f at the end of the ith run of the outer for loop.We prove by induction on i that D e j f (i) = g j for all j ≤ i.Looking at the inner for loop we see that it computes f (i) = f (i−1) + A i (h).For the base case, i = 1, we have D e 1 f (1) = D e 1 (A 1 (h)) = h = g 1 .Now assume the statement holds for i − 1 and we prove it for i.We have: For any j < i, using the induction hypothesis and the compatibility conditions D e j g i = D e i g j we have: Let F be another function F such that D e i F = g i for all i = 1, . . ., k. Therefore D e i (F − f ) = 0 for all i.This happens if and only if F − f has degree zero in each of the variables x 1 , . . ., x k (see Proposition 1(vi)), i.e. it does not depend on those variables.

Remark 1
We treated the case p = 2 separately in Algorithm 1 no only because the formulation becomes simpler, but also because, unlike Algorithm 2, the polynomials for Algorithm 1 do not even need to be given in ANF.They could be given as polynomials which are not in normal form, or even as black box functions (such functions are often considered in attacks on symmetric stream ciphers; they model well the situation when the ANF is too large to be written explicitly).
For k = n we can use Algorithm 1 to obtain by induction an explicit formula for f : e i 1 ,...e i s−1 g i s . ( Keeping in mind that f has to satisfy g i = D e i f , this is equivalent to the following known result (see [1, Proposition 1]): ..e is f .
The complexity of Algorithms 1 and 2 depends on how we represent the polynomial functions.If we represent them in falling factorials ANF, with each polynomial represented as the set of its monomials, the algorithm becomes even simpler, see Algorithm 3.
appearing in its falling factorial ANF, ordered in lexicographical order of (d 1 , . . ., d n ) (or any other admissible monomial order).Then the time complexity of the algorithm is quasilinear in the size of its input.Namely if the total size of the k input polynomials g 1 , . . ., g k is s bits, the algorithm has complexity O(s log k).
Proof For ease of notation, and without loss of generality, let us assume 1 = 1, . . ., k = k.
It can be easily checked, using Proposition 1(ii) that D e i G i = g i .By the construction of G i , we know that all monomials of G i do contain x i .
Consider a term t = x 1 d 1 • • • x n d n which appears in at least one of the polynomials G 1 , . . ., G k .We prove that t appears in all the polynomials G i for which i has the property d i ≥ 1, and moreover it appears with the same coefficient.Without loss of generality, assume that t appears in G 1 and that G 2 , . . ., G s are the other polynomials G i for which i has the property d i ≥ 1.
Let c i be the coefficient of t in G i , for i = 1, . . ., s (at this moment we only know c 1 = 0, the other could be zero).
Let j ∈ {2, . . ., s}.The term On the other hand, D e 1 D e j G j = D e 1 g j and D e 1 D e j G 1 = D e j g 1 , so using the compatibility condition D e j g 1 = D e 1 g j we have that c j = c 1 .
We have therefore proved that taking the union of the sets k i=1 G i maintains the condition that all monomials in the set have different power products.
For the correctness, we need to prove D e i f = g i for i = 1, . . ., k; since we represent polynomials as the sets of their monomials, we will prove equality of sets.Fix an i.Let ct ∈ f For the complexity of the algorithm, note, firstly, that the sorted list of monomials in G i is obtained in linear time from the sorted list of monomials in g i , because all exponents are changed by increasing the exponent of x i by one, which preserves the lexicographical order.Secondly, computing the union of k sets, when each set is a sorted list, can be achieved by an algorithm similar to the merging phase of Mergesort, by log k rounds of merging pairs of lists (ensuring entries which appear in both lists are introduced only once in the merged list), and each merging round having linear complexity in the sizes of the merged lists.
One can verify that they satisfy the compatibility conditions.Applying Algorithm 3 we compute G i = x i g i for i = 1, 2, 3: Finally we collect all the monomials in x 1 g 2 , x 2 g 2 , x 3 g 3 to obtain f : One can verify that they satisfy the compatibility conditions.Using Algorithm 3 we compute: Finally computing the union of G 1 and G 2 we obtain f :

Antiderivatives in arbitrary directions
Next let us consider the general case where we are given the discrete derivatives g 1 , . . ., g k in k ≤ n arbitrary linearly independent directions a 1 , . . ., a k .The following Theorem presents two possible approaches to computing the antiderivative.
Theorem 3 Let a 1 , . . ., a k ∈ F n p with k ≤ n be linearly independent vectors and let g 1 , . . ., g k : a i ,...,a i g i = 0 and D a i g j = D a j g i for all i, j = 1, . . ., k.
Proof (i) This first approach involves a change of coordinates.Obviously Ae i = a i and A −1 a i = e i for all i.Let f be an element of the set on the right hand side of Eq. ( 6), i.e. f (x) = g(A −1 x) for some g ∈ Antiderivative(n, k, (g 1 (Ax), . . ., g k (Ax)), (e 1 , . . ., e k )).Proposition 1(iv) gives for all i hence f does belong to the set on the left side of the Eq. ( 6).Conversely, let f be an element of the set on the left hand side of Eq. ( 6).Let g(x) = f (Ax).Again, Proposition 1(iv) gives for all i hence g ∈ (Antiderivative(n, k, (g 1 (Ax), . . ., g k (Ax)), (e 1 , . . ., e k ))), so f does belong to the set on the right side of the Eq. ( 6).
(ii) This second approach consists of computing the derivatives in the canonical directions from the derivatives in arbitrary directions.Let f be an element of the set on the left hand side of Eq. ( 8).Therefore, D a i f (x) = g i (x) for all i.Using Proposition 2(1) and the fact that e i = n j=1 b i j a j : we have that D e i f = h i with h i as computed in Eq. ( 7).Hence f is indeed in the set on the right side of Eq. ( 8).Any other element f 1 in the set on the right hand side of Eq. (8) differs from f by a constant c ∈ F m p , see Algorithm 2 and Theorem 1.
so f 1 belongs to the set on the left hand side of Eq. (8).
Note that in both approaches in the Theorem above, the complexity of computing the intermediate polynomials g 1 (Ax), . . ., g k (Ax) or h 1 , . . ., h n may be exponential compared to the size of the input polynomials g 1 , . . ., g k .(For the first approach, note that a linear change of coordinates applied to a monomial of degree d can result in up to n d monomials; a similar problem appears when we compute h i in the second approach.)If the size of the output f is exponential compared to the size of the input, any algorithm for computing f in ANF would take exponential time.When the size of f is polynomial in the size of the input, it may still be possible that the computation of intermediate polynomials is exponential and therefore the whole algorithm is exponential.
In [4,Definition 8], an equivalence relation on the class of functions is defined, whereby for a fixed vector space V , two functions f 1 , f 2 are called differentially equivalent with respect to V if D v f 1 = D v f 2 for all v ∈ V .Obviously, is suffices to check whether D a i f 1 = D a i f 2 for some basis a 1 , . . ., a k of V .(see Proposition 2(1)).A characterisation of differential equivalence similar to [4, Proposition 2 and Eq. ( 6)] follows: Corollary 1 Let a 1 , . . ., a k ∈ F n p with k ≤ n be linearly independent vectors.Two functions f 1 , f 2 : F n p → F m p have the same derivatives in the directions a 1 , . . ., a k , i.e.
which does not depend on the variables x 1 , . . ., x k and some invertible matrix A having the first k columns equal to a 1 , . . ., a k .
We then use Theorem 1.

Converting between representations
We will first recall how to convert a function between different representations.For the conversion from univariate to multivariate representation we will work along the lines of [7].
Let g : Fp n → Fp n be represented as a univariate polynomial over Fp n of degree at most p n − 1. Fix a basis α = {α 1 , . . ., α n } for Fp n , viewed as a vector space of dimension n over Fp.The isomorphism ϕ α : Fp n → F n p associates to each element in Fp n its coordinates in the basis α.We then consider n new variables x 1 , . . ., x n representing the coordinates of x in the basis α; note that these new variables are thought of as taking values in Fp, whereas x was taking values in Fp n .Substituting x by α 1 x 1 + . . .+ α n x n and keeping in mind that x p i = x i we convert g(x) into a multivariate polynomial in x 1 , . . ., x n , with coefficients in Fp n .By further applying ϕ α to each coefficient we obtain the multivariate polynomial representation described in Sect.2; we will denote it by ϕ α (g).We can then represent this polynomial in ANF or falling factorial ANF.
To convert the representations in the other direction, let f : F n p → F n p viewed as a polynomial in n variables with coefficients in F n p .First we convert to coefficients to Fp n by applying ϕ −1 α .We then consider a new variable x = α 1 x 1 + . . .+ α n x n .Keeping in mind that x i ∈ Fp means x p i = x i , we have the equations x p i = α . . .
which gives n equations expressing each x i as a linear combination of x, x p , . . ., x p n−1 with coefficients in Fp n .Substituting all x i in f using these equations we obtain a polynomial g in x.One can verify that ϕ α (g) = f ; we denote therefore g = ϕ −1 α ( f ).For computing M −1 α it is convenient to have a basis β which is dual to α; we then have Once we fix a basis of Fp n , differentiation does not depend on the representation, as differentiation only uses the additive group, and (F p n , +) and (F n p , +) are isomorphic.More precisely, we have: Lemma 1 Let α = {α 1 , . . ., α n } be a basis of Fp n as a vector space over Fp.Let f : Fp n → Fp n be a univariate polynomial function and a ∈ Fp n \{0}.Then

Algorithm for antiderivatives of functions over Fp n
We can now give an alternative to the antiderivative algorithms of Suder [4] and of Pasalic et al. [2] for functions over Fp n .Using the changes of representation discussed in Sect.4.1 and Lemma 1 we have: Theorem 4 Let a 1 , . . ., a k ∈ Fp n be linearly independent over Fp and let g 1 , . . ., g k : Fp n → Fp n .There exists a function f : Fp n → Fp n such that D a i f = g i for all i = 1, . . ., k iff the following conditions are satisfied: a i ,...,a i g i = 0 for all i = 1, . . ., k. (ii) D a i g j = D a j g i for all i, j = 1, . . ., k.
If these conditions are satisfied then the set of such functions f can be computed as where the basis α = {α 1 , . . ., α n } is chosen so that α 1 = a 1 , . . ., α k = a k and the function Antiderivative is computed using one of Algorithms 1-3.
Note that although Algorithm 3 is quasilinear, the transformations ϕ α and ϕ −1 α have exponential worst case complexity, so the total complexity of computing the antiderivative for functions given in their univariate representation using Theorem 4 is exponential.
and we would need to apply Theorem 3 in addition to Algorithms 1-3 to compute Antiderivative.
Theorem 6 Let f : F n p → F n p , with p > 2, be a quadratic function in ANF, f = n i=1 b ii x 2 i + 1≤i< j≤n b i j x i x j .We have that f is a PN function if and only if the n × n matrix C = (c i j ) i, j=1,...,n with entries in F n p defined as c ji = c i j = b i j for i < j and c ii = 2b ii has the properties: (i) Each row of C consists of n elements of F n p which are linearly independent over Fp. (ii) Each linear combination (with coefficients in Fp) of the rows in C consists of n vectors in F n p which are linearly independent over Fp (condition (i) is subsumed by (ii) but for convenience we wrote them separately).
Proof (for both Theorems 6 and 5) Here we can use either the APN form or the falling factorial APN form, as x 2 i and x i 2 differ by a polynomial of algebraic degree one, so they are EA-equivalent.For any a = (a 1 , . . ., a n ) ∈ F n p \{0} using Proposition 2(2) and the fact that D e i f is affine (since f is quadratic) we have, for some constant c ∈ F n 2 : Finally, we use the fact that a linear function g : x n is bijective iff c 1 , . . ., c n are linearly independent; for p = 2, g is 2-to-1 if and only if c 1 , . . ., c n span a space of dimension n − 1.
The theorem above can be used for constructing PN functions from scratch, as shown in Example 4 below.It can also be used similarly to the construction in Yu et al. [7], namely we can start with a known quadratic PN function (obtained by some other technique), construct the associated matrix C and then attempt to modify only a few of the entries of C while preserving the property required by Theorem 6.This approach will be the subject of further work.
We also look at examples that use Theorem 5 for constructing APN functions.
Example 5 Let p = 2.For constructing APN functions, due to invariance to EA-equivalence we can assume that the first row (and column) of C consists of vectors in the canonical basis.
For n = 3 we can assume therefore that C = ⎛ ⎝ 0 e 2 e 3 e 2 0 u e 3 u 0 ⎞ ⎠ It is easy to check that for u = e 1 this matrix does satisfy the conditions of Theorem 5, so it corresponds to an APN function, which can be explicitly constructed as f = e 2 x 1 x 2 + e 3 x 1 x 3 + e 1 x 2 x 3 = (x 2 x 3 , x 1 x 2 , x 1 x 3 ).If u is in the vector space generated by e 2 , e 3 then one can easily see that C does not satisfy Theorem 5. Therefore, up to EA-equivalence, the function above is the only APN function in 3 variables.
Note that in Suder [4,Section5], the equivalent of Theorem 5 above reads: "Furthermore, if the linear combinations of these affine derivatives are also 2-to-1, i.e. the columns...span a subspace of dimensions n − 1 over F2 n , then we have a quadratic APN function."While the first part of the statement is correct, the second part can be misleading, as it seems to imply that it is sufficient that each column of C (which, due to symmetry is the same as a each row of C) spans a space of dimension n − 1, in other words, only condition (i) from Theorem 5 need to be satisfied, and not (ii).Condition (i) is necessary, but not sufficient, as the following counterexample demonstrates: If u = e 4 , v = e 1 and w = e 2 then the elements of each row span a space of dimension 3 so condition (i) of Theorem 5 is satisfied.One can also verify that for any sum of two or three rows, the resulting elements also span a space of dimension 3.However, the sum of all four rows consists of the vectors e 2 + e 3 + e 4 , e 1 + e 2 + e 4 , e 2 + e 3 + e 4 , e 1 + e 2 + e 4 which span a space of dimension only 2. So the function f = (x 2 x 4 , x 1 x 2 + x 3 x 4 , x 1 x 3 , x 1 x 4 + x 2 x 3 ) corresponding to this matrix does not satisfy condition (ii) of Theorem 5, therefore it is not APN.Indeed the derivative of f in the direction (1, 1, 1, 1) is not 2-to-1.One choice which does give an APN function is u = e 4 , v = e 3 + e 4 and w = e 1 , corresponding to f = (x 3 x 4 , x 1 x 2 , x 1 x 3 + x 2 x 4 , x 1 x 4 + x 2 x 3 + x 2 x 4 ).

Conclusion
For computing the antiderivative (i.e.retrieving a function from its derivatives) for functions with n bits input and n bits output (or more generally integers modulo a prime p instead of bits), it turns out that representing the function as n multivariate polynomial functions in algebraic normal form over Fp is particularly convenient.For this representation we developed a simple algorithm, which is in general more efficient (quasilinear if the directions of the derivatives are in the canonical basis) compared to previous algorithms which represent the functions as univariate polynomial functions over Fp n and involve solving a system of linear equations of size p n × p n .
The connection pointed out by Suder [4] between antiderivatives and the methods of of Yu et al. [7] and of Weng et al. [5] for constructing new quadratic APN (Almost Perfectly Nonlinear) functions is preserved and simplified by our approach to antiderivatives.Moreover, we develop a similar technique which can be used to construct new quadratic PN (Perfect Nonlinear) functions.

Example 4 ⎠
Let n = 3. Due to invariance under EA-equivalence, we can assume that the first row (and column) of C consists of the canonical basis vectors:For small p we can find all values of u, v and w which satisfy Theorem 6 using an exhaustive computer search.For p = 3 there are 288 solutions; one example solution is u = e 3 , v = e 1 + e 2 , w = e 2 + e 3 .The corresponding quadratic PN function can be constructed explicitly as

Example 6 0 e 2 e 3 e 4 e 2 0 u v e 3 u 0 w e 4
For p = 2 and n = 4 consider the matrix: and each component f i written in its ANF or in the falling factorial ANF, so f i = t c e 1 , . . ., e k )) Input: n, k, g 1 , . . ., g k , e 1 , . . ., e k , with g i : F n p → F m p , deg x i g i ≤ p − 2 and D e i g j = D e j g i for all i, j = 1, . . ., k. Polynomials are represented as the set of monomials that appear in their falling factorial ANF Output: f : F n p → F m p such that D e i f = g i for all i = 1, . . ., k 1: