OXYGENHUB
◆ main openmsx / src / utils / DivModBySame.cc
1.7 KiB · C++ · last changed chore(): initial commit 1 month ago source blame history raw
1
#include "DivModBySame.hh"
2
3
#include "uint128.hh"
4
5
#include <bit>
6
7
namespace openmsx {
8
9
void DivModBySame::setDivisor(uint32_t divisor_)
10
{
11
	//assert(divisor_ < 0x8000000000000000ULL); // when divisor is uint64_t
12
	divisor = divisor_;
13
	if (divisor == 0) {
14
		m = a = s = 0;
15
		return;
16
	}
17
18
	// reduce divisor until it becomes odd
19
	uint32_t n = 0;
20
	uint64_t t = divisor;
21
	while (!(t & 1)) {
22
		t >>= 1;
23
		++n;
24
	}
25
	if (t == 1) {
26
		m = 0xffffffffffffffffULL;
27
		a = m;
28
		s = 0;
29
	} else {
30
		// Generate m, s for algorithm 0. Based on: Granlund, T.; Montgomery,
31
		// P.L.: "Division by Invariant Integers using Multiplication".
32
		// SIGPLAN Notices, Vol. 29, June 1994, page 61.
33
		auto l = narrow<uint32_t>(std::bit_width(t));
34
		uint64_t j = 0xffffffffffffffffULL % t;
35
		uint128 k = (uint128(1) << (64 + l)) / (0xffffffffffffffffULL - j);
36
		uint128 m_low  =  (uint128(1) << (64 + l))      / t;
37
		uint128 m_high = ((uint128(1) << (64 + l)) + k) / t;
38
		while (((m_low >> 1) < (m_high >> 1)) && (l > 0)) {
39
			m_low >>= 1;
40
			m_high >>= 1;
41
			--l;
42
		}
43
		if ((m_high >> 64) == 0) {
44
			m = low64(m_high);
45
			s = l;
46
			a = 0;
47
		} else {
48
			// Generate m, s for algorithm 1. Based on: Magenheimer, D.J.; et al:
49
			// "Integer Multiplication and Division on the HP Precision Architecture".
50
			// IEEE Transactions on Computers, Vol 37, No. 8, August 1988, page 980.
51
			s = narrow<uint32_t>(std::bit_width(t) - 1);
52
			uint128 m_low2 =   (uint128(1) << (64 + s)) / t;
53
			uint64_t r = low64((uint128(1) << (64 + s)) % t);
54
			m = low64(m_low2 + ((r <= (t >> 1)) ? 0 : 1));
55
			a = m;
56
		}
57
		// reduce multiplier to smallest possible
58
		while (!(m & 1)) {
59
			m >>= 1;
60
			a >>= 1;
61
			s--;
62
		}
63
	}
64
	// adjust multiplier for reduction of even divisors
65
	s += n;
66
}
67
68
} // namespace openmsx