diff --git a/Compiler/compilerLib.py b/Compiler/compilerLib.py
index 8c4f41b8d..66febc85c 100644
--- a/Compiler/compilerLib.py
+++ b/Compiler/compilerLib.py
@@ -206,6 +206,12 @@ def build_option_parser(self):
             dest="edabit",
             help="mixing arithmetic and binary computation using edaBits",
         )
+        parser.add_option(
+            "--comparison_rabbit",
+            action="store_true",
+            dest="comparison_rabbit",
+            help="using the rabbit comparison protocol for known prime modulus, instead of truncation",
+        )
         parser.add_option(
             "-Z",
             "--split",
diff --git a/Compiler/library.py b/Compiler/library.py
index 5f7e42e83..d695ee8d2 100644
--- a/Compiler/library.py
+++ b/Compiler/library.py
@@ -131,6 +131,25 @@ def print_ln(s='', *args, **kwargs):
     """
     print_str(str(s) + '\n', *args, **kwargs)
 
+
+def print_without_ln(s='', *args, **kwargs):
+    """ Print line, with optional args for adding variables/registers
+    with ``%s``. By default only player 0 outputs, but the ``-I``
+    command-line option changes that.
+
+    :param s: Python string with same number of ``%s`` as length of :py:obj:`args`
+    :param args: list of public values (regint/cint/int/cfix/cfloat/localint)
+    :param print_secrets: whether to output secret shares
+
+    Example:
+
+    .. code::
+
+        print_ln('a is %s.', a.reveal())
+    """
+    print_str(str(s), *args, **kwargs)
+
+
 def print_both(s, end='\n'):
     """ Print line during compilation and execution. """
     print(s, end=end)
diff --git a/Compiler/non_linear.py b/Compiler/non_linear.py
index 147af9f20..fc0f420a0 100644
--- a/Compiler/non_linear.py
+++ b/Compiler/non_linear.py
@@ -44,8 +44,52 @@ def trunc(self, a, k, m, signed):
         if m == 0:
             return a
         return self._trunc(a, k, m, signed)
+    
+    def LTBits(self, R, x, BIT_SIZE):
+        library.print_ln("in LTBits")
+        R_bits = cint.bit_decompose(R, BIT_SIZE)
+        y = [x[i].bit_xor(R_bits[i]) for i in range(BIT_SIZE)]
+        z = floatingpoint.PreOpL(floatingpoint.or_op, y[::-1])[::-1] + [0]
+        w = [z[i] - z[i + 1] for i in range(BIT_SIZE)]
+
+        return types.sintbit(1) - types.sintbit(sum((R_bits[i] & w[i]) for i in range(BIT_SIZE)))
+
+    def rabbitLTZ(self, x, BIT_SIZE = 64):
+        """
+        s = (x <? 0)
+        BIT_SIZE: bit length of x
+        """
+        length_eda = 64 # BIT_SIZE
+
+        M = P_VALUES[64] # TODO: get program.prime
+        R = 0
+
+        r, r_bits = sint.get_edabit(length_eda, True)
+        masked_a = (x + r).reveal()
+        masked_b = (x + r + M - R).reveal()
+        w = [None, None, None, None]
+
+        library.print_ln("w1, comparing: masked_a=%s edabit=%s", masked_a, r.reveal())
+        w[1] = self.LTBits(masked_a, r_bits, BIT_SIZE)
+
+        library.print_ln("w2, comparing: masked_b=%s edabit=%s", masked_b, r.reveal())
+        w[2] = self.LTBits(masked_b, r_bits, BIT_SIZE)
+
+        library.print_ln("w3, comparing: masked_b=%s with zero", masked_a)
+        w[3] = cint(masked_b < 0)
+
+        result = w[1] - w[2] + w[3]
+
+        library.print_ln("w1=%s w2=%s w3=%s result=%s", w[1].reveal(), w[2].reveal(), w[3], result.reveal())
+        return sint(1 - result)
 
     def ltz(self, a, k):
+        library.print_ln("a=%s k=%s", a.reveal(), k)
+        prog = program.Program.prog
+        if prog.options.comparison_rabbit:
+            return self.rabbitLTZ(a, k)
+        
+        # else, use truncation
         return -self.trunc(a, k, k - 1, True)
 
 class Masking(NonLinear):
@@ -126,6 +170,10 @@ def eqz(self, a, k):
         return 1 - types.sintbit.conv(KORL(self.bit_dec(a, k, k, True)))
 
     def ltz(self, a, k):
+        prog = program.Program.prog
+        if prog.options.comparison_rabbit:
+            return -20
+        
         if k + 1 < self.prime.bit_length():
             # https://dl.acm.org/doi/10.1145/3474123.3486757
             # "negative" values wrap around when doubling, thus becoming odd
diff --git a/Compiler/program.py b/Compiler/program.py
index 5d68948ef..7a32d6731 100644
--- a/Compiler/program.py
+++ b/Compiler/program.py
@@ -69,6 +69,7 @@ class defaults:
     stop = False
     insecure = False
     keep_cisc = False
+    comparison_rabbit = False
 
 
 class Program(object):
@@ -203,11 +204,16 @@ def __init__(self, args, options=defaults, name=None):
             gc.inputbvec,
             gc.reveal,
         ]
+        
         self.use_trunc_pr = False
         """ Setting whether to use special probabilistic truncation. """
+
         self.use_dabit = options.mixed
         """ Setting whether to use daBits for non-linear functionality. """
+
         self._edabit = options.edabit
+        self._comparison_rabbit = options.comparison_rabbit
+
         """ Whether to use the low-level INVPERM instruction (only implemented with the assumption of a semi-honest two-party environment)"""
         self._invperm = options.invperm
         self._split = False
@@ -676,6 +682,19 @@ def use_edabit(self, change=None):
         else:
             self._edabit = change
 
+    def use_comparison_rabbit(self, change=None):
+        """Setting whether to use the rabbit comparison protocol (default: false).
+
+        :param change: change setting if not :py:obj:`None`
+        :returns: setting if :py:obj:`change` is :py:obj:`None`
+        """
+        if change is None:
+            if not self._comparison_rabbit:
+                self.relevant_opts.add("comparison_rabbit")
+            return self._comparison_rabbit
+        else:
+            self._comparison_rabbit = change
+
     def use_invperm(self, change=None):
         """ Set whether to use the low-level INVPERM instruction to inverse a permutation (see sint.inverse_permutation). The INVPERM instruction assumes a semi-honest two-party environment. If false, a general protocol implemented in the high-level language is used.
 
@@ -750,6 +769,8 @@ def options_from_args(self):
             self.always_raw(True)
         if "edabit" in self.args:
             self.use_edabit(True)
+        if "comparison_rabbit" in self.args:
+            self.use_comparison_rabbit(True)
         if "invperm" in self.args:
             self.use_invperm(True)
         if "linear_rounds" in self.args:
diff --git a/Compiler/types.py b/Compiler/types.py
index 5caa21051..5be8e105c 100644
--- a/Compiler/types.py
+++ b/Compiler/types.py
@@ -5107,6 +5107,141 @@ def update(self, other):
         assert self.m == other.m
         self.v.update(other.v)
 
+
+class custom_sfix(sfix):
+
+    # R: clear-text; x: edabit in binary format
+    @vectorize
+    def LTBits(self, R, x, r, BIT_SIZE):
+        R_bits = cint.bit_decompose(R, BIT_SIZE)
+        library.print_ln("\nLTBits: R_bits= ")
+        for i in range(BIT_SIZE):
+            library.print_without_ln("%s", R_bits[i])
+
+        edabit = cint.bit_decompose(r.reveal(), BIT_SIZE)
+        library.print_ln("\nLTBits: edabit= ")
+        for i in range(BIT_SIZE):
+            library.print_without_ln("%s", edabit[i].reveal())
+
+        library.print_ln("\nLTBits from bits: x = ")
+        for i in range(BIT_SIZE):
+            library.print_without_ln("%s", x[i].reveal())
+
+        y = [x[i].bit_xor(R_bits[i]) for i in range(BIT_SIZE)]
+        library.print_ln("\nLTBits: y= ")
+        for i in range(BIT_SIZE):
+            library.print_without_ln("%s", y[i].reveal())
+
+        #z = floatingpoint.PreOpL(floatingpoint.or_op, y[::-1])[::-1] + [0]
+        z = floatingpoint.PreOpL(floatingpoint.or_op, y)
+        library.print_ln("\nLTBits: z= ")
+        for i in range(BIT_SIZE):
+            library.print_without_ln("%s", z[i].reveal())
+
+        z = [0] + z
+        w = [z[i] - z[i - 1] for i in range(BIT_SIZE, 0, -1)]
+        w = w[::-1]
+        library.print_ln("\nLTBits: w= ")
+        for i in range(BIT_SIZE):
+            library.print_without_ln("%s", w[i].reveal())
+
+        s = sum((R_bits[i] & w[i]) for i in range(BIT_SIZE))
+        library.print_ln("\nSum=%s", s.reveal())
+        return_value = 1 - sum((R_bits[i] & w[i]) for i in range(BIT_SIZE))
+        return return_value
+
+    @vectorize
+    def rabbitLTZ(self, x, BIT_SIZE = 64):
+        """
+        s = (c ?< a)
+
+        BIT_SIZE: bit length of a
+        """
+        length_eda = BIT_SIZE
+        library.print_ln("custom sfix: bitsize = %s", BIT_SIZE)
+
+        M = P_VALUES[64]
+        R = (M - 1) // 2  # for field; use 0 for ring
+
+        r, r_bits = sint.get_edabit(length_eda, True)
+        masked_a = (x + r).reveal()
+        masked_b = (x + r + M - R).reveal() # masked_a
+        w = [None, None, None, None]
+
+        w[1] = self.LTBits(masked_a, r_bits, r, BIT_SIZE)
+        library.print_ln("w1, comparing: masked_a=%s edabit=%s w1=%s", masked_a, r.reveal(), w[1].reveal())
+        
+        w[2] = self.LTBits(masked_b, r_bits, r, BIT_SIZE)
+        library.print_ln("w2, comparing: masked_b=%s edabit=%s w2=%s", masked_b, r.reveal(), w[2].reveal())
+
+        w[3] = cint(masked_b < 0)
+        library.print_ln("w3, comparing: masked_b=%s with %s, w3=%s", masked_b, M - R, w[3].reveal())
+
+        result = w[1] - w[2] + w[3]
+
+        library.print_ln("final result = %s and 1- result=%s", result.reveal(), (1-result).reveal())
+        return sint(1 - result)
+
+    @vectorize
+    def rabbitLTS_fix(self, a, b):
+        return 1 - self.rabbitLTS(b, a)
+    
+    @vectorize
+    def rabbitLTS(self, a, b):
+        res = self.rabbitLTZ(a - b)
+        return res
+
+
+    # These are based on the implementation
+    # self.rabbitLTS(a, b) = rabbitLTC(a-b, 0)
+    @vectorize
+    def __lt__(self, other):
+        a = self.v
+        b = other.v
+        result = self.rabbitLTS_fix(a, b)
+        return result
+
+    @vectorize
+    def __le__(self, other):
+        a = self.v
+        b = other.v
+        result = 1 - self.rabbitLTS_fix(b, a)
+        return result
+
+    @vectorize
+    def __gt__(self, other):
+        a = self.v
+        b = other.v
+        result = self.rabbitLTS_fix(b, a)
+        return result
+    
+    @vectorize
+    def __ge__(self, other):
+        a = self.v
+        b = other.v
+        result = 1 - self.rabbitLTS_fix(a, b)
+        return result
+    
+    @vectorize
+    def __eq__(self, other):
+        a = self.v
+        b = other.v
+        result = (1 - self.rabbitLTS_fix(a, b)) * (1 - self.rabbitLTS_fix(b, a))
+        return result
+    
+    @vectorize
+    def __ne__(self, other):
+        a = self.v
+        b = other.v
+        result = 1 - (1 - self.rabbitLTS_fix(a, b)) * (1 - self.rabbitLTS_fix(b, a))
+        return result
+
+    @vectorize
+    def __abs__(self):
+        """ Absolute value. """
+        return (self < custom_sfix(0)).if_else(-self, self)
+
+
 sfix.unreduced_type = unreduced_sfix
 
 sfix.set_precision(16, 31)