Showing posts with label arithmetic. Show all posts
Showing posts with label arithmetic. Show all posts

Introduction to polynomial approximation (WiP)

The goal is to find the best polynomial to approximate a given function $f$ on a given interval $I$.
The best generally meaning the polynomial which minimizes the error (absolute or relative) with the input function. This error is often $L-\infty$ (the maximum of the absolute difference between the function and the approximating polynomial).

Remez:

   The Remez algorithm find solution to linear system of $n+2$ equations $\sum_i{a_i.x_j^i} + (-\epsilon)^i = f(x_j) $ for j in $[0, n+1]$. It starts by setting $x_j$ values to the extrema of the degree $n$ Chebyshev polynomial (mapped to the input interval through affine scaling). Then solve the linear systems which provides a first set of $a_i$ coefficients. Then the extrema of $| P - f $ are found and used as a new set of $x_j$ and the processus is iterated. Theoretically, the process stops when the approximation errors are all equal and of alternating sign, the returned P is the minmax approximation polynomial.
During the iteration the reachable lower bound for the approximation error is the minimum of the absolute value of the extrema the difference $P - f$ (therorem de la Vallée Poussin)

Using Euclidean Lattices

Educational example of polynomial approximation: https://github.com/nibrunie/YourOwnPersonnalRemez

References

