1/*
2 * Copyright 2017, Data61
3 * Commonwealth Scientific and Industrial Research Organisation (CSIRO)
4 * ABN 41 687 119 230.
5 *
6 * This software may be distributed and modified according to the terms of
7 * the BSD 2-Clause license. Note that NO WARRANTY is provided.
8 * See "LICENSE_BSD2.txt" for details.
9 *
10 * @TAG(DATA61_BSD)
11 */
12
13#pragma once
14
15/* macros for doing basic math */
16
17#include <stdint.h>
18
19/* This function can be used to calculate ((a * b) / c) without causing
20 * overflow by doing (a * b). It also allows you to avoid loss of precision
21 * if you attempte to rearrange the algorithm to (a * (c / b))
22 * It is essentially the Ancient Egyption Multiplication
23 * with (a * (b / c)) but we keep (b / c) in quotient+remainder
24 * form to not lose precision.
25 * This function only works with unsigned 64bit types. It can be
26 * trivially abstracted to other sizes. Dealing with signed is possible
27 * but not included here */
28static inline uint64_t muldivu64(uint64_t a, uint64_t b, uint64_t c)
29{
30    /* quotient and remainder of the final calculation */
31    uint64_t quotient = 0;
32    uint64_t remainder = 0;
33    /* running quotient and remainder of the current power of two
34     * bit we are looking at. As we look at the different bits
35     * of a we will increase these */
36    uint64_t cur_quotient = b / c;
37    uint64_t cur_remainder = b % c;
38    /* we will iterate through all the bits of a from least to most
39     * significant, we can stop early once there are no set bits though */
40    while (a) {
41        /* If this bit is set then the power of two is part of the
42         * construction of a */
43        if (a & 1) {
44            /* increase final quotient, taking care of any remainder */
45            quotient += cur_quotient;
46            remainder += cur_remainder;
47            if (remainder >= c) {
48                quotient++;
49                remainder -= c;
50            }
51        }
52        /* Go to the next bit of a. Also means the effective quotient
53         * and remainder of our (b / c) needs increasing */
54        a >>= 1;
55        cur_quotient <<= 1;
56        cur_remainder <<= 1;
57        /* Keep the remainder sensible, otherewise the remainder check
58         * in the if (a & 1) case has to become a loop */
59        if (cur_remainder >= c) {
60            cur_quotient++;
61            cur_remainder -= c;
62        }
63    }
64    return quotient;
65}
66