Long multiplication: Difference between revisions
ATS added ahead of →{{header|AutoHotkey}}
(Long multiplication in BASIC-256) |
(ATS added ahead of →{{header|AutoHotkey}}) |
||
Line 998:
<pre>340282366920938463463374607431768211456</pre>
=={{header|ATS}}==
I perform both binary long multiplication (designed for efficiency) and decimal arithmetic.
The example numbers, though fine for testing decimal arithmetic, are ''very'' bad for testing binary multiplication. They never require a carry. So I have added (for binary multiplication) the example 79 228 162 514 264 337 593 543 950 335 squared equals 6 277 101 735 386 680 763 835 789 423 049 210 091 073 826 769 276 946 612 225. The first number is 96 one-bits.
<syntaxhighlight lang="ats">
(* This is Algorithm 4.3.1M in Volume 2 of Knuth, ‘The Art of Computer
Programming’. *)
#include "share/atspre_staload.hats"
#define NIL list_nil ()
#define :: list_cons
(********************** FOR BINARY ARITHMETIC ***********************)
(* We need to choose a radix for the multiplication, small enough that
intermediate results can be represented, but big for efficiency. To
stay within the POSIX types, I choose 2**32 as my radix. Thus
‘digits’ are stored in uint32 and intermediate results are stored
in uint64.
A number is stored as an array of uint32, with the least
significant uint32 first. *)
extern fn
long_multiplication (* Multiply u and v, giving w. *)
{m, n : int}
(m : size_t m,
n : size_t n,
u : &array (uint32, m),
v : &array (uint32, n),
w : &array (uint32?, m + n) >> array (uint32, m + n))
:<!refwrt> void
%{^
#include <stdint.h>
%}
extern castfn i2u32 : int -<> uint32
extern castfn u32_2i : uint32 -<> int
extern castfn i2u64 : int -<> uint64
extern castfn u32u64 : uint32 -<> uint64
extern castfn u64u32 : uint64 -<> uint32
macdef zero32 = i2u32 0
macdef zero64 = i2u64 0
macdef one32 = i2u32 1
macdef ten32 = i2u32 10
macdef mask32 = $extval (uint32, "UINT32_C (0xFFFFFFFF)")
(* The following implementation is precisely the algorithm suggested
by Knuth, although specialized for b=2**32 and for unsigned
integers of precisely 32 bits. *)
implement
long_multiplication {m, n} (m, n, u, v, w) =
let
(* Establish that the arrays have non-negative lengths. *)
prval () = lemma_array_param u
prval () = lemma_array_param v
(* Knuth initializes only part of the w array. However, if we
initialize ALL of w now, then we will not have to deal with
complicated array views later. *)
val () = array_initize_elt<uint32> (w, m + n, zero32)
(* The following function includes proof of termination. *)
fun
jloop {j : nat | j <= n} .<n - j>.
(u : &array (uint32, m),
v : &array (uint32, n),
w : &array (uint32, m + n),
j : size_t j)
:<!refwrt> void =
if j = n then
()
else if v[j] = zero32 then (* This branch is optional. *)
begin
w[j + m] := zero32;
jloop (u, v, w, succ j)
end
else
let
fun
iloop {i : nat | i <= m} .<m - i>.
(u : &array (uint32, m),
v : &array (uint32, n),
w : &array (uint32, m + n),
i : size_t i,
k : uint64) (* carry *)
:<!refwrt> uint32 = (* carry *)
if i = m then
u64u32 k
else
let
val t = (u32u64 u[i] * u32u64 v[j])
+ u32u64 w[i + j] + k
in
(* The mask here is not actually needed, if uint32
really is treated by the C compiler as 32 bits. *)
w[i + j] := (u64u32 t) land mask32;
iloop (u, v, w, succ i, t >> 32)
end
val k = iloop (u, v, w, i2sz 0, zero64)
in
w[j + m] := k;
jloop (u, v, w, succ j)
end
in
jloop (u, v, w, i2sz 0)
end
fn
big_integer_iseqz (* Is a big integer equal to zero? *)
{m : int}
(m : size_t m,
u : &array (uint32, m))
:<!ref> bool =
let
prval () = lemma_array_param u
fun
loop {n : nat | n <= m} .<n>.
(u : &array (uint32, m),
n : size_t n)
:<!ref> bool =
if n = i2sz 0 then
true
else if u[pred n] = zero32 then
loop (u, pred n)
else
false
in
loop (u, m)
end
(* To print the number in decimal, we need division by 10. So here is
‘short division’: Exercise 4.3.1.16 in Volume 2 of Knuth. *)
fn
short_division
{m : int}
(m : size_t m,
u : &array (uint32, m),
v : uint32,
q : &array (uint32?, m) >> array (uint32, m),
r : &uint32? >> uint32)
:<!refwrt> void =
let
prval () = lemma_array_param u
val () = array_initize_elt<uint32> (q, m, zero32)
val () = r := zero32
fun
loop {i1 : nat | i1 <= m} .<i1>.
(u : &array (uint32, m),
q : &array (uint32, m),
i1 : size_t i1,
r : &uint32)
:<!refwrt> void =
if i1 <> i2sz 0 then
let
val i = pred i1
val tmp = (u32u64 r << 32) lor (u32u64 u[i])
val tmp_q = tmp / u32u64 v and tmp_r = tmp mod (u32u64 v)
in
q[i] := u64u32 tmp_q;
r := u64u32 tmp_r;
loop (u, q, i, r)
end
in
loop (u, q, m, r)
end
fn
fprint_big_integer
{m : int}
(f : FILEref,
m : size_t m,
u : &array (uint32, m))
: void =
let
fun
loop1 (v : &array (uint32, m),
q : &array (uint32, m),
lst : List0 char,
i : uint)
: List0 char =
let
var r : uint32
val () = short_division (m, v, ten32, q, r)
val r = g1ofg0 (u32_2i r)
val () = assertloc ((0 <= r) * (r <= 9))
val digit = int2digit r
in
if big_integer_iseqz (m, q) then
digit :: lst
else if i = 2U then
(* Insert UTF-8 for narrow no-break space U+202F *)
loop1 (q, v, '\xE2' :: '\x80' :: '\xAF' :: digit :: lst, 0U)
else
loop1 (q, v, digit :: lst, succ i)
end
fun
loop2 {n : nat} .<n>.
(lst : list (char, n))
: void =
case+ lst of
| NIL => ()
| hd :: tl => (fprint! (f, hd); loop2 tl)
in
if big_integer_iseqz (m, u) then
fprint! (f, "0")
else
let
val @(pf, pfgc | p) = array_ptr_alloc<uint32> m
val @(qf, qfgc | q) = array_ptr_alloc<uint32> m
val () = array_copy<uint32> (!p, u, m)
val () = array_initize_elt<uint32> (!q, m, zero32)
val () = loop2 (loop1 (!p, !q, NIL, 0U))
val () = array_ptr_free (pf, pfgc | p)
val () = array_ptr_free (qf, qfgc | q)
in
end
end
fn
example_binary (f : FILEref) : void =
let
var u = @[uint32][3] (zero32, zero32, one32)
var v = @[uint32][3] (zero32, zero32, one32)
var w : @[uint32][6]
in
long_multiplication (i2sz 3, i2sz 3, u, v, w);
fprint! (f, "\nBinary long multiplication (b = 2³²)\n\n");
fprint! (f, "u = ");
fprint_big_integer (f, i2sz 3, u);
fprint! (f, "\nv = ");
fprint_big_integer (f, i2sz 3, v);
fprint! (f, "\nu × v = ");
fprint_big_integer (f, i2sz 6, w);
fprint! (f, "\n")
end
fn
test_binary (f : FILEref) : void =
let
var u = @[uint32][3] (mask32, mask32, mask32)
var v = @[uint32][3] (mask32, mask32, mask32)
var w : @[uint32][6]
in
long_multiplication (i2sz 3, i2sz 3, u, v, w);
fprint! (f, "\nThe example numbers specified in the task\n",
"are actually VERY bad for testing binary\n",
"multiplication, because they never need a carry.\n",
"So here is a multiplication full of carries,\n",
"with b = 2³²\n\n");
fprint! (f, "u = ");
fprint_big_integer (f, i2sz 3, u);
fprint! (f, "\nv = ");
fprint_big_integer (f, i2sz 3, v);
fprint! (f, "\nu × v = ");
fprint_big_integer (f, i2sz 6, w);
fprint! (f, "\n")
end
(************** FOR BINARY CODED DECIMAL ARITHMETIC *****************)
(* The following will operate on arrays of BCD digits, with the most
significant digit first. Only the least four bits of a byte will be
considered. This has at least two benefits: any ASCII digit is
treated as its BCD equivalent, and SPACE is treated as zero. *)
extern fn
bcd_multiplication (* Multiply u and v, giving w. *)
{m, n : int}
(m : size_t m,
n : size_t n,
u : &array (char, m),
v : &array (char, n),
w : &array (char?, m + n) >> array (char, m + n))
:<!refwrt> void
fn {}
char2bcd (c : char) :<> intBtwe (0, 9) =
let
val c = char2uchar1 (g1ofg0 c)
val i = g1uint_of_uchar1 c
val i = i mod 16U
val i = i mod 10U (* Guarantees the digit be BCD. *)
in
u2i i
end
extern castfn bcd2char (i : intBtwe (0, 9)) :<> char
(* The following implementation is precisely the algorithm suggested
by Knuth, specialized for b=10. *)
implement
bcd_multiplication {m, n} (m, n, u, v, w) =
let
(* Establish that the arrays have non-negative lengths. *)
prval () = lemma_array_param u
prval () = lemma_array_param v
(* Knuth initializes only part of the w array. However, if we
initialize ALL of w now, then we will not have to deal with
complicated array views later. *)
val () = array_initize_elt<char> (w, m + n, '\0')
(* The following function includes proof of termination. *)
fun
jloop {j : nat | j <= n} .<n - j>.
(u : &array (char, m),
v : &array (char, n),
w : &array (char, m + n),
j : size_t j)
:<!refwrt> void =
if j = n then
()
else if char2bcd v[pred n - j] = 0 then (* Optional branch. *)
begin
w[pred n - j] := '\0';
jloop (u, v, w, succ j)
end
else
let
fun
iloop {i : nat | i <= m} .<m - i>.
(u : &array (char, m),
v : &array (char, n),
w : &array (char, m + n),
i : size_t i,
k : intBtwe (0, 9)) (* carry *)
:<!refwrt> intBtwe (0, 9) = (* carry *)
if i = m then
k
else
let
val ui = char2bcd u[pred m - i]
and vj = char2bcd v[pred n - j]
and wij = char2bcd w[pred (m + n) - (i + j)]
val t = (ui * vj) + wij + k
(* This will prove that 0 <= t *)
prval [ui : int] EQINT () = eqint_make_gint ui
prval [vj : int] EQINT () = eqint_make_gint vj
prval [t : int] EQINT () = eqint_make_gint t
prval () = mul_gte_gte_gte {ui, vj} ()
prval () = prop_verify {0 <= t} ()
(* But I do not feel like proving that t / 10 <= 9. *)
val t_div_10 = t \ndiv 10 and t_mod_10 = t \nmod 10
val () = $effmask_exn assertloc (t_div_10 <= 9)
in
w[pred (m + n) - (i + j)] := bcd2char t_mod_10;
iloop (u, v, w, succ i, t_div_10)
end
val k = iloop (u, v, w, i2sz 0, 0)
in
w[pred n - j] := bcd2char k;
jloop (u, v, w, succ j)
end
in
jloop (u, v, w, i2sz 0)
end
fn
fprint_bcd {m : int}
(f : FILEref,
m : size_t m,
u : &array (char, m))
: void =
let
prval () = lemma_array_param u
fun
skip_zeros {i : nat | i <= m} .<m - i>.
(u : &array (char, m),
i : size_t i)
:<!ref> [i : nat | i <= m] size_t i =
if i = m then
i
else if char2bcd u[i] = 0 then
skip_zeros (u, succ i)
else
i
val [i : int] i = skip_zeros (u, i2sz 0)
fun
loop {j : int | i <= j; j <= m} .<m - j>.
(u : &array (char, m),
j : size_t j)
: void =
if j <> m then
begin
if j <> i && (m - j) mod (i2sz 3) = i2sz 0 then
(* Print UTF-8 for narrow no-break space U+202F *)
fprint! (f, "\xE2\x80\xAF");
fprint! (f, int2digit (char2bcd u[j]));
loop (u, succ j)
end
in
if i = m then
fprint! (f, "0")
else
loop (u, i)
end
fn
string2bcd {n : int}
(s : string n)
: [p : agz]
@(array_v (char, p, n), mfree_gc_v p | ptr p) =
let
val n = strlen s
val @(pf, pfgc | p) = array_ptr_alloc<char> n
implement
array_initize$init<char> (i, x) =
let
val i = g1ofg0 i
prval () = lemma_g1uint_param i
val () = assertloc (i < n)
in
x := s[i]
end
val () = array_initize<char> (!p, n)
in
@(pf, pfgc | p)
end
fn
example_bcd (f : FILEref) : void =
let
val s = g1ofg0 "18446744073709551616"
val m = strlen s
val @(pf_u, pfgc_u | p_u) = string2bcd s
val @(pf_v, pfgc_v | p_v) = string2bcd s
val @(pf_w, pfgc_w | p_w) = array_ptr_alloc<char> (m + m)
macdef u = !p_u
macdef v = !p_v
macdef w = !p_w
in
bcd_multiplication (m, m, u, v, w);
fprint! (f, "\nDecimal long multiplication (b = 10)\n\n");
fprint! (f, "u = ");
fprint_bcd (f, m, u);
fprint! (f, "\nv = ");
fprint_bcd (f, m, v);
fprint! (f, "\nu × v = ");
fprint_bcd (f, m + m, w);
fprint! (f, "\n");
array_ptr_free (pf_u, pfgc_u | p_u);
array_ptr_free (pf_v, pfgc_v | p_v);
array_ptr_free (pf_w, pfgc_w | p_w)
end
(********************************************************************)
implement
main () =
begin
example_binary (stdout_ref);
println! ();
example_bcd (stdout_ref);
println! ();
test_binary (stdout_ref);
println! ();
0
end
</syntaxhighlight>
{{out}}
<pre>$ patscc -std=gnu2x -g -O2 -DATS_MEMALLOC_LIBC long_mult_task.dats && ./a.out
Binary long multiplication (b = 2³²)
u = 18 446 744 073 709 551 616
v = 18 446 744 073 709 551 616
u × v = 340 282 366 920 938 463 463 374 607 431 768 211 456
Decimal long multiplication (b = 10)
u = 18 446 744 073 709 551 616
v = 18 446 744 073 709 551 616
u × v = 340 282 366 920 938 463 463 374 607 431 768 211 456
The example numbers specified in the task
are actually VERY bad for testing binary
multiplication, because they never need a carry.
So here is a multiplication full of carries,
with b = 2³²
u = 79 228 162 514 264 337 593 543 950 335
v = 79 228 162 514 264 337 593 543 950 335
u × v = 6 277 101 735 386 680 763 835 789 423 049 210 091 073 826 769 276 946 612 225
</pre>
=={{header|AutoHotkey}}==
| |||