Line data    Source code

       1             : /* 
       2             :  * math64.h provides the 64 bit unsigned integer add, multiply followed by  modulo operation
       3             :  * The linux/math64.h provides divide and multiply 64 bit integers but:
       4             :  * 1. multiply: mul_u64_u64_shr - only returns 64 bits of the result and has to be called
       5             :  *                     twice to get the complete 128 bits of the result.
       6             :  * 2. Modulo operation of the result of  addition and multiplication of u64 that may result 
       7             :  *                        in integers > 64 bits is not supported
       8             :  * Hence this header to combine add/multiply followed by modulo of u64 integrers
       9             :  * always resulting in u64.
      10             :  * 
      11             :  * Copyright (c) 2016 Cisco and/or its affiliates.
      12             :  * Licensed under the Apache License, Version 2.0 (the "License");
      13             :  * you may not use this file except in compliance with the License.
      14             :  * You may obtain a copy of the License at:
      15             :  *
      16             :  *     http://www.apache.org/licenses/LICENSE-2.0
      17             :  *
      18             :  * Unless required by applicable law or agreed to in writing, software
      19             :  * distributed under the License is distributed on an "AS IS" BASIS,
      20             :  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
      21             :  * See the License for the specific language governing permissions and
      22             :  * limitations under the License.
      23             :  */
      24             : #ifndef include_vnet_math64_h
      25             : #define include_vnet_math64_h
      26             : #include <stdint.h>
      27             : 
      28             : /* 
      29             :  * multiplies and returns result in hi and lo 
      30             :  */
      31           0 : static inline void mul64by64(u64 a, u64 b, u64 * hi, u64 * lo)
      32             : {
      33           0 :     u64 a_lo = (u64) (uint32_t) a;
      34           0 :     u64 a_hi = a >> 32;
      35           0 :     u64 b_lo = (u64) (u32) b;
      36           0 :     u64 b_hi = b >> 32;
      37             : 
      38           0 :     u64 p0 = a_lo * b_lo;
      39           0 :     u64 p1 = a_lo * b_hi;
      40           0 :     u64 p2 = a_hi * b_lo;
      41           0 :     u64 p3 = a_hi * b_hi;
      42             : 
      43           0 :     u32 cy = (u32) (((p0 >> 32) + (u32) p1 + (u32) p2) >> 32);
      44             : 
      45           0 :     *lo = p0 + (p1 << 32) + (p2 << 32);
      46           0 :     *hi = p3 + (p1 >> 32) + (p2 >> 32) + cy;
      47           0 :     return;
      48             : }
      49             : 
      50             : #define TWO64 18446744073709551616.0
      51             : 
      52           0 : static inline u64 mod128by64(u64 x, u64 y, u64 m, double di)
      53             : {
      54             :   u64 q1, q2;
      55             :   u64 p1, p0;
      56             :   double dq;
      57             : 
      58             :   /* calculate quotient first pass 53 bits */
      59           0 :   dq = (TWO64 * (double) x + (double) y) * di;
      60             : 
      61           0 :   if (dq >= TWO64)
      62           0 :     q1 = 0xfffffffffffff800L;
      63             :   else
      64           0 :     q1 = dq;
      65             : 
      66             :   /* q1 * m to compare the product to the dividend.  */
      67           0 :   mul64by64 (q1, m, &p1, &p0);
      68             : 
      69             :   /* Adjust quotient. is it > actual result: */
      70           0 :   if (x < p1 || (x == p1 && y < p0))
      71             :     {
      72             :         /* q1 > quotient.  calculate abs remainder */
      73           0 :         x = p1 - (x + (p0 < y));
      74           0 :         y = p0 - y;
      75             : 
      76             :         /* use the remainder as new dividend to adjust quotient */
      77           0 :         q2 = (u64) ((TWO64 * (double)x + (double)y) * di);
      78           0 :         mul64by64(q2, m, &p1, &p0);
      79             : 
      80           0 :         if (x < p1 || (x == p1 && y <= p0))
      81             :         {
      82           0 :             y = p0 - y;
      83             :         }
      84             :         else
      85             :         {
      86           0 :             y = p0 - y;
      87           0 :             y += m;
      88             :         }
      89             :     }
      90             :     else
      91             :     {
      92           0 :         x = x - (p1 + (y < p0));
      93           0 :         y = y - p0;
      94             : 
      95           0 :         q2 = (u64) ((TWO64 * (double)x + (double)y) * di);
      96           0 :         mul64by64(q2, m, &p1, &p0);
      97             : 
      98           0 :         if (x < p1 || (x == p1 && y < p0))
      99             :         {
     100           0 :             y = y - p0;
     101           0 :             y += m;
     102             :         }
     103             :         else
     104             :         {
     105           0 :             y = y - p0;
     106           0 :             if (y >= m)
     107             :             {
     108           0 :                 y -= m;
     109             :             }
     110             :         }
     111             :     }
     112             : 
     113           0 :     return y;
     114             : }
     115             : 
     116             : /* 
     117             :  * returns a % p
     118             :  */
     119             : static inline u64 mod64by64(u64 a, u64 p, u64 primeinv)
     120             : {
     121             :     return (mod128by64(0, a, p, primeinv));
     122             : }
     123             : 
     124           0 : static inline void add64(u64 a, u64 b, u64 * whi, u64 * wlo)
     125             : {
     126           0 :     *wlo = a + b;
     127           0 :     if (*wlo < a)
     128           0 :         *whi = 1;
     129             : 
     130           0 : }
     131             : 
     132             : /* 
     133             :  * returns (a + b)%p
     134             :  */
     135           0 : static inline u64 add64_mod(u64 a, u64 b, u64 p, double pi)
     136             : {
     137           0 :     u64 shi = 0, slo = 0;
     138             : 
     139           0 :     add64(a, b, &shi, &slo);
     140           0 :     return (mod128by64(shi, slo, p, pi));
     141             : }
     142             : 
     143             : /* 
     144             :  * returns (ab) % p
     145             :  */
     146           0 : static inline u64 mul64_mod(u64 a, u64 b, u64 p, double pi)
     147             : {
     148           0 :     u64 phi = 0, plo = 0;
     149             : 
     150           0 :     mul64by64(a, b, &phi, &plo);
     151           0 :     return (mod128by64(phi, plo, p, pi));
     152             : }
     153             : 
     154             : #endif