Improve PTP model resolution by using Decimal types and wider interna…

…l fns accumulators

Signed-off-by: Alex Forencich <[email protected]>

README.md

@@ -603,6 +603,10 @@ To use this module, import it and connect it to the DUT:
 
 Once the clock is instantiated, it will generate a continuous stream of monotonically increasing PTP timestamps on every clock edge.
 
+Internally, the `PtpClock` module uses 32-bit fractional ns fields for higher frequency resolution.  Only the upper 16 bits are returned in the timestamps, but the full fns value can be accessed with the _ts_tod_fns_ and _ts_rel_fns_ attributes.
+
+All APIs that handle fractional values use the `Decimal` type for maximum precision, as the combination of timestamp range and resolution is usually too much for normal floating point numbers to handle without significant loss of precision.
+
 #### Signals
 
 * `ts_tod`: 96-bit time-of-day timestamp (48 bit seconds, 32 bit ns, 16 bit fractional ns)
@@ -633,24 +637,26 @@ Once the clock is instantiated, it will generate a continuous stream of monotoni
 
 * `set_period(ns, fns)`: set clock period from separate fields
 * `set_drift(num, denom)`: set clock drift from separate fields
-* `set_period_ns(t)`: set clock period in ns (float)
-* `get_period_ns()`: return current clock period in ns (float)
+* `set_period_ns(t)`: set clock period and drift in ns (Decimal)
+* `get_period_ns()`: return current clock period in ns (Decimal)
 * `set_ts_tod(ts_s, ts_ns, ts_fns)`: set 96-bit ToD timestamp from separate fields
 * `set_ts_tod_96(ts)`: set 96-bit ToD timestamp from integer
-* `set_ts_tod_ns(t)`: set 96-bit ToD timestamp from ns (float)
-* `set_ts_tod_s(t)`: set 96-bit ToD timestamp from seconds (float)
+* `set_ts_tod_ns(t)`: set 96-bit ToD timestamp from ns (Decimal)
+* `set_ts_tod_s(t)`: set 96-bit ToD timestamp from seconds (Decimal)
+* `set_ts_tod_sim_time()`: set 96-bit ToD timestamp from sim time
 * `get_ts_tod()`: return current 96-bit ToD timestamp as separate fields
 * `get_ts_tod_96()`: return current 96-bit ToD timestamp as an integer
-* `get_ts_tod_ns()`: return current 96-bit ToD timestamp in ns (float)
-* `get_ts_tod_s()`: return current 96-bit ToD timestamp in seconds (float)
+* `get_ts_tod_ns()`: return current 96-bit ToD timestamp in ns (Decimal)
+* `get_ts_tod_s()`: return current 96-bit ToD timestamp in seconds (Decimal)
 * `set_ts_rel(ts_ns, ts_fns)`: set 64-bit relative timestamp from separate fields
 * `set_ts_rel_64(ts)`: set 64-bit relative timestamp from integer
-* `set_ts_rel_ns(t)`: set 64-bit relative timestamp from ns (float)
-* `set_ts_rel_s(t)`: set 64-bit relative timestamp from seconds (float)
+* `set_ts_rel_ns(t)`: set 64-bit relative timestamp from ns (Decimal)
+* `set_ts_rel_s(t)`: set 64-bit relative timestamp from seconds (Decimal)
+* `set_ts_rel_sim_time()`: set 64-bit relative timestamp from sim time
 * `get_ts_rel()`: return current 64-bit relative timestamp as separate fields
 * `get_ts_rel_64()`: return current 64-bit relative timestamp as an integer
-* `get_ts_rel_ns()`: return current 64-bit relative timestamp in ns (float)
-* `get_ts_rel_s()`: return current 64-bit relative timestamp in seconds (float)
+* `get_ts_rel_ns()`: return current 64-bit relative timestamp in ns (Decimal)
+* `get_ts_rel_s()`: return current 64-bit relative timestamp in seconds (Decimal)
 
 ### PTP clock (sim time)
 
@@ -669,9 +675,11 @@ To use this module, import it and connect it to the DUT:
 
 Once the clock is instantiated, it will generate a continuous stream of monotonically increasing PTP timestamps on every clock edge.
 
+All APIs that handle fractional values use the `Decimal` type for maximum precision, as the combination of timestamp range and resolution is usually too much for normal floating point numbers to handle without significant loss of precision.
+
 #### Signals
 
-* `ts_tod`: 96-bit ToD timestamp (48 bit seconds, 32 bit ns, 16 bit fractional ns)
+* `ts_tod`: 96-bit time-of-day timestamp (48 bit seconds, 32 bit ns, 16 bit fractional ns)
 * `ts_rel`: 64-bit relative timestamp (48 bit ns, 16 bit fractional ns)
 * `pps`: pulse-per-second output, pulsed when ts_tod seconds field increments
 
