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 */
31static 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
52static 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 */
124static inline u64 mod64by64(u64 a, u64 p, u64 primeinv)
125{
126    return (mod128by64(0, a, p, primeinv));
127}
128
129static 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 */
140static 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 */
151static 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
160