From f534c983cd368a85a40a6061e74a93cfb9b4dac4 Mon Sep 17 00:00:00 2001
From: RunasSudo <runassudo@yingtongli.me>
Date: Thu, 14 Jan 2021 22:31:24 +1100
Subject: [PATCH] Implement guarded fixed-point arithmetic

---
 README.md                   |  2 +-
 docs/options.md             |  1 +
 html/index.html             |  1 +
 html/index.js               |  2 +-
 html/worker.js              |  4 ++
 pyRCV2/cli/stv.py           |  5 +-
 pyRCV2/method/meek.py       |  9 +++-
 pyRCV2/numbers/__init__.py  | 24 +++++++++-
 pyRCV2/numbers/base.py      |  4 ++
 pyRCV2/numbers/fixed_js.py  |  5 +-
 pyRCV2/numbers/fixed_py.py  | 11 ++---
 pyRCV2/numbers/gfixed_js.py | 93 +++++++++++++++++++++++++++++++++++++
 pyRCV2/numbers/gfixed_py.py | 92 ++++++++++++++++++++++++++++++++++++
 tests/test_meekm.py         |  4 +-
 tests/test_numbers.py       |  6 +++
 15 files changed, 243 insertions(+), 20 deletions(-)
 create mode 100644 pyRCV2/numbers/gfixed_js.py
 create mode 100644 pyRCV2/numbers/gfixed_py.py

diff --git a/README.md b/README.md
index 47d5d26..b06ccb5 100644
--- a/README.md
+++ b/README.md
@@ -21,7 +21,7 @@ pyRCV2 accepts data in the [BLT file format](http://www.dia.govt.nz/diawebsite.N
 pyRCV2 is highly customisable, including options for:
 
 * different quotas and quota rules (e.g. exact Droop, Hare) or progressively reducing quota
-* calculations using fixed-point arithmetic or exact rational numbers
+* calculations using fixed-point arithmetic, guarded fixed-point ([quasi-exact](http://www.votingmatters.org.uk/ISSUE24/I24P2.pdf)) or exact rational numbers
 * different tie breaking rules (backwards, random, manual) with auditable deterministic random number generation
 * extensible API for other counting methods
 
diff --git a/docs/options.md b/docs/options.md
index 0bb394e..06505d8 100644
--- a/docs/options.md
+++ b/docs/options.md
@@ -18,6 +18,7 @@ This functionality is not available on the Python command line.
 This dropdown allows you to select how numbers (vote totals, etc.) are represented internally in memory. The options are:
 
 * *Fixed*: Numbers are represented as fixed-precision decimals, up to a certain number of decimal places (default: 5).
+* *Fixed (guarded)*: Numbers are represented as fixed-precision decimals with ‘guard digits’ – also known as [‘quasi-exact’ arithmetic](http://www.votingmatters.org.uk/ISSUE24/I24P2.pdf). If *n* decimal places are requested, numbers are represented up to 2*n* decimal places, and two values are considered equal if the absolute difference is less than (10^-*n*)/2.
 * *Rational*: Numbers are represented exactly as fractions, resulting in the elimination of rounding error, but increasing computational complexity when the number of surplus transfers is very large.
 * *Native*: Numbers are represented as native integers or floating-point numbers. This is fast, but not recommended as unexpectedly large rounding errors may be introduced in some circumstances.
 
diff --git a/html/index.html b/html/index.html
index f3aa047..80ea529 100644
--- a/html/index.html
+++ b/html/index.html
@@ -139,6 +139,7 @@
 							<option value="native">Native</option>
 							<option value="rational">Rational</option>
 							<option value="fixed" selected>Fixed</option>
+							<option value="gfixed">Fixed (guarded)</option>
 						</select>
 					</label>
 					<label>
diff --git a/html/index.js b/html/index.js
index 2b832e3..5af1eff 100644
--- a/html/index.js
+++ b/html/index.js
@@ -93,7 +93,7 @@ function changePreset() {
 		document.getElementById('chkBulkElection').checked = true;
 		document.getElementById('chkBulkExclusion').checked = false;
 		document.getElementById('chkDeferSurpluses').checked = false;
-		document.getElementById('selNumbers').value = 'fixed';
+		document.getElementById('selNumbers').value = 'gfixed';
 		document.getElementById('txtDP').value = '9';
 		document.getElementById('txtPPDP').value = '2';
 		document.getElementById('chkRoundQuota').checked = false;
diff --git a/html/worker.js b/html/worker.js
index e382fbb..ab6beaa 100644
--- a/html/worker.js
+++ b/html/worker.js
@@ -33,6 +33,10 @@ onmessage = function(evt) {
 			py.pyRCV2.numbers.set_numclass(py.pyRCV2.numbers.Fixed);
 			py.pyRCV2.numbers.set_dps(evt.data.data.fixedDPs);
 		}
+		if (evt.data.data.numbers === 'gfixed') {
+			py.pyRCV2.numbers.set_numclass(py.pyRCV2.numbers.FixedGuarded);
+			py.pyRCV2.numbers.set_dps(evt.data.data.fixedDPs);
+		}
 		
 		ppDPs = evt.data.data.ppDPs;
 		
diff --git a/pyRCV2/cli/stv.py b/pyRCV2/cli/stv.py
index e97040a..22ec3ea 100644
--- a/pyRCV2/cli/stv.py
+++ b/pyRCV2/cli/stv.py
@@ -34,7 +34,7 @@ def add_parser(subparsers):
 	parser.add_argument('--no-bulk-elect', action='store_true', help='disable bulk election unless absolutely required')
 	parser.add_argument('--bulk-exclude', action='store_true', help='use bulk exclusion')
 	parser.add_argument('--defer-surpluses', action='store_true', help='defer surplus transfers if possible')
-	parser.add_argument('--numbers', '-n', choices=['fixed', 'rational', 'native'], default='fixed', help='numbers mode (default: fixed)')
+	parser.add_argument('--numbers', '-n', choices=['fixed', 'gfixed', 'rational', 'native'], default='fixed', help='numbers mode (default: fixed)')
 	parser.add_argument('--decimals', type=int, default=5, help='decimal places if --numbers fixed (default: 5)')
 	#parser.add_argument('--no-round-quota', action='store_true', help='do not round the quota')
 	parser.add_argument('--round-quota', type=int, help='round quota to specified decimal places')
@@ -102,6 +102,9 @@ def main(args):
 	elif args.numbers == 'fixed':
 		pyRCV2.numbers.set_numclass(pyRCV2.numbers.Fixed)
 		pyRCV2.numbers.set_dps(args.decimals)
+	elif args.numbers == 'gfixed':
+		pyRCV2.numbers.set_numclass(pyRCV2.numbers.FixedGuarded)
+		pyRCV2.numbers.set_dps(args.decimals)
 	
 	with open(args.file, 'r') as f:
 		election = pyRCV2.blt.readBLT(f.read())
diff --git a/pyRCV2/method/meek.py b/pyRCV2/method/meek.py
index f22519d..64a8052 100644
--- a/pyRCV2/method/meek.py
+++ b/pyRCV2/method/meek.py
@@ -14,6 +14,8 @@
 #   You should have received a copy of the GNU Affero General Public License
 #   along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
+DEBUG_MEEK = False
+
 __pragma__ = lambda x: None
 
 from pyRCV2.method.base_stv import BaseSTVCounter, STVException
@@ -170,7 +172,9 @@ class MeekSTVCounter(BaseSTVCounter):
 				# Recompute keep values
 				for candidate, count_card in has_surplus:
 					__pragma__('opov')
-					count_card.keep_value *= self.quota / count_card.votes
+					# Perform in steps to avoid rounding error
+					count_card.keep_value *= self.quota
+					count_card.keep_value /= count_card.votes
 					__pragma__('noopov')
 				
 				# Redistribute votes
@@ -182,6 +186,9 @@ class MeekSTVCounter(BaseSTVCounter):
 				__pragma__('opov')
 				has_surplus = [(c, cc) for c, cc in self.candidates.items() if cc.state == CandidateState.ELECTED and cc.votes / self.quota > self._quota_tolerance]
 				__pragma__('noopov')
+				
+				if DEBUG_MEEK:
+					break
 			
 			if num_iterations == 1:
 				self.logs.append('Surpluses distributed, requiring 1 iteration.')
diff --git a/pyRCV2/numbers/__init__.py b/pyRCV2/numbers/__init__.py
index 4ff3ae9..7ed999b 100644
--- a/pyRCV2/numbers/__init__.py
+++ b/pyRCV2/numbers/__init__.py
@@ -20,17 +20,23 @@ __pragma__('skip')
 is_py = True
 __pragma__('noskip')
 
+# CLASS DEFINITIONS
+
 if is_py:
 	__pragma__('skip')
-	from pyRCV2.numbers.fixed_py import Fixed, set_dps, get_dps
+	from pyRCV2.numbers.fixed_py import Fixed, _fixed_set_dps
+	from pyRCV2.numbers.gfixed_py import FixedGuarded, _gfixed_set_dps
 	from pyRCV2.numbers.native_py import Native
 	from pyRCV2.numbers.rational_py import Rational
 	__pragma__('noskip')
 else: # pragma: no cover
-	from pyRCV2.numbers.fixed_js import Fixed, set_dps, get_dps
+	from pyRCV2.numbers.fixed_js import Fixed, _fixed_set_dps
+	from pyRCV2.numbers.gfixed_js import FixedGuarded, _gfixed_set_dps
 	from pyRCV2.numbers.native_js import Native
 	from pyRCV2.numbers.rational_js import Rational
 
+# GLOBALS
+
 _numclass = Native
 
 def set_numclass(cls):
@@ -39,3 +45,17 @@ def set_numclass(cls):
 
 def Num(val):
 	return _numclass(val)
+
+_dps = 6
+
+def set_dps(dps):
+	global _dps
+	_dps = dps
+	
+	if _numclass is Fixed:
+		_fixed_set_dps(dps)
+	elif _numclass is FixedGuarded:
+		_gfixed_set_dps(dps)
+
+def get_dps():
+	return _dps
diff --git a/pyRCV2/numbers/base.py b/pyRCV2/numbers/base.py
index 1cd7563..9dbe759 100644
--- a/pyRCV2/numbers/base.py
+++ b/pyRCV2/numbers/base.py
@@ -158,6 +158,10 @@ class BaseNum:
 		self.impl = self._truediv_impl(self.impl, other.impl)
 		return self
 	
+	def __idiv__(self, other):
+		# For Transcrypt
+		return self.__itruediv__(other)
+	
 	def __floor__(self):
 		return self.round(0, self.ROUND_DOWN)
 
diff --git a/pyRCV2/numbers/fixed_js.py b/pyRCV2/numbers/fixed_js.py
index f0ccc75..5e1c7c0 100644
--- a/pyRCV2/numbers/fixed_js.py
+++ b/pyRCV2/numbers/fixed_js.py
@@ -18,12 +18,9 @@ from pyRCV2.numbers.base import BaseNum, compatible_types
 
 Big.DP = 6
 
-def set_dps(dps):
+def _fixed_set_dps(dps):
 	Big.DP = dps
 
-def get_dps():
-	return Big.DP
-
 class Fixed(BaseNum):
 	"""
 	Wrapper for big.js (fixed-point arithmetic)
diff --git a/pyRCV2/numbers/fixed_py.py b/pyRCV2/numbers/fixed_py.py
index c3a0ea0..8abde83 100644
--- a/pyRCV2/numbers/fixed_py.py
+++ b/pyRCV2/numbers/fixed_py.py
@@ -18,17 +18,12 @@ from pyRCV2.numbers.base import BasePyNum, compatible_types
 
 import decimal
 
-_dps = 6
-_quantize_exp = decimal.Decimal('10') ** -_dps
+_quantize_exp = decimal.Decimal('10') ** -6
 
-def set_dps(dps):
-	global _dps, _quantize_exp
-	_dps = dps
+def _fixed_set_dps(dps):
+	global _quantize_exp
 	_quantize_exp = decimal.Decimal('10') ** -dps
 
-def get_dps():
-	return _dps
-
 class Fixed(BasePyNum):
 	"""
 	Wrapper for Python Decimal (for fixed-point arithmetic)
diff --git a/pyRCV2/numbers/gfixed_js.py b/pyRCV2/numbers/gfixed_js.py
new file mode 100644
index 0000000..1e3be35
--- /dev/null
+++ b/pyRCV2/numbers/gfixed_js.py
@@ -0,0 +1,93 @@
+#   pyRCV2: Preferential vote counting
+#   Copyright © 2020–2021  Lee Yingtong Li (RunasSudo)
+#
+#   This program is free software: you can redistribute it and/or modify
+#   it under the terms of the GNU Affero General Public License as published by
+#   the Free Software Foundation, either version 3 of the License, or
+#   (at your option) any later version.
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#   GNU Affero General Public License for more details.
+#
+#   You should have received a copy of the GNU Affero General Public License
+#   along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+from pyRCV2.numbers.base import BaseNum, compatible_types
+
+Big.DP = 12
+_pos_epsilon = Big('10').pow(-6).div('2')
+_neg_epsilon = _pos_epsilon.times('-1')
+
+def _gfixed_set_dps(dps):
+	global _pos_epsilon, _neg_epsilon
+	Big.DP = 2 * dps
+	_pos_epsilon = Big('10').pow(-dps).div('2')
+	_neg_epsilon = _pos_epsilon.times('-1')
+
+class FixedGuarded(BaseNum):
+	"""
+	Wrapper for big.js (fixed-point arithmetic)
+	Implements guarded (quasi-exact) fixed-point
+	"""
+	
+	@classmethod
+	def _to_impl(cls, value):
+		"""Implements BaseNum._to_impl"""
+		return Big(value).round(Big.DP)
+	
+	def pp(self, dp):
+		"""Implements BaseNum.pp"""
+		return self.impl.toFixed(dp)
+	
+	@classmethod
+	def _add_impl(cls, i1, i2):
+		"""Implements BaseNum._add_impl"""
+		return i1.plus(i2)
+	@classmethod
+	def _sub_impl(cls, i1, i2):
+		"""Implements BaseNum._sub_impl"""
+		return i1.minus(i2)
+	@classmethod
+	def _mul_impl(cls, i1, i2):
+		"""Implements BaseNum._mul_impl"""
+		return i1.times(i2)
+	@classmethod
+	def _truediv_impl(cls, i1, i2):
+		"""Implements BaseNum._truediv_impl"""
+		return i1.div(i2)
+	
+	@compatible_types
+	def __eq__(self, other):
+		"""Implements BaseNum.__eq__"""
+		diff = self.impl.minus(other.impl)
+		return diff.lt(_pos_epsilon) and diff.gt(_neg_epsilon)
+	@compatible_types
+	def __gt__(self, other):
+		"""Implements BaseNum.__gt__"""
+		diff = self.impl.minus(other.impl)
+		return diff.gt(_pos_epsilon)
+	@compatible_types
+	def __ge__(self, other):
+		"""Implements BaseNum.__ge__"""
+		diff = self.impl.minus(other.impl)
+		return diff.gt(_neg_epsilon)
+	@compatible_types
+	def __lt__(self, other):
+		"""Implements BaseNum.__lt__"""
+		diff = self.impl.minus(other.impl)
+		return diff.lt(_neg_epsilon)
+	@compatible_types
+	def __le__(self, other):
+		"""Implements BaseNum.__le__"""
+		diff = self.impl.minus(other.impl)
+		return diff.lt(_pos_epsilon)
+	
+	def __pow__(self, power):
+		"""Implements BaseNum.__pow__"""
+		return FixedGuarded._from_impl(self.impl.pow(power))
+	
+	def round(self, dps, mode):
+		"""Implements BaseNum.round"""
+		return FixedGuarded._from_impl(self.impl.round(dps, mode))
diff --git a/pyRCV2/numbers/gfixed_py.py b/pyRCV2/numbers/gfixed_py.py
new file mode 100644
index 0000000..49afb08
--- /dev/null
+++ b/pyRCV2/numbers/gfixed_py.py
@@ -0,0 +1,92 @@
+#   pyRCV2: Preferential vote counting
+#   Copyright © 2020–2021  Lee Yingtong Li (RunasSudo)
+#
+#   This program is free software: you can redistribute it and/or modify
+#   it under the terms of the GNU Affero General Public License as published by
+#   the Free Software Foundation, either version 3 of the License, or
+#   (at your option) any later version.
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#   GNU Affero General Public License for more details.
+#
+#   You should have received a copy of the GNU Affero General Public License
+#   along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+from pyRCV2.numbers.base import BasePyNum, compatible_types
+
+import decimal
+
+_quantize_exp = decimal.Decimal('10') ** -12
+_pos_epsilon = decimal.Decimal('10') ** -6 / decimal.Decimal('2')
+_neg_epsilon = -_pos_epsilon
+
+def _gfixed_set_dps(dps):
+	global _quantize_exp
+	_quantize_exp = decimal.Decimal('10') ** -(2 * dps)
+	_pos_epsilon = decimal.Decimal('10') ** -dps / decimal.Decimal('2')
+	_neg_epsilon = -_pos_epsilon
+
+class FixedGuarded(BasePyNum):
+	"""
+	Wrapper for Python Decimal (for fixed-point arithmetic)
+	Implements guarded (quasi-exact) fixed-point
+	"""
+	
+	_py_class = decimal.Decimal # For BasePyNum
+	
+	ROUND_DOWN = decimal.ROUND_DOWN
+	ROUND_HALF_UP = decimal.ROUND_HALF_UP
+	ROUND_HALF_EVEN = decimal.ROUND_HALF_EVEN
+	ROUND_UP = decimal.ROUND_UP
+	
+	@classmethod
+	def _to_impl(cls, value):
+		"""Overrides BasePyNum._to_impl"""
+		return decimal.Decimal(value).quantize(_quantize_exp)
+	
+	@classmethod
+	def _truediv_impl(cls, i1, i2):
+		"""Overrides BasePyNum._truediv_impl"""
+		return (i1 / i2).quantize(_quantize_exp)
+	
+	@compatible_types
+	def __eq__(self, other):
+		"""Overrides BasePyNum.__eq__"""
+		diff = self.impl - other.impl
+		return diff < _pos_epsilon and diff > _neg_epsilon
+	@compatible_types
+	def __ne__(self, other):
+		"""Overrides BasePyNum.__ne__"""
+		return not self.__eq__(other)
+	@compatible_types
+	def __gt__(self, other):
+		"""Overrides BasePyNum.__gt__"""
+		diff = self.impl - other.impl
+		return diff > _pos_epsilon
+	@compatible_types
+	def __ge__(self, other):
+		"""Overrides BasePyNum.__ge__"""
+		diff = self.impl - other.impl
+		return diff > _neg_epsilon
+	@compatible_types
+	def __lt__(self, other):
+		"""Overrides BasePyNum.__lt__"""
+		diff = self.impl - other.impl
+		return diff < _neg_epsilon
+	@compatible_types
+	def __le__(self, other):
+		"""Overrides BasePyNum.__le__"""
+		diff = self.impl - other.impl
+		return diff < _pos_epsilon
+	
+	@compatible_types
+	def __itruediv__(self, other):
+		"""Overrides BasePyNum.__itruediv__"""
+		self.impl = (self.impl / other.impl).quantize(_quantize_exp)
+		return self
+	
+	def round(self, dps, mode):
+		"""Implements BaseNum.round"""
+		return FixedGuarded._from_impl(self.impl.quantize(decimal.Decimal('10') ** -dps, mode))
diff --git a/tests/test_meekm.py b/tests/test_meekm.py
index f4fad4e..e2dd3f6 100644
--- a/tests/test_meekm.py
+++ b/tests/test_meekm.py
@@ -26,8 +26,8 @@ from pyRCV2.ties import TiesBackwards
 def test_meekm():
 	"""Compare count of meekm.blt with model result produced by Hill et al. reference implementation"""
 	
-	pyRCV2.numbers.set_numclass(pyRCV2.numbers.Fixed)
-	pyRCV2.numbers.set_dps(9)
+	pyRCV2.numbers.set_numclass(pyRCV2.numbers.FixedGuarded)
+	pyRCV2.numbers.set_dps(4)
 	
 	with open('tests/data/meekm.blt', 'r') as f:
 		election = pyRCV2.blt.readBLT(f.read())
diff --git a/tests/test_numbers.py b/tests/test_numbers.py
index f2cc405..fb6a35b 100644
--- a/tests/test_numbers.py
+++ b/tests/test_numbers.py
@@ -44,6 +44,7 @@ def maketst(numbers, dps, method, result):
 test_fixed2_add_py, test_fixed2_add_js = maketst('Fixed', 2, '__add__', '356.57')
 test_fixed1_add_py, test_fixed1_add_js = maketst('Fixed', 1, '__add__', '356.6')
 test_fixed0_add_py, test_fixed0_add_js = maketst('Fixed', 0, '__add__', '356')
+test_gfixed2_add_py, test_gfixed2_add_js = maketst('FixedGuarded', 2, '__add__', '356.57')
 test_rational_add_py, test_rational_add_js = maketst('Rational', 0, '__add__', '356.57')
 
 def maketst_round(numbers, dps, num, dps_round, mode_round, result):
@@ -70,6 +71,11 @@ test_fixed_round2_py, test_fixed_round2_js = maketst_round('Fixed', 5, '3141.59'
 test_fixed_round3_py, test_fixed_round3_js = maketst_round('Fixed', 5, '3141.45', 1, 'ROUND_HALF_EVEN', '3141.4')
 test_fixed_round4_py, test_fixed_round4_js = maketst_round('Fixed', 5, '3141.45', 1, 'ROUND_HALF_UP', '3141.5')
 
+test_gfixed_round1_py, test_gfixed_round1_js = maketst_round('FixedGuarded', 5, '3141.59', 1, 'ROUND_DOWN', '3141.5')
+test_gfixed_round2_py, test_gfixed_round2_js = maketst_round('FixedGuarded', 5, '3141.59', 1, 'ROUND_UP', '3141.6')
+test_gfixed_round3_py, test_gfixed_round3_js = maketst_round('FixedGuarded', 5, '3141.45', 1, 'ROUND_HALF_EVEN', '3141.4')
+test_gfixed_round4_py, test_gfixed_round4_js = maketst_round('FixedGuarded', 5, '3141.45', 1, 'ROUND_HALF_UP', '3141.5')
+
 test_native_round1_py, test_native_round1_js = maketst_round('Native', 0, '3141.59', 1, 'ROUND_DOWN', '3141.5')
 test_native_round2_py, test_native_round2_js = maketst_round('Native', 0, '3141.59', 1, 'ROUND_UP', '3141.6')
 test_native_round3_py, _ = maketst_round('Native', 0, '3141.45', 1, 'ROUND_HALF_EVEN', '3141.4')