@@ -694,9 +702,9 @@ Once the clock is instantiated, it will generate a continuous stream of monotoni
 
 * `get_ts_tod()`: return current 96-bit ToD timestamp as separate fields
 * `get_ts_tod_96()`: return current 96-bit ToD timestamp as an integer
-* `get_ts_tod_ns()`: return current 96-bit ToD timestamp in ns (float)
-* `get_ts_tod_s()`: return current 96-bit ToD timestamp in seconds (float)
+* `get_ts_tod_ns()`: return current 96-bit ToD timestamp in ns (Decimal)
+* `get_ts_tod_s()`: return current 96-bit ToD timestamp in seconds (Decimal)
 * `get_ts_rel()`: return current 64-bit relative timestamp as separate fields
 * `get_ts_rel_96()`: return current 64-bit relative timestamp as an integer
-* `get_ts_rel_ns()`: return current 64-bit relative timestamp in ns (float)
-* `get_ts_rel_s()`: return current 64-bit relative timestamp in seconds (float)
+* `get_ts_rel_ns()`: return current 64-bit relative timestamp in ns (Decimal)
+* `get_ts_rel_s()`: return current 64-bit relative timestamp in seconds (Decimal)

cocotbext/eth/ptp.py

@@ -23,7 +23,7 @@
 """
 
 import logging
-import math
+from decimal import Decimal, Context
 from fractions import Fraction
 
 import cocotb
@@ -56,20 +56,22 @@ def __init__(
         self.clock = clock
         self.reset = reset
 
-        self.period_ns = 0
-        self.period_fns = 0
-        self.drift_num = 0
-        self.drift_denom = 0
-        self.drift_cnt = 0
-        self.set_period_ns(period_ns)
-
         self.log.info("PTP clock")
         self.log.info("cocotbext-eth version %s", __version__)
         self.log.info("Copyright (c) 2020 Alex Forencich")
         self.log.info("https://github.com/alexforencich/cocotbext-eth")
 
         super().__init__(*args, **kwargs)
 
+        self.ctx = Context(prec=60)
+
+        self.period_ns = 0
+        self.period_fns = 0
+        self.drift_num = 0
+        self.drift_denom = 0
+        self.drift_cnt = 0
+        self.set_period_ns(period_ns)
+
         self.ts_tod_s = 0
         self.ts_tod_ns = 0
         self.ts_tod_fns = 0
@@ -94,26 +96,29 @@ def __init__(
 
     def set_period(self, ns, fns):
         self.period_ns = int(ns)
-        self.period_fns = int(fns) & 0xffff
+        self.period_fns = int(fns) & 0xffffffff
 
     def set_drift(self, num, denom):
         self.drift_num = int(num)
         self.drift_denom = int(denom)
 
     def set_period_ns(self, t):
-        drift, period = math.modf(t*2**16)
+        t = Decimal(t)
+        period, drift = self.ctx.divmod(Decimal(t) * Decimal(2**32), Decimal(1))
         period = int(period)
-        frac = Fraction(drift).limit_denominator(2**16)
-        self.period_ns = period >> 16
-        self.period_fns = period & 0xffff
-        self.drift_num = frac.numerator
-        self.drift_denom = frac.denominator
+        frac = Fraction(drift).limit_denominator(2**16-1)
+        self.set_period(period >> 32, period & 0xffffffff)
+        self.set_drift(frac.numerator, frac.denominator)
+
+        self.log.info("Set period: %s ns", t)
+        self.log.info("Period: 0x%x ns 0x%08x fns", self.period_ns, self.period_fns)
+        self.log.info("Drift: 0x%04x / 0x%04x fns", self.drift_num, self.drift_denom)
 
     def get_period_ns(self):
-        p = ((self.period_ns << 16) | self.period_fns) / 2**16
+        p = Decimal((self.period_ns << 32) | self.period_fns)
         if self.drift_denom:
-            return p + self.drift_num / self.drift_rate / 2**16
-        return p
+            p += Decimal(self.drift_num) / Decimal(self.drift_denom)
+        return p / Decimal(2**32)
 
     def set_ts_tod(self, ts_s, ts_ns, ts_fns):
         self.ts_tod_s = int(ts_s)
@@ -123,31 +128,37 @@ def set_ts_tod(self, ts_s, ts_ns, ts_fns):
 
     def set_ts_tod_96(self, ts):
         ts = int(ts)
-        self.set_ts_tod(ts >> 48, (ts >> 32) & 0x3fffffff, ts & 0xffff)
+        self.set_ts_tod(ts >> 48, (ts >> 32) & 0x3fffffff, (ts & 0xffff) << 16)
 
     def set_ts_tod_ns(self, t):
-        self.set_ts_tod_s(t*1e-9)
+        ts_s, ts_ns = self.ctx.divmod(Decimal(t), Decimal(1000000000))
+        ts_s = ts_s.scaleb(-9).to_integral_value()
+        ts_ns, ts_fns = self.ctx.divmod(ts_ns, Decimal(1))
+        ts_ns = ts_ns.to_integral_value()
+        ts_fns = (ts_fns * Decimal(2**32)).to_integral_value()
+        self.set_ts_tod(ts_s, ts_ns, ts_fns)
 
     def set_ts_tod_s(self, t):
-        ts_ns, ts_s = math.modf(t)
-        ts_ns *= 1e9
-        ts_fns, ts_ns = math.modf(ts_ns)
-        ts_fns *= 2**16
-        self.set_ts_tod(ts_s, ts_ns, ts_fns)
+        self.set_ts_tod_ns(Decimal(t).scaleb(9, self.ctx))
+
+    def set_ts_tod_sim_time(self):
+        self.set_ts_tod_ns(Decimal(get_sim_time('fs')).scaleb(-6))
 
     def get_ts_tod(self):
         return (self.ts_tod_s, self.ts_tod_ns, self.ts_tod_fns)
 
     def get_ts_tod_96(self):
         ts_s, ts_ns, ts_fns = self.get_ts_tod()
-        return (ts_s << 48) | (ts_ns << 16) | ts_fns
+        return (ts_s << 48) | (ts_ns << 16) | (ts_fns >> 16)
 
     def get_ts_tod_ns(self):
         ts_s, ts_ns, ts_fns = self.get_ts_tod()
-        return ts_s*1e9+ts_ns+ts_fns/2**16
+        ns = Decimal(ts_fns) / Decimal(2**32)
+        ns = self.ctx.add(ns, Decimal(ts_ns))
+        return self.ctx.add(ns, Decimal(ts_s).scaleb(9))
 
     def get_ts_tod_s(self):
-        return self.get_ts_tod_ns()*1e-9
+        return self.get_ts_tod_ns().scaleb(-9, self.ctx)
 
     def set_ts_rel(self, ts_ns, ts_fns):
         self.ts_rel_ns = int(ts_ns)
@@ -159,12 +170,16 @@ def set_ts_rel_64(self, ts):
         self.set_ts_rel(ts >> 16, (ts & 0xffff) << 16)
 
     def set_ts_rel_ns(self, t):
-        ts_fns, ts_ns = math.modf(t)
-        ts_fns *= 2**16
+        ts_ns, ts_fns = self.ctx.divmod(Decimal(t), Decimal(1))
+        ts_ns = ts_ns.to_integral_value()
+        ts_fns = (ts_fns * Decimal(2**32)).to_integral_value()
         self.set_ts_rel(ts_ns, ts_fns)
 
     def set_ts_rel_s(self, t):
-        self.set_ts_rel_ns(t*1e9)
+        self.set_ts_rel_ns(Decimal(t).scaleb(9, self.ctx))
+
+    def set_ts_rel_sim_time(self):
+        self.set_ts_rel_ns(Decimal(get_sim_time('fs')).scaleb(-6))
 
     def get_ts_rel(self):
         return (self.ts_rel_ns, self.ts_rel_fns)
@@ -175,10 +190,10 @@ def get_ts_rel_64(self):
 
     def get_ts_rel_ns(self):
         ts_ns, ts_fns = self.get_ts_rel()
-        return ts_ns + ts_fns/2**16
+        return self.ctx.add(Decimal(ts_fns) / Decimal(2**32), Decimal(ts_ns))
 
     def get_ts_rel_s(self):
-        return self.get_ts_rel()*1e-9
+        return self.get_ts_rel_ns().scaleb(-9, self.ctx)
 
     def _handle_reset(self, state):
         if state:
@@ -219,36 +234,39 @@ async def _run(self):
             if self.pps is not None:
                 self.pps.value = 0
 
-            # increment 96 bit timestamp
-            if self.ts_tod is not None or self.pps is not None:
-                t = ((self.ts_tod_ns << 16) + self.ts_tod_fns) + ((self.period_ns << 16) + self.period_fns)
+            # increment tod bit timestamp
+            self.ts_tod_fns += (self.period_ns << 32) + self.period_fns
 
-                if self.drift_denom and self.drift_cnt == 0:
-                    t += self.drift_num
+            if self.drift_denom and self.drift_cnt == 0:
+                self.ts_tod_fns += self.drift_num
 
-                if t > (1000000000 << 16):
-                    self.ts_tod_s += 1
-                    t -= (1000000000 << 16)
-                    if self.pps is not None:
-                        self.pps.value = 1
+            ns_inc = self.ts_tod_fns >> 32
+            self.ts_tod_fns &= 0xffffffff
 
-                self.ts_tod_fns = t & 0xffff
-                self.ts_tod_ns = t >> 16
+            self.ts_tod_ns += ns_inc
 
-                if self.ts_tod is not None:
-                    self.ts_tod.value = (self.ts_tod_s << 48) | (self.ts_tod_ns << 16) | (self.ts_tod_fns)
+            if self.ts_tod_ns >= 1000000000:
+                self.ts_tod_s += 1
+                self.ts_tod_ns -= 1000000000
+                if self.pps is not None:
+                    self.pps.value = 1
 
-            # increment 64 bit timestamp
-            if self.ts_rel is not None:
-                t = ((self.ts_rel_ns << 16) + self.ts_rel_fns) + ((self.period_ns << 16) + self.period_fns)
+            if self.ts_tod is not None:
+                self.ts_tod.value = (self.ts_tod_s << 48) | (self.ts_tod_ns << 16) | (self.ts_tod_fns >> 16)
 
-                if self.drift_denom and self.drift_cnt == 0:
-                    t += self.drift_num
+            # increment rel bit timestamp
+            self.ts_rel_fns += (self.period_ns << 32) + self.period_fns
 
-                self.ts_rel_fns = t & 0xffff
-                self.ts_rel_ns = t >> 16
+            if self.drift_denom and self.drift_cnt == 0:
+                self.ts_rel_fns += self.drift_num
 
-                self.ts_rel.value = (self.ts_rel_ns << 16) | self.ts_rel_fns
+            ns_inc = self.ts_rel_fns >> 32
+            self.ts_rel_fns &= 0xffffffff
+
+            self.ts_rel_ns = (self.ts_rel_ns + ns_inc) & 0xffffffffffff
+
+            if self.ts_rel is not None:
+                self.ts_rel.value = (self.ts_rel_ns << 16) | (self.ts_rel_fns >> 16)
 
             if self.drift_denom:
                 if self.drift_cnt > 0:
@@ -273,6 +291,8 @@ def __init__(self, ts_tod=None, ts_rel=None, pps=None, clock=None, *args, **kwar
 
         super().__init__(*args, **kwargs)
 
+        self.ctx = Context(prec=60)
+
         self.ts_tod_s = 0
         self.ts_tod_ns = 0
         self.ts_tod_fns = 0
@@ -296,11 +316,16 @@ def get_ts_tod(self):
 
     def get_ts_tod_96(self):
         ts_s, ts_ns, ts_fns = self.get_ts_tod()
-        return (ts_s << 48) | (ts_ns << 16) | ts_fns
+        return (ts_s << 48) | (ts_ns << 16) | (ts_fns >> 16)
 
     def get_ts_tod_ns(self):
         ts_s, ts_ns, ts_fns = self.get_ts_tod()
-        return ts_s*1e9+ts_ns+ts_fns/2**16
+        ns = Decimal(ts_fns) / Decimal(2**32)
+        ns = self.ctx.add(ns, Decimal(ts_ns))
+        return self.ctx.add(ns, Decimal(ts_s).scaleb(9))
+
+    def get_ts_tod_s(self):
+        return self.get_ts_tod_ns().scaleb(-9, self.ctx)
 
     def get_ts_rel(self):
         return (self.ts_rel_ns, self.ts_rel_fns)
@@ -311,23 +336,26 @@ def get_ts_rel_64(self):
 
     def get_ts_rel_ns(self):
         ts_ns, ts_fns = self.get_ts_rel()
-        return ts_ns + ts_fns/2**16
+        return self.ctx.add(Decimal(ts_fns) / Decimal(2**32), Decimal(ts_ns))
 
     def get_ts_rel_s(self):
-        return self.get_ts_rel()*1e-9
+        return self.get_ts_rel_ns().scaleb(-9, self.ctx)
 
     async def _run(self):
         clock_edge_event = RisingEdge(self.clock)
 
         while True:
             await clock_edge_event
 
-            self.ts_rel_fns, self.ts_rel_ns = math.modf(get_sim_time('ns'))
+            ts_ns, ts_fns = self.ctx.divmod(Decimal(get_sim_time('fs')).scaleb(-6), Decimal(1))
+
+            self.ts_rel_ns = int(ts_ns.to_integral_value()) & 0xffffffffffff
+            self.ts_rel_fns = int((ts_fns * Decimal(2**16)).to_integral_value())
 
-            self.ts_rel_ns = int(self.ts_rel_ns)
-            self.ts_rel_fns = int(self.ts_rel_fns*0x10000)
+            ts_s, ts_ns = self.ctx.divmod(ts_ns, Decimal(1000000000))
 
-            self.ts_tod_s, self.ts_tod_ns = divmod(self.ts_rel_ns, 1000000000)
+            self.ts_tod_s = int(ts_s.scaleb(-9).to_integral_value())
+            self.ts_tod_ns = int(ts_ns.to_integral_value())
             self.ts_tod_fns = self.ts_rel_fns
 
             if self.ts_tod is not None: