theory StarForth_Q48_16
  imports "HOL-Library.Word"
begin

(* AND/OR/XOR infix notation moved behind an opt-in bundle at some point after
   2011 -- unbundled by default now to avoid clashing with other uses of the
   same tokens. q48_frac_part below needs it. *)
unbundle bit_operations_syntax

(* =========================================================================
   StarForth_Q48_16 — Q48.16 Fixed-Point Arithmetic (HOL-Word model)

   Mirrors: include/q48_16.h

   This theory uses HOL-Library.Word to model q48_16_t as exactly a 64-bit
   unsigned word — the same wrapping arithmetic as C uint64_t.  This closes
   the "no overflow" domain restriction that the earlier nat model required.

   Format: 64 word with fixed binary point after bit 15.
     Bits  0-15 : fractional part (resolution 1/65536)
     Bits 16-63 : integer part (0 to 2^48 − 1)

   ○ CODE-MUST-MATCH: every operation here matches the corresponding C macro
   or inline function in include/q48_16.h exactly (same shift counts, same
   wrapping behaviour on overflow).
   ======================================================================== *)

type_synonym q48 = "64 word"

(* =========================================================================
   Section 1: Scale constants
   ======================================================================== *)

definition Q48_SCALE :: q48 where "Q48_SCALE = 65536"   \<comment> \<open>2^16 as 64 word\<close>
definition Q48_ONE   :: q48 where "Q48_ONE   = 65536"
definition Q48_HALF  :: q48 where "Q48_HALF  = 32768"

lemma Q48_ONE_eq [simp]: "Q48_ONE = Q48_SCALE"
  by (simp add: Q48_ONE_def Q48_SCALE_def)

lemma Q48_SCALE_nonzero [simp]: "Q48_SCALE \<noteq> 0"
  by (simp add: Q48_SCALE_def)

(* =========================================================================
   Section 2: Conversion operations
   ======================================================================== *)

(* q48_from_u64: integer n → Q48.16 representation (n << 16)
   C: static inline q48_16_t q48_from_u64(uint64_t u) { return u << 16; }
   ○ CODE-MUST-MATCH: shift count = 16, no saturation, wraps on overflow.
   Overflow-free range: n < 2^48.                                            *)
definition q48_from_u64 :: "64 word \<Rightarrow> q48" where
  "q48_from_u64 n = push_bit 16 n"

(* q48_to_u64: Q48.16 value → truncated integer (q >> 16)
   C: static inline uint64_t q48_to_u64(q48_16_t q) { return q >> 16; }
   ○ CODE-MUST-MATCH: logical right shift by 16, fractional bits discarded.  *)
definition q48_to_u64 :: "q48 \<Rightarrow> 64 word" where
  "q48_to_u64 q = drop_bit 16 q"

(* Round-trip: exact when n fits in the 48-bit integer part.
   Proof: drop_bit 16 (push_bit 16 n) = n AND mask 48 = n (when n < 2^48). *)
lemma q48_round_trip:
  assumes "unat n < 2 ^ 48"
  shows "q48_to_u64 (q48_from_u64 n) = n"
proof -
  have "unat n * 65536 < 2 ^ 64" using assms by simp
  hence h: "unat (n * 65536) = unat n * 65536"
    using unat_mult_lem[of n "65536::q48"] by simp
  show ?thesis
    unfolding q48_to_u64_def q48_from_u64_def
    by (simp add: word_unat_eq_iff push_bit_eq_mult drop_bit_eq_div
                  unat_div_distrib h)
qed

lemma q48_from_u64_zero [simp]: "q48_from_u64 0 = 0"
  by (simp add: q48_from_u64_def)

lemma q48_to_u64_zero [simp]: "q48_to_u64 0 = 0"
  by (simp add: q48_to_u64_def)

