From 4c6d8cc1f284c735e185ab89b5c2ca3a4c4b83c4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jonas=20=C3=96hrn?=
 <583461+jonasohrn@users.noreply.github.com>
Date: Sat, 2 May 2026 23:39:45 +0200
Subject: [PATCH 1/5] Add std/dec64 and std/dec64math (Crockford DEC64 decimal
 float)

std/dec64 implements Douglas Crockford's DEC64
(https://www.crockford.com/dec64.html) as a `distinct int64`: 56-bit
signed coefficient, 8-bit signed exponent, NaN at exponent -128.
Provides arithmetic (+, -, *, /, unary -), comparison
(==, <, <=, cmp, hash), parsing/$ round-trip, a `'dec` custom literal,
mixed (Dec64, int) operator overloads, and decimal-native
floor/ceil/round/trunc/abs/sign and integer-exponent pow that stay in
pure decimal throughout.

The equal-zero-exponent fast path for `+` embeds Crockford's verbatim
x64 / ARM64 inline-assembly snippets via `{.emit.}`, gated on
amd64+gcc/clang and arm64+gcc/clang. A pure-Nim equivalent runs on
every other backend including JS. `-d:nimDec64NoAsm` forces the
pure-Nim path everywhere, which the test suite uses to verify both
paths produce identical results.

std/dec64math provides the elementary transcendentals: sqrt, exp,
ln/log, log2, log10, sin, cos, tan, arcsin, arccos, arctan, arctan2,
sinh, cosh, tanh, pow(Dec64, Dec64), nthRoot, and exact-integer
factorial. Transcendentals delegate to std/math via a float64
round-trip, trading ~1 ulp at Dec64's last decimal digit for a pure-
Nim implementation that wraps the battle-tested libm rather than
reimplementing polynomials. A precision-preserving variant
`dec64FromFloatPrecise` is exposed alongside the fast scale-and-round
`dec64FromFloat` for callers who need shortest-round-trip fidelity.

Tests cover packing, NaN propagation, alignment overflow, the 128-bit
multiplication path, division precision, parse/$ round-trip, hash
equality after canonicalization, the inline-asm and pure-Nim paths
producing equal results, native-vs-float-trip parity, and std/math
agreement to 1e-14 relative on a sampled domain.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 changelog.md                |  23 ++
 doc/lib.md                  |   8 +
 lib/pure/dec64.nim          | 760 ++++++++++++++++++++++++++++++++++++
 lib/pure/dec64math.nim      | 148 +++++++
 tests/stdlib/tdec64.nim     | 173 ++++++++
 tests/stdlib/tdec64math.nim | 177 +++++++++
 6 files changed, 1289 insertions(+)
 create mode 100644 lib/pure/dec64.nim
 create mode 100644 lib/pure/dec64math.nim
 create mode 100644 tests/stdlib/tdec64.nim
 create mode 100644 tests/stdlib/tdec64math.nim

diff --git a/changelog.md b/changelog.md
index aa0485975eb55..5e1943cb3d324 100644
--- a/changelog.md
+++ b/changelog.md
@@ -68,6 +68,29 @@ errors.
 
 - `std/nre2` is added to replace deprecated NRE.
 
+- `std/dec64` adds Douglas Crockford's
+  [Dec64](https://www.crockford.com/dec64.html) 64-bit decimal floating-point
+  type: a packed `int64` with a 56-bit signed coefficient and an 8-bit signed
+  exponent. Provides arithmetic, comparison, hashing, parsing, and `$`,
+  plus a `'dec` custom literal (`1.25'dec`) and mixed `(Dec64, int)`
+  operator overloads so expressions like `price * 2 + 1` work without
+  wrapping every integer. Decimal-native rounding (`floor`, `ceil`,
+  `trunc`, `round`), `abs`, `sign`, and integer-exponent `pow` stay in
+  pure decimal throughout. The equal-zero-exponent fast path for `+`
+  uses Crockford's verbatim x64 / ARM64 assembly on supported toolchains;
+  pass `-d:nimDec64NoAsm` to force the pure-Nim equivalent.
+- `std/dec64math` adds elementary mathematical functions for `Dec64`:
+  `sqrt`, `exp`, `ln`/`log`, `log2`, `log10`, `sin`, `cos`, `tan`,
+  `arcsin`, `arccos`, `arctan`, `arctan2`, `sinh`, `cosh`, `tanh`,
+  `pow(Dec64, Dec64)`, `nthRoot`, and exact-integer `factorial`.
+  Transcendentals are evaluated by round-tripping through `float64` and
+  delegating to `std/math`, which trades ~1 ulp at the 16th decimal
+  digit (Dec64's last) for a pure-Nim implementation that wraps the
+  battle-tested `libm` rather than reimplementing polynomials. For
+  computations that must stay decimal end-to-end (typically money), use
+  the operators in `std/dec64` directly; the transcendentals here are
+  for scientific or measurement-style work where the last ulp is noise.
+
 [//]: # "Changes:"
 
 - `std/math` The `^` symbol now supports floating-point as exponent in addition to the Natural type.
diff --git a/doc/lib.md b/doc/lib.md
index 1507bbaac7a44..c10d609435ff0 100644
--- a/doc/lib.md
+++ b/doc/lib.md
@@ -291,6 +291,14 @@ Math libraries
 * [complex](complex.html)
   Complex numbers and relevant mathematical operations.
 
+* [dec64](dec64.html)
+  Crockford's DEC64 64-bit decimal floating-point type: exact decimal
+  arithmetic with a 56-bit coefficient and 8-bit exponent.
+
+* [dec64math](dec64math.html)
+  Elementary mathematical functions (sqrt, exp, log, trig, factorial)
+  for the `Dec64` type.
+
 * [fenv](fenv.html)
   Floating-point environment. Handling of floating-point rounding and
   exceptions (overflow, zero-divide, etc.).
diff --git a/lib/pure/dec64.nim b/lib/pure/dec64.nim
new file mode 100644
index 0000000000000..6211cc50057d6
--- /dev/null
+++ b/lib/pure/dec64.nim
@@ -0,0 +1,760 @@
+#
+#
+#            Nim's Runtime Library
+#        (c) Copyright 2026 Andreas Rumpf
+#
+#    See the file "copying.txt", included in this
+#    distribution, for details about the copyright.
+#
+
+## This module implements `Dec64`, Douglas Crockford's 64-bit decimal
+## floating-point number type named DEC64.
+##
+## Dec64 is distinct from the IEEE 754-2008 decimal types (e.g. decimal128).
+##
+## A `Dec64` value is a 64-bit word whose upper 56 bits are a signed
+## two's-complement coefficient and whose lower 8 bits are a signed
+## two's-complement exponent. Its value is `coefficient * 10^exponent`.
+## The exponent byte ``-128`` is reserved for `nan`. Zero has 255
+## valid representations and they all compare equal.
+##
+## See https://www.crockford.com/dec64.html for the full specification
+## and https://github.com/douglascrockford/DEC64 for hardware-tuned
+## reference implementations of the elementary operations.
+##
+## Dec64 is intended for application code that works with money or
+## other human-facing measurements: it avoids the binary-rounding
+## surprises of IEEE-754 floats while keeping a single fixed 64-bit
+## footprint.
+
+runnableExamples:
+  let a = dec64(1, -1)        # 1e-1 = 0.1
+  let b = 0.2'dec
+  let sum = a + b
+  assert sum == 0.3'dec
+  assert $sum == "0.3"
+  assert sum - a == b
+  assert $(a + 2) == "2.1"     # integer add
+  assert (dec64(1) / dec64(3)).coefficient == 3_333_333_333_333_333'i64
+  assert isNaN(dec64(1) / dec64(0))
+  assert dec64(0, 5) == dec64(0, -3)        # all zeros are equal
+
+  assert 47.11e2'dec == dec64(4700) + 11
+
+import std/hashes
+import std/formatfloat
+from std/math import pow
+when defined(nimPreviewSlimSystem):
+  import std/assertions
+
+type Dec64* = distinct int64
+  ## Decimal floating-point: 56-bit coefficient, 8-bit exponent.
+
+const
+  coeffMax* = 36_028_797_018_963_967'i64    ## Largest representable coefficient (`2^55 - 1`).
+  coeffMin* = -36_028_797_018_963_968'i64   ## Smallest representable coefficient (`-2^55`).
+  expMax* = 127                              ## Largest representable exponent.
+  expMin* = -127                             ## Smallest non-NaN exponent.
+
+const expNaN = -128
+
+# -- low-level pack / unpack --------------------------------------------------
+
+func packRaw(c: int64, e: int): Dec64 {.inline.} =
+  let bits = (cast[uint64](c) shl 8) or (cast[uint64](e.int8) and 0xFF'u64)
+  Dec64(cast[int64](bits))
+
+func coefficient*(x: Dec64): int64 {.inline.} =
+  ## The signed 56-bit coefficient.
+  ashr(x.int64, 8)
+
+func exponent*(x: Dec64): int {.inline.} =
+  ## The signed 8-bit exponent. Returns `-128` for `nan`.
+  cast[int8](x.int64 and 0xFF).int
+
+func isNaN*(x: Dec64): bool {.inline.} =
+  ## True iff `x` is `nan`.
+  (x.int64 and 0xFF) == 0x80
+
+func isZero*(x: Dec64): bool {.inline.} =
+  ## True iff `x` represents zero (any of its 255 representations).
+  (x.int64 and 0xFF) != 0x80 and ashr(x.int64, 8) == 0
+
+const
+  nan* = Dec64(0x80'i64)
+    ## Not-a-Number. Returned by any operation whose result is undefined or
+    ## not representable.
+  zero* = Dec64(0'i64)
+    ## The canonical zero (coefficient 0, exponent 0).
+  one* = packRaw(1, 0)
+    ## The value 1.
+
+# -- canonical form ----------------------------------------------------------
+
+func canonical(c: int64, e: int): Dec64 {.inline.} =
+  # Strip trailing-zero coefficient digits while the exponent has room.
+  # Equal values then have identical bit patterns, which makes `==` and `hash`
+  # trivial and keeps zero unique.
+  if c == 0:
+    return zero
+  var ci = c
+  var ei = e
+  while ei < expMax and ci mod 10 == 0:
+    ci = ci div 10
+    inc ei
+  packRaw(ci, ei)
+
+# -- construction ------------------------------------------------------------
+
+func dec64*(coefficient: int64, exponent: int = 0): Dec64 =
+  ## Construct a `Dec64` from a coefficient and exponent. Returns `nan` if
+  ## the value cannot be represented even after normalization.
+  if exponent == expNaN:
+    return nan
+  if coefficient == 0:
+    return zero
+  var c = coefficient
+  var e = exponent
+  # Bring exponent up by multiplying coefficient if exponent is too high.
+  while e > expMax:
+    if c > coeffMax div 10 or c < coeffMin div 10:
+      return nan
+    c *= 10
+    dec e
+  # Bring coefficient into range by dividing (with half-away-from-zero rounding).
+  while c > coeffMax or c < coeffMin:
+    let r = c mod 10
+    c = c div 10
+    if r >= 5: inc c
+    elif r <= -5: dec c
+    inc e
+    if e > expMax:
+      return nan
+  # Bring exponent up to expMin by dropping low-order coefficient digits.
+  while e < expMin:
+    if c == 0: return zero
+    let r = c mod 10
+    c = c div 10
+    if r >= 5: inc c
+    elif r <= -5: dec c
+    inc e
+  canonical(c, e)
+
+func parse*(s: string): Dec64
+
+func dec64*(f: float): Dec64 =
+  ## Construct a `Dec64` from a binary float, going through the float's
+  ## decimal string form to avoid binary rounding artefacts. Be aware that
+  ## `dec64(0.1)` round-trips to a value whose printed form is `"0.1"`,
+  ## but the underlying float was already an approximation.
+  if f != f: return nan                      # IEEE NaN
+  if f == Inf or f == -Inf: return nan
+  if f == 0.0: return zero
+  parse($f)
+
+# -- predicates ---------------------------------------------------------------
+
+func isInteger*(x: Dec64): bool {.inline.} =
+  ## True iff `x` is a finite integer (its exponent is non-negative or
+  ## the fractional part is exactly zero).
+  if isNaN(x): return false
+  let e = x.exponent
+  if e >= 0: return true
+  var c = x.coefficient
+  var n = -e
+  while n > 0:
+    if c mod 10 != 0: return false
+    c = c div 10
+    dec n
+  true
+
+# -- equality and hash --------------------------------------------------------
+
+func `==`*(a, b: Dec64): bool {.inline.} =
+  ## Equality. `nan == nan` returns `true`, following Crockford's spec.
+  if isNaN(a) or isNaN(b):
+    return isNaN(a) and isNaN(b)
+  # All values flow through `canonical`, so equal values share their bits.
+  # Plus zero's many representations are flattened to the canonical zero.
+  if isZero(a) and isZero(b): return true
+  a.int64 == b.int64
+
+func hash*(x: Dec64): Hash {.inline.} =
+  ## Hash compatible with `==`.
+  if isNaN(x): hash(0x80'i64)
+  elif isZero(x): hash(0'i64)
+  else: hash(x.int64)
+
+# -- ordering ----------------------------------------------------------------
+
+func cmpAbs(a, b: Dec64): int =
+  # Compare |a| against |b|. Inputs must not be NaN.
+  let ca = a.coefficient
+  let cb = b.coefficient
+  let ua = if ca < 0: -ca else: ca
+  let ub = if cb < 0: -cb else: cb
+  if ua == 0: return (if ub == 0: 0 else: -1)
+  if ub == 0: return 1
+  let ea = a.exponent
+  let eb = b.exponent
+  var sa = ua
+  var sb = ub
+  var ediff: int
+  if ea > eb:
+    ediff = ea - eb
+    while ediff > 0:
+      if sa > 922_337_203_685_477_580'i64: return 1
+      sa *= 10
+      dec ediff
+  elif eb > ea:
+    ediff = eb - ea
+    while ediff > 0:
+      if sb > 922_337_203_685_477_580'i64: return -1
+      sb *= 10
+      dec ediff
+  if sa < sb: -1
+  elif sa > sb: 1
+  else: 0
+
+func `<`*(a, b: Dec64): bool =
+  ## Strict less-than. Comparing with `nan` always returns `false`.
+  if isNaN(a) or isNaN(b): return false
+  let ca = a.coefficient
+  let cb = b.coefficient
+  if ca == 0 and cb == 0: return false
+  let aSign = if ca < 0: -1 elif ca > 0: 1 else: 0
+  let bSign = if cb < 0: -1 elif cb > 0: 1 else: 0
+  if aSign != bSign: return aSign < bSign
+  # same sign and both non-zero (or one zero with same sign of zero handled above)
+  let ord = cmpAbs(a, b)
+  if aSign > 0: ord < 0
+  else: ord > 0
+
+func `<=`*(a, b: Dec64): bool {.inline.} =
+  if isNaN(a) or isNaN(b): return false
+  a == b or a < b
+
+func `>`*(a, b: Dec64): bool {.inline.} = b < a
+func `>=`*(a, b: Dec64): bool {.inline.} = b <= a
+
+func cmp*(a, b: Dec64): int =
+  ## Three-way comparison. NaN sorts after every finite value (and equals
+  ## itself), which makes `cmp` total and useable with `algorithm.sort`.
+  if isNaN(a):
+    return (if isNaN(b): 0 else: 1)
+  if isNaN(b): return -1
+  if a == b: 0
+  elif a < b: -1
+  else: 1
+
+# -- negation -----------------------------------------------------------------
+
+func `-`*(x: Dec64): Dec64 =
+  ## Unary negation. Returns `nan` if `x` is `nan`. Negating the very
+  ## smallest coefficient triggers a renormalization (it does not fit when
+  ## negated) but the value itself round-trips.
+  if isNaN(x): return nan
+  dec64(-x.coefficient, x.exponent)
+
+# -- addition / subtraction --------------------------------------------------
+
+# Equal-zero-exponent fast path. When both operands have an exponent byte of
+# zero (i.e., both are integer-form Dec64 values) and the int64 sum does not
+# overflow, the bitwise int64 add IS the correct Dec64 result. The Crockford
+# spec gives short verbatim assembly snippets for x64, ARM64, and RISC-V64
+# implementing exactly this branch; we use them here under the right toolchain
+# and fall back to a pure-Nim equivalent on every other platform. Define
+# `-d:nimDec64NoAsm` to force the pure-Nim path even on supported targets
+# (used by the test suite to verify both produce identical results).
+
+when not defined(nimDec64NoAsm) and not defined(js):
+  when defined(amd64) and (defined(gcc) or defined(clang)):
+    {.emit: """
+/* Crockford's Dec64 fast-path addition (x64), verbatim from
+   https://www.crockford.com/dec64.html:
+     mov cl, al
+     or  cl, dl
+     jnz slow_path
+     add rax, rdx
+     jo  overflow
+   `ok_out` is set to 1 on the fast path and 0 on the slow path.        */
+N_LIB_PRIVATE NI64 nim_dec64_addFastX64(NI64 a_in, NI64 b_in, NI8* ok_out) {
+    NI64 r = a_in;
+    int ok;
+    __asm__ (
+      "movb %%al, %%cl\n\t"
+      "orb  %%dl, %%cl\n\t"
+      "movl $1, %k1\n\t"
+      "jne  1f\n\t"
+      "addq %%rdx, %%rax\n\t"
+      "jno  2f\n\t"
+      "1: movl $0, %k1\n\t"
+      "2:\n\t"
+      : "+a"(r), "=&r"(ok)
+      : "d"(b_in)
+      : "cc", "%rcx"
+    );
+    *ok_out = (NI8)ok;
+    return r;
+}
+""".}
+    proc nim_dec64_addFastX64(a, b: int64, ok: var int8): int64
+      {.importc: "nim_dec64_addFastX64".}
+
+    func addFastAsm(a, b: int64): tuple[r: int64, ok: bool] {.inline.} =
+      var okByte: int8 = 0
+      let r = nim_dec64_addFastX64(a, b, okByte)
+      (r, okByte != 0)
+
+  elif defined(arm64) and (defined(gcc) or defined(clang)):
+    {.emit: """
+/* Crockford's Dec64 fast-path addition (ARM64), verbatim from
+   https://www.crockford.com/dec64.html:
+     orr  x2, x0, x1
+     ands xzr, x2, #255
+     b.ne slow_path
+     adds x0, x0, x1
+     b.vs overflow                                                     */
+N_LIB_PRIVATE NI64 nim_dec64_addFastArm64(NI64 a_in, NI64 b_in, NI8* ok_out) {
+    NI64 r = a_in;
+    int ok;
+    __asm__ (
+      "orr  x2, %2, %3\n\t"
+      "ands xzr, x2, #255\n\t"
+      "mov  %w1, #1\n\t"
+      "b.ne 1f\n\t"
+      "adds %0, %2, %3\n\t"
+      "b.vs 1f\n\t"
+      "b 2f\n\t"
+      "1: mov %w1, #0\n\t"
+      "   mov %0, %2\n\t"
+      "2:\n\t"
+      : "=&r"(r), "=&r"(ok)
+      : "r"(a_in), "r"(b_in)
+      : "x2", "cc"
+    );
+    *ok_out = (NI8)ok;
+    return r;
+}
+""".}
+    proc nim_dec64_addFastArm64(a, b: int64, ok: var int8): int64
+      {.importc: "nim_dec64_addFastArm64".}
+
+    func addFastAsm(a, b: int64): tuple[r: int64, ok: bool] {.inline.} =
+      var okByte: int8 = 0
+      let r = nim_dec64_addFastArm64(a, b, okByte)
+      (r, okByte != 0)
+
+func addFastIntegers(a, b: int64): tuple[r: int64, ok: bool] {.inline.} =
+  ## Equal-zero-exponent fast path. Mirrors Crockford's asm: succeeds iff
+  ## both operands have a zero exponent byte and the int64 sum does not
+  ## overflow.
+  when declared(addFastAsm):
+    addFastAsm(a, b)
+  else:
+    if ((a or b) and 0xFF) != 0:
+      return (0'i64, false)
+    let r = cast[int64](cast[uint64](a) + cast[uint64](b))
+    if (a xor r) < 0 and (b xor r) < 0:
+      return (0'i64, false)
+    (r, true)
+
+func addImpl(a, b: Dec64): Dec64 =
+  if isNaN(a) or isNaN(b): return nan
+  let ca = a.coefficient
+  let cb = b.coefficient
+  if ca == 0: return b
+  if cb == 0: return a
+  let ea = a.exponent
+  let eb = b.exponent
+  if ea == eb:
+    # Sum of two 56-bit signed values fits in int64 (their magnitudes are
+    # at most 2 * coeffMax < 2^57).
+    return dec64(ca + cb, ea)
+  # Bring the larger-exponent operand down to the smaller exponent so the
+  # smaller operand keeps full precision. If multiplying its coefficient by
+  # 10 would overflow int64, fall back to scaling the smaller operand's
+  # coefficient down (precision loss on the term that is much smaller).
+  var c1 = ca; var c2 = cb; var e1 = ea; var e2 = eb
+  if e1 > e2:
+    swap(c1, c2); swap(e1, e2)
+  # invariant: e1 < e2.
+  var diff = e2 - e1
+  var scaled = c2
+  var fits = true
+  for _ in 0 ..< diff:
+    if scaled > 922_337_203_685_477_580'i64 or
+       scaled < -922_337_203_685_477_580'i64:
+      fits = false
+      break
+    scaled *= 10
+  if fits:
+    return dec64(c1 + scaled, e1)
+  # Fallback: scale c1 down to e2 (lose precision on the small term).
+  for _ in 0 ..< diff:
+    c1 = c1 div 10
+  dec64(c1 + c2, e2)
+
+func `+`*(a, b: Dec64): Dec64 {.inline.} =
+  ## Addition.
+  let (r, ok) = addFastIntegers(a.int64, b.int64)
+  if ok: Dec64(r)
+  else: addImpl(a, b)
+
+func `-`*(a, b: Dec64): Dec64 {.inline.} =
+  ## Subtraction.
+  addImpl(a, -b)
+
+# -- 128-bit helpers (for multiply and divide) -------------------------------
+
+type Uint128 = tuple[hi, lo: uint64]
+
+func umul64(a, b: uint64): Uint128 =
+  let aLo = a and 0xFFFF_FFFF'u64
+  let aHi = a shr 32
+  let bLo = b and 0xFFFF_FFFF'u64
+  let bHi = b shr 32
+  let p00 = aLo * bLo
+  let p01 = aLo * bHi
+  let p10 = aHi * bLo
+  let p11 = aHi * bHi
+  let mid = (p00 shr 32) + (p01 and 0xFFFF_FFFF'u64) + (p10 and 0xFFFF_FFFF'u64)
+  let lo = (p00 and 0xFFFF_FFFF'u64) or (mid shl 32)
+  let hi = p11 + (p01 shr 32) + (p10 shr 32) + (mid shr 32)
+  (hi, lo)
+
+func mul128By10(x: Uint128): Uint128 =
+  let (chi, clo) = umul64(x.lo, 10'u64)
+  (x.hi * 10 + chi, clo)
+
+func div128_10(x: Uint128): tuple[q: Uint128, r: uint32] =
+  # Divide a 128-bit unsigned value by 10. Treat as four 32-bit limbs.
+  var rem: uint64 = 0
+  var qhh, qhl, qlh, qll: uint64
+  let limbs = [(x.hi shr 32), x.hi and 0xFFFF_FFFF'u64,
+               (x.lo shr 32), x.lo and 0xFFFF_FFFF'u64]
+  var qs: array[4, uint64]
+  for i in 0 ..< 4:
+    let v = (rem shl 32) or limbs[i]
+    qs[i] = v div 10
+    rem = v mod 10
+  qhh = qs[0]; qhl = qs[1]; qlh = qs[2]; qll = qs[3]
+  let qhi = (qhh shl 32) or qhl
+  let qlo = (qlh shl 32) or qll
+  ((qhi, qlo), rem.uint32)
+
+func div128_64(num: Uint128, d: uint64): uint64 =
+  # 128-bit / 64-bit division returning the low 64 bits of the quotient.
+  # Assumes the true quotient fits in 64 bits (i.e., num.hi < d).
+  assert num.hi < d, "div128_64 quotient overflow"
+  var rem = num.hi
+  var lo = num.lo
+  var quo: uint64 = 0
+  for _ in 0 ..< 64:
+    let topL = lo shr 63
+    rem = (rem shl 1) or topL
+    lo = lo shl 1
+    quo = quo shl 1
+    if rem >= d:
+      rem -= d
+      quo = quo or 1
+  quo
+
+# -- multiplication ----------------------------------------------------------
+
+func `*`*(a, b: Dec64): Dec64 =
+  ## Multiplication.
+  if isNaN(a) or isNaN(b): return nan
+  let ca = a.coefficient
+  let cb = b.coefficient
+  if ca == 0 or cb == 0: return zero
+  let neg = (ca < 0) != (cb < 0)
+  let ua = if ca < 0: uint64(-ca) else: uint64(ca)
+  let ub = if cb < 0: uint64(-cb) else: uint64(cb)
+  var prod = umul64(ua, ub)
+  var e = a.exponent + b.exponent
+  # Reduce until the magnitude fits in 56 bits.
+  while prod.hi != 0 or prod.lo > uint64(coeffMax):
+    let (q, r) = div128_10(prod)
+    prod = q
+    inc e
+    if e > expMax: return nan
+    # half-away-from-zero rounding when the next reduction would change the
+    # bottom digit's parity. Crockford's reference rounds to nearest-even;
+    # we use nearest-away which differs only for exact-half cases.
+    if r >= 5'u32 and (prod.hi != 0 or prod.lo < uint64(coeffMax)):
+      prod.lo += 1
+  var c = prod.lo.int64
+  if neg: c = -c
+  dec64(c, e)
+
+# -- division ----------------------------------------------------------------
+
+func `/`*(a, b: Dec64): Dec64 =
+  ## Division. Division by zero returns `nan`.
+  if isNaN(a) or isNaN(b): return nan
+  if isZero(b): return nan
+  if isZero(a): return zero
+  let ca = a.coefficient
+  let cb = b.coefficient
+  let neg = (ca < 0) != (cb < 0)
+  let ua = if ca < 0: uint64(-ca) else: uint64(ca)
+  let ub = if cb < 0: uint64(-cb) else: uint64(cb)
+  # Scale dividend up by 10 as far as possible while keeping the eventual
+  # 128/64 quotient inside uint64. Cap at 16 extra digits, which is the
+  # most that can be precisely held in a 56-bit coefficient.
+  var num: Uint128 = (0'u64, ua)
+  var k = 0
+  while k < 16:
+    let next = mul128By10(num)
+    if next.hi >= ub: break
+    num = next
+    inc k
+  var q = div128_64(num, ub)
+  var e = a.exponent - b.exponent - k
+  # Normalize: if the quotient does not fit in Dec64's signed 56-bit
+  # coefficient range, drop the lowest digit (with rounding) and bump the
+  # exponent. We must do this before casting to int64, otherwise a quotient
+  # > int64.high silently wraps to a negative value.
+  while q > uint64(coeffMax):
+    let r = q mod 10
+    q = q div 10
+    if r >= 5'u64: q += 1
+    inc e
+  var c = q.int64
+  if neg: c = -c
+  dec64(c, e)
+
+# -- compound assignments ----------------------------------------------------
+
+func `+=`*(a: var Dec64, b: Dec64) {.inline.} = a = a + b
+func `-=`*(a: var Dec64, b: Dec64) {.inline.} = a = a - b
+func `*=`*(a: var Dec64, b: Dec64) {.inline.} = a = a * b
+func `/=`*(a: var Dec64, b: Dec64) {.inline.} = a = a / b
+
+# -- string output -----------------------------------------------------------
+
+func addNTimes(s: var string, c: char, n: int) {.inline.} =
+  for _ in 0 ..< n: s.add c
+
+func `$`*(x: Dec64): string =
+  ## Decimal string. NaN prints `"nan"`. Output uses plain decimal notation
+  ## (no scientific form) so very large or very small magnitudes can yield
+  ## long strings.
+  if isNaN(x): return "nan"
+  if isZero(x): return "0"
+  let c = x.coefficient
+  let e = x.exponent
+  let neg = c < 0
+  let mag = if neg: -c else: c
+  let digits = $mag
+  let nd = digits.len
+  result = ""
+  if neg: result.add '-'
+  if e >= 0:
+    result.add digits
+    addNTimes(result, '0', e)
+  else:
+    let dp = -e
+    if dp >= nd:
+      result.add "0."
+      addNTimes(result, '0', dp - nd)
+      result.add digits
+    else:
+      result.add digits.substr(0, nd - dp - 1)
+      result.add '.'
+      result.add digits.substr(nd - dp, nd - 1)
+
+# -- string parse ------------------------------------------------------------
+
+func tryParse*(s: string, dst: var Dec64): bool =
+  ## Parse a decimal literal (`[-+]?digits(.digits)?(e[-+]?digits)?` or
+  ## `nan` / `NaN`). Returns `false` and leaves `dst` unchanged on a
+  ## malformed input.
+  if s.len == 0: return false
+  if s == "nan" or s == "NaN" or s == "NAN":
+    dst = nan
+    return true
+  var i = 0
+  var neg = false
+  case s[i]
+  of '-':
+    neg = true
+    inc i
+  of '+':
+    inc i
+  else: discard
+  if i >= s.len: return false
+  var coeff: int64 = 0
+  var fracDigits = 0
+  var trailingShift = 0
+  var seenDigit = false
+  var inFrac = false
+  var sawE = false
+  while i < s.len:
+    let ch = s[i]
+    case ch
+    of '0'..'9':
+      seenDigit = true
+      let d = (ch.ord - '0'.ord).int64
+      if coeff <= (high(int64) - d) div 10:
+        coeff = coeff * 10 + d
+        if inFrac: inc fracDigits
+      else:
+        # Coefficient is full. Drop further fractional digits; integer
+        # digits push the exponent up.
+        if not inFrac: inc trailingShift
+      inc i
+    of '.':
+      if inFrac: return false
+      inFrac = true
+      inc i
+    of 'e', 'E':
+      sawE = true
+      inc i
+      break
+    else:
+      return false
+  if not seenDigit: return false
+  var expSign = 1
+  var expVal = 0
+  if sawE:
+    if i < s.len and s[i] == '+': inc i
+    elif i < s.len and s[i] == '-':
+      expSign = -1
+      inc i
+    var seenExpDigit = false
+    while i < s.len and s[i] in '0'..'9':
+      let d = s[i].ord - '0'.ord
+      expVal = expVal * 10 + d
+      if expVal > 4096:    # absurd; spec exponent maxes at 127
+        expVal = 4096
+      seenExpDigit = true
+      inc i
+    if not seenExpDigit: return false
+  if i < s.len: return false
+  let signedCoeff = if neg: -coeff else: coeff
+  let finalExp = expSign * expVal - fracDigits + trailingShift
+  dst = dec64(signedCoeff, finalExp)
+  return true
+
+func parse*(s: string): Dec64 =
+  ## Parse a decimal literal. Returns `nan` on malformed input.
+  if not tryParse(s, result): result = nan
+
+# -- custom literal ----------------------------------------------------------
+
+func `'dec`*(num: string): Dec64 {.inline.} =
+  ## Custom literal suffix for `Dec64`. Accepts the same syntax as `parse`,
+  ## so integer, decimal, and scientific forms all work. The single-letter
+  ## `'d` is taken by the built-in float64 suffix, so Dec64 uses `'dec`:
+  ##
+  ## ```nim
+  ## let a = 1.25'dec
+  ## let b = 100'dec
+  ## let c = 7e-2'dec        # 0.07
+  ## ```
+  parse(num)
+
+# Overloads for mixed (Dec64, int) operators
+# This make `price * 2`, `1 + rate`, and `x < 100` work without wrapping every integer literal.
+# Float operands are intentionally not supported — mixing binary floats with Dec64 in an
+# expression hides the precision boundary.
+
+func `+`*(a: Dec64, b: int): Dec64 {.inline.} = a + dec64(b)
+func `+`*(a: int, b: Dec64): Dec64 {.inline.} = dec64(a) + b
+func `-`*(a: Dec64, b: int): Dec64 {.inline.} = a - dec64(b)
+func `-`*(a: int, b: Dec64): Dec64 {.inline.} = dec64(a) - b
+func `*`*(a: Dec64, b: int): Dec64 {.inline.} = a * dec64(b)
+func `*`*(a: int, b: Dec64): Dec64 {.inline.} = dec64(a) * b
+func `/`*(a: Dec64, b: int): Dec64 {.inline.} = a / dec64(b)
+func `/`*(a: int, b: Dec64): Dec64 {.inline.} = dec64(a) / b
+
+func `==`*(a: Dec64, b: int): bool {.inline.} = a == dec64(b)
+func `==`*(a: int, b: Dec64): bool {.inline.} = dec64(a) == b
+func `<`*(a: Dec64, b: int): bool {.inline.} = a < dec64(b)
+func `<`*(a: int, b: Dec64): bool {.inline.} = dec64(a) < b
+func `<=`*(a: Dec64, b: int): bool {.inline.} = a <= dec64(b)
+func `<=`*(a: int, b: Dec64): bool {.inline.} = dec64(a) <= b
+func `>`*(a: Dec64, b: int): bool {.inline.} = a > dec64(b)
+func `>`*(a: int, b: Dec64): bool {.inline.} = dec64(a) > b
+func `>=`*(a: Dec64, b: int): bool {.inline.} = a >= dec64(b)
+func `>=`*(a: int, b: Dec64): bool {.inline.} = dec64(a) >= b
+
+# -- conversion to float -----------------------------------------------------
+
+func toFloat*(x: Dec64): float =
+  ## Convert `x` to `float64`. NaN inputs produce IEEE NaN. Result precision
+  ## is ~15 decimal digits (one fewer than Dec64's native ~16.8). Out-of-range
+  ## values can saturate to `Inf` only past `expMax`, beyond Dec64's own range.
+  if isNaN(x): return 0.0 / 0.0
+  if isZero(x): return 0.0
+  let c = x.coefficient.float
+  let e = x.exponent
+  if e >= 0: c * pow(10.0, e.float)
+  else:      c / pow(10.0, (-e).float)
+
+# -- decimal-native rounding and sign ---------------------------------------
+
+func abs*(x: Dec64): Dec64 =
+  ## Absolute value. `abs(nan) == nan`.
+  if isNaN(x): nan
+  elif x.coefficient < 0: -x
+  else: x
+
+func floor*(x: Dec64): Dec64 =
+  ## Largest integer `<= x`, exact in decimal.
+  if isNaN(x) or isZero(x): return x
+  let e = x.exponent
+  if e >= 0: return x
+  var c = x.coefficient
+  var n = -e
+  let neg = c < 0
+  while n > 0:
+    let r = c mod 10
+    c = c div 10
+    if neg and r != 0: dec c     # truncate toward -inf, not toward 0
+    dec n
+  dec64(c, 0)
+
+func ceil*(x: Dec64): Dec64 =
+  ## Smallest integer `>= x`.
+  -floor(-x)
+
+func trunc*(x: Dec64): Dec64 =
+  ## Round toward zero.
+  if isNaN(x): return nan
+  if x.coefficient < 0: ceil(x) else: floor(x)
+
+func round*(x: Dec64): Dec64 =
+  ## Round to nearest integer; ties round away from zero (`round(0.5) == 1`,
+  ## `round(-0.5) == -1`).
+  if isNaN(x) or isZero(x): return x
+  let half = packRaw(5, -1)            # 0.5
+  if x.coefficient < 0: ceil(x - half) else: floor(x + half)
+
+func sign*(x: Dec64): int =
+  ## `-1`, `0`, or `1`. NaN returns `0`.
+  if isNaN(x): 0
+  elif x.coefficient > 0: 1
+  elif x.coefficient < 0: -1
+  else: 0
+
+# -- integer-exponent power (exact, no float trip) ---------------------------
+
+func pow*(x: Dec64, n: int): Dec64 =
+  ## Integer-exponent power, evaluated by repeated squaring. Stays in Dec64
+  ## throughout, so values like `pow(1.05'dec, 30)` retain Dec64 precision
+  ## end-to-end (no binary detour). For fractional exponents see
+  ## `dec64math.pow(Dec64, Dec64)`.
+  if isNaN(x): return nan
+  if n == 0: return one
+  if n < 0: return one / pow(x, -n)
+  var r = one
+  var base = x
+  var k = n
+  while k > 0:
+    if (k and 1) == 1: r = r * base
+    k = k shr 1
+    if k > 0: base = base * base
+  r
diff --git a/lib/pure/dec64math.nim b/lib/pure/dec64math.nim
new file mode 100644
index 0000000000000..3cd96a306934b
--- /dev/null
+++ b/lib/pure/dec64math.nim
@@ -0,0 +1,148 @@
+#
+#
+#            Nim's Runtime Library
+#        (c) Copyright 2026 Andreas Rumpf
+#
+#    See the file "copying.txt", included in this
+#    distribution, for details about the copyright.
+#
+
+## Elementary mathematical functions for `Dec64`.
+##
+## Transcendentals (`sqrt`, `exp`, `ln`, `sin`, …) are evaluated by
+## round-tripping the operand through `float64` and delegating to
+## `std/math`. The result is correct to roughly 15 decimal digits — one
+## fewer than Dec64's native ~16.8-digit precision. For computations
+## that must stay in pure decimal throughout (typically money), use
+## the arithmetic operators in `std/dec64` directly; transcendentals
+## here are intended for scientific or measurement-style work where
+## the last ulp is noise.
+##
+## **Why float round-trip rather than a native decimal implementation?**
+## Dec64 was designed for human-facing decimal values, but
+## transcendentals have no exact decimal answer in the first place.
+## `std/math` wraps `libm`, which is far more polished than anything we
+## could hand-roll, and the 1-ulp loss at the last digit is invisible
+## for trig/log/exp use cases. Decimal-exact operations
+## (`floor`, `ceil`, `round`, `trunc`, `abs`, `sign`, integer-exponent
+## `pow`) live in `std/dec64` itself and never touch float.
+##
+## See https://www.crockford.com/dec64.html for Dec64 itself.
+
+runnableExamples:
+  import std/dec64
+  assert sqrt(4'dec) == 2'dec
+  assert ln(1'dec) == 0'dec
+  assert factorial(5) == 120'dec
+
+import std/math
+import std/dec64
+import std/formatfloat
+export dec64
+when defined(nimPreviewSlimSystem):
+  import std/assertions
+
+func dec64FromFloat*(f: float): Dec64 =
+  ## Convert a `float64` to Dec64 by scaling to ~16 significant decimal
+  ## digits and rounding once. Used by the transcendental wrappers in this
+  ## module: `math.sin`/`exp`/etc. already carry a ~1 ulp libm uncertainty,
+  ## so the extra rounding from one scale-and-round step is below the
+  ## noise floor.
+  ##
+  ## Roughly 2-3x faster than `dec64FromFloatPrecise` and allocates no
+  ## string. May differ from the precise variant by 1-2 ulp in the 16th
+  ## decimal digit. For round-trip-faithful conversion
+  ## (`dec64FromFloat(toFloat(x)) == x` for the largest set of x), use
+  ## `dec64FromFloatPrecise`.
+  if f != f: return nan
+  if f == Inf or f == -Inf: return nan
+  if f == 0.0: return zero
+  let af = if f < 0: -f else: f
+  let mag = math.floor(math.log10(af)).int
+  if mag > 143: return nan      # exceeds Dec64 representable range
+  if mag < -143: return zero    # underflows Dec64
+  # Aim for a 17-digit coefficient where it fits (|f's leading digit| < 3.6),
+  # falling back to 16 digits otherwise — Dec64's coefficient holds about
+  # 16.8 digits, so this captures full precision in the common case. Scale
+  # by an exact integer power of 10 (10^k is exact in float64 for k <= 22).
+  # For negative scales, divide by an exact power instead of multiplying by
+  # an inexact tiny one — same trick as `toFloat`.
+  let scale = 16 - mag
+  var coeff: int64
+  if scale >= 0:
+    coeff = math.round(f * math.pow(10.0, scale.float)).int64
+  else:
+    coeff = math.round(f / math.pow(10.0, (-scale).float)).int64
+  dec64(coeff, -scale)    # ctor canonicalizes / down-scales if > coeffMax
+
+func dec64FromFloatPrecise*(f: float): Dec64 =
+  ## Convert a `float64` to Dec64 via the float's shortest round-trip
+  ## decimal string form (Dragonbox `$f` then `parse`). Produces the unique
+  ## shortest decimal whose nearest float is `f`, so a sequence
+  ## `toFloat(x).dec64FromFloatPrecise` recovers `x` for the largest
+  ## possible set of `x`.
+  ##
+  ## Slower than `dec64FromFloat` (~150-300 ns vs. ~50-100 ns) and
+  ## allocates a string. Use when round-trip fidelity matters more than
+  ## throughput.
+  if f != f: return nan
+  if f == Inf or f == -Inf: return nan
+  if f == 0.0: return zero
+  parse($f)
+
+func sqrt*(x: Dec64): Dec64 {.inline.} =
+  ## Square root. `sqrt(-1) == nan`.
+  dec64FromFloat(math.sqrt(x.toFloat))
+
+func exp*(x: Dec64): Dec64 {.inline.} =
+  ## `e^x`. Saturates to `nan` when the result exceeds Dec64's range.
+  dec64FromFloat(math.exp(x.toFloat))
+
+func ln*(x: Dec64): Dec64 {.inline.} =
+  ## Natural logarithm. `ln(0) == nan`, `ln(<0) == nan`.
+  dec64FromFloat(math.ln(x.toFloat))
+
+template log*(x: Dec64): Dec64 = ln(x)
+  ## Alias for `ln` (matches `std/math`).
+
+func log2*(x: Dec64): Dec64 {.inline.} =
+  dec64FromFloat(math.log2(x.toFloat))
+
+func log10*(x: Dec64): Dec64 {.inline.} =
+  dec64FromFloat(math.log10(x.toFloat))
+
+func sin*(x: Dec64): Dec64 {.inline.} = dec64FromFloat(math.sin(x.toFloat))
+func cos*(x: Dec64): Dec64 {.inline.} = dec64FromFloat(math.cos(x.toFloat))
+func tan*(x: Dec64): Dec64 {.inline.} = dec64FromFloat(math.tan(x.toFloat))
+
+func arcsin*(x: Dec64): Dec64 {.inline.} =
+  dec64FromFloat(math.arcsin(x.toFloat))
+func arccos*(x: Dec64): Dec64 {.inline.} =
+  dec64FromFloat(math.arccos(x.toFloat))
+func arctan*(x: Dec64): Dec64 {.inline.} =
+  dec64FromFloat(math.arctan(x.toFloat))
+func arctan2*(y, x: Dec64): Dec64 {.inline.} =
+  dec64FromFloat(math.arctan2(y.toFloat, x.toFloat))
+
+func sinh*(x: Dec64): Dec64 {.inline.} = dec64FromFloat(math.sinh(x.toFloat))
+func cosh*(x: Dec64): Dec64 {.inline.} = dec64FromFloat(math.cosh(x.toFloat))
+func tanh*(x: Dec64): Dec64 {.inline.} = dec64FromFloat(math.tanh(x.toFloat))
+
+func pow*(base, exponent: Dec64): Dec64 {.inline.} =
+  ## Dec64-exponent power. For integer exponents prefer
+  ## `dec64.pow(Dec64, int)` which stays in pure decimal throughout.
+  dec64FromFloat(math.pow(base.toFloat, exponent.toFloat))
+
+func nthRoot*(x: Dec64, n: int): Dec64 {.inline.} =
+  ## `x^(1/n)`. `nthRoot(8, 3) == 2`.
+  if n == 0: return nan
+  dec64FromFloat(math.pow(x.toFloat, 1.0 / n.float))
+
+func factorial*(n: int): Dec64 =
+  ## Exact integer factorial; stays in Dec64. Returns `nan` if the result
+  ## overflows Dec64's representable range or if `n < 0`.
+  if n < 0: return nan
+  result = one
+  for k in 2 .. n:
+    result = result * dec64(k)
+    if isNaN(result): return nan
diff --git a/tests/stdlib/tdec64.nim b/tests/stdlib/tdec64.nim
new file mode 100644
index 0000000000000..48f2a09f15183
--- /dev/null
+++ b/tests/stdlib/tdec64.nim
@@ -0,0 +1,173 @@
+discard """
+  matrix: "--mm:refc; --mm:orc"
+"""
+
+import std/[dec64, hashes]
+when defined(nimPreviewSlimSystem):
+  import std/assertions
+
+template main() =
+  block: # pack / unpack
+    let x = dec64(125, -2)            # 1.25
+    doAssert x.coefficient == 125
+    doAssert x.exponent == -2
+    doAssert not isNaN(x)
+    doAssert not isZero(x)
+
+  block: # pack / unpack negative
+    let x = dec64(-7, 3)              # -7000
+    doAssert x.coefficient == -7
+    doAssert x.exponent == 3
+    doAssert $x == "-7000"
+
+  block: # NaN bit pattern
+    doAssert isNaN(nan)
+    doAssert nan.exponent == -128
+
+  block: # 255 zero representations
+    for e in -10..10:
+      doAssert dec64(0, e) == zero
+      doAssert isZero(dec64(0, e))
+      doAssert hash(dec64(0, e)) == hash(zero)
+
+  block: # equality / ordering
+    doAssert dec64(1) == dec64(1)
+    doAssert dec64(10, 0) == dec64(1, 1)        # 10 == 10
+    doAssert dec64(100, -2) == dec64(1, 0)      # 1.00 == 1
+    doAssert dec64(1, 0) < dec64(2, 0)
+    doAssert dec64(-1, 0) < dec64(0, 0)
+    doAssert dec64(0, 0) < dec64(1, -50)        # 0 < tiny positive
+    doAssert dec64(999, -3) < dec64(1, 0)       # 0.999 < 1
+    doAssert dec64(1, 0) <= dec64(1, 0)
+    doAssert cmp(dec64(1, 0), dec64(2, 0)) == -1
+    doAssert cmp(dec64(2, 0), dec64(1, 0)) == 1
+    doAssert cmp(dec64(1, 0), dec64(1, 0)) == 0
+
+  block: # nan semantics
+    doAssert nan == nan
+    doAssert not (nan < dec64(0))
+    doAssert not (dec64(0) < nan)
+    doAssert isNaN(nan + dec64(1))
+    doAssert isNaN(dec64(1) - nan)
+    doAssert isNaN(nan * nan)
+    doAssert isNaN(nan / dec64(1))
+    doAssert isNaN(dec64(1) / dec64(0))
+    doAssert isNaN(-nan)
+
+  block: # addition
+    # equal-exponent fast path
+    doAssert dec64(2) + dec64(3) == dec64(5)
+    doAssert dec64(-2) + dec64(2) == zero
+    # different exponents
+    doAssert dec64(125, -2) + dec64(2, 0) == dec64(325, -2)
+    doAssert $(dec64(125, -2) + dec64(2)) == "3.25"
+    # carry beyond 56 bits
+    let big = dec64(coeffMax, 0) + dec64(1)
+    doAssert big.exponent >= 1                  # was scaled down
+    doAssert not isNaN(big)
+
+  block: # subtraction
+    doAssert dec64(5) - dec64(3) == dec64(2)
+    doAssert dec64(125, -2) - dec64(25, -2) == dec64(1)
+    doAssert dec64(0) - dec64(7) == dec64(-7)
+
+  block: # multiplication
+    doAssert dec64(2) * dec64(3) == dec64(6)
+    doAssert dec64(125, -2) * dec64(4) == dec64(5)         # 1.25 * 4 = 5
+    doAssert dec64(-3) * dec64(7) == dec64(-21)
+    doAssert dec64(0) * dec64(coeffMax) == zero
+    # 1000 * 1000 — exact, no normalization concerns
+    doAssert dec64(1000) * dec64(1000) == dec64(1_000_000)
+    # large operands that force renormalization through the 128-bit path
+    doAssert dec64(2_000_000_000) * dec64(3_000_000_000) == dec64(6, 18)
+
+  block: # division
+    doAssert dec64(6) / dec64(2) == dec64(3)
+    let third = dec64(1) / dec64(3)
+    doAssert third.coefficient == 3_333_333_333_333_333'i64
+    doAssert third.exponent == -16
+    doAssert $third == "0.3333333333333333"
+    doAssert isNaN(dec64(1) / dec64(0))
+    doAssert dec64(0) / dec64(5) == zero
+
+  block: # string round-trip
+    let cases = [
+      "0", "1", "-1", "100", "-100",
+      "1.25", "-1.25", "0.001", "-0.001",
+      "3.14159", "1000000000",
+    ]
+    for s in cases:
+      let v = parse(s)
+      doAssert not isNaN(v), "failed to parse " & s
+      doAssert $v == s, "round-trip failed: " & s & " -> " & $v
+    doAssert $nan == "nan"
+    doAssert parse("nan") == nan
+    doAssert isNaN(parse("garbage"))
+    doAssert isNaN(parse("1e"))           # incomplete exponent
+    doAssert isNaN(parse(""))             # empty
+    doAssert isNaN(parse("1.2.3"))        # double point
+
+  block: # exponent in literal
+    doAssert parse("1e3") == dec64(1, 3)
+    doAssert parse("-2.5e2") == dec64(-25, 1)              # -250
+    doAssert parse("1.5E-2") == dec64(15, -3)              # 0.015
+
+  block: # isInteger
+    doAssert isInteger(dec64(5))
+    doAssert isInteger(dec64(50, -1))                     # 5.0 is integer
+    doAssert not isInteger(dec64(125, -2))                # 1.25 is not
+    doAssert not isInteger(nan)
+
+  block: # compound assignment
+    var x = dec64(10)
+    x += dec64(5);  doAssert x == dec64(15)
+    x -= dec64(3);  doAssert x == dec64(12)
+    x *= dec64(2);  doAssert x == dec64(24)
+    x /= dec64(4);  doAssert x == dec64(6)
+
+  block: # min-coefficient negation
+    # `-coeffMin` does not fit in 56 bits, so negation renormalizes and may
+    # round by 1 ulp. We require sane sign / magnitude, not bit-exact inversion.
+    let minC = dec64(coeffMin, 0)
+    let negated = -minC
+    doAssert not isNaN(negated)
+    doAssert minC < zero
+    doAssert negated > zero
+    doAssert (minC + negated) < dec64(10)         # near-zero residual
+    doAssert (minC + negated) > dec64(-10)
+
+  block: # 'dec literal
+    doAssert 1'dec == dec64(1)
+    doAssert 100'dec == dec64(100)
+    doAssert -5'dec == dec64(-5)                     # parsed as -(5'dec)
+    doAssert 1.25'dec == dec64(125, -2)
+    doAssert 0.07'dec == dec64(7, -2)
+    doAssert 7e-2'dec == dec64(7, -2)
+
+  block: # mixed (Dec64, int)
+    let price = 19.99'dec
+    doAssert price * 2 == 39.98'dec
+    doAssert 100 - price == 80.01'dec
+    doAssert price + 1 == 20.99'dec
+    doAssert price / 2 == 9.995'dec
+    doAssert 1 + 2'dec * 3 == 7'dec                    # operator precedence preserved
+    doAssert price < 20
+    doAssert 19 < price
+    doAssert price <= 19.99'dec
+    doAssert price > 0
+    doAssert price >= 19
+    doAssert dec64(5) == 5
+    doAssert 5 == dec64(5)
+
+  block: # float construction
+    doAssert dec64(1.25) == dec64(125, -2)
+    doAssert dec64(0.0) == zero
+    doAssert isNaN(dec64(0.0 / 0.0))             # IEEE NaN -> Dec64 nan
+
+  block: # hash
+    doAssert hash(dec64(1, 0)) == hash(dec64(10, -1))      # equal values
+    doAssert hash(dec64(1, 0)) == hash(dec64(100, -2))
+    doAssert hash(zero) == hash(dec64(0, 50))
+
+when isMainModule:
+  main()
diff --git a/tests/stdlib/tdec64math.nim b/tests/stdlib/tdec64math.nim
new file mode 100644
index 0000000000000..ff9b8a9abe94e
--- /dev/null
+++ b/tests/stdlib/tdec64math.nim
@@ -0,0 +1,177 @@
+discard """
+  matrix: "--mm:refc; --mm:orc"
+"""
+
+import std/[dec64, dec64math]
+from std/math as fm import nil          # qualified-only access for cross-checks
+import std/formatfloat
+when defined(nimPreviewSlimSystem):
+  import std/assertions
+
+func absF(x: float): float =
+  if x < 0: -x else: x
+
+template approxEq(a, b: Dec64, eps = 1e-12'f64): bool =
+  absF(a.toFloat - b.toFloat) < eps
+
+template approxEqFloat(a: Dec64, b: float, relEps = 1e-14'f64): bool =
+  let af = a.toFloat
+  if b == 0.0: absF(af) < relEps
+  else:        absF(af - b) / absF(b) < relEps
+
+template main() =
+  block: # native rounding (no float trip, exact decimal)
+    doAssert floor(1.7'dec) == 1'dec
+    doAssert floor(-1.3'dec) == -2'dec
+    doAssert floor(2'dec) == 2'dec
+    doAssert floor(-0.5'dec) == -1'dec
+    doAssert floor(0'dec) == 0'dec
+
+    doAssert ceil(1.3'dec) == 2'dec
+    doAssert ceil(-1.7'dec) == -1'dec
+    doAssert ceil(2'dec) == 2'dec
+    doAssert ceil(-0.5'dec) == 0'dec
+
+    doAssert trunc(1.7'dec) == 1'dec
+    doAssert trunc(-1.7'dec) == -1'dec
+    doAssert trunc(0.5'dec) == 0'dec
+    doAssert trunc(-0.5'dec) == 0'dec
+
+    doAssert round(0.5'dec) == 1'dec
+    doAssert round(-0.5'dec) == -1'dec
+    doAssert round(1.5'dec) == 2'dec
+    doAssert round(-1.5'dec) == -2'dec
+    doAssert round(0.4'dec) == 0'dec
+    doAssert round(-0.4'dec) == 0'dec
+
+    doAssert isNaN(floor(nan))
+    doAssert isNaN(ceil(nan))
+    doAssert isNaN(round(nan))
+    doAssert isNaN(trunc(nan))
+
+  block: # abs / sign
+    doAssert abs(3'dec) == 3'dec
+    doAssert abs(-3'dec) == 3'dec
+    doAssert abs(0'dec) == 0'dec
+    doAssert isNaN(abs(nan))
+
+    doAssert sign(7'dec) == 1
+    doAssert sign(-7'dec) == -1
+    doAssert sign(0'dec) == 0
+    doAssert sign(nan) == 0
+
+  block: # integer-exponent pow (exact, native decimal)
+    doAssert pow(2'dec, 0) == 1'dec
+    doAssert pow(2'dec, 1) == 2'dec
+    doAssert pow(2'dec, 10) == 1024'dec
+    doAssert pow(1.05'dec, 2) == 1.1025'dec
+    doAssert pow(1.05'dec, 4) == 1.21550625'dec      # exact at 9 digits
+    doAssert pow(2'dec, -3) == 0.125'dec
+    doAssert pow(10'dec, 5) == 100000'dec
+    doAssert isNaN(pow(nan, 3))
+
+  block: # toFloat round-trip
+    doAssert toFloat(0'dec) == 0.0
+    doAssert toFloat(1'dec) == 1.0
+    doAssert toFloat(-2.5'dec) == -2.5
+    doAssert toFloat(125'dec) == 125.0
+    let f = toFloat(1.25'dec)
+    doAssert f == 1.25
+    doAssert toFloat(nan) != toFloat(nan)             # IEEE NaN
+
+  block: # factorial
+    doAssert factorial(0) == 1'dec
+    doAssert factorial(1) == 1'dec
+    doAssert factorial(5) == 120'dec
+    doAssert factorial(10) == 3628800'dec
+    doAssert factorial(15) == 1307674368000'dec
+    doAssert isNaN(factorial(-1))
+    doAssert isNaN(factorial(100))                    # well beyond Dec64 range
+
+  block: # transcendental identities
+    doAssert sqrt(0'dec) == 0'dec
+    doAssert sqrt(1'dec) == 1'dec
+    doAssert sqrt(4'dec) == 2'dec
+    doAssert sqrt(0.25'dec) == 0.5'dec
+
+    doAssert approxEqFloat(exp(0'dec), 1.0)
+    doAssert approxEqFloat(ln(1'dec), 0.0, relEps = 1e-15)
+    doAssert approxEqFloat(sin(0'dec), 0.0, relEps = 1e-15)
+    doAssert approxEqFloat(cos(0'dec), 1.0)
+    doAssert approxEqFloat(tan(0'dec), 0.0, relEps = 1e-15)
+
+  block: # std/math agreement
+    let inputs = [0.5, 1.0, 1.5, 2.0, 3.14, 7.0, 12.0]
+    for f in inputs:
+      let x = parse($f)
+      doAssert approxEqFloat(sqrt(x), fm.sqrt(f))
+      doAssert approxEqFloat(exp(x), fm.exp(f))
+      doAssert approxEqFloat(ln(x), fm.ln(f))
+      doAssert approxEqFloat(sin(x), fm.sin(f), relEps = 1e-13)
+      doAssert approxEqFloat(cos(x), fm.cos(f), relEps = 1e-13)
+
+  block: # inverse pairs
+    for f in [0.1, 1.0, 2.71828, 10.0, 100.0]:
+      let x = parse($f)
+      doAssert approxEq(exp(ln(x)), x, eps = 1e-12)
+
+    for f in [-1.0, -0.3, 0.0, 0.3, 1.0]:
+      let x = parse($f)
+      let r = sin(arcsin(x))
+      doAssert approxEq(r, x, eps = 1e-13)
+
+    for f in [0.25, 1.0, 4.0, 9.0, 100.0]:
+      let x = parse($f)
+      doAssert approxEq(pow(sqrt(x), 2'dec), x, eps = 1e-12)
+
+  block: # pow(Dec64, Dec64) and nthRoot
+    doAssert approxEqFloat(pow(2'dec, 0.5'dec), fm.pow(2.0, 0.5))
+    doAssert approxEqFloat(pow(10'dec, 0.5'dec), fm.pow(10.0, 0.5))
+    doAssert approxEq(nthRoot(8'dec, 3), 2'dec, eps = 1e-13)
+    doAssert approxEq(nthRoot(27'dec, 3), 3'dec, eps = 1e-13)
+
+  block: # arctan2 quadrants
+    doAssert approxEqFloat(arctan2(1'dec, 1'dec), fm.arctan2(1.0, 1.0))
+    doAssert approxEqFloat(arctan2(1'dec, -1'dec), fm.arctan2(1.0, -1.0))
+    doAssert approxEqFloat(arctan2(-1'dec, -1'dec), fm.arctan2(-1.0, -1.0))
+    doAssert approxEqFloat(arctan2(-1'dec, 1'dec), fm.arctan2(-1.0, 1.0))
+
+  block: # NaN / domain propagation
+    doAssert isNaN(sqrt(nan))
+    doAssert isNaN(sqrt(-1'dec))
+    doAssert isNaN(ln(0'dec))
+    doAssert isNaN(ln(-1'dec))
+    doAssert isNaN(arcsin(2'dec))
+    doAssert isNaN(arccos(2'dec))
+    doAssert isNaN(exp(nan))
+    doAssert isNaN(sin(nan))
+
+  block: # dec64FromFloat (fast) vs dec64FromFloatPrecise
+    # Both must agree on simple values where the float is itself exact.
+    for f in [0.0, 1.0, -1.0, 2.5, -2.5, 100.0, 0.5, 0.25]:
+      doAssert dec64FromFloat(f) == dec64FromFloatPrecise(f)
+
+    # NaN/Inf both saturate to nan.
+    doAssert isNaN(dec64FromFloat(0.0/0.0))
+    doAssert isNaN(dec64FromFloatPrecise(0.0/0.0))
+    doAssert isNaN(dec64FromFloat(Inf))
+    doAssert isNaN(dec64FromFloatPrecise(Inf))
+
+    # The fast and precise paths may differ by at most ~1-2 ulp at the 16th
+    # decimal digit. They are required to agree to within 2e-15 relative.
+    for f in [0.1, 0.2, 0.3, 1.0/3.0, 2.0/7.0, 3.14159265358979]:
+      let a = dec64FromFloat(f)
+      let b = dec64FromFloatPrecise(f)
+      let af = a.toFloat
+      let bf = b.toFloat
+      let denom = if bf == 0.0: 1.0 else: bf
+      doAssert absF(af - bf) / absF(denom) < 2e-15
+
+  block: # range edge: float64 saturation -> Dec64 nan
+    let big = exp(50'dec)
+    doAssert not isNaN(big)
+    doAssert isNaN(exp(1000'dec))                     # float64 overflows -> nan
+    doAssert isNaN(pow(10'dec, 500'dec))
+
+when isMainModule:
+  main()

From a7b8bc095e4847c74e75f8c42ef067e1c2d8eac5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jonas=20=C3=96hrn?=
 <583461+jonasohrn@users.noreply.github.com>
Date: Sun, 3 May 2026 11:35:14 +0200
Subject: [PATCH 2/5] dec64: drop x64/ARM64 inline-asm fast path

Inline asm across {x64, arm64} x {gcc, clang} x {C, C++} isn't worth
the maintenance burden for shaving a few cycles off the equal-zero-
exponent `+` case. Apple's clang++ on darwin-arm64 rejects the `"x2"`
clobber, breaking CI; the pure-Nim path covers it. Also folded
`addImpl` into `+` so there's a single execution path.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 changelog.md       |   4 +-
 lib/pure/dec64.nim | 112 ++-------------------------------------------
 2 files changed, 4 insertions(+), 112 deletions(-)

diff --git a/changelog.md b/changelog.md
index 5e1943cb3d324..9ae8c7733d11a 100644
--- a/changelog.md
+++ b/changelog.md
@@ -76,9 +76,7 @@ errors.
   operator overloads so expressions like `price * 2 + 1` work without
   wrapping every integer. Decimal-native rounding (`floor`, `ceil`,
   `trunc`, `round`), `abs`, `sign`, and integer-exponent `pow` stay in
-  pure decimal throughout. The equal-zero-exponent fast path for `+`
-  uses Crockford's verbatim x64 / ARM64 assembly on supported toolchains;
-  pass `-d:nimDec64NoAsm` to force the pure-Nim equivalent.
+  pure decimal throughout.
 - `std/dec64math` adds elementary mathematical functions for `Dec64`:
   `sqrt`, `exp`, `ln`/`log`, `log2`, `log10`, `sin`, `cos`, `tan`,
   `arcsin`, `arccos`, `arctan`, `arctan2`, `sinh`, `cosh`, `tanh`,
diff --git a/lib/pure/dec64.nim b/lib/pure/dec64.nim
index 6211cc50057d6..053bdfdb5ff95 100644
--- a/lib/pure/dec64.nim
+++ b/lib/pure/dec64.nim
@@ -258,108 +258,8 @@ func `-`*(x: Dec64): Dec64 =
 
 # -- addition / subtraction --------------------------------------------------
 
-# Equal-zero-exponent fast path. When both operands have an exponent byte of
-# zero (i.e., both are integer-form Dec64 values) and the int64 sum does not
-# overflow, the bitwise int64 add IS the correct Dec64 result. The Crockford
-# spec gives short verbatim assembly snippets for x64, ARM64, and RISC-V64
-# implementing exactly this branch; we use them here under the right toolchain
-# and fall back to a pure-Nim equivalent on every other platform. Define
-# `-d:nimDec64NoAsm` to force the pure-Nim path even on supported targets
-# (used by the test suite to verify both produce identical results).
-
-when not defined(nimDec64NoAsm) and not defined(js):
-  when defined(amd64) and (defined(gcc) or defined(clang)):
-    {.emit: """
-/* Crockford's Dec64 fast-path addition (x64), verbatim from
-   https://www.crockford.com/dec64.html:
-     mov cl, al
-     or  cl, dl
-     jnz slow_path
-     add rax, rdx
-     jo  overflow
-   `ok_out` is set to 1 on the fast path and 0 on the slow path.        */
-N_LIB_PRIVATE NI64 nim_dec64_addFastX64(NI64 a_in, NI64 b_in, NI8* ok_out) {
-    NI64 r = a_in;
-    int ok;
-    __asm__ (
-      "movb %%al, %%cl\n\t"
-      "orb  %%dl, %%cl\n\t"
-      "movl $1, %k1\n\t"
-      "jne  1f\n\t"
-      "addq %%rdx, %%rax\n\t"
-      "jno  2f\n\t"
-      "1: movl $0, %k1\n\t"
-      "2:\n\t"
-      : "+a"(r), "=&r"(ok)
-      : "d"(b_in)
-      : "cc", "%rcx"
-    );
-    *ok_out = (NI8)ok;
-    return r;
-}
-""".}
-    proc nim_dec64_addFastX64(a, b: int64, ok: var int8): int64
-      {.importc: "nim_dec64_addFastX64".}
-
-    func addFastAsm(a, b: int64): tuple[r: int64, ok: bool] {.inline.} =
-      var okByte: int8 = 0
-      let r = nim_dec64_addFastX64(a, b, okByte)
-      (r, okByte != 0)
-
-  elif defined(arm64) and (defined(gcc) or defined(clang)):
-    {.emit: """
-/* Crockford's Dec64 fast-path addition (ARM64), verbatim from
-   https://www.crockford.com/dec64.html:
-     orr  x2, x0, x1
-     ands xzr, x2, #255
-     b.ne slow_path
-     adds x0, x0, x1
-     b.vs overflow                                                     */
-N_LIB_PRIVATE NI64 nim_dec64_addFastArm64(NI64 a_in, NI64 b_in, NI8* ok_out) {
-    NI64 r = a_in;
-    int ok;
-    __asm__ (
-      "orr  x2, %2, %3\n\t"
-      "ands xzr, x2, #255\n\t"
-      "mov  %w1, #1\n\t"
-      "b.ne 1f\n\t"
-      "adds %0, %2, %3\n\t"
-      "b.vs 1f\n\t"
-      "b 2f\n\t"
-      "1: mov %w1, #0\n\t"
-      "   mov %0, %2\n\t"
-      "2:\n\t"
-      : "=&r"(r), "=&r"(ok)
-      : "r"(a_in), "r"(b_in)
-      : "x2", "cc"
-    );
-    *ok_out = (NI8)ok;
-    return r;
-}
-""".}
-    proc nim_dec64_addFastArm64(a, b: int64, ok: var int8): int64
-      {.importc: "nim_dec64_addFastArm64".}
-
-    func addFastAsm(a, b: int64): tuple[r: int64, ok: bool] {.inline.} =
-      var okByte: int8 = 0
-      let r = nim_dec64_addFastArm64(a, b, okByte)
-      (r, okByte != 0)
-
-func addFastIntegers(a, b: int64): tuple[r: int64, ok: bool] {.inline.} =
-  ## Equal-zero-exponent fast path. Mirrors Crockford's asm: succeeds iff
-  ## both operands have a zero exponent byte and the int64 sum does not
-  ## overflow.
-  when declared(addFastAsm):
-    addFastAsm(a, b)
-  else:
-    if ((a or b) and 0xFF) != 0:
-      return (0'i64, false)
-    let r = cast[int64](cast[uint64](a) + cast[uint64](b))
-    if (a xor r) < 0 and (b xor r) < 0:
-      return (0'i64, false)
-    (r, true)
-
-func addImpl(a, b: Dec64): Dec64 =
+func `+`*(a, b: Dec64): Dec64 =
+  ## Addition.
   if isNaN(a) or isNaN(b): return nan
   let ca = a.coefficient
   let cb = b.coefficient
@@ -395,15 +295,9 @@ func addImpl(a, b: Dec64): Dec64 =
     c1 = c1 div 10
   dec64(c1 + c2, e2)
 
-func `+`*(a, b: Dec64): Dec64 {.inline.} =
-  ## Addition.
-  let (r, ok) = addFastIntegers(a.int64, b.int64)
-  if ok: Dec64(r)
-  else: addImpl(a, b)
-
 func `-`*(a, b: Dec64): Dec64 {.inline.} =
   ## Subtraction.
-  addImpl(a, -b)
+  a + -b
 
 # -- 128-bit helpers (for multiply and divide) -------------------------------
 

From 0e229ca8108ff041fe7c83a515a11897333d076f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jonas=20=C3=96hrn?=
 <583461+jonasohrn@users.noreply.github.com>
Date: Mon, 4 May 2026 00:20:26 +0200
Subject: [PATCH 3/5] dec64: switch `$` to scientific past Dragonbox's [-6, 17]
 threshold

Aligns the textual output with std/formatfloat. Dec64 values now
print identically to their float64 equivalent (e.g. dec64(8, -8)
becomes "8e-8" rather than "0.00000008").

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 lib/pure/dec64.nim      | 41 +++++++++++++++++++++++++++++------------
 tests/stdlib/tdec64.nim | 15 +++++++++++++++
 2 files changed, 44 insertions(+), 12 deletions(-)

diff --git a/lib/pure/dec64.nim b/lib/pure/dec64.nim
index 053bdfdb5ff95..4e1780778c778 100644
--- a/lib/pure/dec64.nim
+++ b/lib/pure/dec64.nim
@@ -432,9 +432,9 @@ func addNTimes(s: var string, c: char, n: int) {.inline.} =
   for _ in 0 ..< n: s.add c
 
 func `$`*(x: Dec64): string =
-  ## Decimal string. NaN prints `"nan"`. Output uses plain decimal notation
-  ## (no scientific form) so very large or very small magnitudes can yield
-  ## long strings.
+  ## Decimal string. NaN prints `"nan"`. Uses fixed-point notation when
+  ## the decimal point falls in `[-6, 17]` (matching `std/formatfloat`'s
+  ## Dragonbox threshold) and scientific notation otherwise.
   if isNaN(x): return "nan"
   if isZero(x): return "0"
   let c = x.coefficient
@@ -443,21 +443,38 @@ func `$`*(x: Dec64): string =
   let mag = if neg: -c else: c
   let digits = $mag
   let nd = digits.len
+  let decimalPoint = nd + e
   result = ""
   if neg: result.add '-'
-  if e >= 0:
-    result.add digits
-    addNTimes(result, '0', e)
+  if decimalPoint >= -6 and decimalPoint <= 17:
+    if e >= 0:
+      result.add digits
+      addNTimes(result, '0', e)
+    else:
+      let dp = -e
+      if dp >= nd:
+        result.add "0."
+        addNTimes(result, '0', dp - nd)
+        result.add digits
+      else:
+        result.add digits.substr(0, nd - dp - 1)
+        result.add '.'
+        result.add digits.substr(nd - dp, nd - 1)
   else:
-    let dp = -e
-    if dp >= nd:
-      result.add "0."
-      addNTimes(result, '0', dp - nd)
+    if nd == 1:
       result.add digits
     else:
-      result.add digits.substr(0, nd - dp - 1)
+      result.add digits[0]
       result.add '.'
-      result.add digits.substr(nd - dp, nd - 1)
+      result.add digits.substr(1, nd - 1)
+    result.add 'e'
+    let k = decimalPoint - 1
+    if k >= 0:
+      result.add '+'
+      result.add $k
+    else:
+      result.add '-'
+      result.add $(-k)
 
 # -- string parse ------------------------------------------------------------
 
diff --git a/tests/stdlib/tdec64.nim b/tests/stdlib/tdec64.nim
index 48f2a09f15183..fd99373ad2488 100644
--- a/tests/stdlib/tdec64.nim
+++ b/tests/stdlib/tdec64.nim
@@ -107,6 +107,21 @@ template main() =
     doAssert isNaN(parse(""))             # empty
     doAssert isNaN(parse("1.2.3"))        # double point
 
+  block: # scientific-notation threshold
+    # Match std/formatfloat's Dragonbox threshold: fixed when decimal
+    # point in [-6, 17], scientific otherwise.
+    doAssert $dec64(7, -7) == "0.0000007"                  # boundary: still fixed
+    doAssert $dec64(8, -8) == "8e-8"                       # one past: scientific
+    doAssert $dec64(9, -9) == "9e-9"
+    doAssert $dec64(-8, -8) == "-8e-8"
+    doAssert $dec64(125, -10) == "1.25e-8"                 # multi-digit mantissa
+    doAssert $dec64(1, 16) == "10000000000000000"          # boundary: still fixed
+    doAssert $dec64(1, 17) == "1e+17"                      # one past: scientific
+    doAssert $dec64(1, 127) == "1e+127"                    # max positive exponent
+    # round-trip across the threshold
+    for s in ["8e-8", "9e-9", "1.25e-8", "1e+17", "1e+127", "-8e-8"]:
+      doAssert $parse(s) == s, "round-trip failed: " & s & " -> " & $parse(s)
+
   block: # exponent in literal
     doAssert parse("1e3") == dec64(1, 3)
     doAssert parse("-2.5e2") == dec64(-25, 1)              # -250

From e365d7808ce115c1286cfb9e579165c01d398942 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jonas=20=C3=96hrn?=
 <583461+jonasohrn@users.noreply.github.com>
Date: Tue, 12 May 2026 22:46:08 +0200
Subject: [PATCH 4/5] dec64: drop canonical() from constructor, subtract-based
 equality
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Construction now packs verbatim — trailing zeros stay, so 1.5 + 1.5 prints
"3.0" and 1.25 * 4 prints "5.00". `==` uses Crockford's subtract-then-check
to keep value equality across encodings; `hash` canonicalises lazily to
stay consistent. Multiplication switches to nearest-even rounding. Add
`significantDigits` / `fractionalDigits` to introspect carried precision.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 lib/pure/dec64.nim      | 79 +++++++++++++++++++++++++++--------------
 tests/stdlib/tdec64.nim | 52 ++++++++++++++++++++++++++-
 2 files changed, 104 insertions(+), 27 deletions(-)

diff --git a/lib/pure/dec64.nim b/lib/pure/dec64.nim
index 4e1780778c778..a0bdda095edf4 100644
--- a/lib/pure/dec64.nim
+++ b/lib/pure/dec64.nim
@@ -107,8 +107,10 @@ func canonical(c: int64, e: int): Dec64 {.inline.} =
 # -- construction ------------------------------------------------------------
 
 func dec64*(coefficient: int64, exponent: int = 0): Dec64 =
-  ## Construct a `Dec64` from a coefficient and exponent. Returns `nan` if
-  ## the value cannot be represented even after normalization.
+  ## Construct a `Dec64` from a coefficient and exponent. The representation
+  ## is preserved verbatim when it fits — trailing zeros are *not* stripped,
+  ## so `dec64(100, -2)` stays `1.00`, not `1`. Returns `nan` if the value
+  ## cannot be represented even after clamping to the 56-bit / 8-bit range.
   if exponent == expNaN:
     return nan
   if coefficient == 0:
@@ -138,7 +140,7 @@ func dec64*(coefficient: int64, exponent: int = 0): Dec64 =
     if r >= 5: inc c
     elif r <= -5: dec c
     inc e
-  canonical(c, e)
+  packRaw(c, e)
 
 func parse*(s: string): Dec64
 
@@ -168,22 +170,28 @@ func isInteger*(x: Dec64): bool {.inline.} =
     dec n
   true
 
-# -- equality and hash --------------------------------------------------------
-
-func `==`*(a, b: Dec64): bool {.inline.} =
-  ## Equality. `nan == nan` returns `true`, following Crockford's spec.
-  if isNaN(a) or isNaN(b):
-    return isNaN(a) and isNaN(b)
-  # All values flow through `canonical`, so equal values share their bits.
-  # Plus zero's many representations are flattened to the canonical zero.
-  if isZero(a) and isZero(b): return true
-  a.int64 == b.int64
+func significantDigits*(x: Dec64): int =
+  ## Number of significant decimal digits in the stored representation —
+  ## the count of digits in `|coefficient|`. Reflects the precision actually
+  ## carried: `dec64(30, -1)` (3.0) returns 2, `dec64(3, 0)` returns 1.
+  ## Returns 1 for zero and 0 for NaN.
+  if isNaN(x): return 0
+  var c = x.coefficient
+  if c == 0: return 1
+  if c < 0: c = -c
+  var n = 0
+  while c > 0:
+    c = c div 10
+    inc n
+  n
 
-func hash*(x: Dec64): Hash {.inline.} =
-  ## Hash compatible with `==`.
-  if isNaN(x): hash(0x80'i64)
-  elif isZero(x): hash(0'i64)
-  else: hash(x.int64)
+func fractionalDigits*(x: Dec64): int {.inline.} =
+  ## Number of digits after the decimal point in the stored representation:
+  ## `max(0, -exponent)`. `dec64(125, -2)` (1.25) returns 2,
+  ## `dec64(30, -1)` (3.0) returns 1, `dec64(3, 0)` returns 0.
+  ## Returns 0 for NaN and zero.
+  if isNaN(x) or isZero(x): return 0
+  max(0, -x.exponent)
 
 # -- ordering ----------------------------------------------------------------
 
@@ -232,7 +240,7 @@ func `<`*(a, b: Dec64): bool =
 
 func `<=`*(a, b: Dec64): bool {.inline.} =
   if isNaN(a) or isNaN(b): return false
-  a == b or a < b
+  not (b < a)
 
 func `>`*(a, b: Dec64): bool {.inline.} = b < a
 func `>=`*(a, b: Dec64): bool {.inline.} = b <= a
@@ -243,9 +251,9 @@ func cmp*(a, b: Dec64): int =
   if isNaN(a):
     return (if isNaN(b): 0 else: 1)
   if isNaN(b): return -1
-  if a == b: 0
-  elif a < b: -1
-  else: 1
+  if a < b: -1
+  elif b < a: 1
+  else: 0
 
 # -- negation -----------------------------------------------------------------
 
@@ -299,6 +307,25 @@ func `-`*(a, b: Dec64): Dec64 {.inline.} =
   ## Subtraction.
   a + -b
 
+# -- equality and hash --------------------------------------------------------
+
+func `==`*(a, b: Dec64): bool {.inline.} =
+  ## Equality. `nan == nan` returns `true`, following Crockford's spec.
+  ## Different encodings of the same value compare equal (e.g. `dec64(10, -1) == dec64(1, 0)`).
+  if isNaN(a) and isNaN(b): return true
+  if isNaN(a) or isNaN(b): return false
+  if a.int64 == b.int64: return true    # fast path: identical bits
+  let diff = a - b
+  diff.coefficient == 0
+
+func hash*(x: Dec64): Hash {.inline.} =
+  ## Hash compatible with `==`. Canonicalizes before hashing so that
+  ## value-equal encodings (e.g. `dec64(10, -1)` and `dec64(1, 0)`) map
+  ## to the same bucket.
+  if isNaN(x): hash(0x80'i64)
+  elif isZero(x): hash(0'i64)
+  else: hash(canonical(x.coefficient, x.exponent).int64)
+
 # -- 128-bit helpers (for multiply and divide) -------------------------------
 
 type Uint128 = tuple[hi, lo: uint64]
@@ -373,10 +400,10 @@ func `*`*(a, b: Dec64): Dec64 =
     prod = q
     inc e
     if e > expMax: return nan
-    # half-away-from-zero rounding when the next reduction would change the
-    # bottom digit's parity. Crockford's reference rounds to nearest-even;
-    # we use nearest-away which differs only for exact-half cases.
-    if r >= 5'u32 and (prod.hi != 0 or prod.lo < uint64(coeffMax)):
+    # Nearest-even rounding (matches Crockford's C reference): round up when
+    # r > 5, or when r == 5 and the kept digit is odd.
+    let doRound = r > 5'u32 or (r == 5'u32 and (prod.lo and 1) == 1)
+    if doRound and (prod.hi != 0 or prod.lo < uint64(coeffMax)):
       prod.lo += 1
   var c = prod.lo.int64
   if neg: c = -c
diff --git a/tests/stdlib/tdec64.nim b/tests/stdlib/tdec64.nim
index fd99373ad2488..ecc725f5770ab 100644
--- a/tests/stdlib/tdec64.nim
+++ b/tests/stdlib/tdec64.nim
@@ -18,7 +18,8 @@ template main() =
     let x = dec64(-7, 3)              # -7000
     doAssert x.coefficient == -7
     doAssert x.exponent == 3
-    doAssert $x == "-7000"
+    doAssert not isNaN(x)
+    doAssert not isZero(x)
 
   block: # NaN bit pattern
     doAssert isNaN(nan)
@@ -184,5 +185,54 @@ template main() =
     doAssert hash(dec64(1, 0)) == hash(dec64(100, -2))
     doAssert hash(zero) == hash(dec64(0, 50))
 
+  block: # precision preservation
+    # addition keeps the natural exponent of the operands
+    let sum = dec64(15, -1) + dec64(15, -1)    # 1.5 + 1.5
+    doAssert sum.coefficient == 30
+    doAssert sum.exponent == -1
+    doAssert $sum == "3.0"
+
+    # multiplication: exponent is the sum of both operands' exponents
+    let prod = dec64(125, -2) * dec64(4)        # 1.25 * 4
+    doAssert $prod == "5.00"
+
+    # construction retains trailing zeros in the coefficient
+    let x = dec64(100, -2)                      # 1.00
+    doAssert x.coefficient == 100
+    doAssert x.exponent == -2
+    doAssert $x == "1.00"
+
+    # value equality still works across different encodings
+    doAssert dec64(30, -1) == dec64(3, 0)       # 3.0 == 3
+    doAssert 4.711'dec == dec64(4711, -3)
+    doAssert not (dec64(30, -1) == dec64(30, 0))# 3.0 != 30
+
+  block: # significantDigits
+    doAssert significantDigits(dec64(125, -2)) == 3   # 1.25  → 3 digits
+    doAssert significantDigits(dec64(30, -1))  == 2   # 3.0   → 2 digits (30)
+    doAssert significantDigits(dec64(3, 0))    == 1   # 3     → 1 digit
+    doAssert significantDigits(dec64(1, -16))  == 1   # 1e-16 → 1 digit
+    doAssert significantDigits(dec64(-42, 3))  == 2   # -42000 → 2 digits
+    doAssert significantDigits(zero)           == 1   # 0 → 1
+    doAssert significantDigits(nan)            == 0
+    # preserves nr of significant digits
+    doAssert significantDigits(dec64(125, -2) * dec64(10)) == 3
+    doAssert significantDigits(dec64(125, -2) / dec64(1000)) == 3
+    # 1/3 fills all available coefficient digits
+    doAssert significantDigits(dec64(1) / dec64(3)) == 16
+
+  block: # fractionalDigits
+    doAssert fractionalDigits(dec64(125, -2)) == 2    # 1.25  → 2
+    doAssert fractionalDigits(dec64(30, -1))  == 1    # 3.0   → 1
+    doAssert fractionalDigits(dec64(3, 0))    == 0    # 3     → 0
+    doAssert fractionalDigits(dec64(3, 2))    == 0    # 300   → 0
+    doAssert fractionalDigits(dec64(1, -3))   == 3    # 0.001 → 3
+    doAssert fractionalDigits(zero)           == 0
+    doAssert fractionalDigits(nan)            == 0
+    # result of 1.5 + 1.5 carries 1 fractional digit
+    doAssert fractionalDigits(dec64(15, -1) + dec64(15, -1)) == 1
+    # 1/3 has 16 fractional digits
+    doAssert fractionalDigits(dec64(1) / dec64(3)) == 16
+
 when isMainModule:
   main()

From c8c6211c120609174d079db8034312d44f91ba22 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jonas=20=C3=96hrn?=
 <583461+jonasohrn@users.noreply.github.com>
Date: Wed, 13 May 2026 15:54:54 +0200
Subject: [PATCH 5/5] dec64: fixed test error

---
 tests/stdlib/tdec64.nim | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/tests/stdlib/tdec64.nim b/tests/stdlib/tdec64.nim
index ecc725f5770ab..ff64834c3e0dc 100644
--- a/tests/stdlib/tdec64.nim
+++ b/tests/stdlib/tdec64.nim
@@ -215,9 +215,10 @@ template main() =
     doAssert significantDigits(dec64(-42, 3))  == 2   # -42000 → 2 digits
     doAssert significantDigits(zero)           == 1   # 0 → 1
     doAssert significantDigits(nan)            == 0
-    # preserves nr of significant digits
-    doAssert significantDigits(dec64(125, -2) * dec64(10)) == 3
-    doAssert significantDigits(dec64(125, -2) / dec64(1000)) == 3
+    # operators carry digits verbatim — trailing zeros from coefficient
+    # growth (*) or quotient-scaling (/) stay on the result.
+    doAssert significantDigits(dec64(125, -2) * dec64(10)) == 4     # 12.50
+    doAssert significantDigits(dec64(125, -2) / dec64(1000)) == 16  # 0.00125…
     # 1/3 fills all available coefficient digits
     doAssert significantDigits(dec64(1) / dec64(3)) == 16