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