-- Module:PythonScript/Librarys/Cmath
-- 複雜(Complex)的數學(Math)庫
local MathCore = {}
-- 內含
-- 1. math庫
-- 2. cmath庫
local MathCoreMain = require("Module:PythonScript/Runtime/Cmath")
local Runtime = require("Module:PythonScript/Runtime")
local PyClass = require("Module:PythonScript/Class")
local Cmath = MathCoreMain.Cmath
local Matrix = MathCoreMain.Matrix
local BigInt = MathCoreMain.BigInt
local function raise_err(err_type, msg)
require("Module:PythonScript/Runtime").raise(err_type, msg)
end
local CmathLib = {}
local function is_inf(x) return x == math.huge or x == -math.huge end
local function is_nan(x) return x ~= x end
local function is_decimal_obj(x)
return type(x) == "table" and x.__type == "instance" and x.__dict and x.__dict._is_decimal == true
end
local function is_fraction_obj(x)
return type(x) == "table" and x.__type == "instance" and x.__dict and x.__dict._is_fraction == true
end
local function unwrap_num(x)
if type(x) == "number" then return x end
if type(x) == "boolean" then return x and 1 or 0 end
if type(x) == "table" then
if x.__type == "int" then return x.value or tonumber(BigInt.to_string(x.sign, x.limbs))
elseif x.__type == "instance" and x.__dict and x.__dict._is_decimal then
return require("Module:PythonScript/Librarys/Decimal").to_number(x)
elseif x.__type == "instance" then
local ok, flt_m = pcall(PyClass.get_attr, x, "__float__")
if ok and flt_m and flt_m.__type == "bound_method" then
return require("Module:PythonScript/Evaluator").invoke(flt_m.__func, {}, {}, flt_m.__self)
end
if x.__builtin_value ~= nil then return unwrap_num(x.__builtin_value) end
end
end
if type(x) == "string" then return tonumber(x) end
raise_err("TypeError", "must be real number")
end
-- 判斷是否為 2D 矩陣 (走矩陣泛函 expm, sinm 等)
local function is_2d_matrix(x)
return type(x) == "table" and (x.__type == "list" or x.__type == "tuple")
and #x.elements > 0 and type(x.elements[1]) == "table"
and (x.elements[1].__type == "list" or x.elements[1].__type == "tuple")
end
-- 判斷是否為任意序列 (1D 陣列走逐元素映射)
local function is_any_seq(x)
return type(x) == "table" and (x.__type == "list" or x.__type == "tuple")
end
-- =========================================================================
-- 無損 BigInt 自動升降級引擎 (閥值：MAX_SAFE_INT = 9007199254740990)
-- =========================================================================
local MAX_SAFE_INT = 9007199254740990
-- 自動探測底層 MathCoreMain.BigInt 的基底 (預設 10^7)
local BI_BASE = tonumber(BigInt.to_string and BigInt.to_string(1, { 0, 1 })) or 10000000
local function is_bigint_obj(x)
return type(x) == "table" and x.__type == "int"
end
local function ul_trim(limbs)
local n = #limbs
while n > 1 and limbs[n] == 0 do
limbs[n] = nil
n = n - 1
end
if n == 0 then limbs[1] = 0 end
return limbs
end
local function ul_is_zero(limbs)
return #limbs == 1 and limbs[1] == 0
end
local function ul_cmp(a, b)
if #a ~= #b then return (#a < #b) and -1 or 1 end
for i = #a, 1, -1 do
if a[i] ~= b[i] then return (a[i] < b[i]) and -1 or 1 end
end
return 0
end
local function ul_add(a, b)
local res = {}
local carry = 0
local n = math.max(#a, #b)
for i = 1, n do
local s = (a[i] or 0) + (b[i] or 0) + carry
if s >= BI_BASE then
res[i] = s - BI_BASE
carry = 1
else
res[i] = s
carry = 0
end
end
if carry > 0 then res[n + 1] = carry end
return res
end
-- 假設 a >= b
local function ul_sub(a, b)
local res = {}
local borrow = 0
for i = 1, #a do
local d = a[i] - (b[i] or 0) - borrow
if d < 0 then
d = d + BI_BASE
borrow = 1
else
borrow = 0
end
res[i] = d
end
return ul_trim(res)
end
local function ul_mul(a, b)
if ul_is_zero(a) or ul_is_zero(b) then return { 0 } end
local res = {}
for i = 1, #a + #b do res[i] = 0 end
for i = 1, #a do
local ai = a[i]
if ai ~= 0 then
local carry = 0
for j = 1, #b do
local idx = i + j - 1
local prod = ai * b[j] + res[idx] + carry
res[idx] = prod % BI_BASE
carry = math.floor(prod / BI_BASE)
end
local idx = i + #b
while carry > 0 do
local s = (res[idx] or 0) + carry
res[idx] = s % BI_BASE
carry = math.floor(s / BI_BASE)
idx = idx + 1
end
end
end
return ul_trim(res)
end
-- 無號長整數除法：回傳商 q, 餘數 r (支援任意大數除數)
local function ul_divmod(a, b)
if ul_is_zero(b) then raise_err("ZeroDivisionError", "integer division or modulo by zero") end
local cmp = ul_cmp(a, b)
if cmp < 0 then return { 0 }, { unpack(a) } end
if cmp == 0 then return { 1 }, { 0 } end
-- 單一 limb 快速通道
if #b == 1 then
local d = b[1]
local q = {}
local rem = 0
for i = #a, 1, -1 do
local cur = rem * BI_BASE + a[i]
q[i] = math.floor(cur / d)
rem = cur % d
end
return ul_trim(q), { rem }
end
-- 多 limb 二分搜尋長除法
local q = {}
local rem = { 0 }
for i = #a, 1, -1 do
if not ul_is_zero(rem) then
table.insert(rem, 1, a[i])
else
rem[1] = a[i]
end
ul_trim(rem)
if ul_cmp(rem, b) < 0 then
q[i] = 0
else
local lo, hi = 1, BI_BASE - 1
local best = 0
while lo <= hi do
local mid = math.floor((lo + hi) * 0.5)
local prod = ul_mul(b, { mid })
local c = ul_cmp(prod, rem)
if c <= 0 then
best = mid
lo = mid + 1
else
hi = mid - 1
end
end
q[i] = best
if best > 0 then
rem = ul_sub(rem, ul_mul(b, { best }))
end
end
end
return ul_trim(q), ul_trim(rem)
end
local function make_bi(sign, limbs)
limbs = ul_trim(limbs)
if ul_is_zero(limbs) then sign = 1 end
return { __type = "int", sign = sign, limbs = limbs }
end
local function bi_add(A, B)
if A.sign == B.sign then
return make_bi(A.sign, ul_add(A.limbs, B.limbs))
end
local c = ul_cmp(A.limbs, B.limbs)
if c == 0 then return make_bi(1, { 0 })
elseif c > 0 then return make_bi(A.sign, ul_sub(A.limbs, B.limbs))
else return make_bi(B.sign, ul_sub(B.limbs, A.limbs)) end
end
local function bi_sub(A, B)
return bi_add(A, { __type = "int", sign = -B.sign, limbs = B.limbs })
end
local function bi_mul(A, B)
return make_bi(A.sign * B.sign, ul_mul(A.limbs, B.limbs))
end
local function bi_divExact(A, B)
local q, _ = ul_divmod(A.limbs, B.limbs)
return make_bi(A.sign * B.sign, q)
end
local function bi_mod_abs(A, B)
local _, r = ul_divmod(A.limbs, B.limbs)
return make_bi(1, r)
end
local function bi_pow_uint(base_bi, exp_bi)
local res, _ = make_bi(1, { 1 }), {}
local b = make_bi(base_bi.sign, { unpack(base_bi.limbs) })
local e_limbs = { unpack(exp_bi.limbs) }
while not ul_is_zero(e_limbs) do
if e_limbs[1] % 2 == 1 then
res = bi_mul(res, b)
end
b = bi_mul(b, b)
e_limbs, _ = ul_divmod(e_limbs, { 2 })
end
return res
end
-- 將精确小整數 (|n| <= 2^53) 轉為 BigInt
local function small_int_to_bi(n)
if n == 0 then return make_bi(1, { 0 }) end
local sign = (n < 0) and -1 or 1
local v = math.abs(n)
local limbs = {}
while v > 0 do
table.insert(limbs, v % BI_BASE)
v = math.floor(v / BI_BASE)
end
return make_bi(sign, limbs)
end
-- 將任意 IEEE 754 整數浮點數 (含 > 2^53 甚至 1e300) 無損轉為 BigInt
local function float_int_to_bi(x)
if x == 0 then return make_bi(1, { 0 }) end
local sign = (x < 0) and -1 or 1
local ax = math.abs(x)
if ax <= 9007199254740991 then
return small_int_to_bi(x)
end
local m, e = math.frexp(ax)
local mant = math.ldexp(m, 53) -- 精確落在 [2^52, 2^53 - 1]
local b = small_int_to_bi(mant)
local shift = e - 53
if shift > 0 then
local p2 = bi_pow_uint(small_int_to_bi(2), small_int_to_bi(shift))
b = bi_mul(b, p2)
end
b.sign = sign
return b
end
-- 將 BigInt 依閥值 9007199254740990 決定回傳 Lua number 或保持 BigInt
local MAX_SAFE_BI = small_int_to_bi(MAX_SAFE_INT)
local function bi_to_py_or_num(B_obj)
if ul_cmp(B_obj.limbs, MAX_SAFE_BI.limbs) <= 0 then
local val = 0
for i = #B_obj.limbs, 1, -1 do
val = val * BI_BASE + B_obj.limbs[i]
end
return B_obj.sign * val
end
-- 💡 關鍵修復：將原始表轉換為字串，並呼叫 Python int 類別重新包裝！
local s = BigInt.to_string(B_obj.sign, B_obj.limbs)
return Runtime.builtins.int.__call({ s })
end
-- 判斷輸入是否為整數 (含 BigInt 或誤差 <= 1e-17 的浮點數)
-- 回傳: is_int (boolean), num_val (number 或 nil), bi_val (BigInt)
local function inspect_int_arg(v)
if is_bigint_obj(v) then
local norm_bi = make_bi(v.sign or 1, v.limbs or { 0 })
local maybe_num = bi_to_py_or_num(norm_bi)
if type(maybe_num) == "number" then
return true, maybe_num, norm_bi
else
return true, (norm_bi.sign * math.huge), norm_bi
end
end
local x = unwrap_num(v)
if is_nan(x) or is_inf(x) then return false, x, nil end
local nearest = math.floor(x + 0.5)
if math.abs(x - nearest) <= 1e-17 then
return true, nearest, float_int_to_bi(nearest)
end
return false, x, nil
end
-- 偵測輸入是否為實數整數 (支援 int, BigInt, float 整數，以及虛部為 0 的 complex 如 3+0j)
local function inspect_complex_int_arg(v)
if v == nil or v == Runtime.builtins.None or is_any_seq(v) then
return false, nil, nil
end
if is_bigint_obj(v) then
return inspect_int_arg(v)
end
if type(v) == "table" and v.__type == "complex" then
local r, i, k, jk = Cmath.get_parts(v)
if math.abs(i) <= 1e-17 and math.abs(k) <= 1e-17 and math.abs(jk) <= 1e-17 then
return inspect_int_arg(r)
end
return false, nil, nil
end
local ok, is_i, num_v, bi_v = pcall(inspect_int_arg, v)
if ok and is_i then return true, num_v, bi_v end
return false, nil, nil
end
-- 計算正整數組合 C(N_bi, K_num) 並自動處理 BigInt 溢位
local function comb_pos_exact(n_is_huge, n_num, n_bi, k_num)
if k_num == 0 then return 1 end
if not n_is_huge then
if k_num > n_num then return 0 end
if k_num == n_num then return 1 end
k_num = math.min(k_num, n_num - k_num)
if k_num == 0 then return 1 end
local c = 1
local is_big = false
local c_bi = nil
local function num_gcd(a, b)
while b ~= 0 do a, b = b, a % b end
return a
end
for i = 1, k_num do
local num = n_num - i + 1
local den = i
if not is_big then
local g = num_gcd(num, den)
num = num / g
den = den / g
c = c / den -- 必定整除
if c * num > MAX_SAFE_INT then
is_big = true
c_bi = bi_mul(float_int_to_bi(c), float_int_to_bi(num))
else
c = c * num
end
else
c_bi = bi_divExact(bi_mul(c_bi, float_int_to_bi(num)), float_int_to_bi(den))
end
end
return is_big and bi_to_py_or_num(c_bi) or c
else
-- N 本身已超過 9007199254740990
local c_bi = make_bi(1, { 1 })
local one_bi = make_bi(1, { 1 })
local cur_n = n_bi
for i = 1, k_num do
c_bi = bi_divExact(bi_mul(c_bi, cur_n), float_int_to_bi(i))
cur_n = bi_sub(cur_n, one_bi)
end
return bi_to_py_or_num(c_bi)
end
end
-- Python cmath 標準庫 (支援複數、四元數純量、1D 向量逐元素映射與 2D 矩陣泛函微積分)
function CmathLib.init()
local get_parts = Cmath.get_parts
-- 修補內建 abs() 以支援矩陣 Frobenius 範數與四元數模長
local B = Runtime.builtins
if B and B.abs and not B._matrix_abs_patched then
B._matrix_abs_patched = true
local orig_abs = B.abs
B.abs = function(args)
local x = args and args[1]
if is_any_seq(x) then return Matrix.frobenius_norm(x, false)
elseif type(x) == "table" and x.__type == "complex" then return Cmath.abs(x) end
return orig_abs(args)
end
end
local function to_q4(val)
local a, b, c, d = get_parts(val)
return { a, b, c, d }
end
local function from_q4(q, force_complex)
local a = Cmath.clean_zero(q[1])
local b = (math.abs(q[2]) < 1e-12) and 0 or Cmath.clean_zero(q[2])
local c = (math.abs(q[3]) < 1e-12) and 0 or Cmath.clean_zero(q[3])
local d = (math.abs(q[4]) < 1e-12) and 0 or Cmath.clean_zero(q[4])
if not force_complex and b == 0 and c == 0 and d == 0 then return a end
return Cmath.make_complex(a, b, c, d)
end
-- 四元數純量底層算子
local function q_add(p, q) return { p[1]+q[1], p[2]+q[2], p[3]+q[3], p[4]+q[4] } end
local function q_sub(p, q) return { p[1]-q[1], p[2]-q[2], p[3]-q[3], p[4]-q[4] } end
local function q_neg(p) return { -p[1], -p[2], -p[3], -p[4] } end
local function q_scale(p, s) return { p[1]*s, p[2]*s, p[3]*s, p[4]*s } end
local function q_abs_sq(p) return p[1]*p[1] + p[2]*p[2] + p[3]*p[3] + p[4]*p[4] end
local function q_abs(p) return math.sqrt(q_abs_sq(p)) end
local function q_vec_norm(p) return math.sqrt(p[2]*p[2] + p[3]*p[3] + p[4]*p[4]) end
local function q_mul(p, q)
local a, b, c, d = p[1], p[2], p[3], p[4]
local t, x, y, z = q[1], q[2], q[3], q[4]
return {
a*t - b*x - c*y - d*z,
a*x + b*t + c*z - d*y,
a*y - b*z + c*t + d*x,
a*z + b*y - c*x + d*t
}
end
local function q_inv(p)
local n2 = q_abs_sq(p)
if n2 == 0 then raise_err("ZeroDivisionError", "complex/quaternion division by zero") end
return { p[1]/n2, -p[2]/n2, -p[3]/n2, -p[4]/n2 }
end
local function q_rdiv(p, q) return q_mul(p, q_inv(q)) end
local function q_vec_unit(p, default_j)
local vn = q_vec_norm(p)
if vn == 0 then return default_j and { 0, 1, 0, 0 } or { 0, 0, 0, 0 }, 0 end
return { 0, p[2]/vn, p[3]/vn, p[4]/vn }, vn
end
local function q_arg(p)
local vn = q_vec_norm(p)
if p[1] == 0 and vn == 0 then return 0.0 end
return math.atan2(vn, p[1])
end
local function q_exp(p)
local ea = math.exp(p[1])
local u, vn = q_vec_unit(p, false)
if vn == 0 then return { ea, 0, 0, 0 } end
local cv, sv = ea * math.cos(vn), ea * math.sin(vn)
return { cv, u[2]*sv, u[3]*sv, u[4]*sv }
end
local function q_log(p)
local mag = q_abs(p)
if mag == 0 then raise_err("ValueError", "math domain error") end
local u, vn = q_vec_unit(p, true)
local theta = math.atan2(vn, p[1])
return { math.log(mag), u[2]*theta, u[3]*theta, u[4]*theta }
end
local function q_sqrt(p)
local mag = q_abs(p)
if mag == 0 then return { 0, 0, 0, 0 } end
local u, vn = q_vec_unit(p, true)
if vn == 0 and p[1] >= 0 then return { math.sqrt(p[1]), 0, 0, 0 } end
local wr = math.sqrt(math.max(0, (mag + p[1]) * 0.5))
local wv = math.sqrt(math.max(0, (mag - p[1]) * 0.5))
return { wr, u[2]*wv, u[3]*wv, u[4]*wv }
end
local function r_sinh(x) if x == 0 then return 0.0 end; local e = math.exp(x); return (e - 1.0/e)*0.5 end
local function r_cosh(x) if x == 0 then return 1.0 end; local e = math.exp(x); return (e + 1.0/e)*0.5 end
local function q_sin(p)
local a = p[1]; local u, vn = q_vec_unit(p, false)
local sc, cs = math.sin(a) * r_cosh(vn), math.cos(a) * r_sinh(vn)
return { sc, u[2]*cs, u[3]*cs, u[4]*cs }
end
local function q_cos(p)
local a = p[1]; local u, vn = q_vec_unit(p, false)
local cc, ss = math.cos(a) * r_cosh(vn), -math.sin(a) * r_sinh(vn)
return { cc, u[2]*ss, u[3]*ss, u[4]*ss }
end
local function q_tan(p) return q_rdiv(q_sin(p), q_cos(p)) end
local function q_sinh(p)
local a = p[1]; local u, vn = q_vec_unit(p, false)
local sc, cs = r_sinh(a) * math.cos(vn), r_cosh(a) * math.sin(vn)
return { sc, u[2]*cs, u[3]*cs, u[4]*cs }
end
local function q_cosh(p)
local a = p[1]; local u, vn = q_vec_unit(p, false)
local cc, ss = r_cosh(a) * math.cos(vn), r_sinh(a) * math.sin(vn)
return { cc, u[2]*ss, u[3]*ss, u[4]*ss }
end
local function q_tanh(p) return q_rdiv(q_sinh(p), q_cosh(p)) end
local function q_asinh(p)
local p2 = q_mul(p, p)
return q_log(q_add(p, q_sqrt({ p2[1] + 1, p2[2], p2[3], p2[4] })))
end
local function q_acosh(p)
local sp1 = q_sqrt({ p[1] + 1, p[2], p[3], p[4] })
local sm1 = q_sqrt({ p[1] - 1, p[2], p[3], p[4] })
return q_log(q_add(p, q_mul(sp1, sm1)))
end
local function q_atanh(p)
local one_plus_p = { 1 + p[1], p[2], p[3], p[4] }
local one_minus_p = { 1 - p[1], -p[2], -p[3], -p[4] }
return q_scale(q_sub(q_log(one_plus_p), q_log(one_minus_p)), 0.5)
end
local function q_asin(p) local u = q_vec_unit(p, true); return q_mul(q_neg(u), q_asinh(q_mul(p, u))) end
local function q_acos(p) local u = q_vec_unit(p, true); return q_mul(q_neg(u), q_acosh(p)) end
local function q_atan(p) local u = q_vec_unit(p, true); return q_mul(q_neg(u), q_atanh(q_mul(p, u))) end
local cmath_api = {}
cmath_api.pi = math.pi; cmath_api.e = math.exp(1); cmath_api.tau = math.pi * 2
cmath_api.inf = math.huge; cmath_api.nan = 0/0
cmath_api.infj = Cmath.make_complex(0, math.huge, 0, 0); cmath_api.nanj = Cmath.make_complex(0, 0/0, 0, 0)
cmath_api.infk = Cmath.make_complex(0, 0, math.huge, 0); cmath_api.nank = Cmath.make_complex(0, 0, 0/0, 0)
cmath_api.infjk = Cmath.make_complex(0, 0, 0, math.huge); cmath_api.nanjk = Cmath.make_complex(0, 0, 0, 0/0)
cmath_api.j = Cmath.make_complex(0, 1, 0, 0)
cmath_api.k = Cmath.make_complex(0, 0, 1, 0)
cmath_api.jk = Cmath.make_complex(0, 0, 0, 1)
cmath_api.quaternion = function(args)
return Cmath.make_complex(
args[1] and unwrap_num(args[1]) or 0,
args[2] and unwrap_num(args[2]) or 0,
args[3] and unwrap_num(args[3]) or 0,
args[4] and unwrap_num(args[4]) or 0
)
end
cmath_api.sgn = function(args)
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], Cmath.sgn) end
return Cmath.sgn(args[1])
end
cmath_api.arg = function(args)
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], Cmath.arg) end
return Cmath.arg(args[1])
end
cmath_api.dot = function(args) return Cmath.dot(args[1], args[2]) end
cmath_api.cross = function(args) return Cmath.cross(args[1], args[2]) end
cmath_api.conjugate = function(args) return Cmath.conj(args[1]) end
cmath_api.scalar = function(args) return Cmath.scalar(args[1]) end
cmath_api.vector = function(args) return Cmath.vector(args[1]) end
cmath_api.phase = function(args)
local function scalar_phase(z)
local a, b, k_val, d_val = get_parts(z)
if k_val == 0 and d_val == 0 then return math.atan2(b, a) end
return q_arg({ a, b, k_val, d_val })
end
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], scalar_phase) end
return scalar_phase(args[1])
end
cmath_api.polar = function(args)
local q = to_q4(args[1])
local ph = (q[3] == 0 and q[4] == 0) and math.atan2(q[2], q[1]) or q_arg(q)
return { __type = "tuple", elements = { q_abs(q), ph } }
end
cmath_api.rect = function(args)
local mag, ph = unwrap_num(args[1]), unwrap_num(args[2])
return Cmath.make_complex(mag * math.cos(ph), mag * math.sin(ph), 0, 0)
end
-- 泛函與逐元素自動分流器 (2D 走矩陣泛函，1D 走逐元素，純量走四元數算子)
local function make_unary_dispatcher(mat_fn, q_fn)
local function scalar_op(x) return from_q4(q_fn(to_q4(x)), true) end
return function(args)
local a1 = args[1]
if is_2d_matrix(a1) then return mat_fn(a1, false)
elseif is_any_seq(a1) then return Matrix.unary_op(a1, scalar_op) end
return scalar_op(a1)
end
end
cmath_api.gamma = function(args)
local function scalar_cgamma(v)
local r, i = Cmath.extract_complex_2d(v, "cmath.gamma")
local gr, gi = Cmath.c2_gamma(r, i)
return from_q4({ gr, gi, 0, 0 }, true)
end
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], scalar_cgamma) end
return scalar_cgamma(args[1])
end
cmath_api.lgamma = function(args)
local function scalar_clgamma(v)
local r, i = Cmath.extract_complex_2d(v, "cmath.lgamma")
local lgr, lgi = Cmath.c2_lgamma(r, i)
return from_q4({ lgr, lgi, 0, 0 }, true)
end
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], scalar_clgamma) end
return scalar_clgamma(args[1])
end
cmath_api.factorial = function(args)
local py_int = Runtime.builtins.int
local Evaluator = require("Module:PythonScript/Evaluator")
local function scalar_cfact(v)
local r, i, k, jk = get_parts(v)
-- 🛡️ 拒絕四元數不交換律的複雜操作
if (k and k ~= 0) or (jk and jk ~= 0) then
raise_err("NotImplementedError", "cmath.factorial() does not support Quaternions")
end
local nearest_r = math.floor(r + 0.5)
-- 判斷是否為純實數整數 (虛部趨近於 0，且實部與整數誤差小於 1e-17)
if math.abs(i) <= 1e-17 and math.abs(r - nearest_r) <= 1e-17 then
if nearest_r < 0 then
raise_err("ValueError", "factorial() not defined for negative integers")
end
local res = 1
local is_big = false
for idx = 1, nearest_r do
if not is_big and res * idx > 9007199254740990 then
res = py_int.__call({ tostring(math.floor(res)) })
is_big = true
end
if is_big then
local i_obj = py_int.__call({ tostring(idx) })
local mul_m = PyClass.get_attr(res, "__mul__")
res = Evaluator.invoke(mul_m.__func, {i_obj}, {}, mul_m.__self)
else
res = res * idx
end
end
return res
else
-- 複數或非整數，使用高精度複變 Gamma 函數計算
return cmath_api.gamma({ Cmath.make_complex(r + 1.0, i, 0, 0) })
end
end
-- 同樣支援陣列廣播
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], scalar_cfact) end
return scalar_cfact(args[1])
end
cmath_api.pow = function(args)
local a1, a2, a3 = args[1], args[2], args[3]
local has_mod = (a3 ~= nil and a3 ~= Runtime.builtins.None)
if not is_any_seq(a1) and not is_any_seq(a2) then
local ok1, x_num, x_bi = inspect_complex_int_arg(a1)
local ok2, y_num, y_bi = inspect_complex_int_arg(a2)
if has_mod then
local ok3, _, z_bi = inspect_complex_int_arg(a3)
if not (ok1 and ok2 and ok3) then
raise_err("TypeError", "pow() 3rd argument not allowed unless all arguments are integers")
end
return BigInt.powmod(x_bi, y_bi, z_bi)
end
if ok1 and ok2 then
if ul_is_zero(y_bi.limbs) then return 1 end
if y_bi.sign > 0 then
local need_big = (math.abs(x_num) > MAX_SAFE_INT) or (y_num > MAX_SAFE_INT)
if not need_big then
local est = math.abs(x_num) ^ y_num
if est > MAX_SAFE_INT or is_inf(est) then need_big = true end
end
if need_big then
return bi_to_py_or_num(bi_pow_uint(x_bi, y_bi))
end
return x_num ^ y_num
end
end
elseif has_mod then
raise_err("TypeError", "pow() 3rd argument not allowed unless all arguments are integers")
end
return Cmath.pow(a1, a2)
end
cmath_api.exp = make_unary_dispatcher(Matrix.expm, q_exp)
cmath_api.sqrt = make_unary_dispatcher(Matrix.sqrtm, q_sqrt)
cmath_api.sin = make_unary_dispatcher(Matrix.sinm, q_sin)
cmath_api.cos = make_unary_dispatcher(Matrix.cosm, q_cos)
cmath_api.tan = make_unary_dispatcher(Matrix.tanm, q_tan)
cmath_api.sinh = make_unary_dispatcher(Matrix.sinhm, q_sinh)
cmath_api.cosh = make_unary_dispatcher(Matrix.coshm, q_cosh)
cmath_api.tanh = make_unary_dispatcher(Matrix.tanhm, q_tanh)
cmath_api.asin = make_unary_dispatcher(Matrix.asinm, q_asin)
cmath_api.acos = make_unary_dispatcher(Matrix.acosm, q_acos)
cmath_api.atan = make_unary_dispatcher(Matrix.atanm, q_atan)
cmath_api.asinh = make_unary_dispatcher(Matrix.asinhm, q_asinh)
cmath_api.acosh = make_unary_dispatcher(Matrix.acoshm, q_acosh)
cmath_api.atanh = make_unary_dispatcher(Matrix.atanhm, q_atanh)
cmath_api.log = function(args)
local z, base = args[1], args[2]
if is_2d_matrix(z) then return Matrix.logm(z, base, false) end
local function scalar_log(v)
local lz = q_log(to_q4(v))
if base ~= nil and not is_2d_matrix(base) then
return from_q4(q_rdiv(lz, q_log(to_q4(base))), true)
end
return from_q4(lz, true)
end
if is_any_seq(z) then return Matrix.unary_op(z, scalar_log) end
local lz = q_log(to_q4(z))
if is_2d_matrix(base) then return Matrix.rdiv(from_q4(lz, true), Matrix.logm(base, nil, false)) end
return scalar_log(z)
end
cmath_api.log10 = function(args)
local a1 = args[1]
if is_2d_matrix(a1) then return Matrix.logm(a1, 10, false) end
local function scalar_log10(v) return from_q4(q_scale(q_log(to_q4(v)), 1 / math.log(10)), true) end
if is_any_seq(a1) then return Matrix.unary_op(a1, scalar_log10) end
return scalar_log10(a1)
end
cmath_api.isfinite = function(args)
local function fn_fin(a)
local r, i, k_val, d_val = get_parts(a)
return not (is_inf(r) or is_inf(i) or is_inf(k_val) or is_inf(d_val) or is_nan(r) or is_nan(i) or is_nan(k_val) or is_nan(d_val))
end
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], fn_fin) end
return fn_fin(args[1])
end
cmath_api.isinf = function(args)
local function fn_inf(a)
local r, i, k_val, d_val = get_parts(a)
return is_inf(r) or is_inf(i) or is_inf(k_val) or is_inf(d_val)
end
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], fn_inf) end
return fn_inf(args[1])
end
cmath_api.isnan = function(args)
local function fn_nan(a)
local r, i, k_val, d_val = get_parts(a)
return is_nan(r) or is_nan(i) or is_nan(k_val) or is_nan(d_val)
end
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], fn_nan) end
return fn_nan(args[1])
end
cmath_api.isclose = function(args, kwargs)
local rel_tol = (kwargs and kwargs.rel_tol) or 1e-09
local abs_tol = (kwargs and kwargs.abs_tol) or 0.0
local function check_close(x, y)
if is_any_seq(x) and is_any_seq(y) then
if #x.elements ~= #y.elements then return false end
for i = 1, #x.elements do
if not check_close(x.elements[i], y.elements[i]) then return false end
end
return true
end
local qa, qb = to_q4(x), to_q4(y)
local diff_mag = q_abs(q_sub(qa, qb))
local max_mag = math.max(q_abs(qa), q_abs(qb))
return diff_mag <= math.max(rel_tol * max_mag, abs_tol)
end
return check_close(args[1], args[2])
end
local inst = PyClass.make_instance(PyClass.make_class("module", {}, {}))
inst.__dict["__name__"] = "cmath"
for k, v in pairs(cmath_api) do
if type(v) == "function" then inst.__dict[k] = { __type = "builtin_method", __call = v } else inst.__dict[k] = v end
end
for k, v in pairs(MathCoreMain.linalg) do
if type(v) == "function" then inst.__dict[k] = { __type = "builtin_method", __call = v } else inst.__dict[k] = v end
end
return inst
end
MathCore.CmathLib = CmathLib
-- Python math 標準庫 (支援純量、1D 偽 NumPy 逐元素向量化與 2D 矩陣泛函微積分)
local MathLib = {}
function MathLib.init()
-- 修補內建 abs()
local Builtins = Runtime.builtins
if Builtins and Builtins.abs and not Builtins._matrix_abs_patched then
Builtins._matrix_abs_patched = true
local orig_abs = Builtins.abs
Builtins.abs = function(args)
local x = args and args[1]
if is_any_seq(x) then return Matrix.frobenius_norm(x, false)
elseif type(x) == "table" and x.__type == "complex" then return Cmath.abs(x) end
return orig_abs(args)
end
end
local math_api = {}
math_api.pi = math.pi; math_api.e = math.exp(1); math_api.tau = math.pi * 2
math_api.inf = math.huge; math_api.nan = 0/0
-- 自動分流器：2D 矩陣 -> 矩陣泛函；1D 陣列 -> 逐元素運算；純量 -> 實數運算
local function make_real_dispatcher(mat_fn, real_fn)
return function(args)
local a1 = args[1]
if is_2d_matrix(a1) then return mat_fn(a1, true)
elseif is_any_seq(a1) then return Matrix.unary_op(a1, real_fn) end
return real_fn(a1)
end
end
math_api.exp = make_real_dispatcher(Matrix.expm, function(v) return math.exp(unwrap_num(v)) end)
math_api.expm1 = make_real_dispatcher(Matrix.expm1, function(v) return math.exp(unwrap_num(v)) - 1 end)
math_api.log1p = make_real_dispatcher(Matrix.log1p, function(v) return math.log(1 + unwrap_num(v)) end)
math_api.sqrt = make_real_dispatcher(Matrix.sqrtm, function(v)
local x = unwrap_num(v); if x < 0 then raise_err("ValueError", "math domain error") end
return math.sqrt(x)
end)
math_api.cbrt = make_real_dispatcher(Matrix.cbrtm, function(v)
local x = unwrap_num(v); return x < 0 and -math.pow(-x, 1/3) or math.pow(x, 1/3)
end)
math_api.gamma = function(args)
local function scalar_gamma(v)
local x = unwrap_num(v)
if x == -math.huge then raise_err("ValueError", "math domain error") end
if x == math.huge then return math.huge end
local gr, _ = Cmath.c2_gamma(x, 0.0)
if is_inf(gr) then raise_err("OverflowError", "math range error") end
return Cmath.clean_zero(gr)
end
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], scalar_gamma) end
return scalar_gamma(args[1])
end
math_api.lgamma = function(args)
local function scalar_lgamma(v)
local res = Cmath.r_lgamma(unwrap_num(v))
return Cmath.clean_zero(res)
end
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], scalar_lgamma) end
return scalar_lgamma(args[1])
end
math_api.sumprod = function(args)
local p, q = args[1], args[2]
if not (is_any_seq(p) and is_any_seq(q)) then
raise_err("TypeError", "sumprod() arguments must be sequences")
end
if #p.elements ~= #q.elements then
raise_err("ValueError", "sumprod() inputs do not have the same length")
end
local total = 0.0
for i = 1, #p.elements do
total = total + unwrap_num(p.elements[i]) * unwrap_num(q.elements[i])
end
return total
end
math_api.pow = function(args)
local a1, a2, a3 = args[1], args[2], args[3]
local has_mod = (a3 ~= nil and a3 ~= Runtime.builtins.None)
if is_any_seq(a1) or is_any_seq(a2) then
if has_mod then
raise_err("TypeError", "pow() 3rd argument not allowed unless all arguments are integers")
end
return Matrix.pow(a1, a2)
end
local ok1, x_num, x_bi = inspect_int_arg(a1)
local ok2, y_num, y_bi = inspect_int_arg(a2)
if has_mod then
local ok3, _, z_bi = inspect_int_arg(a3)
if not (ok1 and ok2 and ok3) then
raise_err("TypeError", "pow() 3rd argument not allowed unless all arguments are integers")
end
return BigInt.powmod(x_bi, y_bi, z_bi)
end
if ok1 and ok2 and y_bi.sign > 0 and not ul_is_zero(y_bi.limbs) then
local need_big = (math.abs(x_num) > MAX_SAFE_INT) or (y_num > MAX_SAFE_INT)
if not need_big then
local est = math.abs(x_num) ^ y_num
if est > MAX_SAFE_INT or is_inf(est) then need_big = true end
end
if need_big then
return bi_to_py_or_num(bi_pow_uint(x_bi, y_bi))
end
elseif ok1 and ok2 and ul_is_zero(y_bi.limbs) then
return 1
end
return unwrap_num(a1) ^ unwrap_num(a2)
end
math_api.log = function(args)
local a1, a2 = args[1], args[2]
if is_2d_matrix(a1) then return Matrix.logm(a1, a2, true) end
if is_2d_matrix(a2) then
local x = unwrap_num(a1); if x <= 0 then raise_err("ValueError", "math domain error") end
return Matrix.rdiv(math.log(x), Matrix.logm(a2, nil, true))
end
local function scalar_log(v)
if is_decimal_obj(v) or is_decimal_obj(a2) or (type(v) == "table" and v.__type == "int") or (type(a2) == "table" and a2.__type == "int") then
return require("Module:PythonScript/Librarys/Decimal").math_log(v, a2)
end
local x = unwrap_num(v)
if x <= 0 then raise_err("ValueError", "math domain error") end
if a2 ~= nil then return math.log(x) / math.log(unwrap_num(a2)) end
return math.log(x)
end
if is_any_seq(a1) then return Matrix.unary_op(a1, scalar_log) end
return scalar_log(a1)
end
math_api.log10 = function(args)
local a1 = args[1]
if is_2d_matrix(a1) then return Matrix.logm(a1, 10, true) end
local function scalar_log10(v)
if is_decimal_obj(v) or (type(v) == "table" and v.__type == "int") then
return require("Module:PythonScript/Librarys/Decimal").math_log10(v)
end
local x = unwrap_num(v); if x <= 0 then raise_err("ValueError", "math domain error") end
return math.log10(x)
end
if is_any_seq(a1) then return Matrix.unary_op(a1, scalar_log10) end
return scalar_log10(a1)
end
math_api.log2 = function(args)
local a1 = args[1]
if is_2d_matrix(a1) then return Matrix.logm(a1, 2, true) end
local function scalar_log2(v)
if is_decimal_obj(v) or (type(v) == "table" and v.__type == "int") then
return require("Module:PythonScript/Librarys/Decimal").math_log2(v)
end
local x = unwrap_num(v); if x <= 0 then raise_err("ValueError", "math domain error") end
return math.log(x) / math.log(2)
end
if is_any_seq(a1) then return Matrix.unary_op(a1, scalar_log2) end
return scalar_log2(a1)
end
math_api.sin = make_real_dispatcher(Matrix.sinm, function(v) return math.sin(unwrap_num(v)) end)
math_api.cos = make_real_dispatcher(Matrix.cosm, function(v) return math.cos(unwrap_num(v)) end)
math_api.tan = make_real_dispatcher(Matrix.tanm, function(v) return math.tan(unwrap_num(v)) end)
math_api.sinh = make_real_dispatcher(Matrix.sinhm, function(v) local x = unwrap_num(v); return (math.exp(x) - math.exp(-x)) / 2 end)
math_api.cosh = make_real_dispatcher(Matrix.coshm, function(v) local x = unwrap_num(v); return (math.exp(x) + math.exp(-x)) / 2 end)
math_api.tanh = make_real_dispatcher(Matrix.tanhm, function(v) local e = math.exp(2 * unwrap_num(v)); return (e - 1) / (e + 1) end)
math_api.asin = make_real_dispatcher(Matrix.asinm, function(v) return math.asin(unwrap_num(v)) end)
math_api.acos = make_real_dispatcher(Matrix.acosm, function(v) return math.acos(unwrap_num(v)) end)
math_api.atan = make_real_dispatcher(Matrix.atanm, function(v) return math.atan(unwrap_num(v)) end)
math_api.asinh = make_real_dispatcher(Matrix.asinhm, function(v) local x = unwrap_num(v); return math.log(x + math.sqrt(x*x + 1)) end)
math_api.acosh = make_real_dispatcher(Matrix.acoshm, function(v) local x = unwrap_num(v); return math.log(x + math.sqrt(x*x - 1)) end)
math_api.atanh = make_real_dispatcher(Matrix.atanhm, function(v) local x = unwrap_num(v); return 0.5 * math.log((1 + x) / (1 - x)) end)
-- 純量降維與逐元素映射
math_api.fabs = function(args)
if is_any_seq(args[1]) then return Matrix.frobenius_norm(args[1], true) end
return math.abs(unwrap_num(args[1]))
end
math_api.fsum = function(args)
local sum = 0
local function accum(obj)
if is_any_seq(obj) then for _, v in ipairs(obj.elements) do accum(v) end
else sum = sum + unwrap_num(obj) end
end
accum(args[1]); return sum
end
math_api.prod = function(args)
local p = 1
local function accum(obj)
if is_any_seq(obj) then for _, v in ipairs(obj.elements) do accum(v) end
else p = p * unwrap_num(obj) end
end
accum(args[1]); return p
end
math_api.dist = function(args)
local p, q = args[1], args[2]
local diff = Matrix.elementwise_op(p, q, "-")
return Matrix.frobenius_norm(diff, true)
end
math_api.hypot = function(args)
local sum = 0
for _, v in ipairs(args) do
local num = is_any_seq(v) and Matrix.frobenius_norm(v, true) or unwrap_num(v)
sum = sum + num * num
end
return math.sqrt(sum)
end
local function wrap_int_output(val)
if is_bigint_obj(val) then return bi_to_py_or_num(val) end
if type(val) == "number" and not is_nan(val) and not is_inf(val) and math.abs(val) > MAX_SAFE_INT then
return float_int_to_bi(val)
end
return val
end
local function scalar_ceil(x)
if is_bigint_obj(x) then return bi_to_py_or_num(x) end
if is_decimal_obj(x) then return wrap_int_output(require("Module:PythonScript/Librarys/Decimal").math_ceil(x)) end
if is_fraction_obj(x) then return wrap_int_output(require("Module:PythonScript/Librarys/Fractions").math_ceil(x)) end
return wrap_int_output(math.ceil(unwrap_num(x)))
end
local function scalar_floor(x)
if is_bigint_obj(x) then return bi_to_py_or_num(x) end
if is_decimal_obj(x) then return wrap_int_output(require("Module:PythonScript/Librarys/Decimal").math_floor(x)) end
if is_fraction_obj(x) then return wrap_int_output(require("Module:PythonScript/Librarys/Fractions").math_floor(x)) end
return wrap_int_output(math.floor(unwrap_num(x)))
end
local function scalar_trunc(x)
if is_bigint_obj(x) then return bi_to_py_or_num(x) end
if is_decimal_obj(x) then return wrap_int_output(require("Module:PythonScript/Librarys/Decimal").math_trunc(x)) end
if is_fraction_obj(x) then return wrap_int_output(require("Module:PythonScript/Librarys/Fractions").math_trunc(x)) end
local n = unwrap_num(x)
return wrap_int_output(n > 0 and math.floor(n) or math.ceil(n))
end
local function scalar_isqrt(x)
if is_decimal_obj(x) or is_bigint_obj(x) then
return wrap_int_output(require("Module:PythonScript/Librarys/Decimal").math_isqrt(x))
end
local n = unwrap_num(x)
if n < 0 then raise_err("ValueError", "isqrt() argument must be nonnegative") end
if n > MAX_SAFE_INT then
local b = float_int_to_bi(math.floor(n))
return wrap_int_output(require("Module:PythonScript/Librarys/Decimal").math_isqrt(b))
end
return math.floor(math.sqrt(n))
end
local function scalar_erf(x)
local n = unwrap_num(x); local sign = n < 0 and -1 or 1; n = math.abs(n)
local t = 1.0 / (1.0 + 0.3275911 * n)
local y = 1.0 - (((((1.061405429 * t - 1.453152027) * t) + 1.421413741) * t - 0.284496736) * t + 0.254829592) * t * math.exp(-n * n)
return sign * y
end
math_api.ceil = function(args) if is_any_seq(args[1]) then return Matrix.unary_op(args[1], scalar_ceil) end; return scalar_ceil(args[1]) end
math_api.floor = function(args) if is_any_seq(args[1]) then return Matrix.unary_op(args[1], scalar_floor) end; return scalar_floor(args[1]) end
math_api.trunc = function(args) if is_any_seq(args[1]) then return Matrix.unary_op(args[1], scalar_trunc) end; return scalar_trunc(args[1]) end
math_api.isqrt = function(args) if is_any_seq(args[1]) then return Matrix.unary_op(args[1], scalar_isqrt) end; return scalar_isqrt(args[1]) end
math_api.degrees = function(args) if is_any_seq(args[1]) then return Matrix.unary_op(args[1], function(x) return math.deg(unwrap_num(x)) end) end; return math.deg(unwrap_num(args[1])) end
math_api.radians = function(args) if is_any_seq(args[1]) then return Matrix.unary_op(args[1], function(x) return math.rad(unwrap_num(x)) end) end; return math.rad(unwrap_num(args[1])) end
math_api.erf = function(args) if is_any_seq(args[1]) then return Matrix.unary_op(args[1], scalar_erf) end; return scalar_erf(args[1]) end
math_api.erfc = function(args) if is_any_seq(args[1]) then return Matrix.unary_op(args[1], function(x) return 1.0 - scalar_erf(x) end) end; return 1.0 - scalar_erf(args[1]) end
math_api.isinf = function(args)
local function fn_inf(a) if is_decimal_obj(a) then return a.__dict.special == "Inf" end; local x = unwrap_num(a); return x == math.huge or x == -math.huge end
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], fn_inf) end
return fn_inf(args[1])
end
math_api.isnan = function(args)
local function fn_nan(a) if is_decimal_obj(a) then return a.__dict.special == "NaN" or a.__dict.special == "sNaN" end; local x = unwrap_num(a); return x ~= x end
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], fn_nan) end
return fn_nan(args[1])
end
math_api.isfinite = function(args)
local function fn_fin(a) if is_decimal_obj(a) then return a.__dict.special == nil end; local x = unwrap_num(a); return x ~= math.huge and x ~= -math.huge and x == x end
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], fn_fin) end
return fn_fin(args[1])
end
math_api.atan2 = function(args) return math.atan2(unwrap_num(args[1]), unwrap_num(args[2])) end
math_api.fmod = function(args) return math.fmod(unwrap_num(args[1]), unwrap_num(args[2])) end
math_api.frexp = function(args) local m, e = math.frexp(unwrap_num(args[1])); return {__type="tuple", elements={m, e}} end
math_api.ldexp = function(args) return math.ldexp(unwrap_num(args[1]), unwrap_num(args[2])) end
math_api.modf = function(args) local f, i = math.modf(unwrap_num(args[1])); return {__type="tuple", elements={f, i}} end
math_api.copysign = function(args) local x, y = unwrap_num(args[1]), unwrap_num(args[2]); return (y < 0) and -math.abs(x) or math.abs(x) end
math_api.fma = function(args)
local ok1, x_num, x_bi = inspect_int_arg(args[1])
local ok2, y_num, y_bi = inspect_int_arg(args[2])
local ok3, z_num, z_bi = inspect_int_arg(args[3])
if ok1 and ok2 and ok3 then
local prod_est = math.abs(x_num) * math.abs(y_num)
if prod_est > MAX_SAFE_INT or math.abs(z_num) > MAX_SAFE_INT or (prod_est + math.abs(z_num) > MAX_SAFE_INT) then
local res_bi = bi_add(bi_mul(x_bi, y_bi), z_bi)
return bi_to_py_or_num(res_bi)
end
return x_num * y_num + z_num
end
return unwrap_num(args[1]) * unwrap_num(args[2]) + unwrap_num(args[3])
end
math_api.factorial = function(args)
local function scalar_fact(v)
local ok_i, nearest, _ = inspect_int_arg(v)
if ok_i then
if nearest < 0 then raise_err("ValueError", "factorial() not defined for negative values") end
local res = 1
local is_big = false
local res_bi = nil
for idx = 2, nearest do
if not is_big then
if res * idx > MAX_SAFE_INT then
is_big = true
res_bi = bi_mul(small_int_to_bi(res), small_int_to_bi(idx))
else
res = res * idx
end
else
res_bi = bi_mul(res_bi, small_int_to_bi(idx))
end
end
return is_big and bi_to_py_or_num(res_bi) or res
else
return math_api.gamma({ unwrap_num(v) + 1.0 })
end
end
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], scalar_fact) end
return scalar_fact(args[1])
end
math_api.exp2 = function(args)
local function scalar_exp2(v)
local ok_i, x_num, x_bi = inspect_int_arg(v)
if ok_i and x_num >= 53 then
return bi_to_py_or_num(bi_pow_uint(small_int_to_bi(2), x_bi))
end
return math.pow(2, unwrap_num(v))
end
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], scalar_exp2) end
return scalar_exp2(args[1])
end
math_api.ulp = function(args)
local function scalar_ulp(v)
local x = unwrap_num(v)
if is_nan(x) then return 0/0 end
if is_inf(x) then return math.huge end
x = math.abs(x)
if x == 0 then return 5e-324 end -- IEEE 754 雙精度最小非正規數 (2^-1074)
local _, e = math.frexp(x)
local ulp_val = math.ldexp(1, e - 53)
-- 防護：ULP 不可小於最小非正規數
if ulp_val < 5e-324 then return 5e-324 end
return ulp_val
end
if is_any_seq(args[1]) then return Matrix.unary_op(args[1], scalar_ulp) end
return scalar_ulp(args[1])
end
math_api.nextafter = function(args)
local x, y = unwrap_num(args[1]), unwrap_num(args[2])
if is_nan(x) or is_nan(y) then return 0/0 end
if x == y then return y end
if x == 0 then return y > 0 and 5e-324 or -5e-324 end
local abs_x = math.abs(x)
local m, e = math.frexp(abs_x)
-- 判斷是向外擴張 (step up) 還是向 0 縮小 (step down)
local step_up = (x > 0 and y > x) or (x < 0 and y < x)
-- 如果向 0 縮小且 x 剛好是 2 的冪次 (尾數 m == 0.5)，則間距要縮小一半 (e - 54)
local exp_offset = step_up and 53 or ((m == 0.5 and abs_x > 5e-324) and 54 or 53)
local diff = math.ldexp(1, e - exp_offset)
if diff < 5e-324 then diff = 5e-324 end
local res = x + (x > 0 and (step_up and diff or -diff) or (step_up and -diff or diff))
if is_inf(res) then return x > 0 and math.huge or -math.huge end
if math.abs(res) < 5e-324 and res ~= 0 then return 0 end
return res
end
math_api.remainder = function(args)
local x, y = unwrap_num(args[1]), unwrap_num(args[2])
if is_nan(x) or is_nan(y) then return 0/0 end
if y == 0 or is_inf(x) then raise_err("ValueError", "math domain error") end
if is_inf(y) then return x end
local q = x / y
local f = math.floor(q)
local c = math.ceil(q)
local n
-- IEEE 754 捨入到最接近的整數；若剛好在正中間 (x.5)，則捨入到偶數
if c - q == q - f then
n = (f % 2 == 0) and f or c
else
n = math.floor(q + 0.5)
end
return x - n * y
end
local function bi_gcd_core(A, B)
local a = make_bi(1, { unpack(A.limbs) })
local b = make_bi(1, { unpack(B.limbs) })
while not ul_is_zero(b.limbs) do
a, b = b, bi_mod_abs(a, b)
end
return a
end
math_api.gcd = function(args)
if #args == 0 then return 0 end
local g_bi = make_bi(1, { 0 })
for i = 1, #args do
local ok_i, _, x_bi = inspect_int_arg(args[i])
if not ok_i then raise_err("TypeError", "gcd() arguments must be integers") end
g_bi = bi_gcd_core(g_bi, x_bi)
end
return bi_to_py_or_num(g_bi)
end
math_api.lcm = function(args)
if #args == 0 then return 1 end
local l_bi = make_bi(1, { 1 })
for i = 1, #args do
local ok_i, _, x_bi = inspect_int_arg(args[i])
if not ok_i then raise_err("TypeError", "lcm() arguments must be integers") end
if ul_is_zero(x_bi.limbs) then return 0 end
local g_bi = bi_gcd_core(l_bi, x_bi)
local q_bi = bi_divExact(l_bi, g_bi)
l_bi = bi_mul(q_bi, make_bi(1, x_bi.limbs))
end
return bi_to_py_or_num(l_bi)
end
local function apply_neg1_power(exp_is_int, exp_num, exp_bi, val)
if exp_is_int then
local is_odd = (exp_bi.limbs[1] % 2 == 1)
if not is_odd then return val end
if is_bigint_obj(val) then
return make_bi(-val.sign, val.limbs)
else
return -val
end
else
-- 非整數次方：(-1)^exp = cos(pi * exp) + j * sin(pi * exp)
local phase = math.pi * exp_num
local sign_c = Cmath.make_complex(Cmath.clean_zero(math.cos(phase)), Cmath.clean_zero(math.sin(phase)), 0, 0)
return Cmath.mul(sign_c, val)
end
end
-- 安全計算 Gamma(a) / (Gamma(b) * Gamma(c))，自動處理分母極點 (1 / inf = 0) 與大數對數防溢位
local function gamma_comb_ratio(a, b, c)
local function is_nonpos_int(z)
return z <= 0 and math.abs(z - math.floor(z + 0.5)) <= 1e-15
end
-- 若分母落在 Gamma 極點 (非正整數) 且分子不在極點，則極限值為 0
if not is_nonpos_int(a) and (is_nonpos_int(b) or is_nonpos_int(c)) then
return 0.0
end
local ga = Cmath.c2_gamma(a, 0.0)
local gb = Cmath.c2_gamma(b, 0.0)
local gc = Cmath.c2_gamma(c, 0.0)
local denom = gb * gc
if not is_inf(ga) and not is_inf(denom) and denom ~= 0 then
return Cmath.clean_zero(ga / denom)
end
-- 大數防溢位：改用 lgamma 相減並追蹤正負號
local function gamma_sign(z)
if z > 0 then return 1.0 end
return (math.floor(z) % 2 == 0) and -1.0 or 1.0
end
local sgn = gamma_sign(a) * gamma_sign(b) * gamma_sign(c)
local lg = Cmath.r_lgamma(a) - Cmath.r_lgamma(b) - Cmath.r_lgamma(c)
return Cmath.clean_zero(sgn * math.exp(lg))
end
math_api.comb = function(args)
local n_ok, n_num, n_bi = inspect_int_arg(args[1])
local k_ok, k_num, k_bi = inspect_int_arg(args[2])
local both_int = n_ok and k_ok
-- Case 1: k == 0 -> 1
if k_ok and ul_is_zero(k_bi.limbs) then
return 1
end
-- Case 2: (0 <= n < k) or (k < 0 <= n) -> 0
if (n_num >= 0 and n_num < k_num) or (k_num < 0 and 0 <= n_num) then
return 0
end
-- Case 3: (n < 0) and (k > 0)
if n_num < 0 and k_num > 0 then
if both_int then
-- ((-1)^k) * C(|n| + k - 1, k)
local abs_n_bi = make_bi(1, n_bi.limbs)
local top_bi = bi_sub(bi_add(abs_n_bi, k_bi), make_bi(1, { 1 }))
local top_maybe_num = bi_to_py_or_num(top_bi)
local c_val = comb_pos_exact(type(top_maybe_num) ~= "number", type(top_maybe_num) == "number" and top_maybe_num or math.huge, top_bi, k_num)
return apply_neg1_power(true, k_num, k_bi, c_val)
else
-- ((-1)^k) * Gamma(|n| + k) / (Gamma(k + 1) * Gamma(|n|))
local abs_n = math.abs(n_num)
local g_val = gamma_comb_ratio(abs_n + k_num, k_num + 1.0, abs_n)
return apply_neg1_power(k_ok, k_num, k_bi, g_val)
end
end
-- Case 4: (n < 0) and (k < 0)
if n_num < 0 and k_num < 0 then
if both_int then
-- ((-1)^(n+k)) * C(|k| - 1, |n| - 1)
local abs_k_m1_bi = bi_sub(make_bi(1, k_bi.limbs), make_bi(1, { 1 }))
local abs_n_m1_bi = bi_sub(make_bi(1, n_bi.limbs), make_bi(1, { 1 }))
if ul_cmp(abs_k_m1_bi.limbs, abs_n_m1_bi.limbs) < 0 then
return 0
end
local top_maybe_num = bi_to_py_or_num(abs_k_m1_bi)
local bot_maybe_num = bi_to_py_or_num(abs_n_m1_bi)
local c_val = comb_pos_exact(type(top_maybe_num) ~= "number", type(top_maybe_num) == "number" and top_maybe_num or math.huge, abs_k_m1_bi, bot_maybe_num)
local nk_bi = bi_add(n_bi, k_bi)
return apply_neg1_power(true, n_num + k_num, nk_bi, c_val)
else
-- ((-1)^(n+k)) * Gamma(|k|) / (Gamma(|n|) * Gamma(|k| - |n| + 1))
local abs_k, abs_n = math.abs(k_num), math.abs(n_num)
local g_val = gamma_comb_ratio(abs_k, abs_n, abs_k - abs_n + 1.0)
local nk_ok, nk_num, nk_bi = inspect_int_arg(n_num + k_num)
return apply_neg1_power(nk_ok, nk_num, nk_bi, g_val)
end
end
-- Case 5: n 或 k 有一項不是整數 (且 n >= k > 0)
if not both_int then
return gamma_comb_ratio(n_num + 1.0, k_num + 1.0, n_num - k_num + 1.0)
end
-- Case 6: n 與 k 皆為正整數且 n >= k > 0
return comb_pos_exact(math.abs(n_num) > MAX_SAFE_INT, n_num, n_bi, k_num)
end
math_api.perm = function(args)
local n_ok, n_num, n_bi = inspect_int_arg(args[1])
local k_arg = (args[2] ~= nil and args[2] ~= Runtime.builtins.None) and args[2] or args[1]
local k_ok, k_num, _ = inspect_int_arg(k_arg)
if n_ok and k_ok then
if n_num < 0 or k_num < 0 then raise_err("ValueError", "perm() not defined for negative values") end
if k_num > n_num then return 0 end
if k_num == 0 then return 1 end
if math.abs(n_num) <= MAX_SAFE_INT then
local p = 1
local is_big = false
local p_bi = nil
for i = 1, k_num do
local factor = n_num - i + 1
if not is_big then
if p * factor > MAX_SAFE_INT then
is_big = true
p_bi = bi_mul(float_int_to_bi(p), float_int_to_bi(factor))
else
p = p * factor
end
else
p_bi = bi_mul(p_bi, float_int_to_bi(factor))
end
end
return is_big and bi_to_py_or_num(p_bi) or p
else
local p_bi = make_bi(1, { 1 })
local one_bi = make_bi(1, { 1 })
local cur_n = n_bi
for _ = 1, k_num do
p_bi = bi_mul(p_bi, cur_n)
cur_n = bi_sub(cur_n, one_bi)
end
return bi_to_py_or_num(p_bi)
end
end
-- 非整數時使用 Gamma(n + 1) / Gamma(n - k + 1)
return gamma_comb_ratio(n_num + 1.0, 1.0, n_num - k_num + 1.0)
end
math_api.isclose = function(args, kwargs)
local rel_tol = (kwargs and kwargs.rel_tol) or 1e-09; local abs_tol = (kwargs and kwargs.abs_tol) or 0.0
local function check_close(x, y)
if is_any_seq(x) and is_any_seq(y) then
if #x.elements ~= #y.elements then return false end
for i = 1, #x.elements do
if not check_close(x.elements[i], y.elements[i]) then return false end
end
return true
end
local a, b = unwrap_num(x), unwrap_num(y)
return math.abs(a - b) <= math.max(rel_tol * math.max(math.abs(a), math.abs(b)), abs_tol)
end
return check_close(args[1], args[2])
end
local inst = PyClass.make_instance(PyClass.make_class("module", {}, {}))
inst.__dict["__name__"] = "math"
for k, v in pairs(math_api) do
if type(v) == "function" then inst.__dict[k] = { __type = "builtin_method", __call = v } else inst.__dict[k] = v end
end
for k, v in pairs(MathCoreMain.linalg) do
if type(v) == "function" then inst.__dict[k] = { __type = "builtin_method", __call = v } else inst.__dict[k] = v end
end
return inst
end
MathCore.MathLib = MathLib
return MathCore
