-- Copyright 1993-1998, by the Cecil Project
-- Department of Computer Science and Engineering, University of Washington
-- See the LICENSE file for license information.

--DOC Big integers are an arbitrary-precision representation of integers.

module BigInt {

template object big_int isa integer;
  -- a collection of digits in range [0..big_int_base)
  private field digits(@:big_int):indexed[int];
  private field sign(@:big_int):int;  -- 1, 0, or -1

private method new_big_int(sign:int, digits:indexed[int]):big_int {
    assert(digits.non_empty => { digits.first != 0 },
	   "leading zeros on big_int");
    -- do some canonicalization; no leading zeros, 0 sign <=> no digits
    concrete object isa big_int {
	sign := if(digits.is_empty, { 0 }, { sign }),
	digits := if(sign = 0, { [] }, { digits }) } }

private concrete object zero_big_int isa big_int { digits := [], sign := 0 };


-- this is the greatest power of 10 that is somewhat smaller than sqrt(max_int)
-- (we desire a power of 10 to speed up decimal-base printing)
private let big_int_base:int :=
  as_small_int(10 ** (sqrt(max_int / 2).succ _log_base 10).pred);


-- arithmetic

method =(l@:big_int, r@:big_int):bool {
    l.sign = r.sign & { l.digits = r.digits } }

method < (l@:big_int, r@:big_int):bool {
    compare(l, r, { true }, { false }, { false }) }
method <=(l@:big_int, r@:big_int):bool {
    compare(l, r, { true }, { true }, { false }) }
method >=(l@:big_int, r@:big_int):bool {
    compare(l, r, { false }, { true }, { true }) }
method > (l@:big_int, r@:big_int):bool {
    compare(l, r, { false }, { false }, { true }) }

method compare(l@:big_int, r@:big_int, 
	       if_less:&():`T, if_equal:&():`T, if_greater:&():`T
	       ):T (** recursive_customization **) {
    compare(l.sign, r.sign, if_less, {
	compare_digits(l.digits, r.digits,
		       { -- abs(l) < abs(r)
			   if(l.sign > 0, if_less, if_greater) },
		       if_equal, -- l = r
		       { -- abs(l) > abs(r)
			   if(l.sign > 0, if_greater, if_less) })
    }, if_greater) }


method -(l@:big_int):big_int {
    new_big_int(-l.sign, l.digits) }


method +(l@:big_int, r@:big_int):big_int {
    -- check quickie cases
    if(l.sign = 0, { ^ r });
    if(r.sign = 0, { ^ l });

    -- if same sign, then add digits and preserve sign
    if(l.sign = r.sign, {
	^ new_big_int(l.sign, add_digits(l.digits, r.digits));
    });

    -- opposite sign: subtract smaller from bigger, choosing sign of bigger
    compare_digits(l.digits, r.digits,
		   { -- r bigger than l
		       new_big_int(r.sign, sub_digits(r.digits, l.digits)) },
		   { -- same digits: cancel out
		       zero_big_int },
		   { -- l bigger than r
		       new_big_int(l.sign, sub_digits(l.digits, r.digits)) }) }

method -(l@:big_int, r@:big_int):big_int {
    l + -r }

method *(l@:big_int, r@:big_int):big_int {
    -- check quickie cases
    if(l.sign = 0, { ^ l });
    if(r.sign = 0, { ^ r });
    -- do the work
    new_big_int(l.sign * r.sign, mul_digits(l.digits, r.digits)) }

method /(l@:big_int, r@:big_int):big_int {
    -- check quickie cases
    if(l.sign = 0, { ^ l });
    if(r.sign = 0, { error("division by zero") });
    -- do the work
    new_big_int(l.sign * r.sign, div_digits(l.digits, r.digits)) }

-- returns abs(l) % abs(r), with the sign of l;
-- the identity ((X / Y) * Y + X % Y) = X holds, i.e. X % Y = X - ((X / Y) * Y)
method %(l@:big_int, r@:big_int):big_int {
    l - (l / r) * r }
method %(l@:big_int, r@:int):int {
    as_small_int(l % r.as_big_int) }


method <<(l@:big_int, r@:integer):integer {
    l * 2**r }

method >>(l@:big_int, r@:integer):integer {
    l / 2**r }


-- Low-level helper functions

private method add_digits(l@:indexed[int], r@:indexed[int]):indexed[int] {
    let len:int := l.length;
    let r_len:int := r.length;

    if(len < r_len, { ^ add_digits(r, l) });
    -- assert: l is >= r in length

    -- the number of implicit zeros before r's digits:
    let r_pad:int := len - r_len;

    let result:array[int] := new_array[int](len, 0);

    let var carry:int := 0;
    for(len.pred, 0, -1, &(i:int){
        let var sum:int := l.fetch(i) + r.fetch(i-r_pad, { 0 }) + carry;
	if(sum < big_int_base, {
	    carry := 0;
	}, {
	    sum := sum - big_int_base;
	    carry := 1;
	});
	result!i := sum;
    });

    if(carry > 0, {
	result.add_first(carry);
    });

    result }

private method sub_digits(l@:indexed[int], r@:indexed[int]):indexed[int] {
    let len:int := l.length;
    let r_len:int := r.length;

    assert(len >= r_len, "should be subtracting bigger from smaller");

    -- the number of implicit zeros before r's digits:
    let r_pad:int := len - r_len;

    let result:array[int] := new_array[int](len, 0);

    let var borrow:int := 0;
    for(len.pred, 0, -1, &(i:int){
        let var diff:int := l.fetch(i) - r.fetch(i-r_pad, { 0 }) - borrow;
	if(diff >= 0, {
	    borrow := 0;
	}, {
	    diff := diff + big_int_base;
	    borrow := 1;
	});
	result!i := diff;
    });

    assert(borrow = 0, "should have no borrow left");

    -- remove leading zeros
    while({ result.first = 0 }, { result.remove_first; });

    result }

private method mul_digits(l@:indexed[int], r@:indexed[int]):indexed[int] {
    let l_len:int := l.length;
    let r_len:int := r.length;

    let result:array[int] := new_array[int](l_len + r_len, 0);

    -- for each digit of op1, from least significant to most significant:
    for(l_len.pred, 0, -1, &(i:int){	
	-- compute 0th digit position of this row of multiply
	let pos:int := l_len + r_len - i - 1;
	let var carry:int := 0;
	-- for each digit of op2, from least significant to most significant:
	for(r_len.pred, 0, -1, &(j:int){
	    let prev:int := result!(i+j+1); -- fetch previous pass's value
	    let var digit:int := (l!i) * (r!j) + prev + carry;
	    assert(digit >= 0, "digit overflowed!");
	    if(digit < big_int_base, {
		carry := 0;
	    }, {
		carry := digit / big_int_base;
		digit := digit - carry * big_int_base; -- faster than d = c % b
		assert(digit >= 0, "subtracting carry went too far!");
		assert(digit < big_int_base,
		       "subtracting carry didn't go far enough!");
		assert(carry < big_int_base, "carry too big");
	    });
	    result!(i+j+1) := digit; -- store computed digit
	});
	assert(carry >= 0, "carry is negative!");
	result!i := carry; -- store final carried value
    });

    -- remove leading zeros
    while({ result.first = 0 }, { result.remove_first; });

    result }

private method div_digits(top@:indexed[int],
	                  bottom@:indexed[int]):indexed[int] {
    -- check for base cases
    compare_digits(top, bottom,
		   { -- divisor bigger; returns zero
		       ^ [] },
		   { -- divisor and dividend are the same; returns 1
		       ^ [1] },
		   { -- divisor smaller; must do the division
		             });

    -- extract most significant digits
    let var top_d:int := top.first;
    let bottom_d:int := bottom.first;

    -- initially, peel off top's prefix that's same length as bottom
    let var top_prefix_len:int := bottom.length;
    
    -- compute relevant prefix of top
    let var top_prefix:indexed[int] := view_subrange(top, 0, top_prefix_len-1);

    -- check if this is enough
    compare_digits(bottom, top_prefix, {}, {}, {
        -- need to grab another digit of top, and use 2 most significant digits
	top_prefix_len := top_prefix_len + 1;
	top_prefix := view_subrange(top, 0, top_prefix_len-1);
	top_d := top_d * big_int_base + top.second;
    });

    -- approximate first digit of quotient by dividing most sig digits
    -- (real quotient digit will be q or something smaller)
    let var q:int := top_d / bottom_d;
    assert(q > 0, "should be a non-zero quotient digit");

    -- compute q * bottom; compare to top's prefix
    let var partial_product:indexed[int] := mul_digits([q], bottom);

    -- search for a q that works (binary search would be faster...)
    loop_exit(&(exit:&():none){
	compare_digits(partial_product, top_prefix, exit, exit, {
	    -- q was too big: retry with q-1
	    q := q-1;
	    assert(q > 0, "should be a non-zero quotient digit");
	    partial_product := mul_digits([q], bottom);
	});
    });

    -- compute partial remainder, left over after first quotient digit
    let partial_remainder :=
      if(top_prefix = partial_product,
         { [] },
         { sub_digits(top_prefix, partial_product) });

    -- construct full remainder as concatenation of partial remainder and
    -- rest of top
    let top_suffix := view_subrange(top, top_prefix_len);
    let var remainder:indexed[int] := partial_remainder || top_suffix;

    -- remove leading zeros of remainder
    let var first_non_zero:int := 0;
    while({ first_non_zero < remainder.length &
	    { remainder!first_non_zero = 0 } }, {
	first_non_zero := first_non_zero.succ;
    });
    if(first_non_zero > 0, {
	remainder :=
	  remainder.view_subrange(first_non_zero);
    });

    -- recursively divide rest of number
    let result:indexed[int] := div_digits(remainder, bottom);

    -- prepend first quotient digit and any needed zeros, and return
    [q] || new_i_vector[int](top_suffix.length - result.length, 0) || result }

private method compare_digits(l@:indexed[int], r@:indexed[int],
			      if_less:&():`T, if_equal:&():`T,
			      if_greater:&():`T):T {
    -- quick check of rough sizes
    compare(l.length, r.length,
	    { ^ eval(if_less) }, {}, { ^ eval(if_greater) });

    -- compare digits left to right
    l.length.do(&(i:int){
	compare(l!i, r!i,
		{ ^ eval(if_less) }, {}, { ^ eval(if_greater) });
    });

    eval(if_equal) }


-- Coercions

method as_small_int(l@:big_int, if_overflow:&():`T):int|T {
    let ov:&():none := { ^ if(l = min_int, { min_int }, if_overflow) };
    let var result:int := 0;
    l.digits.do(&(d:int){
	result := add_ov(mul_ov(result, big_int_base, ov), d, ov);
    });
    result * l.sign }

signature as_big_int(integer):big_int;
implementation as_big_int(l@:big_int):big_int { l }
implementation as_big_int(l@:int):big_int {
    if(l = 0, { ^ zero_big_int });
    -- check for this special case, since we can't do abs(min_int)
    -- w/o overflowing
    if(l = min_int, { ^ as_big_int(l + 1) - as_big_int(1) });
    let l_p:int := abs(l);
    -- check another easy case: one digit
    if(l_p < big_int_base, { ^ new_big_int(sign(l), [l_p]) });
    -- the general case: some number of digits
    let digits:array[int] := new_array[int](3);  -- max 3 digits
    -- compute digits right to left
    let var n:int := l_p;
    while({n != 0}, {
	digits.add_first(n % big_int_base);
	n := n / big_int_base;
    });
    new_big_int(sign(l), digits) }

method as_single_float(l@:big_int):single_float {
    let var f:single_float := 0.0;
    let fbase:single_float := big_int_base.as_single_float;
    l.digits.do(&(d:int){ f := f * fbase + d.as_single_float; });
    if(sign(l) < 0, { -f }, { f }) }

method as_double_float(l@:big_int):double_float {
    let var f:double_float := 0.0.as_double_float (-- 0.0d --);
    let fbase:double_float := big_int_base.as_double_float;
    l.digits.do(&(d:int){ f := f * fbase + d.as_double_float; });
    if(sign(l) < 0, { -f }, { f }) }

-- for faster printing

-- number of digits in base-10 print string for each digit in big_int
private let digit_length:int := big_int_base.pred.print_string.length;

method print_string_base(x@:big_int, base:int):string {
    if(x.sign = 0, { ^ "0" });
    if(base != 10, { ^ resend });
    -- base 10 is an easy case to do,
    -- since digits are stored as some power of 10
    let var s:string := if(x.sign < 0, {"-"}, {""});
    let var first:bool := true;
    x.digits.do(&(d:int){
	s := s || if(first,
		     { first := false; d.print_string },
		     { pad_left(d.print_string, digit_length, '0') });
    });
    s }

} -- end module