diff options
Diffstat (limited to 'kernel/uint63_amd64_63.ml')
| -rw-r--r-- | kernel/uint63_amd64_63.ml | 90 |
1 files changed, 45 insertions, 45 deletions
diff --git a/kernel/uint63_amd64_63.ml b/kernel/uint63_amd64_63.ml index 20b2f58496..5c4028e1c8 100644 --- a/kernel/uint63_amd64_63.ml +++ b/kernel/uint63_amd64_63.ml @@ -82,6 +82,8 @@ let div (x : int) (y : int) = let rem (x : int) (y : int) = if y = 0 then 0 else Int64.to_int (Int64.rem (to_uint64 x) (to_uint64 y)) +let diveucl x y = (div x y, rem x y) + let addmuldiv p x y = l_or (l_sl x p) (l_sr y (uint_size - p)) @@ -94,55 +96,32 @@ let le (x : int) (y : int) = (x lxor 0x4000000000000000) <= (y lxor 0x4000000000000000) [@@ocaml.inline always] -(* A few helper functions on 128 bits *) -let lt128 xh xl yh yl = - lt xh yh || (xh = yh && lt xl yl) - -let le128 xh xl yh yl = - lt xh yh || (xh = yh && le xl yl) - (* division of two numbers by one *) -(* precondition: y <> 0 *) -(* outputs: q % 2^63, r s.t. x = q * y + r, r < y *) +(* precondition: xh < y *) +(* outputs: q, r s.t. x = q * y + r, r < y *) let div21 xh xl y = - let maskh = ref 0 in - let maskl = ref 1 in - let dh = ref 0 in - let dl = ref y in - let cmp = ref true in - (* n = ref 0 *) - (* loop invariant: mask = 2^n, d = mask * y, (2 * d <= x -> cmp), n >= 0 *) - while !dh >= 0 && !cmp do (* dh >= 0 tests that dh highest bit is zero *) - (* We don't use addmuldiv below to avoid checks on 1 *) - dh := (!dh lsl 1) lor (!dl lsr (uint_size - 1)); - dl := !dl lsl 1; - maskh := (!maskh lsl 1) lor (!maskl lsr (uint_size - 1)); - maskl := !maskl lsl 1; - (* incr n *) - cmp := lt128 !dh !dl xh xl; - done; (* mask = 2^n, d = 2^n * y, 2 * d > x *) - let remh = ref xh in - let reml = ref xl in - (* quotienth = ref 0 *) - let quotientl = ref 0 in - (* loop invariant: x = quotient * y + rem, y * 2^(n+1) > r, - mask = floor(2^n), d = mask * y, n >= -1 *) - while !maskh lor !maskl <> 0 do - if le128 !dh !dl !remh !reml then begin (* if rem >= d, add one bit and subtract d *) - (* quotienth := !quotienth lor !maskh *) - quotientl := !quotientl lor !maskl; - remh := if lt !reml !dl then !remh - !dh - 1 else !remh - !dh; - reml := !reml - !dl; - end; - maskl := (!maskl lsr 1) lor (!maskh lsl (uint_size - 1)); - maskh := !maskh lsr 1; - dl := (!dl lsr 1) lor (!dh lsl (uint_size - 1)); - dh := !dh lsr 1; - (* decr n *) + (* nh might temporarily grow as large as 2*y - 1 in the loop body, + so we store it as a 64-bit unsigned integer *) + let nh = ref xh in + let nl = ref xl in + let q = ref 0 in + for _i = 0 to 62 do + (* invariants: 0 <= nh < y, nl = (xl*2^i) % 2^63, + (q*y + nh) * 2^(63-i) + (xl % 2^(63-i)) = (xh%y) * 2^63 + xl *) + nh := Int64.logor (Int64.shift_left !nh 1) (Int64.of_int (!nl lsr 62)); + nl := !nl lsl 1; + q := !q lsl 1; + (* TODO: use "Int64.unsigned_compare !nh y >= 0", + once OCaml 4.08 becomes the minimal required version *) + if Int64.compare !nh 0L < 0 || Int64.compare !nh y >= 0 then + begin q := !q lor 1; nh := Int64.sub !nh y; end done; - !quotientl, !reml + !q, Int64.to_int !nh -let div21 xh xl y = if y = 0 then 0, 0 else div21 xh xl y +let div21 xh xl y = + let xh = to_uint64 xh in + let y = to_uint64 y in + if Int64.compare y xh <= 0 then 0, 0 else div21 xh xl y (* exact multiplication *) (* TODO: check that none of these additions could be a logical or *) @@ -225,3 +204,24 @@ let tail0 x = let is_uint63 t = Obj.is_int t [@@ocaml.inline always] + +(* Arithmetic with explicit carries *) + +(* Analog of Numbers.Abstract.Cyclic.carry *) +type 'a carry = C0 of 'a | C1 of 'a + +let addc x y = + let r = x + y in + if lt r x then C1 r else C0 r + +let addcarryc x y = + let r = x + y + 1 in + if le r x then C1 r else C0 r + +let subc x y = + let r = x - y in + if le y x then C0 r else C1 r + +let subcarryc x y = + let r = x - y - 1 in + if lt y x then C0 r else C1 r |