plugins/ioam-plugin/ioam/lib-pot/math64.h

   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 static inline void mul64by64(u64 a, u64 b, u64 * hi, u64 * lo)
  32 {
  33     u64 a_lo = (u64) (uint32_t) a;
  34     u64 a_hi = a >> 32;
  35     u64 b_lo = (u64) (u32) b;
  36     u64 b_hi = b >> 32;
  37
  38     u64 p0 = a_lo * b_lo;
  39     u64 p1 = a_lo * b_hi;
  40     u64 p2 = a_hi * b_lo;
  41     u64 p3 = a_hi * b_hi;
  42
  43     u32 cy = (u32) (((p0 >> 32) + (u32) p1 + (u32) p2) >> 32);
  44
  45     *lo = p0 + (p1 << 32) + (p2 << 32);
  46     *hi = p3 + (p1 >> 32) + (p2 >> 32) + cy;
  47     return;
  48 }
  49
  50 #define TWO64 18446744073709551616.0
  51
  52 static inline u64 mod128by64(u64 x, u64 y, u64 m, double di)
  53 {
  54     u64 q1, q2, q;
  55     u64 p1, p0;
  56     double dq;
  57
  58     /* calculate quotient first pass 53 bits */
  59     dq = (TWO64 * (double)x + (double)y) * di;
  60
  61     if (dq >= TWO64)
  62         q1 = 0xfffffffffffff800L;
  63     else
  64         q1 = dq;
  65
  66     /* q1 * m to compare the product to the dividend.  */
  67     mul64by64(q1, m, &p1, &p0);
  68
  69     /* Adjust quotient. is it > actual result: */
  70     if (x < p1 || (x == p1 && y < p0))
  71     {
  72         /* q1 > quotient.  calculate abs remainder */
  73         x = p1 - (x + (p0 < y));
  74         y = p0 - y;
  75
  76         /* use the remainder as new dividend to adjust quotient */
  77         q2 = (u64) ((TWO64 * (double)x + (double)y) * di);
  78         mul64by64(q2, m, &p1, &p0);
  79
  80         q = q1 - q2;
  81         if (x < p1 || (x == p1 && y <= p0))
  82         {
  83             y = p0 - y;
  84         }
  85         else
  86         {
  87             y = p0 - y;
  88             y += m;
  89             q--;
  90         }
  91     }
  92     else
  93     {
  94         x = x - (p1 + (y < p0));
  95         y = y - p0;
  96
  97         q2 = (u64) ((TWO64 * (double)x + (double)y) * di);
  98         mul64by64(q2, m, &p1, &p0);
  99
 100         q = q1 + q2;
 101         if (x < p1 || (x == p1 && y < p0))
 102         {
 103             y = y - p0;
 104             y += m;
 105             q--;
 106         }
 107         else
 108         {
 109             y = y - p0;
 110             if (y >= m)
 111             {
 112                 y -= m;
 113                 q++;
 114             }
 115         }
 116     }
 117
 118     return y;
 119 }
 120
 121 /*
 122  * returns a % p
 123  */
 124 static inline u64 mod64by64(u64 a, u64 p, u64 primeinv)
 125 {
 126     return (mod128by64(0, a, p, primeinv));
 127 }
 128
 129 static inline void add64(u64 a, u64 b, u64 * whi, u64 * wlo)
 130 {
 131     *wlo = a + b;
 132     if (*wlo < a)
 133         *whi = 1;
 134
 135 }
 136
 137 /*
 138  * returns (a + b)%p
 139  */
 140 static inline u64 add64_mod(u64 a, u64 b, u64 p, double pi)
 141 {
 142     u64 shi = 0, slo = 0;
 143
 144     add64(a, b, &shi, &slo);
 145     return (mod128by64(shi, slo, p, pi));
 146 }
 147
 148 /*
 149  * returns (ab) % p
 150  */
 151 static inline u64 mul64_mod(u64 a, u64 b, u64 p, double pi)
 152 {
 153     u64 phi = 0, plo = 0;
 154
 155     mul64by64(a, b, &phi, &plo);
 156     return (mod128by64(phi, plo, p, pi));
 157 }
 158
 159 #endif