mirror of
https://github.com/odin-lang/Odin.git
synced 2025-12-29 17:34:34 +00:00
540 lines
13 KiB
Odin
540 lines
13 KiB
Odin
package math_big
|
|
|
|
import "base:builtin"
|
|
import "base:intrinsics"
|
|
import "core:math"
|
|
|
|
Rat :: struct {
|
|
a, b: Int,
|
|
}
|
|
|
|
rat_set_f64 :: proc(dst: ^Rat, f: f64, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(dst)
|
|
context.allocator = allocator
|
|
|
|
EXP_MASK :: 1<<11 - 1
|
|
|
|
bits := transmute(u64)f
|
|
mantissa := bits & (1<<52 - 1)
|
|
exp := int((bits>>52) & EXP_MASK)
|
|
|
|
int_set_from_integer(&dst.b, 1) or_return
|
|
|
|
switch exp {
|
|
case EXP_MASK:
|
|
dst.a.flags += {.Inf}
|
|
return
|
|
case 0:
|
|
exp -= 1022
|
|
case:
|
|
mantissa |= 1<<52
|
|
exp -= 1023
|
|
}
|
|
|
|
shift := 52 - exp
|
|
|
|
for mantissa&1 == 0 && shift > 0 {
|
|
mantissa >>= 1
|
|
shift -= 1
|
|
}
|
|
|
|
int_set_from_integer(&dst.a, mantissa) or_return
|
|
dst.a.sign = .Negative if f < 0 else .Zero_or_Positive
|
|
|
|
if shift > 0 {
|
|
internal_int_shl(&dst.b, &dst.b, shift) or_return
|
|
} else {
|
|
internal_int_shl(&dst.a, &dst.a, -shift) or_return
|
|
}
|
|
|
|
return internal_rat_norm(dst)
|
|
}
|
|
|
|
rat_set_f32 :: proc(dst: ^Rat, f: f32, allocator := context.allocator) -> (err: Error) {
|
|
return rat_set_f64(dst, f64(f), allocator)
|
|
}
|
|
rat_set_f16 :: proc(dst: ^Rat, f: f16, allocator := context.allocator) -> (err: Error) {
|
|
return rat_set_f64(dst, f64(f), allocator)
|
|
}
|
|
|
|
|
|
rat_set_frac :: proc{rat_set_frac_digit, rat_set_frac_int}
|
|
|
|
rat_set_frac_digit :: proc(dst: ^Rat, a, b: DIGIT, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(dst)
|
|
if b == 0 {
|
|
return .Division_by_Zero
|
|
}
|
|
context.allocator = allocator
|
|
internal_set(&dst.a, a) or_return
|
|
internal_set(&dst.b, b) or_return
|
|
return internal_rat_norm(dst)
|
|
}
|
|
|
|
rat_set_frac_int :: proc(dst: ^Rat, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(dst)
|
|
assert_if_nil(a, b)
|
|
if internal_is_zero(b) {
|
|
return .Division_by_Zero
|
|
}
|
|
context.allocator = allocator
|
|
internal_set(&dst.a, a) or_return
|
|
internal_set(&dst.b, b) or_return
|
|
return internal_rat_norm(dst)
|
|
}
|
|
|
|
rat_set_int :: proc(dst: ^Rat, a: ^Int, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(dst)
|
|
assert_if_nil(a)
|
|
context.allocator = allocator
|
|
internal_set(&dst.a, a) or_return
|
|
internal_set(&dst.b, 1) or_return
|
|
return
|
|
}
|
|
|
|
rat_set_digit :: proc(dst: ^Rat, a: DIGIT, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(dst)
|
|
context.allocator = allocator
|
|
internal_set(&dst.a, a) or_return
|
|
internal_set(&dst.b, 1) or_return
|
|
return
|
|
}
|
|
|
|
rat_set_rat :: proc(dst, x: ^Rat, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(dst, x)
|
|
context.allocator = allocator
|
|
internal_set(&dst.a, &x.a) or_return
|
|
internal_set(&dst.b, &x.b) or_return
|
|
return
|
|
}
|
|
|
|
rat_set_u64 :: proc(dst: ^Rat, x: u64, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(dst)
|
|
context.allocator = allocator
|
|
internal_set(&dst.a, x) or_return
|
|
internal_set(&dst.b, 1) or_return
|
|
return
|
|
}
|
|
rat_set_i64 :: proc(dst: ^Rat, x: i64, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(dst)
|
|
context.allocator = allocator
|
|
internal_set(&dst.a, x) or_return
|
|
internal_set(&dst.b, 1) or_return
|
|
return
|
|
}
|
|
|
|
rat_copy :: proc(dst, src: ^Rat, minimize := false, allocator := context.allocator) -> (err: Error) {
|
|
if (dst == src) { return nil }
|
|
|
|
assert_if_nil(dst, src)
|
|
context.allocator = allocator
|
|
int_copy(&dst.a, &src.a, minimize, allocator) or_return
|
|
int_copy(&dst.b, &src.b, minimize, allocator) or_return
|
|
internal_rat_norm(dst) or_return
|
|
return nil
|
|
}
|
|
|
|
internal_rat_destroy :: proc(rationals: ..^Rat) {
|
|
rationals := rationals
|
|
|
|
for &z in rationals {
|
|
internal_int_destroy(&z.a, &z.b)
|
|
}
|
|
}
|
|
|
|
internal_rat_norm :: proc(z: ^Rat, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(z)
|
|
context.allocator = allocator
|
|
switch {
|
|
case internal_is_zero(&z.a):
|
|
z.a.sign = .Zero_or_Positive
|
|
fallthrough
|
|
case internal_is_zero(&z.b):
|
|
int_set_from_integer(&z.b, 1) or_return
|
|
case:
|
|
sign := z.a.sign
|
|
z.a.sign = .Zero_or_Positive
|
|
z.b.sign = .Zero_or_Positive
|
|
|
|
f := &Int{}
|
|
internal_int_gcd(f, &z.a, &z.b) or_return
|
|
if !internal_int_equals_digit(f, 1) {
|
|
f.sign = .Zero_or_Positive
|
|
internal_int_div(&z.a, &z.a, f) or_return
|
|
internal_int_div(&z.b, &z.b, f) or_return
|
|
}
|
|
z.a.sign = sign
|
|
}
|
|
return
|
|
}
|
|
|
|
rat_swap :: proc(a, b: ^Rat) {
|
|
assert_if_nil(a, b)
|
|
#force_inline internal_swap(a, b)
|
|
}
|
|
|
|
internal_rat_swap :: #force_inline proc(a, b: ^Rat) {
|
|
internal_int_swap(&a.a, &b.a)
|
|
internal_int_swap(&a.b, &b.b)
|
|
}
|
|
|
|
rat_sign :: proc(z: ^Rat) -> Sign {
|
|
if z == nil {
|
|
return .Zero_or_Positive
|
|
}
|
|
return z.a.sign
|
|
}
|
|
|
|
rat_is_int :: proc(z: ^Rat) -> bool {
|
|
assert_if_nil(z)
|
|
return internal_is_zero(&z.a) || internal_int_equals_digit(&z.b, 1)
|
|
}
|
|
|
|
rat_is_zero :: proc(z: ^Rat) -> bool {
|
|
return internal_rat_is_zero(z)
|
|
}
|
|
internal_rat_is_zero :: #force_inline proc(z: ^Rat) -> bool {
|
|
assert_if_nil(z)
|
|
return internal_is_zero(&z.a)
|
|
}
|
|
|
|
internal_int_mul_denom :: proc(dst, x, y: ^Int, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(dst, x, y)
|
|
context.allocator = allocator
|
|
switch {
|
|
case internal_is_zero(x) && internal_is_zero(y):
|
|
return internal_set(dst, 1)
|
|
case internal_is_zero(x):
|
|
return internal_set(dst, y)
|
|
case internal_is_zero(y):
|
|
return internal_set(dst, x)
|
|
}
|
|
return int_mul(dst, x, y)
|
|
}
|
|
|
|
internal_int_scale_denom :: proc(dst, x, y: ^Int, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(dst, x, y)
|
|
if internal_is_zero(y) {
|
|
return internal_set(dst, x)
|
|
}
|
|
int_mul(dst, x, y) or_return
|
|
dst.sign = x.sign
|
|
return
|
|
}
|
|
|
|
|
|
rat_add_rat :: proc(dst, x, y: ^Rat, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(dst, x, y)
|
|
context.allocator = allocator
|
|
|
|
a1, a2: Int
|
|
defer internal_destroy(&a1, &a2)
|
|
|
|
internal_int_scale_denom(&a1, &x.a, &y.b) or_return
|
|
internal_int_scale_denom(&a2, &y.a, &x.b) or_return
|
|
int_add(&dst.a, &a1, &a2) or_return
|
|
internal_int_mul_denom(&dst.b, &x.b, &y.b) or_return
|
|
return internal_rat_norm(dst)
|
|
}
|
|
|
|
rat_sub_rat :: proc(dst, x, y: ^Rat, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(dst, x, y)
|
|
context.allocator = allocator
|
|
|
|
a1, a2 := &Int{}, &Int{}
|
|
defer internal_destroy(a1, a2)
|
|
|
|
internal_int_scale_denom(a1, &x.a, &y.b) or_return
|
|
internal_int_scale_denom(a2, &y.a, &x.b) or_return
|
|
int_sub(&dst.a, a1, a2) or_return
|
|
internal_int_mul_denom(&dst.b, &x.b, &y.b) or_return
|
|
return internal_rat_norm(dst)
|
|
}
|
|
|
|
rat_mul_rat :: proc(dst, x, y: ^Rat, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(dst, x, y)
|
|
context.allocator = allocator
|
|
|
|
if x == y {
|
|
internal_sqr(&dst.a, &x.a) or_return
|
|
if internal_is_zero(&x.b) {
|
|
internal_set(&dst.b, 1) or_return
|
|
} else {
|
|
internal_sqr(&dst.a, &x.b) or_return
|
|
}
|
|
return
|
|
}
|
|
|
|
int_mul(&dst.a, &x.a, &y.a) or_return
|
|
internal_int_mul_denom(&dst.b, &x.b, &y.b) or_return
|
|
return internal_rat_norm(dst)
|
|
}
|
|
|
|
rat_div_rat :: proc(dst, x, y: ^Rat, allocator := context.allocator) -> (err: Error) {
|
|
if internal_rat_is_zero(y) {
|
|
return .Division_by_Zero
|
|
}
|
|
context.allocator = allocator
|
|
|
|
a, b := &Int{}, &Int{}
|
|
defer internal_destroy(a, b)
|
|
|
|
internal_int_scale_denom(a, &x.a, &y.b) or_return
|
|
internal_int_scale_denom(b, &y.a, &x.b) or_return
|
|
internal_set(&dst.a, a) or_return
|
|
internal_set(&dst.b, b) or_return
|
|
internal_int_abs(&dst.a, &dst.a)
|
|
internal_int_abs(&dst.b, &dst.b)
|
|
dst.a.sign = .Negative if a.sign != b.sign else .Zero_or_Positive
|
|
return internal_rat_norm(dst)
|
|
}
|
|
|
|
|
|
rat_abs :: proc(dst, x: ^Rat, allocator := context.allocator) -> (err: Error) {
|
|
rat_set_rat(dst, x, allocator) or_return
|
|
internal_abs(&dst.a, &dst.a, allocator) or_return
|
|
return
|
|
}
|
|
rat_neg :: proc(dst, x: ^Rat, allocator := context.allocator) -> (err: Error) {
|
|
rat_set_rat(dst, x, allocator) or_return
|
|
internal_neg(&dst.a, &dst.a, allocator) or_return
|
|
return
|
|
}
|
|
|
|
|
|
rat_is_positive :: proc(z: ^Rat, allocator := context.allocator) -> (ok: bool, err: Error) {
|
|
assert_if_nil(z)
|
|
a := int_is_positive(&z.a, allocator) or_return
|
|
b := int_is_positive(&z.b, allocator) or_return
|
|
return !(a ~ b), nil
|
|
}
|
|
rat_is_negative :: proc(z: ^Rat, allocator := context.allocator) -> (ok: bool, err: Error) {
|
|
assert_if_nil(z)
|
|
a := int_is_positive(&z.a, allocator) or_return
|
|
b := int_is_positive(&z.b, allocator) or_return
|
|
return (a ~ b), nil
|
|
}
|
|
|
|
rat_is_even :: proc(z: ^Rat, allocator := context.allocator) -> (ok: bool, err: Error) {
|
|
assert_if_nil(z)
|
|
if rat_is_int(z) {
|
|
return int_is_even(&z.a, allocator)
|
|
}
|
|
return false, nil
|
|
}
|
|
rat_is_odd :: proc(z: ^Rat, allocator := context.allocator) -> (ok: bool, err: Error) {
|
|
assert_if_nil(z)
|
|
if rat_is_int(z) {
|
|
return int_is_odd(&z.a, allocator)
|
|
}
|
|
return false, nil
|
|
}
|
|
|
|
rat_to_f16 :: proc(z: ^Rat, allocator := context.allocator) -> (f: f16, exact: bool, err: Error) {
|
|
assert_if_nil(z)
|
|
return internal_rat_to_float(f16, z, allocator)
|
|
}
|
|
rat_to_f32 :: proc(z: ^Rat, allocator := context.allocator) -> (f: f32, exact: bool, err: Error) {
|
|
assert_if_nil(z)
|
|
return internal_rat_to_float(f32, z, allocator)
|
|
}
|
|
rat_to_f64 :: proc(z: ^Rat, allocator := context.allocator) -> (f: f64, exact: bool, err: Error) {
|
|
assert_if_nil(z)
|
|
return internal_rat_to_float(f64, z, allocator)
|
|
}
|
|
|
|
internal_rat_to_float :: proc($T: typeid, z: ^Rat, allocator := context.allocator) -> (f: T, exact: bool, err: Error) where intrinsics.type_is_float(T) {
|
|
FSIZE :: 8*size_of(T)
|
|
when FSIZE == 16 {
|
|
MSIZE :: 10
|
|
} else when FSIZE == 32 {
|
|
MSIZE :: 23
|
|
} else when FSIZE == 64 {
|
|
MSIZE :: 52
|
|
} else {
|
|
#panic("unsupported float type")
|
|
}
|
|
|
|
MSIZE1 :: MSIZE+1
|
|
MSIZE2 :: MSIZE+2
|
|
|
|
ESIZE :: FSIZE - MSIZE1
|
|
EBIAS :: 1<<(ESIZE-1) - 1
|
|
EMIN :: 1 - EBIAS
|
|
EMAX :: EBIAS
|
|
|
|
assert_if_nil(z)
|
|
a, b := &z.a, &z.b
|
|
|
|
context.allocator = allocator
|
|
|
|
alen := internal_count_bits(a)
|
|
if alen == 0 {
|
|
return 0, true, nil
|
|
}
|
|
blen := internal_count_bits(b)
|
|
if blen == 0 {
|
|
return T(math.nan_f64()), false, .Division_by_Zero
|
|
}
|
|
|
|
has_sign := a.sign != b.sign
|
|
defer if has_sign {
|
|
f = -builtin.abs(f)
|
|
}
|
|
|
|
exp := alen - blen
|
|
a2, b2 := &Int{}, &Int{}
|
|
defer internal_destroy(a2, b2)
|
|
internal_int_abs(a2, a) or_return
|
|
internal_int_abs(b2, b) or_return
|
|
|
|
if shift := MSIZE2 - exp; shift > 0 {
|
|
internal_int_shl(a2, a2, shift) or_return
|
|
} else if shift < 0 {
|
|
internal_int_shl(b2, b2, -shift) or_return
|
|
}
|
|
|
|
q, r := &Int{}, &Int{}
|
|
defer internal_destroy(q, r)
|
|
|
|
internal_int_divmod(q, r, a2, b2) or_return
|
|
|
|
has_rem := !internal_is_zero(r)
|
|
mantissa := internal_int_get_u64(q) or_return
|
|
|
|
if mantissa>>MSIZE2 == 1 {
|
|
if mantissa&1 == 1 {
|
|
has_rem = true
|
|
}
|
|
mantissa >>= 1
|
|
exp += 1
|
|
}
|
|
|
|
assert(mantissa>>MSIZE1 == 1, "invalid bit result")
|
|
|
|
|
|
if EMIN-MSIZE <= exp && exp <= EMIN {
|
|
shift := uint(EMIN - (exp - 1))
|
|
lost_bits := mantissa & (1<<shift - 1)
|
|
has_rem ||= lost_bits != 0
|
|
mantissa >>= shift
|
|
exp = 2 - EBIAS // exp + shift
|
|
}
|
|
|
|
|
|
exact = !has_rem
|
|
if mantissa&1 != 0 {
|
|
exact = false
|
|
if has_rem || mantissa&2 != 0 {
|
|
mantissa += 1
|
|
if mantissa >= 1<<MSIZE2 {
|
|
mantissa >>= 1
|
|
exp += 1
|
|
}
|
|
}
|
|
}
|
|
|
|
mantissa >>= 1
|
|
|
|
f = T(math.ldexp(f64(mantissa), exp-MSIZE1))
|
|
if math.is_inf(f, 0) {
|
|
exact = false
|
|
}
|
|
return
|
|
}
|
|
|
|
|
|
rat_compare :: proc(x, y: ^Rat, allocator := context.allocator) -> (comparison: int, error: Error) {
|
|
assert_if_nil(x, y)
|
|
context.allocator = allocator
|
|
|
|
a, b: Int
|
|
internal_init_multi(&a, &b) or_return
|
|
defer internal_destroy(&a, &b)
|
|
internal_int_scale_denom(&a, &x.a, &y.b) or_return
|
|
internal_int_scale_denom(&b, &y.a, &x.b) or_return
|
|
return int_compare(&a, &b)
|
|
}
|
|
|
|
|
|
|
|
rat_add_int :: proc(dst, x: ^Rat, y: ^Int, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(dst, x)
|
|
assert_if_nil(y)
|
|
|
|
z: Rat
|
|
rat_set_int(&z, y, allocator) or_return
|
|
defer internal_destroy(&z)
|
|
return rat_add_rat(dst, x, &z, allocator)
|
|
}
|
|
|
|
rat_sub_int :: proc(dst, x: ^Rat, y: ^Int, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(dst, x)
|
|
assert_if_nil(y)
|
|
|
|
z: Rat
|
|
rat_set_int(&z, y, allocator) or_return
|
|
defer internal_destroy(&z)
|
|
return rat_sub_rat(dst, x, &z, allocator)
|
|
}
|
|
|
|
rat_mul_int :: proc(dst, x: ^Rat, y: ^Int, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(dst, x)
|
|
assert_if_nil(y)
|
|
|
|
z: Rat
|
|
rat_set_int(&z, y, allocator) or_return
|
|
defer internal_destroy(&z)
|
|
return rat_mul_rat(dst, x, &z, allocator)
|
|
}
|
|
|
|
rat_div_int :: proc(dst, x: ^Rat, y: ^Int, allocator := context.allocator) -> (err: Error) {
|
|
if internal_is_zero(y) {
|
|
return .Division_by_Zero
|
|
}
|
|
z: Rat
|
|
rat_set_int(&z, y, allocator) or_return
|
|
defer internal_destroy(&z)
|
|
return rat_div_rat(dst, x, &z, allocator)
|
|
}
|
|
|
|
|
|
int_add_rat :: proc(dst: ^Rat, x: ^Int, y: ^Rat, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(x)
|
|
assert_if_nil(dst, y)
|
|
|
|
w: Rat
|
|
rat_set_int(&w, x, allocator) or_return
|
|
defer internal_destroy(&w)
|
|
return rat_add_rat(dst, &w, y, allocator)
|
|
}
|
|
|
|
int_sub_rat :: proc(dst: ^Rat, x: ^Int, y: ^Rat, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(x)
|
|
assert_if_nil(dst, y)
|
|
|
|
w: Rat
|
|
rat_set_int(&w, x, allocator) or_return
|
|
defer internal_destroy(&w)
|
|
return rat_sub_rat(dst, &w, y, allocator)
|
|
}
|
|
|
|
int_mul_rat :: proc(dst: ^Rat, x: ^Int, y: ^Rat, allocator := context.allocator) -> (err: Error) {
|
|
assert_if_nil(x)
|
|
assert_if_nil(dst, y)
|
|
|
|
w: Rat
|
|
rat_set_int(&w, x, allocator) or_return
|
|
defer internal_destroy(&w)
|
|
return rat_mul_rat(dst, &w, y, allocator)
|
|
}
|
|
|
|
int_div_rat :: proc(dst: ^Rat, x: ^Int, y: ^Rat, allocator := context.allocator) -> (err: Error) {
|
|
if internal_is_zero(y) {
|
|
return .Division_by_Zero
|
|
}
|
|
w: Rat
|
|
rat_set_int(&w, x, allocator) or_return
|
|
defer internal_destroy(&w)
|
|
return rat_div_rat(dst, &w, y, allocator)
|
|
} |