rowl1-bigint.rlc

;
; rowl - 1st generation
; Copyright (C) 2012 nineties
;
; $Id: rowl1-bigint.rlc 2013-03-05 01:12:44 nineties $
;

(import "rlvm-compile")
(import "rowl1-types")

(compile `object `(

(import "rowl1-util")
(import "rowl1-node")
(import "rowl1-module")
(import "rowl1-compile")
(import "rowl1-assemble")
(import "rowl1-error")

(extern object current_loc)

(export fun pp_bint (ochan n) (
    (output_string ochan (bint_to_s n))
    ))

(export fun bint_digit (u i) (
    (return (array_get int (field_get u 0) i))
    ))

(export fun bint_set_digit (u i d) (
    (array_set int (field_get u 0) i d)
    ))

(export fun bint_set_ndigit (u len) (
    (field_set u 1 len)
    ))

(export fun bint_is_zero (v) (
    (if (& v 1) (return (== v 1)))
    (return (&& (== (field_get v 1) 1) (== (array_get int (field_get v 0) 0) 0)))
    ))

(export fun int_to_bint (v) (
    (var b (allocate_bint 1))
    (if (>= v 0)
        (array_set int (field_get b 0) 0 v)
        (do
            (array_set int (field_get b 0) 0 (- v))
            (field_set b 3 1)
        ))
    (return b)
    ))

(export fun uint_to_bint (v) (
    (var b (allocate_bint 1))
    (array_set int (field_get b 0) 0 v)
    (return b)
    ))

; number of digits in base 2
(export fun bint_ndigits2 (v) (
    (var len (field_get v 1))
    (return (+ (* 32 (- len 1)) (+ (bsr (array_get int (field_get v 0) (- len 1))) 1)))
    ))

(export fun bigneg (u) (
    (field_set u 3 (- 1 (field_get u 3)))
    (return u)
    ))

; assumption: 0 <= u.ndigit - v.ndigit <= 2
; returns quotient and `u' will be overwritten by remainder
(fun bigdiv_sub (u v) (
    (if (== u v) (do
        (clear_bint u) ; remainder == 0
        (return 1)
        ))
    (var ulen (field_get u 1))
    (var vlen (field_get v 1))
    (var udigits (field_get u 0))
    (var vdigits (field_get v 0))
    (var q 0)
    (if (== ulen vlen)
        (= q (udiv (array_get int udigits (- ulen 1))
                (array_get int vdigits (- vlen 1))))
        (do
            (var u1 (array_get int udigits (- ulen 1)))
            (var u0 (array_get int udigits (- ulen 2)))
            (var v0 (array_get int vdigits (- vlen 1)))
            (= q (ludiv u1 u0 v0))
        ))

    (var qv (to_bint_copy v))
    (bigmul qv (uint_to_bint q))
    (while (< (bigcmp u qv) 0) (do
        (-= q 1)
        (bigadd u v)
        ))
    (bigsub u qv)
    (return q)
    ))

; u = u/v and returns u%v (signs are ignored)
(fun bigdiv_nat (u v) (
    (if (bint_is_zero v)
        (throw (division_by_zero current_loc))
        )
    (var c (bigcmp u v))
    (if (== c 0) (do
        (clear_bint u)
        (array_set int (field_get u 0) 0 1)
        (return (uint_to_bint 0))
        ))
    
    (if (< c 0) (do
        (var r (to_bint_copy u))
        (var udigits (field_get u 0))
        (memset udigits 0 (* 4 (field_get u 1)))
        (array_set int udigits 0 0)
        (field_set u 1 1)
        (return r)
        ))

    ; most significant word of v
    (var vms (array_get int (field_get v 0) (- (field_get v 1) 1)))
    (var shift 0)
    (while (! (& vms 0x80000000)) (do
        (+= shift 1)
        (<<= vms 1)
        ))
    (bigshl u shift)
    (bigshl v shift)
    (var udigit (field_get u 0))
    (var vdigit (field_get v 0))
    (var ulen   (field_get u 1))
    (var vlen   (field_get v 1))

    (var q (allocate_bint (+ (- ulen vlen) 2)))
    (var qdigit (field_get q 0))

    (var r (allocate_bint (+ vlen 1)))
    (field_set r 1 (+ vlen 1))
    (var rdigit (field_get r 0))
    (var rlen (+ vlen 1))
    (var p (+ (- ulen vlen) 1))
    (memcpy rdigit (+ udigit (* 4 p)) (* 4 (- vlen 1)))
    (array_set int rdigit vlen 0)

    (array_set int qdigit p (bigdiv_sub r v))
    (while (> p 0) (do
        (-= p 1)
        (= rdigit (field_get r 0))
        (memcpy (+ rdigit 4) rdigit (* 4 vlen))
        (array_set int rdigit 0 (array_get int udigit p))
        (field_set r 1 (+ (field_get r 1) 1))
        (array_set int qdigit p (bigdiv_sub r v))
        ))


    ; set length of q
    (var i (+ (- ulen vlen) 1))
    (while (== (array_get int qdigit i) 0) (do
        (-= i 1)
        ))
    (field_set q 1 (+ i 1))

    ; copy q to u
    (field_set u 0 (field_get q 0))
    (field_set u 1 (field_get q 1))
    (field_set u 2 (field_get q 2))
    (field_set u 3 (field_get q 3))

    (bigshr r shift)
    (bigshr v shift)
    (return r)
    ))

; u = u/v and returns u%v
(export fun bigdiv (u v) (
    (var usign (field_get u 3))
    (var vsign (field_get v 3))
    (if (&& (== usign 0) (== vsign 0)) (return (bigdiv_nat u v)))
    (if (&& (== usign 1) (== vsign 0)) (do
        ; let q, r are numbers such that
        ; |u| = qv + r    (0 <= r < v)
        ; if r == 0
        ;     quotient  => -q
        ;     remainder => 0
        ; else
        ;     quotient  => -(q+1)
        ;     remainder => v - r
        (bigneg u)
        (var r (bigdiv_nat u v))
        (if (bint_is_zero r)
            (do
                (bigneg u)
                (return (uint_to_bint 0))
            )
            (do
                (bigadd u (uint_to_bint 1))
                (bigneg u)
                (bigneg r)
                (bigadd r v)
                (return r)
            ))
        ))
    (if (&& (== usign 0) (== vsign 1)) (do
        ; let q, r are numbers such that
        ; u = q|v| + r  (0 <= r < |v|)
        ; if r == 0
        ;     quotient  => -q
        ;     remainder => 0
        ; else
        ;     quotient  => -(q+1)
        ;     remainder => r + v
        (bigneg v)
        (var r (bigdiv_nat u v))
        (bigneg v)
        (if (bint_is_zero r)
            (do
                (bigneg u)
                (return (uint_to_bint 0))
            )
            (do
                (bigadd u (uint_to_bint 1))
                (bigneg u)
                (bigadd r v)
                (return r)
            ))
        ))
    (if (&& (== usign 1) (== vsign 1)) (do
        (bigneg u)
        (bigneg v)
        (var r (bigdiv_nat u v))
        (bigneg v)
        (bigneg r)
        (return r)
        ))
    ))

(export fun bint_even (u) (
    (if (& (bint_digit u 0) 1)
        (return @FALSE)
        (return @TRUE)
        )
    ))

(export fun bint_to_s (u) (
    (= u (to_bint_copy u))

    (if (bint_is_zero u) (return (string "0")))

    (var udigit (field_get u 0))
    (var ulen (field_get u 1))
    (var usign (field_get u 3))
    (field_set u 3 0)
    ; convert base 2^32 to 10^9
    ; number of digits in base 10^9 is
    ;    ( ulen * log(2^32)/log(10^9) ).ceil
    ; =  ( ulen * 1.070... ).ceil
    ; <= ulen + ( ulen * 3/32 ). ceil
    ; = ulen + ( (ulen * 3 + 31) / 32 ).floor
    (var ary  (allocate_iarray (+ ulen (>> (+ (* 3 ulen) 31) 5))))
    (var base (uint_to_bint 1000000000))
    (var i 0)
    (while (! (bint_is_zero u)) (do
        (var r (bigdiv u base))
        (array_set int ary i (bint_digit r 0))
        (+= i 1)
        ))

    ; convert base 10^9 integer to string
    ; length of the string is
    ;    ( ulen * log(2^32)/log(10) ).ceil
    ; =  ( ulen * 9.632... ).ceil
    ; <= 9*ulen + ( ulen * 3/4 ).ceil
    ; = 9*ulen  + ( (ulen * 3 + 3)/4 ).floor
    (var slen (+ (* 9 ulen) (/ (+ 3 (* 3 ulen)) 4)))
    (if usign (+= slen 1)) ; for '-'
    (var str (allocate_string slen))

    (var p 0) ; position for string
    (if usign (do
        (array_set char str 0 '-')
        (+= p 1)
        ))
    (-= i 1)

    ; skip leading zeros
    (while (== (array_get int ary i) 0) (-= i 1))

    ; skip leading zeros in most significant word
    (var w (array_get int ary i))
    (var b 100000000)
    (while (== (/ w b) 0) (do
        (/= b 10)
        ))
    (while (> b 0) (do
        (array_set char str p (+ (umod (/ w b) 10) '0'))
        (+= p 1)
        (/= b 10)
        ))

    ; fill remaining digits
    (while (> i 0) (do
        (-= i 1)
        (= w (array_get int ary i))
        (= b 100000000)
        (while (> b 0) (do
            (array_set char str p (+ (umod (/ w b) 10) '0'))
            (+= p 1)
            (/= b 10)
            ))
        ))
    ; fill '\0'
    (array_set char str p 0)

    ; recover sign of u
    (field_set u 3 usign)
    (return str)
    ))

(export fun parse_bint (str) (
    (var slen (strlen str))
    (var sign 0)
    (var p 0)
    (if (== (array_get char str 0) '-') (do
        (= sign 1)
        (+= p 1)
        ))
    ; convert base 10 integer to base 2^32
    ; number of digits for bint is
    ;    ( (slen-p+1) * log(10)/log(2^32) ).ceil
    ; =  ( (slen-p+1) * 0.1038... ).ceil
    ; <= ( (slen-p+1) * 1/8 ).ceil
    ; =  ( (slen-p)/8 + 1).floor
    (var ulen (+ (>> (- slen p) 3) 1))
    (var u (allocate_bint ulen))
    (var d (allocate_bint 1))
    (var b (uint_to_bint 1000000000))

    ; construct u by Horner's method with base 10^9
    (var r (umod (- slen p) 9))
    (while (< p slen) (do
        (var w 0)
        (while (> r 0) (do
            (var c (array_get char str p))
            (if (|| (< c '0') (< '9' c))
                (throw (invalid_format current_loc (string str)))
                )
            (= w (+ (* 10 w) (- c '0')))
            (-= r 1)
            (+= p 1)
            ))
        (bint_set_digit d 0 w)
        (bigmul u b)
        (bigadd u d)
        (= r 9)
        ))
    (field_set u 3 sign)
    (return u)
    ))

(export fun parse_bint_binary (str) (
    (var slen (strlen str))
    (var p 0)

    ; skip "0b"
    (+= str 2)
    (-= slen 2)

    (var ulen (/ (+ slen 31) 32))
    (var u (allocate_bint ulen))

    (var j 0)
    (var c 0)
    (var w 0)
    (rfor i 0 slen (do
        (if (== (array_get char str i) '1') (|= w (<< 1 c)))
        (+= c 1)
        (if (== c 32) (do
            (bint_set_digit u j w)
            (= c 0)
            (= w 0)
            (+= j 1)
            ))
        ))
    (bint_set_digit u j w)
    (bint_set_ndigit u ulen)
    (return u)
    ))

(export fun parse_bint_octal (str) (
    (var slen (strlen str))
    (var p 0)

    ; skip "0o"
    (+= str 2)
    (-= slen 2)

    (var ulen (/ (+ (* 3 slen) 31) 32))
    (var u (allocate_bint ulen))

    (var c 0)
    (var w 0)
    (for i 0 slen (do
        (<<= w 3)
        (|= w (- (array_get char str i) '0'))
        (+= c 3)
        (if (== c 30) (do
            (bigshl u 30)
            (bint_set_digit u 0 (| (bint_digit u 0) w))
            (= c 0)
            (= w 0)
            ))
        ))
    (if (> c 0) (do
        (bigshl u c)
        (bint_set_digit u 0 (| (bint_digit u 0) w))
        ))
    (return u)
    ))

(export fun parse_bint_hex (str) (
    (var slen (strlen str))
    (var p 0)

    ; skip "0x"
    (+= str 2)
    (-= slen 2)

    (var ulen (/ (+ (* 4 slen) 31) 32))
    (var u (allocate_bint ulen))

    (var j 0)
    (var c 0)
    (var w 0)
    (rfor i 0 slen (do
        (var d (array_get char str i))
        (if (>= d 'a') (= d (+ (- d 'a') 10))
        (if (>= d 'A') (= d (+ (- d 'A') 10))
            (-= d '0')
            ))
        (|= w (<< d c))
        (+= c 4)
        (if (== c 32) (do
            (bint_set_digit u j w)
            (= c 0)
            (= w 0)
            (+= j 1)
            ))
        ))
    (bint_set_digit u j w)
    (bint_set_ndigit u ulen)
    (return u)
    ))

(fun bint_hash (n) (
    (var h 0)
    (var digit (field_get n 0))
    (var len   (field_get n 1))
    (for i 0 len (do
        (= h (+ (* 13 h) (array_get int digit i)))
        ))
    (if (== (field_get n 3) 1)
        (= h (+ (* 13 h) 1))
        )
    (return h)
    ))

(fun bint_to_i (str) (
    (return (norm (parse_bint str)))
    ))

(fun bint_to_f (b) (
    (= b (copy_bint b))
    (var f (allocate_float))
    (var n (bint_ndigits2 b))
    (var e (- n 1))
    (= b (bint_shiftR b (box (- n 53))))
    (field_set f 0 (bint_digit b 0))
    (var s (field_get b 3))
    (field_set f 1
        (| (<< s 31)
        (| (<< (& (+ e 1023) 0x7ff) 20)
           (& (bint_digit b 1) 0xfffff))))
    (return f)
    ))

(export fun norm (n) (
    (if (& n 1) (return n))
    (if (!= (field_get n 1) 1) (return n))
    (var d (bint_digit n 0))
    (if (u> d @FIXNUM_MAX) (return n))
    (if (field_get n 3)
        (return (box (- d)))
        (return (box d))
        )
    ))

(export fun to_bint (n) (
    (if (& n 1)
        (return (int_to_bint (unbox n)))
        (return n)
        )
    ))
(export fun to_bint_copy (n) (
    (if (& n 1)
        (return (int_to_bint (unbox n)))
        (return (copy_bint n))
        )
    ))

(export fun bint_uplus (n) (
    (return (norm n))
    ))
(export fun bint_uminus (n) (
    (return (bigneg (to_bint_copy n)))
    ))
(export fun bint_add (a b) (
    (= a (to_bint_copy a))
    (return (norm (bigadd a (to_bint b))))
    ))
(export fun bint_sub (a b) (
    (= a (to_bint_copy a))
    (return (norm (bigsub a (to_bint b))))
    ))
(export fun bint_mul (a b) (
    (= a (to_bint_copy a))
    (return (norm (bigmul a (to_bint b))))
    ))
(export fun bint_div (a b) (
    (= a (to_bint_copy a))
    (bigdiv a (to_bint b))
    (return (norm a))
    ))
(export fun bint_mod (a b) (
    (= a (to_bint_copy a))
    (return (norm (bigdiv a (to_bint b))))
    ))
(export fun bint_shiftL (a b) (
    (= a (to_bint_copy a))
    (= b (norm b))
    (if (== (node_type b) @BigIntE)
        (throw (invalid_argument current_loc (string "shift length is too long") b))
        )
    (= b (unbox b))
    (if (> b 0)
        (return (norm (bigshl a b)))
    (if (== b 0)
        (return (norm a))
        (return (norm (bigshr a (- b))))
        ))
    ))
(export fun bint_shiftR (a b) (
    (= a (to_bint_copy a))
    (= b (norm b))
    (if (== (node_type b) @BigIntE)
        (throw (invalid_argument current_loc (string "shift length is too long") b))
        )
    (= b (unbox b))
    (if (> b 0)
        (return (norm (bigshr a b)))
    (if (== b 0)
        (return (norm a))
        (return (norm (bigshl a (- b))))
        ))
    ))
(export fun bint_lt (a b) (
    (= a (to_bint a))
    (= b (to_bint b))
    (if (< (bigcmp a b) 0) (return @C_TRUE) (return @C_FALSE))
    ))
(export fun bint_gt (a b) (
    (= a (to_bint a))
    (= b (to_bint b))
    (if (> (bigcmp a b) 0) (return @C_TRUE) (return @C_FALSE))
    ))
(export fun bint_le (a b) (
    (= a (to_bint a))
    (= b (to_bint b))
    (if (<= (bigcmp a b) 0) (return @C_TRUE) (return @C_FALSE))
    ))
(export fun bint_ge (a b) (
    (= a (to_bint a))
    (= b (to_bint b))
    (if (>= (bigcmp a b) 0) (return @C_TRUE) (return @C_FALSE))
    ))
(export fun bint_eq (a b) (
    (= a (to_bint a))
    (= b (to_bint b))
    (if (== (bigcmp a b) 0) (return @C_TRUE) (return @C_FALSE))
    ))

; comput a^n by binary method
(export fun bint_pow (a n) (
    (= n (norm n))
    (if (== (node_type n) @BigIntE)
        (throw (invalid_argument current_loc (string "exponent is too big") n))
        )
    (if (== n 1) (return (box 1))) ; 1 == box(0)
    (if (< n 0) (throw (invalid_argument current_loc (string "positive integer is required") n)))

    (var v (int_to_bint 1))
    (= a (to_bint a))
    (= n (unbox n))
    (var m 0x40000000)
    (while (! (& n m)) (>>= m 1))

    (while (> m 0) (do
        (bigmul v v)
        (if (& n m) (bigmul v a))
        (>>= m 1)
        ))
    (return (norm v))
    ))

(export fun int_mul2 (a b) (
    (if (== b 1) (return b))    ; 1 == box(0)
    (if (== a 1) (return a))    ; 1 == box(0)

    (= a (unbox a))
    (= b (unbox b))
    (var c (* a b))
    (if (&&
            (&& (<= c @FIXNUM_MAX) (>= c @FIXNUM_MIN))
            (== (/ c b) a)
        )
        (return (box c))
        )
    (return (bigmul (int_to_bint a) (int_to_bint b)))
    ))
(export fun int_add2 (a b) (
    (var c (+ (unbox a) (unbox b)))
    (if (&& (<= c @FIXNUM_MAX) (>= c @FIXNUM_MIN))
        (return (box c))
        (return (int_to_bint c))
        )
    ))
(export fun int_sub2 (a b) (
    (var c (- (unbox a) (unbox b)))
    (if (&& (<= c @FIXNUM_MAX) (>= c @FIXNUM_MIN))
        (return (box c))
        (return (int_to_bint c))
        )
    ))

(fun ib_coerce (a b) ((return (list2 (int_to_bint (unbox a)) b))))
(fun bi_coerce (a b) ((return (list2 a (int_to_bint (unbox b))))))
(fun fb_coerce (a b) ((return (list2 a (bint_to_f b)))))
(fun bf_coerce (a b) ((return (list2 (bint_to_f a) b))))

(extern fun float_infinity)
(extern fun float_uminus)

(export fun setup_bigint (std) (
    (var mod (create_module std (to_sym "BigInt")))

    (add_module_variable mod (to_sym "INTEGER_MIN") @C_FALSE (float_uminus (float_infinity)))
    (add_module_variable mod (to_sym "INTEGER_MAX") @C_FALSE (float_infinity))

    (add_function1 mod (to_sym "uplus") bintT bint_uplus 0)
    (add_function1 mod (to_sym "uminus") bintT bint_uminus 0)
    (add_function1 mod (to_sym "to_s") bintT bint_to_s 0)
    (add_function1 mod (to_sym "hash") bintT bint_hash 0)
    (add_function1 mod (to_sym "to_i") stringT bint_to_i 0)
    (add_function1 mod (to_sym "to_f") bintT bint_to_f 0)

    (add_function2 mod (to_sym "add") bintT bintT bint_add 0)
    (add_function2 mod (to_sym "sub") bintT bintT bint_sub 0)
    (add_function2 mod (to_sym "mul") bintT bintT bint_mul 0)
    (add_function2 mod (to_sym "pow") bintT bintT bint_pow 0)
    (add_function2 mod (to_sym "div") bintT bintT bint_div 0)
    (add_function2 mod (to_sym "mod") bintT bintT bint_mod 0)
    (add_function2 mod (to_sym "lt") bintT bintT bint_lt 0)
    (add_function2 mod (to_sym "gt") bintT bintT bint_gt 0)
    (add_function2 mod (to_sym "le") bintT bintT bint_le 0)
    (add_function2 mod (to_sym "ge") bintT bintT bint_ge 0)
    (add_function2 mod (to_sym "equal") bintT bintT bint_eq 0)
    (add_function2 mod (to_sym "shiftL") bintT intT bint_shiftL 0)
    (add_function2 mod (to_sym "shiftR") bintT intT bint_shiftR 0)
    

    ; overwrite operations between int and int
    (add_function2 mod (to_sym "add") intT intT int_add2 0)
    (add_function2 mod (to_sym "sub") intT intT int_sub2 0)
    (add_function2 mod (to_sym "mul") intT intT int_mul2 0)
    (add_function2 mod (to_sym "pow") intT intT bint_pow 0)
    (add_function2 mod (to_sym "shiftL") intT intT bint_shiftL 0)
    (add_function2 mod (to_sym "shiftR") intT intT bint_shiftR 0)

    (add_function2 mod (to_sym "coerce") intT bintT ib_coerce 0)
    (add_function2 mod (to_sym "coerce") bintT intT bi_coerce 0)
    (add_function2 mod (to_sym "coerce") floatT bintT fb_coerce 0)
    (add_function2 mod (to_sym "coerce") bintT floatT bf_coerce 0)
    ))

    ))