(* CORRECTED 2026-08-13: the original statement had no upper bound and is
   false as such -- push_bit 16 wraps mod 2^64, so e.g. a=1, b=2^48 satisfies
   unat a \<le> unat b (1 \<le> 2^48) while push_bit 16 a = 65536 and push_bit 16 b
   wraps to 0, breaking the conclusion. Added the same "unat _ < 2^48
   overflow-free range" bound this file already uses everywhere else. *)
lemma q48_from_u64_mono:
  assumes "unat a \<le> unat b" and "unat b < 2 ^ 48"
  shows "unat (q48_from_u64 a) \<le> unat (q48_from_u64 b)"
proof -
  have ha: "unat a < 2 ^ 48" using assms by simp
  have hb64: "unat b * 65536 < 2 ^ 64" using assms(2) by simp
  have ha64: "unat a * 65536 < 2 ^ 64" using ha by simp
  have eb: "unat (b * 65536) = unat b * 65536"
    using unat_mult_lem[of b "65536::q48"] hb64 by simp
  have ea: "unat (a * 65536) = unat a * 65536"
    using unat_mult_lem[of a "65536::q48"] ha64 by simp
  show ?thesis
    unfolding q48_from_u64_def
    by (simp add: push_bit_eq_mult ea eb assms(1))
qed

(* =========================================================================
   Section 3: Arithmetic operations
   ======================================================================== *)

(* q48_add: pointwise addition, wrapping on overflow (matches C + on uint64_t)
   C: static inline q48_16_t q48_add(q48_16_t a, q48_16_t b) { return a+b; } *)
definition q48_add :: "q48 \<Rightarrow> q48 \<Rightarrow> q48" where
  "q48_add a b = a + b"

lemma q48_add_comm: "q48_add a b = q48_add b a"
  by (simp add: q48_add_def add.commute)

lemma q48_add_assoc: "q48_add (q48_add a b) c = q48_add a (q48_add b c)"
  by (simp add: q48_add_def add.assoc)

lemma q48_add_zero_right [simp]: "q48_add a 0 = a"
  by (simp add: q48_add_def)

lemma q48_add_zero_left [simp]: "q48_add 0 a = a"
  by (simp add: q48_add_def)

(* q48_sub: pointwise subtraction, wrapping on underflow (matches C - on
   uint64_t). Added for src/word_source/q48_words.c's Q.- (StarForth_Q48_
   Words.thy).
   C: static inline q48_16_t q48_sub(q48_16_t a, q48_16_t b) { return a-b; } *)
definition q48_sub :: "q48 \<Rightarrow> q48 \<Rightarrow> q48" where
  "q48_sub a b = a - b"

lemma q48_sub_zero_right [simp]: "q48_sub a 0 = a"
  by (simp add: q48_sub_def)

lemma q48_sub_self [simp]: "q48_sub a a = 0"
  by (simp add: q48_sub_def)

lemma q48_add_sub_cancel [simp]: "q48_sub (q48_add a b) b = a"
  by (simp add: q48_sub_def q48_add_def)

