1206497Sluigi/* 2206552Sluigi * Copyright 2017, Data61 3206552Sluigi * Commonwealth Scientific and Industrial Research Organisation (CSIRO) 4206497Sluigi * ABN 41 687 119 230. 5206497Sluigi * 6206497Sluigi * This software may be distributed and modified according to the terms of 7206497Sluigi * the BSD 2-Clause license. Note that NO WARRANTY is provided. 8206497Sluigi * See "LICENSE_BSD2.txt" for details. 9206497Sluigi * 10206497Sluigi * @TAG(DATA61_BSD) 11206497Sluigi */ 12206497Sluigi 13206497Sluigi#pragma once 14206497Sluigi 15206497Sluigi/* macros for doing basic math */ 16206497Sluigi 17206497Sluigi#include <stdint.h> 18206497Sluigi 19206497Sluigi/* This function can be used to calculate ((a * b) / c) without causing 20206497Sluigi * overflow by doing (a * b). It also allows you to avoid loss of precision 21206497Sluigi * if you attempte to rearrange the algorithm to (a * (c / b)) 22206497Sluigi * It is essentially the Ancient Egyption Multiplication 23206497Sluigi * with (a * (b / c)) but we keep (b / c) in quotient+remainder 24206497Sluigi * form to not lose precision. 25206497Sluigi * This function only works with unsigned 64bit types. It can be 26206497Sluigi * trivially abstracted to other sizes. Dealing with signed is possible 27206497Sluigi * but not included here */ 28206497Sluigistatic inline uint64_t muldivu64(uint64_t a, uint64_t b, uint64_t c) 29206497Sluigi{ 30206497Sluigi /* quotient and remainder of the final calculation */ 31206497Sluigi uint64_t quotient = 0; 32206497Sluigi uint64_t remainder = 0; 33206497Sluigi /* running quotient and remainder of the current power of two 34206497Sluigi * bit we are looking at. As we look at the different bits 35206497Sluigi * of a we will increase these */ 36206497Sluigi uint64_t cur_quotient = b / c; 37206497Sluigi uint64_t cur_remainder = b % c; 38206497Sluigi /* we will iterate through all the bits of a from least to most 39206497Sluigi * significant, we can stop early once there are no set bits though */ 40206497Sluigi while (a) { 41206497Sluigi /* If this bit is set then the power of two is part of the 42206497Sluigi * construction of a */ 43206497Sluigi if (a & 1) { 44206497Sluigi /* increase final quotient, taking care of any remainder */ 45206497Sluigi quotient += cur_quotient; 46206497Sluigi remainder += cur_remainder; 47206497Sluigi if (remainder >= c) { 48206497Sluigi quotient++; 49206497Sluigi remainder -= c; 50206497Sluigi } 51206497Sluigi } 52206497Sluigi /* Go to the next bit of a. Also means the effective quotient 53206497Sluigi * and remainder of our (b / c) needs increasing */ 54206497Sluigi a >>= 1; 55206497Sluigi cur_quotient <<= 1; 56206497Sluigi cur_remainder <<= 1; 57206497Sluigi /* Keep the remainder sensible, otherewise the remainder check 58206497Sluigi * in the if (a & 1) case has to become a loop */ 59206497Sluigi if (cur_remainder >= c) { 60206497Sluigi cur_quotient++; 61206497Sluigi cur_remainder -= c; 62206497Sluigi } 63206497Sluigi } 64206497Sluigi return quotient; 65206497Sluigi} 66206497Sluigi