- Efficient L-inf polynomial approximations (https://hal.archives-ouvertes.fr/inria-00119513v1)

Determining numerical bound for floating-point reciprocal overflow


This article is both an example of pythonsollya use and a small study of floating-point property of the reciprocal function.

Context and objectives

Our goal is to find a bound $b$ such that when computing $\frac{1}{x}$ for any $x \lt b $  an IEEE overflow exception is raised (and the value corresponding to an overflow in the given rounding mode is returned). $b^*$ designs the exact value and $b$ the floating-point number of precision $p$ which verify the above property. Let us note that $b$ is not necessarily equal to $b^*$ rounded in a pre-defined rounding mode.

We will use the following notations: $\circ$ is the rounding operation which can be specialized by a rounding mode, $RNE$ is rounding to nearest (tie-to-even), $RU$ is rounding up (to $+\infty$) and $RD$ is rounding downward (to $-\infty)$, $next(x)$ is the minimal floating-point number strictly greater than $x$ and $prev(x)$ is the minimal floating-point number strictly lower than $x$.

Let us first define an overflow bound:

Defintion 1: $b$, (respectively $b^*$) is defined as the floating-point (resp. exact) overflow bound of the reciprocal function if it is the minimal floating-point (resp. real) number which verifies: $ \circ(\frac{1}{b})$ does not overflow.

Section 7.4.0 of the IEEE-754 standard defines overflow as: "The overflow exception shall be signaled if and only if the destination format’s largest finite number is exceeded in magnitude by what would have been the rounded floating-point result (...) were the exponent range unbounded. " This means that $b$ shoud depend on the rounding mode.

Let us first consider $b_{RNE}$ the bound associated to the round to nearest (tie-to-even) mode and defined as the maximal $x$ such that: $ \frac{1}{x} \ge B_{RNE} = \omega + \frac{1}{2} ulp(\omega) $.
$B_{RNE}$ is the mid-point between the largest floating-point number $\omega$ and the first value defined with an unbounded exponent range which will be rounded to $+\infty$ as per Section 7.4.
Note that the inequality does not need to be strict since $\omega$ has an odd mantissa ($1 ... 1$) if $\frac{1}{x} = B_{RNE}$ it should be rounded to a value greater than $\omega$ because of the tie-to-even rule, thus it will overflow.

$p$ is the format precision (e.g. $p=24$ for binary32 format).
$ b^*_{RNE} = \frac{1}{B_{RNE}} = \frac{1}{\omega + 2^{emax - p}} $
$ b^*_{RU} = \frac{1}{B_{RU}} = \frac{1}{\omega} $

Evaluating bound values

Let us compute (approximation to) $b^*_{RNE}$ and $b^*_{RU}$ using pythonsollya:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# importing pythonsollya
>>> import sollya
# forcing display setting to hexadecimal value
# (easier to read)
>>> sollya.settings.display = hexadecimal
# building 2.0 as a Sollya Object
>>> S2 = sollya.SollyaObject(2)
# building binary32's omega value (greatest floating-point number)
>>> omega = S2**127 * (S2**24 - 1) * S2**-23
>>> omega
0x1.fffffep127
>>> B_RU = omega
# omega + 1/2 ulp(omega)
>>> B_RNE = omega  + S2**103
# rounding 1 / B_RU with 100-bit precision
>>> sollya.round(1 / B_RU, 100, sollya.RD)
0x1.000001000001000001000001p-128
# rounding B_RNE to 100-bit precision
>>> sollya.round(1 / B_RNE, 100, sollya.RD)
0x1.0000008000004000002p-128
# rounding 1 / B_RNE to binary32 precision
>>> sollya.round(1 / B_RNE, sollya.binary32, sollya.RN)
0x1p-128
# rounding 1 / B_RU to binary32 precision
>>> sollya.round(1 / B_RU, sollya.binary32, sollya.RN)
0x1p-128

The first interestings values are given lines 16 and 19: correctly rounded approximation on 100 bits of $b_{RU} = \frac{1}{B_{RU}}$ and $b_{RNE} = \frac{1}{B_{RNE}}$. Let us first notice that those two values are greater than they counterpart rounded to binary32 (visible on lines 22 and 25).
Let us call $c_- = 2^{-128}$ and $c_+ = next(c_-) = 2^{-128} + 2^{-149} $ , we have $c_- \lt b_{RU} \lt c_+ $ and $c_- \lt b_{RNE}\lt c_+ $ which means $ \frac{1}{c_-} \gt \frac{1}{b_{RNE}} \gt \frac{1}{c_+} $ and $ \frac{1}{c_-} \gt \frac{1}{b_{RU}} \gt \frac{1}{c_+} $ and $\forall x , x \le c, \frac{1}{x} \gt \frac{1}{b_{RNE}}$

Claim 1: $ \forall mode, \circ_{mode}(\frac{1}{c_+}) $ is a normal value.
Claim 2: $ \forall mode, \circ_{mode}(\frac{1}{c_-}) $ overflows.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
>>> c_minus = S2**-128
>>> c_plus = S2**-128 + S2**-149
>>> c_minus
0x1p-128
>>> c_plus
0x1.000008p-128
>>> for rnd_mode in [sollya.RN, sollya.RU, sollya.RD, sollya.RZ]:
...     print sollya.round(1 / c_minus, sollya.binary32, rnd_mode), sollya.round(1 / c_plus, sollya.binary32, rnd_mode)
infty 0x1.fffffp127
infty 0x1.fffff2p127
0x1.fffffep127 0x1.fffffp127
0x1.fffffep127 0x1.fffffp127

$c_+ $ is greater than both $b^*_{RU}$ and $b^*_{RNE}$, and $c_-$ is lower (strictly) than both of them. Thus, as there is no floating point number between $c_-$ and $c_+$ we can chose $c_-$ as the overflow bound (included) or $c_+$(excluded). There is no need to distinguish several values of $c$ based on the rounding mode. $\forall \ x \ge c_+$, $\frac{1}{x}$ will not overflow and $\forall \ x \le c_-$, $\frac{1}{x}$ will overflow in any rounding mode.

Let us now consider the remaining rounding modes (rounding down and towards zero). Without loss of generality we can reduce the study to positive overflow and thus to a single rounding mode, let us say rounding down.
$B_{RD} = next(\omega) = 2^{128} $ (assuming unbounded exponent range)


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
>>> B_RD = S2**128
>>> b_RD = S2**-128
>>> sollya.round(1 / B_RD, sollya.binary32, sollya.RD)
0x1p-128
>>> sollya.round(1 / B_RD, sollya.binary32, sollya.RU)
0x1p-128
>>> sollya.round(1 / c_minus, sollya.binary32, sollya.RD)
0x1.fffffep127
>>> sollya.round(1 / c_plus, sollya.binary32, sollya.RD)
0x1.fffffp127

$\frac{1}{B_{RD}}$ is exactly equal to the floating point number $c_-$ and $\frac{1}{c_-}$ will indeed lead to an overflow even in rounded down mode. This overflow will return the $\omega$ value as output.
$\circ_{RD}(\frac{1}{c_+}) = \omega - 2^{104} - 2^{105} - 2^{106} $ which is a normal number strictly lower than omega. Thus we get the same behavior as for $RNE$ and $RU$: $c_-$ is the overflow bound (included) and $c_+$ is the excluded overflow bound $\square$.

References:

  • Sollya: a tool for manipulating arbitrary precision numbers, computing approximation polynomial ... (Sollya's website)
  • pythonsollya: a python wrapper for sollya (pythonsollya's gitlab)
  • Quick introduction to IEEE-754 floating-point arithmetic standard (web)


New instruction for Fast GCM multiplication


New Operation/Instruction for Fast GCM multiplication


In this article we will study how we could compute GCM multiplication with only 3 instructions (and two exclusive 64-bit or) by extending a very interesting white paper by Intel:  Intel ® Carry-Less Multiplication Instruction and its Usage for computing the GCM Mode ([4]).
We will first introduce GCM, then survey the basic of the carry less multiplication and the Karatsuba algorithm for multiplication before detailing the modulo reduction method of [4]. Then we will introduce our new operation and evaluate its hardware implementation.

Galois Counter Mode

Galois-Counter Mode or GCM ([1]) is a block cipher mode of operation which uses a multiplication in a finite field (Galois Field or $GF$) as basis for the authentication procedure. The Galois Field used is $ GF(2)[X] \!/\!_{X^{128} + X^7 + X ^2 + 1}$. The fact that the irreducible polynomial used to build the field has a very low number of non-zero terms is very useful to implement a fast modulo reduction.

Carry Less Multiplication

In the following most binary values encode polynomial over $GF(2)[X]$. For example a 64-bit value $A$ encodes $A=\sum_{i=0}^{63}a_i . X^i$ where a_i is the i-th bit of A (and a coefficient in $GF(2)$). The addition of two polynomials $A$ and $B$ is simply their binary xor: $A + B = \sum_{i=0}^{63}(a_i + b_i \mod 2) = \sum_{i=0}^{63}(a_i \oplus b_i ) = A \oplus B $.
A multiplication of polynomials over $GF(2)[X] is called a Carry Less Multiplication because there are no carry propagation between monomial term of distinct degrees.
The following listing shows a possible implementation of a Carry-Less multiplication of two 64-bit polynomials, giving a 128-bit result. As we operate on polynomials, the accumulate operation is a XOR which contrary to an addition does not propagate a carry.

__int128 poly_mul_64x64_128(uint64_t a, uint64_t b)
{
   int i;
   __int128 result = 0;
   for (i = 0; i < 64; ++i) {
       // result +=  b . X^i if a_i is non-zero
       result ^= ((a >> i) & 1) ? (b << i) : 0;
   }

}

Recently, in x86 architecture,  the carry-less multiplication operation is implemented by a specific instruction: PCLMULQDQ ([2]).

Assuming the availability of a $64-bit \times 64-bit \rightarrow 128-bit $ carry-less operation, a full GCM reduction can be performed by first computing the full $128 \times 128$ CLM using 4 $64 \times 64$ multiply operations and then performing a $mod G$ reduction on the $128$ most significant bits of the multiplication result.

Fast GCM Multiplication

A very interesting white paper by intel [4] sums up a lot of good ideas on how to implement a fast GCM multiplication.
   The proposal is two-fold:
  • Use Karatsuba algorithm [5] to accelerate the multiplication
  • Use the structure of the field to accelerate the reduction part 
This white paper contains other interesting idea (folding ...) that will not be detailed in this article.

Karatsuba multiplication

We wish to compute $M = A \times B $, where $A$ and $B$ are 128-bit polynomial (degree 127 in $GF(2)[X]$). Thus M is at of degree 255. Let us split A and B into high and low part:
$$ A = A_{hi}. X^{64} + A_{lo}, B = B_{hi}. X^{64} + B_{lo}$$
$$ M = A_{hi} \times B_{hi} . X^{128} +  (A_{hi} \times B_{lo} + A_{lo} \times B_{hi}) . X^{64} + A_{lo} \times B_{lo} $$
Assuming $ C =  A_{hi} \times B_{hi}$, and $ D = A_{lo} \times B_{lo}$ and $E = (A_{hi} + A_{lo}) $ and $F = B_ {hi} + B_{lo}$,   $ (A_{hi} \times B_{lo} + A_{lo} \times B_{hi}) $ can be writen as $$  (A_{hi} \times B_{lo} + A_{lo} \times B_{hi}) = E \times F - C - D $$.
As we are working in $GF(2)[X]$, the operation $+$ is equivalent to $-$: both are binary $xor$.

We went from using 4 half multiplications in the standard algorithm to only 3 with Karatsuba (while appending a few additions). Additions are considered cheap compared to multiplications (and they generally are) so it should not be an issue.

Fast Reduction

[4] provides the following  algorithm (Algorithm 3) to compute, $R=M_{hi} \times X^t (\mod G)$, assuming an input $M_{hi}$ and $G$ an irreducible polynomial of degree $t$, $s=max(degree( M_{hi}))+1$ (for GCM, $s=128$):

  1. Preprocessing: we compute $g$ and $q$, $g$ is the polynomial built by the $t-1$ least coefficient of $G$ and $q$ is the quotient of the division of $X^{t+s}$ by $G$ where $s$ is the degree of the quotient of the division of $X^t \times M$ by $G$
  2. Compute $A = M_{hi} \times q$, 
  3. Multiply the $s$ most significant terms of $A$ by $g$
  4. Output the $t$ least significant terms of the Step 3.
Let us apply this algorithm to the case of GCM: $G = X^{128} + X^7 + X^2 + X + 1$, $g=X^7 + X^2 + X + 1$. We want to compute $M = A \times B $, where $A$ and $B$ are polynomials of degree 127. 
Let us first compute $q$ such that $X^{256}=q \times G + r $ with $degree(r) < degree(G)$, it is easy to verify that $X^{256}=G^2 \oplus x^{14}+x^4+x^2+1 $, thus $q = G$. 

We split M in 4 part of 64 term each: from most significant to least signifiant $M_3, M_2, M_1, M_0$, $M_{hi} =[ M_3 : M_2] = M_3 . X^{64} \oplus M_2 $.

$$ M [G] = M_1 . X^{64} \oplus M_0 \oplus ((M_3 X^{64} \oplus M_2) . X^{128} [G]) $$

To get $((M_3 X^{64} \oplus M_2) . X^{128} [G])$, we use the algorithm from [4]:
$$ A =  (M_3 X^{64} \oplus M_2) \times q $$
As we only need the 128 most significan terms of $A$, $A_{hi}$ we can compute them as follows:

\begin{equation}
\begin{split}
A_{hi} & = [M_3 : M_2] \oplus [M_3 : M_2 ] >> (128 - 7) \oplus [M_3 : M_2 ] >> (128 - 2) \\
             &\oplus [M_3 : M_2 ] >> (128 - 1) \oplus  [M_3 : M_2 ] >> (128 - 0)   \\
             & = [M_3 : M_2] \oplus [0 : M_3 ] >> 57 \oplus [0 : M_3 ] >> 62 \oplus [0 : M_3] >> 63 \\
      & = [M_3 : (M_2 \oplus M_3 >> 57 \oplus M_3 >> 62 \oplus M_3 >> 63)]
\end{split}
\end{equation}

We now compute $R_{lo}$ from $ R = A_{hi} \times g $
\begin{equation}
\begin{split}
R &= A_{hi} \times g \\
 &=  A_{hi} \times (X^7 + X ^2 + X + 1) \\
&= [A_{hi,1} : A_{hi,0}] << 7 \oplus [A_{hi,1} : A_{hi,0}] << 2 \oplus [A_{hi,1} : A_{hi,0}] << 1 \oplus [A_{hi,1} : A_{hi,0}]
\end{split}
\end{equation}

Computing $A_{hi}$ requires 3 64-bit right shifts and 3 64-bit exclusive-or operations.
Computing $R$ from $A_{hi}$ is a little more expensive: it requires 3 128-bit shifts and 3 128-bit exclusive or operations. Overall a hardware implementation of this reduction will require around 576 standard xor gates (static shifts being free).

Merging operation for GCM multiplication

We are now going to merge the operations described above and the Karatsuba algorithm. The idea is to suggest a "generic" operation that could be implemented as processor instruction to accelerate GCM multiplication with more efficiency than the expanded version described in [4].
  1. Compute $O_1 = (A_{lo} \times B_{lo} \oplus (A_{lo} \times B_{lo} . X^{64})) [G] $
  2. Compute $O_2 = (A_{hi} \times B_{hi} . X^{128} \oplus (A_{hi} \times B_{hi} . X^{64})) [G] \oplus O_1 $
  3. Compute $(O_3 = ((A_{hi} \oplus A_{lo}) \times (B_{hi} \oplus B_{lo}) . X^{64}) [G] \oplus O_2 $
  4. Output the final result $O_3 $ 
Each step consists in:
  • One $64-bit \times 64-bit$ carry-less multiplication
  • One 128-bit modulo $G$ reduction
  • A few exclusive-or operations
The product operands for the third step ($O_3$) may be computed in parallel to step 1 and 2, they consist in two 64-bit exclusive or operations.
We can factorize the previous steps into a new primitive: GCMRED.

$GCMRED(A, B, \alpha, \beta, C)$, with $A$ and $B$ 64-bit operand, $C$ a 128-bit operand and $\alpha, \beta$ coefficients in $\{0, 1, X^{64}, X^{128}\}$  implements the following operations:
$$ GCMRED(A, B, \alpha, \beta, C) = (A \times B) . (\alpha + \beta)  [G] \oplus C $$
We can then compact a full GCM multiplication into 3 consecutives application of $GCMRED$:

\begin{equation*}
\begin{split}
 GCM(A_{full}, B_{full}) &= GCMRED(A_{hi} \oplus A_{lo}, B_{hi} \oplus B_{lo}, X^{64}, 0, \\
                                               & \ \ \  \ \ \  GCMRED(A_{hi} , B_{hi} , X^{64}, X^{128}, \\
                                               &  \ \ \ \ \ \ \ \ \      GCMRED(A_{lo} , B_{lo} , 1, X^{64}, 0 )))
\end{split}
\end{equation*}

As you can see we do not need the full range of $\alpha$ and $\beta$: only 3 variants are required, namely: $(X^{64}, 0)$, $(X^{64}, X^{128})$ and $(1, X^{64})$. Thus only those 3 may be encoded as possible instructions.

Complexity of GCMRED

A possible implementation of GCMRED is illustrated by the figure below. There are 4 64-bit inputs: $A$, $B$ ,$C_{hi}$ and $C_{lo}$ and 2 64-bit output corresponding to the high and low part of the output. 

The desisgn consists in :
  • A $64 \times 64$ Carry-Less Multiplication
  • A 128-bit GCM modulo reducer
  • 6 64-bit exclusive or operators
  • 4 64-bit and operators (to zeroify signals)
  • 2 64-bit 3-input multiplexers 
A naive CLM implementations would require around 4096 $and$ gates (with two 1-bit inputs and one 1-bit output) and about the same number of $xor$ gates (give or take a few). The GCM modulo reducer as described in a previous section requires about 576 $xor$ gates. 
   If we consider the silicon area of a xor gates is about twice those of a nand gate (generally used as reference for silicon design comparison),  the silicon area of an $and$ gate is $ \frac{3}{2}$ times a nand gate and the silicon area of a 3-input multiplexer (with one input to zero) to be 5 times a nan gate we can evaluate the overcost of GCMRED compare to a simple CLM multiplication to be roughly around $25\%$.
   This delta may be under estimated as the area to implement a carry less multiply can be further reduced from the naive implementation by using the Karatsuba algorithm recursively to obtain a more efficient design.

Endianess and GCM

In the AES version of Galois Counter Mode, input are assumed to be received in reverse order (least significant bit correspond to the term $X^{127}$ in a 128-bit value) which implies both inputs and outputs must be permuted before and after standard polynomial computations to ensure correct results. Such a permutation is not performed by GCMRED. The cost of supporting it may be limited to introducing a level of inputs multiplexers and a level of output multiplexers to select between standard I/Os and permuted ones (allowing only two permutations).

Conclusion

This small technical report suggests a small addition to the PCLMULQDQ design with which a 25% area increase will bring down the cost of computing a GCM multiplication to five operations (with a critical path length of 3). 

Future works include: synthesising the design and a standard CLM using a standard silicon technology node for better comparison (with critical path and latency analysis), comparing more accurately a full GCM implementation using the new GCMRED instruction operation with Intel's best implementation (including folding and factorized bit reversing).

Updated on January 3rd, 2018
Thank you to Julien LM and Hugues de LSG and Arnaud O for pointing many mistakes in the first versions.
Second update on August 27th, 2018, fixing typo in GCMRED description

References:

Short note: consequence of arithmetic overflow

This is not a full length article but simply a short note.

I recently came across a very interesting post on OpenSSL blog which detailed how a "bad" use of pointer arithmetic lead to a possible low-level security breach in the library. I will let the interested reader go through the original post to understand all the technical details and I will focus on a short explanation. In C90 pointer arithmetic is only defined/specified for malloc-ed pointers with operands less than the malloc-ed size, everything else is undefined behavior.
  It means that in the following example, len must be less or equal to SIZE for the test to work as expected, which makes it pretty much completely useless since this code is supposed to test that specific condition.


/* implied */
char* p = malloc(SIZE);

// buffer overflow detection
if (p + len < limit) {
  // manage overflow case
}

    The reason is as follows: if len is superior to SIZE, then p + len may overflow from the memory bound and wrap, thus resulting in a lesser value which will then be wrongly compared as lesser than limit. In the case of OpenSSL it could leads to decryption or encryption outside the assumed area (at the least).
    This is a good example showing that arithmetic overflow (and pointer arithmetic) must be considered carefully, especially when security is a concern,

References

  1. https://www.openssl.org/blog/blog/2016/06/27/undefined-pointer-arithmetic/