(* q48_abs: C compares the raw uint64 bit pattern against 0x8000...0 (the
   sign bit threshold) rather than using a signed comparison operator, but
   this is bit-for-bit the same test as `cell_abs`'s `n <s 0`
   (StarForth_Base.thy) -- both are exactly "is the sign bit set". q48_abs
   is therefore not modelled as a separate definition here: StarForth_Q48_
   Words.thy's Q.ABS uses `cell_abs` directly under the `q48` type synonym. *)

(* q48_mul: (a * b) >> 16, wrapping.
   C: q48_16_t q48_mul(q48_16_t a, q48_16_t b) { return ((__uint128_t)a*b) >> 16; }
   ⚠ HUMAN-REVIEW: The C implementation uses __uint128_t for the intermediate
   product to avoid overflow before shifting.  The HOL model uses 64-word
   multiplication (wrapping), which matches C only when a*b < 2^64 before
   the shift.  Verify that the physics loops stay within this range.         *)
definition q48_mul :: "q48 \<Rightarrow> q48 \<Rightarrow> q48" where
  "q48_mul a b = drop_bit 16 (a * b)"

lemma q48_mul_comm: "q48_mul a b = q48_mul b a"
  by (simp add: q48_mul_def mult.commute)

lemma q48_mul_zero_right [simp]: "q48_mul a 0 = 0"
  by (simp add: q48_mul_def)

lemma q48_mul_zero_left [simp]: "q48_mul 0 a = 0"
  by (simp add: q48_mul_def)

(* q48_mul(a, Q48_ONE) = a when a < 2^48 (overflow-free range) *)
lemma q48_mul_one_right:
  assumes "unat a < 2 ^ 48"
  shows "q48_mul a Q48_ONE = a"
proof -
  have "unat a * 65536 < 2 ^ 64" using assms by simp
  hence h: "unat (a * 65536) = unat a * 65536"
    using unat_mult_lem[of a "65536::q48"] by simp
  show ?thesis
    unfolding q48_mul_def Q48_ONE_def Q48_SCALE_def
    by (simp add: word_unat_eq_iff drop_bit_eq_div unat_div_distrib h)
qed

(* q48_div: (a << 16) / b
   C: return ((__uint128_t)a << 16) / b;
   Same intermediate-precision note as q48_mul applies.                      *)
definition q48_div :: "q48 \<Rightarrow> q48 \<Rightarrow> q48" where
  "q48_div a b = (if b = 0 then 0 else push_bit 16 a div b)"

lemma q48_div_zero_denom [simp]: "q48_div a 0 = 0"
  by (simp add: q48_div_def)

(* CORRECTED 2026-08-13: the original statement had no bound and is false
   as such -- push_bit 16 wraps mod 2^64, so e.g. a=2^48 gives push_bit 16 a
   = 0, so q48_div a Q48_ONE = 0 \<noteq> a. Added the same overflow-free bound
   used throughout this file. *)
lemma q48_div_one:
  assumes "unat a < 2 ^ 48"
  shows "q48_div a Q48_ONE = a"
proof -
  have "unat a * 65536 < 2 ^ 64" using assms by simp
  hence h: "unat (a * 65536) = unat a * 65536"
    using unat_mult_lem[of a "65536::q48"] by simp
  show ?thesis
    unfolding q48_div_def Q48_ONE_def Q48_SCALE_def
    by (simp add: word_unat_eq_iff push_bit_eq_mult unat_div_distrib h)
qed

(* =========================================================================
   Section 4: Accuracy ratio in Q48.16 (used by Loop #4 / Loop #5)
   ======================================================================== *)

(* Prefetch accuracy: hits / total, represented in Q48.16.
   Arguments are natural numbers (counters); result is a 64 word.            *)
definition q48_accuracy :: "nat \<Rightarrow> nat \<Rightarrow> q48" where
  "q48_accuracy hits tot =
     (if tot = 0 then 0
      else word_of_nat ((hits * 65536) div tot))"

lemma q48_accuracy_zero_total [simp]: "q48_accuracy hits 0 = 0"
  by (simp add: q48_accuracy_def)

lemma q48_accuracy_upper_bound:
  assumes "hits \<le> tot"
  shows "unat (q48_accuracy hits tot) \<le> 65536"
proof (cases "tot = 0")
  case True thus ?thesis by simp
next
  case False
  have "(hits * 65536) div tot \<le> (tot * 65536) div tot"
    using assms by (intro div_le_mono) simp
  also have "\<dots> = 65536"
    using False by simp
  finally have "(hits * 65536) div tot \<le> 65536" .
  thus ?thesis
    by (simp add: q48_accuracy_def False unat_of_nat)
qed

(* =========================================================================
   Section 5: Bit-level properties (using HOL-Word bit operations)
   ======================================================================== *)

(* The integer part of a Q48.16 value is its upper 48 bits. *)
definition q48_int_part :: "q48 \<Rightarrow> 64 word" where
  "q48_int_part q = drop_bit 16 q"

(* The fractional part is the lower 16 bits. *)
definition q48_frac_part :: "q48 \<Rightarrow> 64 word" where
  "q48_frac_part q = q AND mask 16"

lemma q48_decompose:
  "push_bit 16 (q48_int_part q) + q48_frac_part q = q"
  unfolding q48_int_part_def q48_frac_part_def
proof -
  have disj: "push_bit 16 (drop_bit 16 q) AND (q AND mask 16) = 0"
    by (rule bit_word_eqI) (auto simp: bit_simps)
  have "push_bit 16 (drop_bit 16 q) + (q AND mask 16)
      = push_bit 16 (drop_bit 16 q) OR (q AND mask 16)"
    by (rule disjunctive_add_eq_or) (rule disj)
  also have "\<dots> = q"
    by (rule bit_word_eqI) (auto simp: bit_simps)
  finally show "push_bit 16 (drop_bit 16 q) + (q AND mask 16) = q" .
qed

end