-- Module:PythonScript/Runtime/Cmath
-- 整合核心：BigInt 大數引擎、Cmath 四元數/複數引擎、Matrix 廣義超複數泛函線性代數引擎、math 與 cmath 標準庫
local MathCore = {}
-- 延遲載入的例外拋出引擎 (避免與 Runtime 產生循環依賴)
local function raise_err(err_type, msg)
require("Module:PythonScript/Runtime").raise(err_type, msg)
end
-- ==========================================
-- [1] 大數運算引擎 (BigInt Engine)
-- ==========================================
local BASE = 10000000
local SAFE_LIMIT = 9007199254740990
local BigInt = {}
local function strip_limbs(limbs)
while #limbs > 1 and limbs[#limbs] == 0 do table.remove(limbs) end
if #limbs > 10000 then
raise_err("ValueError", "Integer exceeds maximum length limit (10000 limbs)")
end
return limbs
end
local function add_limbs(a, b)
local res = {}; local carry = 0; local n = math.max(#a, #b)
for i = 1, n do
local sum = (a[i] or 0) + (b[i] or 0) + carry
carry = math.floor(sum / BASE)
table.insert(res, sum % BASE)
end
if carry > 0 then table.insert(res, carry) end
return strip_limbs(res)
end
local function sub_limbs(a, b)
local res = {}; local borrow = 0
for i = 1, #a do
local diff = a[i] - (b[i] or 0) - borrow
if diff < 0 then diff = diff + BASE; borrow = 1 else borrow = 0 end
table.insert(res, diff)
end
return strip_limbs(res)
end
local function cmp_limbs(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 mul_limbs(a, b)
if math.min(#a, #b) <= 16 then
local res = {}; for i = 1, #a + #b do res[i] = 0 end
for i = 1, #a do
local carry = 0
for j = 1, #b do
local sum = res[i+j-1] + a[i] * b[j] + carry
carry = math.floor(sum / BASE); res[i+j-1] = sum % BASE
end
res[i+#b] = res[i+#b] + carry
end
return strip_limbs(res)
end
local n = math.max(#a, #b)
local m = math.floor(n / 2)
local low1, high1, low2, high2 = {}, {}, {}, {}
for i = 1, m do low1[i] = a[i] or 0; low2[i] = b[i] or 0 end
for i = m + 1, n do high1[i - m] = a[i] or 0; high2[i - m] = b[i] or 0 end
low1, high1 = strip_limbs(low1), strip_limbs(high1)
low2, high2 = strip_limbs(low2), strip_limbs(high2)
local z0 = mul_limbs(low1, low2)
local z2 = mul_limbs(high1, high2)
local z1 = mul_limbs(add_limbs(low1, high1), add_limbs(low2, high2))
z1 = sub_limbs(sub_limbs(z1, z2), z0)
local res = {}; for i = 1, n * 2 do res[i] = 0 end
for i = 1, #z0 do res[i] = z0[i] end
for i = 1, #z1 do res[i+m] = res[i+m] + z1[i] end
for i = 1, #z2 do res[i+m*2] = res[i+m*2] + z2[i] end
local carry = 0
for i = 1, #res do
local sum = res[i] + carry
carry = math.floor(sum / BASE); res[i] = sum % BASE
end
while carry > 0 do table.insert(res, carry % BASE); carry = math.floor(carry / BASE) end
return strip_limbs(res)
end
local function divmod_limbs(a, b)
if #b == 1 and b[1] == 0 then raise_err("ZeroDivisionError", "division by zero") end
local cmp = cmp_limbs(a, b)
if cmp < 0 then return {0}, { unpack(a) } end
if cmp == 0 then return {1}, {0} end
local q, r = {}, {}
for i = 1, #a do r[i] = a[i] end
if #b == 1 then
local divisor = b[1]; local carry = 0
for i = #r, 1, -1 do
local cur = carry * BASE + r[i]
q[i] = math.floor(cur / divisor); carry = cur % divisor
end
return strip_limbs(q), strip_limbs({carry})
end
-- Knuth Algorithm D 標準正規化：保證 #norm_b == #b 不溢位且最高位 d >= BASE / 2
local scale = math.floor(BASE / (b[#b] + 1))
local norm_b = (scale > 1) and mul_limbs(b, {scale}) or { unpack(b) }
local norm_a = (scale > 1) and mul_limbs(r, {scale}) or r
local n_b = #norm_b
local m_len = #norm_a - n_b
table.insert(norm_a, 0)
local d = norm_b[n_b]
local v_next = norm_b[n_b - 1] or 0
for i = m_len + 1, 1, -1 do
local u2 = norm_a[i + n_b] or 0
local u1 = norm_a[i + n_b - 1] or 0
local u0 = norm_a[i + n_b - 2] or 0
local num_hi = u2 * BASE + u1
local q_guess = math.floor(num_hi / d)
if q_guess >= BASE then q_guess = BASE - 1 end
local r_hat = num_hi - q_guess * d
-- Knuth Step D3 次高位精確修正：消除絕大多數高估誤差
while r_hat < BASE and (q_guess * v_next) > (r_hat * BASE + u0) do
q_guess = q_guess - 1
r_hat = r_hat + d
end
if q_guess > 0 then
local carry = 0
for j = 1, n_b do
local sub_val = norm_b[j] * q_guess + carry
carry = math.floor(sub_val / BASE)
local diff = (norm_a[i+j-1] or 0) - (sub_val % BASE)
if diff < 0 then diff = diff + BASE; carry = carry + 1 end
norm_a[i+j-1] = diff
end
local diff_last = (norm_a[i+n_b] or 0) - carry
while diff_last < 0 do
q_guess = q_guess - 1; carry = 0
for j = 1, n_b do
local sum = (norm_a[i+j-1] or 0) + norm_b[j] + carry
carry = math.floor(sum / BASE); norm_a[i+j-1] = sum % BASE
end
diff_last = diff_last + carry
end
norm_a[i+n_b] = diff_last
end
q[i] = q_guess
end
local rem = (scale > 1) and divmod_limbs(strip_limbs(norm_a), {scale}) or strip_limbs(norm_a)
return strip_limbs(q), rem
end
function BigInt.to_limbs(val)
if type(val) == "table" and val.__type == "int" and val.is_big then
return val.limbs
elseif type(val) == "number" then
if val == 0 then return {0} end; local abs_val = math.abs(val); local limbs = {}
while abs_val > 0 do table.insert(limbs, abs_val % BASE); abs_val = math.floor(abs_val / BASE) end
return limbs
elseif type(val) == "string" then
local s = val:gsub("^-", ""); if s == "" or s == "0" then return {0} end
if #s > 70000 then
raise_err("ValueError", "Integer string exceeds maximum length limit (70000 digits)")
end
local limbs = {}
for i = #s, 1, -7 do local start_idx = math.max(1, i - 6); table.insert(limbs, tonumber(s:sub(start_idx, i))) end
return strip_limbs(limbs)
end; return {0}
end
function BigInt.to_string(sign, limbs)
if #limbs == 1 and limbs[1] == 0 then return "0" end
local res = {}
table.insert(res, tostring(limbs[#limbs]))
for i = #limbs - 1, 1, -1 do table.insert(res, string.format("%07d", limbs[i])) end
return (sign < 0 and "-" or "") .. table.concat(res)
end
local function extract_int_sign_limbs(val)
if type(val) == "boolean" then val = val and 1 or 0 end
if type(val) == "table" and val.__type == "int" then
local l = strip_limbs({ unpack(val.limbs or {0}) })
local s = (#l == 1 and l[1] == 0) and 1 or (val.sign or 1)
return s, l
elseif type(val) == "number" then
if val ~= val or val == math.huge or val == -math.huge or math.floor(val) ~= val then
raise_err("TypeError", "pow() 3rd argument not allowed unless all arguments are integers")
end
if val == 0 then return 1, {0} end
local s = (val < 0) and -1 or 1
local av = math.abs(val)
if av <= SAFE_LIMIT then
return s, BigInt.to_limbs(av)
end
local m, e = math.frexp(av)
local b_limbs = BigInt.to_limbs(math.ldexp(m, 53))
local shift = e - 53
local two_pow = {2}
while shift > 0 do
if shift % 2 == 1 then b_limbs = mul_limbs(b_limbs, two_pow) end
shift = math.floor(shift * 0.5)
if shift > 0 then two_pow = mul_limbs(two_pow, two_pow) end
end
return s, strip_limbs(b_limbs)
end
raise_err("TypeError", "pow() 3rd argument not allowed unless all arguments are integers")
end
-- 大數快速模冪演算法 (Binary Exponentiation with Modulo + 負指數擴展歐幾里得模反元素)
function BigInt.powmod(base, exp, mod)
local b_sign, b_limbs = extract_int_sign_limbs(base)
local e_sign, e_limbs = extract_int_sign_limbs(exp)
local m_sign, m_limbs = extract_int_sign_limbs(mod)
if #m_limbs == 1 and m_limbs[1] == 0 then
raise_err("ValueError", "pow() 3rd argument cannot be 0")
end
if #m_limbs == 1 and m_limbs[1] == 1 then
return 0
end
-- 極速通道：當 |mod| <= 94906265 (平方不超過 2^53-1) 且指數為非負小整數時，直接以原生 53-bit 整數運算
if e_sign > 0 and #m_limbs <= 2 and #e_limbs <= 2 then
local m_val = (m_limbs[2] or 0) * BASE + m_limbs[1]
if m_val <= 94906265 then
local _, r_limbs = divmod_limbs(b_limbs, m_limbs)
local b_val = (r_limbs[2] or 0) * BASE + r_limbs[1]
if b_sign < 0 and b_val ~= 0 then b_val = m_val - b_val end
local e_val = (e_limbs[2] or 0) * BASE + e_limbs[1]
local res = 1
while e_val > 0 do
if e_val % 2 == 1 then res = (res * b_val) % m_val end
e_val = math.floor(e_val * 0.5)
if e_val > 0 then b_val = (b_val * b_val) % m_val end
end
if m_sign < 0 and res > 0 then res = res - m_val end
return res
end
end
-- 1. 先將底數正規化至 [0, |mod|) 區間
local _, rem = divmod_limbs(b_limbs, m_limbs)
if b_sign < 0 and not (#rem == 1 and rem[1] == 0) then
b_limbs = sub_limbs(m_limbs, rem)
else
b_limbs = rem
end
-- 2. 若指數為負數 (Python 3.8+)，透過擴展歐幾里得演算法求模反元素 base^-1 (mod |mod|)
if e_sign < 0 then
if #b_limbs == 1 and b_limbs[1] == 0 then
raise_err("ValueError", "base is not invertible for the given modulus")
end
local old_r, r = m_limbs, b_limbs
local old_t_s, old_t_l = 1, {0}
local t_s, t_l = 1, {1}
while not (#r == 1 and r[1] == 0) do
local q, rem_r = divmod_limbs(old_r, r)
old_r, r = r, rem_r
local qt_l = mul_limbs(q, t_l)
local qt_s = (#qt_l == 1 and qt_l[1] == 0) and 1 or -t_s
local nt_s, nt_l
if old_t_s == qt_s then
nt_s = old_t_s
nt_l = add_limbs(old_t_l, qt_l)
else
local cmp = cmp_limbs(old_t_l, qt_l)
if cmp == 0 then nt_s, nt_l = 1, {0}
elseif cmp > 0 then nt_s, nt_l = old_t_s, sub_limbs(old_t_l, qt_l)
else nt_s, nt_l = qt_s, sub_limbs(qt_l, old_t_l) end
end
old_t_s, old_t_l = t_s, t_l
t_s, t_l = nt_s, nt_l
end
if not (#old_r == 1 and old_r[1] == 1) then
raise_err("ValueError", "base is not invertible for the given modulus")
end
if old_t_s < 0 and not (#old_t_l == 1 and old_t_l[1] == 0) then
b_limbs = sub_limbs(m_limbs, old_t_l)
else
b_limbs = old_t_l
end
end
-- 3. 大數二進位平方取模迴圈 (支援任意千位數以上之超大指數)
local res_limbs = {1}
local cur_base = b_limbs
local cur_exp = e_limbs
while not (#cur_exp == 1 and cur_exp[1] == 0) do
if cur_exp[1] % 2 == 1 then
local prod = mul_limbs(res_limbs, cur_base)
local _, r_mod = divmod_limbs(prod, m_limbs)
res_limbs = r_mod
end
local next_exp = {}
local carry = 0
for i = #cur_exp, 1, -1 do
local cur = carry * BASE + cur_exp[i]
next_exp[i] = math.floor(cur * 0.5)
carry = cur % 2
end
cur_exp = strip_limbs(next_exp)
if not (#cur_exp == 1 and cur_exp[1] == 0) then
local sq = mul_limbs(cur_base, cur_base)
local _, r_sq = divmod_limbs(sq, m_limbs)
cur_base = r_sq
end
end
-- 4. 對齊 Python 負模數餘數正負號規範
local res_sign = 1
if m_sign < 0 and not (#res_limbs == 1 and res_limbs[1] == 0) then
res_limbs = sub_limbs(m_limbs, res_limbs)
res_sign = -1
end
if #res_limbs == 1 and res_limbs[1] == 0 then return 0 end
local num_str = BigInt.to_string(res_sign, res_limbs)
local num = tonumber(num_str)
if num and math.abs(num) <= SAFE_LIMIT then return num end
return { __type = "int", is_big = true, value = num, sign = res_sign, limbs = res_limbs }
end
MathCore.BigInt = {
BASE = BASE,
SAFE_LIMIT = SAFE_LIMIT,
add_limbs = add_limbs, sub_limbs = sub_limbs,
mul_limbs = mul_limbs, divmod_limbs = divmod_limbs, cmp_limbs = cmp_limbs,
to_limbs = BigInt.to_limbs, to_string = BigInt.to_string,
powmod = BigInt.powmod
}
-- ==========================================
-- [2] Python `int` 的專屬擴充方法
-- ==========================================
MathCore.IntMethods = {
bit_length = { __type = "builtin_method", __call = function(args)
local self = args[1]
local v = self.__dict.value
if not self.__dict.is_big and v then
if v == 0 then return 0 end
return math.floor(math.log(math.abs(v)) / math.log(2)) + 1
end
local s = BigInt.to_string(1, self.__dict.limbs)
return math.floor(#s * 3.321928094887362)
end},
bit_count = { __type = "builtin_method", __call = function(args)
local self = args[1]; local v = self.__dict.value
if not self.__dict.is_big and v then
local count = 0; v = math.abs(v)
while v > 0 do count = count + (v % 2); v = math.floor(v / 2) end
return count
end
raise_err("NotImplementedError", "bit_count for BigInt requires full base2 conversion")
end},
to_bytes = { __type = "builtin_method", __call = function(args, kwargs)
local self, length = args[1], args[2] or 1
local byteorder = args[3] or (kwargs and kwargs.byteorder) or "big"
return "<bytes: " .. length .. " bytes, " .. byteorder .. ">"
end},
from_bytes = { __type = "classmethod", __func = { __type = "builtin_method", __call = function()
raise_err("NotImplementedError", "from_bytes requires raw byte parsing")
end}},
as_integer_ratio = { __type = "builtin_method", __call = function(args)
return { __type = "tuple", elements = {args[1], { __type = "int", is_big = false, value = 1, sign = 1 }} }
end},
is_integer = { __type = "builtin_method", __call = function() return true 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 is_matrix_obj(x)
return type(x) == "table" and (x.__type == "list" or x.__type == "tuple") and (
x.is_matrix == true or (#x.elements > 0 and type(x.elements[1]) == "table" and (x.elements[1].__type == "list" or x.elements[1].__type == "tuple"))
)
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(MathCore.BigInt.to_string(x.sign, x.limbs))
elseif x.__type == "complex" then
raise_err("TypeError", "must be real number, not " .. (((x.k and x.k ~= 0) or (x.jk and x.jk ~= 0)) and "quaternion" or "complex"))
elseif is_decimal_obj(x) then
return require("Module:PythonScript/Librarys/Decimal").to_number(x)
elseif x.__type == "instance" then
local PyClass = require("Module:PythonScript/Class")
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
raise_err("TypeError", "must be real number, not " .. (x.__class and x.__class.__name or x.__type or type(x)))
end
if type(x) == "string" then return tonumber(x) end
raise_err("TypeError", "must be real number")
end
-- ==========================================
-- [3] 超複數 / 四元數引擎 (Cmath & Quaternion Engine)
-- Python 表示法：p = a + bj + ck + djk (對應數學 i -> 1j, j -> 1k, k -> 1jk)
-- 內部極輕量結構：q = {a, b, c, d} (4元素陣列，0 雜湊表開銷)
-- ==========================================
local Cmath = {}
local Matrix = {}
function Cmath.make_complex(real, imag, k, jk)
return {
__type = "complex",
real = real or 0,
imag = imag or 0,
k = k or 0,
jk = jk or 0
}
end
Cmath.make_quat = Cmath.make_complex
function Cmath.is_quaternion(x)
return type(x) == "table" and x.__type == "complex" and ((x.k and x.k ~= 0) or (x.jk and x.jk ~= 0))
end
local function get_parts(val)
if type(val) == "table" and val.__type == "complex" then
return val.real or 0, val.imag or 0, val.k or 0, val.jk or 0
end
local n = unwrap_num(val)
if type(n) == "number" then return n, 0, 0, 0 end
raise_err("TypeError", "unsupported operand type(s) for complex/quaternion math")
end
Cmath.get_parts = get_parts
local function to_q4(val, forbid_hyper)
if type(val) == "number" then return { val, 0, 0, 0 } end
if type(val) == "boolean" then return { val and 1 or 0, 0, 0, 0 } end
if type(val) == "table" and val.__type == "complex" then
if forbid_hyper then
raise_err("TypeError", "must be real number, not " .. (Cmath.is_quaternion(val) and "quaternion" or "complex"))
end
return { val.real or 0, val.imag or 0, val.k or 0, val.jk or 0 }
end
return { unwrap_num(val), 0, 0, 0 }
end
local function clean_zero(v)
if type(v) ~= "number" then return v end
if v == 0 or math.abs(v) < 1e-14 then return 0 end
local rounded = math.floor(v + 0.5)
if math.abs(v - rounded) < 1e-13 then return rounded end
return v
end
local function from_q4(q, force_complex)
local a = clean_zero(q[1])
local b = (math.abs(q[2]) < 1e-12) and 0 or clean_zero(q[2])
local c = (math.abs(q[3]) < 1e-12) and 0 or clean_zero(q[3])
local d = (math.abs(q[4]) < 1e-12) and 0 or 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
-- 極輕量 4D 四元數底層算子 {a, b, c, d}
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_conj(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_dot(p, q) return p[1]*q[1] + p[2]*q[2] + p[3]*q[3] + p[4]*q[4] 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
-- 漢密爾頓積 (Hamilton Product): p * q
-- (a*t-b*x-c*y-d*z) + (a*x+b*t+c*z-d*y)j + (a*y-b*z+c*t+d*x)k + (a*z+b*y-c*x+d*t)jk
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
-- 四元數外積: p @ q
-- (t*b-a*x+c*z-d*y)j + (t*c-a*y+d*x-b*z)k + (t*d-a*z+b*y-x*c)jk
local function q_cross(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 {
0,
t*b - a*x + c*z - d*y,
t*c - a*y + d*x - b*z,
t*d - a*z + b*y - x*c
}
end
-- 乘法逆元: p^-1 = conj(p) / dot(p, p)
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
-- 右除 p / q = p * q^-1 與 左除 p \ q = p^-1 * q
local function q_rdiv(p, q) return q_mul(p, q_inv(q)) end
local function q_ldiv(p, q) return q_mul(q_inv(p), q) end
-- 單位號函數 sgn(p) = p / |p|
local function q_sgn(p)
local m = q_abs(p)
if m == 0 then return { 0, 0, 0, 0 } end
return { p[1]/m, p[2]/m, p[3]/m, p[4]/m }
end
-- 向量部單位純虛四元數 u = sgn(vec(p))，若 |vec(p)| == 0 則回傳 u (預設可選 0 或 1j)
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
-- 幅角 arg(p) = arccos(scalar(p) / abs(p)) = atan2(|v|, a)
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
-- exp(p) = exp(a) * (cos(|v|) + sgn(v)*sin(|v|))
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 = ea * math.cos(vn)
local sv = ea * math.sin(vn)
return { cv, u[2]*sv, u[3]*sv, u[4]*sv }
end
-- ln(p) = ln(|p|) + sgn(v)*arg(p)
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
-- sqrt(p) = sqrt(|p|) * (cos(arg(p)/2) + sgn(v)*sin(arg(p)/2))
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
-- p ** q = exp(q * ln(p)) (由左至右約定；整數次方支援精確快速冪)
local function q_pow(p, q)
if q_abs_sq(p) == 0 then
if q_abs_sq(q) == 0 then return { 1, 0, 0, 0 }
elseif q[1] > 0 and q[2] == 0 and q[3] == 0 and q[4] == 0 then return { 0, 0, 0, 0 }
else raise_err("ZeroDivisionError", "0.0 cannot be raised to a negative or complex power") end
end
if q[2] == 0 and q[3] == 0 and q[4] == 0 then
local exp_r = q[1]
if exp_r == -1 then return q_inv(p) end
if exp_r == 0 then return { 1, 0, 0, 0 } end
if math.floor(exp_r) == exp_r and math.abs(exp_r) <= 1024 then
local k = math.floor(math.abs(exp_r))
local base = (exp_r < 0) and q_inv(p) or { p[1], p[2], p[3], p[4] }
local res = { 1, 0, 0, 0 }
while k > 0 do
if k % 2 == 1 then res = q_mul(res, base) end
k = math.floor(k / 2)
if k > 0 then base = q_mul(base, base) end
end
return res
end
end
return q_exp(q_mul(q, q_log(p)))
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
-- sin(p) = sin(a)cosh(|v|) + cos(a)sgn(v)sinh(|v|)
local function q_sin(p)
local a = p[1]
local u, vn = q_vec_unit(p, false)
local sc = math.sin(a) * r_cosh(vn)
local cs = math.cos(a) * r_sinh(vn)
return { sc, u[2]*cs, u[3]*cs, u[4]*cs }
end
-- cos(p) = cos(a)cosh(|v|) - sin(a)sgn(v)sinh(|v|)
local function q_cos(p)
local a = p[1]
local u, vn = q_vec_unit(p, false)
local cc = math.cos(a) * r_cosh(vn)
local ss = -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
-- sinh(p) = sinh(a)cos(|v|) + cosh(a)sgn(v)sin(|v|)
local function q_sinh(p)
local a = p[1]
local u, vn = q_vec_unit(p, false)
local sc = r_sinh(a) * math.cos(vn)
local cs = r_cosh(a) * math.sin(vn)
return { sc, u[2]*cs, u[3]*cs, u[4]*cs }
end
-- cosh(p) = cosh(a)cos(|v|) + sinh(a)sgn(v)sin(|v|)
local function q_cosh(p)
local a = p[1]
local u, vn = q_vec_unit(p, false)
local cc = r_cosh(a) * math.cos(vn)
local ss = 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
-- asinh(p) = ln(p + sqrt(p^2 + 1))
local function q_asinh(p)
local p2 = q_mul(p, p)
local p2_plus_1 = { p2[1] + 1, p2[2], p2[3], p2[4] }
return q_log(q_add(p, q_sqrt(p2_plus_1)))
end
-- acosh(p) = ln(p + sqrt(p + 1)*sqrt(p - 1)) (消除負實軸分支切割誤差，等價於 ln(p + sqrt(p^2 - 1)))
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
-- atanh(p) = (ln(1 + p) - ln(1 - p)) / 2
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
-- asin(p) = -sgn(v) * asinh(p * sgn(v))
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
-- acos(p) = -sgn(v) * acosh(p)
local function q_acos(p)
local u = q_vec_unit(p, true)
return q_mul(q_neg(u), q_acosh(p))
end
-- atan(p) = -sgn(v) * atanh(p * sgn(v))
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
-- =========================================================================
-- 高精度複變與實數 Gamma / Lgamma 引擎 (Godfrey Lanczos g = 607/128, N = 15)
-- =========================================================================
local LANCZOS_G = 4.7421875 -- 607 / 128
local LANCZOS_P = {
0.99999999999999709182,
57.156235665862923517,
-59.597960355475491248,
14.136097974741747174,
-0.49191381609762019978,
0.33994649984811888699e-4,
0.46523628927048575665e-4,
-0.98374475304879564677e-4,
0.15808870322491248884e-3,
-0.21026444172410488319e-3,
0.21743961811521264320e-3,
-0.16431810653676389022e-3,
0.84418223983852743293e-4,
-0.26190838401581408670e-4,
0.36899182659531622704e-5
}
local SQRT_2PI = 2.506628274631000502415765284811
local LN_SQRT_2PI = 0.918938533204672741780329736406
local function is_nan(x) return x ~= x end
-- 🛡️ 四元數防護網：確保輸入僅為實數或標準複數 (a + bj)
local function extract_complex_2d(val, func_name)
local r, i, k, jk = get_parts(val)
if (k and k ~= 0) or (jk and jk ~= 0) then
raise_err("NotImplementedError", func_name .. "() does not support Quaternions; input must be real or complex")
end
return r, i
end
-- 複數二維基礎算子 (專供 Lanczos 高速內迴圈使用)
local function c2_add(ar, ai, br, bi) return ar + br, ai + bi end
local function c2_sub(ar, ai, br, bi) return ar - br, ai - bi end
local function c2_mul(ar, ai, br, bi) return ar*br - ai*bi, ar*bi + ai*br end
local function c2_div(ar, ai, br, bi)
local denom = br*br + bi*bi
if denom == 0 then raise_err("ZeroDivisionError", "complex division by zero in gamma") end
return (ar*br + ai*bi) / denom, (ai*br - ar*bi) / denom
end
local function c2_exp(ar, ai)
local ea = math.exp(ar)
return ea * math.cos(ai), ea * math.sin(ai)
end
local function c2_log(ar, ai)
return 0.5 * math.log(ar*ar + ai*ai), math.atan2(ai, ar)
end
local function c2_pow(ar, ai, br, bi)
local lr, li = c2_log(ar, ai)
local pr, pi = c2_mul(br, bi, lr, li)
return c2_exp(pr, pi)
end
-- 高精度計算 sin(pi * z)，先對實部做模數歸約以防止大數精度流失
local function c2_sin_pi(zr, zi)
local n = math.floor(zr + 0.5)
local dr = zr - n
local sign = (n % 2 == 0) and 1.0 or -1.0
local pr, pi = math.pi * dr, math.pi * zi
return sign * math.sin(pr) * r_cosh(pi), sign * math.cos(pr) * r_sinh(pi)
end
-- 計算 Lanczos 級數和 A_g(z - 1)
local function lanczos_sum(zr, zi)
local xm1_r, xm1_i = zr - 1.0, zi
local sr, si = LANCZOS_P[1], 0.0
for idx = 2, #LANCZOS_P do
local dr, di = xm1_r + (idx - 1), xm1_i
local tr, ti = c2_div(LANCZOS_P[idx], 0.0, dr, di)
sr, si = sr + tr, si + ti
end
return sr, si
end
-- 複變 Gamma 核心：回傳 (real, imag)
local function c2_gamma(zr, zi)
if is_nan(zr) or is_nan(zi) then return 0/0, 0/0 end
-- 極點檢查：非正整數 0, -1, -2, ...
if zi == 0 and zr <= 0 and zr == math.floor(zr) then
raise_err("ValueError", "math domain error")
end
-- 正整數快速通道 (1 <= z <= 23 完全無浮點捨入誤差)
if zi == 0 and zr >= 1 and zr <= 23 and zr == math.floor(zr) then
local prod = 1.0
for k = 2, zr - 1 do prod = prod * k end
return prod, 0.0
end
-- 左半平面 (Re(z) < 0.5)：套用歐拉反射公式 Gamma(z) = pi / (sin(pi * z) * Gamma(1 - z))
if zr < 0.5 then
local g1mr, g1mi = c2_gamma(1.0 - zr, -zi)
local spr, spi = c2_sin_pi(zr, zi)
local den_r, den_i = c2_mul(spr, spi, g1mr, g1mi)
return c2_div(math.pi, 0.0, den_r, den_i)
end
local tr, ti = zr + LANCZOS_G - 0.5, zi
local sr, si = lanczos_sum(zr, zi)
sr, si = c2_mul(SQRT_2PI, 0.0, sr, si)
-- 半指數拆分防溢位：t^((z-0.5)/2) * exp(-t) * t^((z-0.5)/2) * (sqrt(2*pi) * A_g(z))
local half_zr, half_zi = (zr - 0.5) * 0.5, zi * 0.5
local p_half_r, p_half_i = c2_pow(tr, ti, half_zr, half_zi)
local exp_mt_r, exp_mt_i = c2_exp(-tr, -ti)
local res_r, res_i = c2_mul(p_half_r, p_half_i, exp_mt_r, exp_mt_i)
res_r, res_i = c2_mul(res_r, res_i, sr, si)
res_r, res_i = c2_mul(res_r, res_i, p_half_r, p_half_i)
return res_r, res_i
end
-- 複變 Log-Gamma 核心：回傳 (real, imag)
local function c2_lgamma(zr, zi)
if is_nan(zr) or is_nan(zi) then return 0/0, 0/0 end
if zi == 0 and zr <= 0 and zr == math.floor(zr) then
raise_err("ValueError", "math domain error")
end
if zi == 0 and (zr == 1.0 or zr == 2.0) then return 0.0, 0.0 end
if zr < 0.5 then
local lg1mr, lg1mi = c2_lgamma(1.0 - zr, -zi)
local spr, spi = c2_sin_pi(zr, zi)
local lsin_r, lsin_i = c2_log(spr, spi)
return math.log(math.pi) - lsin_r - lg1mr, -lsin_i - lg1mi
end
local tr, ti = zr + LANCZOS_G - 0.5, zi
local lntr, lnti = c2_log(tr, ti)
local term1_r, term1_i = c2_mul(zr - 0.5, zi, lntr, lnti)
local sr, si = lanczos_sum(zr, zi)
local lnsr, lnsi = c2_log(sr, si)
return term1_r - tr + LN_SQRT_2PI + lnsr, term1_i - ti + lnsi
end
-- 實數專用 lgamma(x) = ln(|Gamma(x)|)，完全對齊 CPython math.lgamma
local function r_lgamma(x)
if is_nan(x) then return 0/0 end
if x == math.huge or x == -math.huge then
if x < 0 then raise_err("ValueError", "math domain error") end
return math.huge
end
if x <= 0 and x == math.floor(x) then raise_err("ValueError", "math domain error") end
if x == 1.0 or x == 2.0 then return 0.0 end
if x < 0.5 then
local rem = x - math.floor(x)
local sin_p = math.abs(math.sin(math.pi * rem))
if sin_p == 0 then raise_err("ValueError", "math domain error") end
return math.log(math.pi) - math.log(sin_p) - r_lgamma(1.0 - x)
end
local lgr, _ = c2_lgamma(x, 0.0)
return lgr
end
-- ==========================================
-- Cmath 對外運算介面 (自動路由純量與矩陣)
-- ==========================================
Cmath.clean_zero = clean_zero
Cmath.c2_gamma = c2_gamma
Cmath.c2_lgamma = c2_lgamma
Cmath.r_lgamma = r_lgamma
Cmath.extract_complex_2d = extract_complex_2d
function Cmath.add(c1, c2)
if is_matrix_obj(c1) or is_matrix_obj(c2) then return Matrix.elementwise_op(c1, c2, "+") end
return from_q4(q_add(to_q4(c1, false), to_q4(c2, false)), true)
end
function Cmath.sub(c1, c2)
if is_matrix_obj(c1) or is_matrix_obj(c2) then return Matrix.elementwise_op(c1, c2, "-") end
return from_q4(q_sub(to_q4(c1, false), to_q4(c2, false)), true)
end
function Cmath.mul(c1, c2)
if is_matrix_obj(c1) or is_matrix_obj(c2) then return Matrix.elementwise_op(c1, c2, "*") end
return from_q4(q_mul(to_q4(c1, false), to_q4(c2, false)), true)
end
function Cmath.div(c1, c2)
if is_matrix_obj(c1) or is_matrix_obj(c2) then return Matrix.rdiv(c1, c2) end
return from_q4(q_rdiv(to_q4(c1, false), to_q4(c2, false)), true)
end
function Cmath.ldiv(c1, c2)
if is_matrix_obj(c1) or is_matrix_obj(c2) then return Matrix.ldiv(c1, c2) end
return from_q4(q_ldiv(to_q4(c1, false), to_q4(c2, false)), true)
end
function Cmath.pow(c1, c2)
if is_matrix_obj(c1) or is_matrix_obj(c2) then return Matrix.pow(c1, c2) end
return from_q4(q_pow(to_q4(c1, false), to_q4(c2, false)), true)
end
-- 四元數外積 p @ q
function Cmath.cross(c1, c2)
return from_q4(q_cross(to_q4(c1, false), to_q4(c2, false)), true)
end
Cmath.matmul = Cmath.cross
function Cmath.dot(c1, c2)
return clean_zero(q_dot(to_q4(c1, false), to_q4(c2, false)))
end
function Cmath.conj(c)
return from_q4(q_conj(to_q4(c, false)), true)
end
function Cmath.scalar(c)
local a = get_parts(c)
return clean_zero(a)
end
-- 雙模態 vector：作為屬性讀取時是四元數 0+bj+ck+djk；作為方法呼叫 x.vector() 時傳回 3D 向量 [b, c, d]
function Cmath.vector(c)
local _, b, k_val, jk_val = get_parts(c)
local vec_quat = Cmath.make_complex(0, b, k_val, jk_val)
vec_quat.__call = function()
return { __type = "list", elements = { clean_zero(b), clean_zero(k_val), clean_zero(jk_val) }, is_matrix = true, shape = {3} }
end
return vec_quat
end
function Cmath.sgn(c)
return from_q4(q_sgn(to_q4(c, false)), true)
end
function Cmath.arg(c)
return clean_zero(q_arg(to_q4(c, false)))
end
function Cmath.neg(c) return from_q4(q_neg(to_q4(c, false)), true) end
function Cmath.pos(c) return from_q4(to_q4(c, false), true) end
function Cmath.abs(c) return clean_zero(q_abs(to_q4(c, false))) end
function Cmath.eq(c1, c2)
local s1, a1, b1, k1, d1 = pcall(get_parts, c1)
local s2, a2, b2, k2, d2 = pcall(get_parts, c2)
if not s1 or not s2 then return false end
return a1 == a2 and b1 == b2 and k1 == k2 and d1 == d2
end
function Cmath.to_str(c)
local a, b, k_val, d_val = get_parts(c)
local function fmt(n)
local s = tostring(n)
if s:match("%.0$") then return s:sub(1, -3) end
return s
end
local function signed_fmt(n, unit)
return (n >= 0 and "+" or "") .. fmt(n) .. unit
end
-- 若 k 與 jk 為 0，維持標準 Python complex 字串格式
if k_val == 0 and d_val == 0 then
if a == 0 then return fmt(b) .. "j" end
return "(" .. fmt(a) .. signed_fmt(b, "j") .. ")"
end
-- 單一純虛基元簡寫 (如 2k 或 3jk)
if a == 0 and b == 0 and d_val == 0 then return fmt(k_val) .. "k" end
if a == 0 and b == 0 and k_val == 0 then return fmt(d_val) .. "jk" end
if a == 0 then
return "(" .. fmt(b) .. "j" .. signed_fmt(k_val, "k") .. signed_fmt(d_val, "jk") .. ")"
end
return "(" .. fmt(a) .. signed_fmt(b, "j") .. signed_fmt(k_val, "k") .. signed_fmt(d_val, "jk") .. ")"
end
-- 供 Evaluator / Class 直接查詢 complex / quaternion 的屬性與方法
function Cmath.get_attr(c, attr)
local a, b, k_val, d_val = get_parts(c)
if attr == "real" or attr == "scalar" then return clean_zero(a)
elseif attr == "imag" then return clean_zero(b)
elseif attr == "k" then return clean_zero(k_val)
elseif attr == "jk" then return clean_zero(d_val)
elseif attr == "vector" then return Cmath.vector(c)
elseif attr == "conjugate" then
return { __type = "builtin_method", __name = "conjugate", __call = function() return Cmath.conj(c) end }
elseif attr == "sgn" then
return { __type = "builtin_method", __name = "sgn", __call = function() return Cmath.sgn(c) end }
elseif attr == "arg" then
return { __type = "builtin_method", __name = "arg", __call = function() return Cmath.arg(c) end }
end
return nil
end
MathCore.Cmath = Cmath
-- ==========================================
-- [4] 線性代數與廣義超複數矩陣泛函引擎 (Matrix Functional Calculus)
-- ==========================================
local eval_expr_cached = nil
local dummy_scope_cached = nil
local function py_exact_binop(op, a, b)
if type(a) == "number" and type(b) == "number" then
if op == "+" then return a + b
elseif op == "-" then return a - b
elseif op == "*" then return a * b
elseif op == "/" then
if b == 0 then raise_err("ZeroDivisionError", "division by zero") end
return a / b
elseif op == "//" then
if b == 0 then raise_err("ZeroDivisionError", "integer division or modulo by zero") end
return math.floor(a / b)
elseif op == "%" then
if b == 0 then raise_err("ZeroDivisionError", "integer division or modulo by zero") end
return a % b
elseif op == "**" then return a ^ b end
end
if not eval_expr_cached then
eval_expr_cached = require("Module:PythonScript/Evaluator/Core").getCoreFunction()
dummy_scope_cached = require("Module:PythonScript/Evaluator/Utils").create_scope(nil)
end
return eval_expr_cached({
type = "BinaryOp", op = op,
left = { type = "DirectValue", value = a },
right = { type = "DirectValue", value = b }
}, dummy_scope_cached)
end
local function has_exact_algebraic_elements(obj)
if type(obj) ~= "table" then return false end
if obj.__type == "list" or obj.__type == "tuple" then
for _, el in ipairs(obj.elements) do
if has_exact_algebraic_elements(el) then return true end
end
return false
end
if obj.__type == "complex" then return false end
return true -- int (BigInt), Decimal, Fraction, etc.
end
function Matrix.get_shape(obj)
if type(obj) ~= "table" or (obj.__type ~= "list" and obj.__type ~= "tuple") then return {} end
if #obj.elements == 0 then return {0} end
local sub_shape = Matrix.get_shape(obj.elements[1])
if not sub_shape then return nil end
for i = 2, #obj.elements do
local s = Matrix.get_shape(obj.elements[i])
if not s or #s ~= #sub_shape then return nil end
for j = 1, #s do if s[j] ~= sub_shape[j] then return nil end end
end
local shape = {#obj.elements}
for i = 1, #sub_shape do table.insert(shape, sub_shape[i]) end
return shape
end
function Matrix.to_lua_2d(obj)
local shape = Matrix.get_shape(obj)
if not shape or #shape ~= 2 or shape[1] == 0 or shape[2] == 0 then
raise_err("ValueError", "expected a non-empty 2D matrix")
end
local m, n = shape[1], shape[2]
local mat = {}
for i = 1, m do
mat[i] = {}
local row = obj.elements[i].elements
for j = 1, n do mat[i][j] = unwrap_num(row[j]) end
end
return mat, m, n
end
function Matrix.from_lua_2d(mat, m, n)
local rows = {}
for i = 1, m do
local r = {}
for j = 1, n do r[j] = clean_zero(mat[i][j]) end
rows[i] = { __type = "list", elements = r, is_matrix = true, shape = {n} }
end
return { __type = "list", elements = rows, is_matrix = true, shape = {m, n} }
end
-- 轉為廣義超複數 2D 表 q2d (每個元素為 {a, b, c, d})
function Matrix.to_q2d(obj, forbid_hyper)
local shape = Matrix.get_shape(obj)
if not shape or #shape ~= 2 or shape[1] == 0 or shape[2] == 0 then
raise_err("ValueError", "expected a non-empty 2D matrix")
end
local m, n = shape[1], shape[2]
local mat = {}
for i = 1, m do
mat[i] = {}
local row = obj.elements[i].elements
for j = 1, n do mat[i][j] = to_q4(row[j], forbid_hyper) end
end
return mat, m, n
end
function Matrix.from_q2d(mat, m, n)
local rows = {}
for i = 1, m do
local r = {}
for j = 1, n do
r[j] = from_q4(mat[i][j], false)
end
rows[i] = { __type = "list", elements = r, is_matrix = true, shape = {n} }
end
return { __type = "list", elements = rows, is_matrix = true, shape = {m, n} }
end
function Matrix.as_matrix(obj)
local shape = Matrix.get_shape(obj)
if not shape or #shape == 0 then raise_err("ValueError", "invalid nested list for matrix") end
return Matrix.unary_op(obj, function(x) return x end)
end
function Matrix.elementwise_op(a, b, op)
local is_list_a = (type(a) == "table" and (a.__type == "list" or a.__type == "tuple"))
local is_list_b = (type(b) == "table" and (b.__type == "list" or b.__type == "tuple"))
if is_list_a and is_list_b then
local shapeA = Matrix.get_shape(a)
local shapeB = Matrix.get_shape(b)
if not shapeA or not shapeB then raise_err("ValueError", "invalid shapes for matrix operation") end
if #shapeA == #shapeB then
for i = 1, #shapeA do
if shapeA[i] ~= shapeB[i] then
raise_err("ValueError", "operands could not be broadcast together with shapes")
end
end
local res = { __type = "list", elements = {}, is_matrix = true, shape = shapeA }
for i = 1, #a.elements do
table.insert(res.elements, Matrix.elementwise_op(a.elements[i], b.elements[i], op))
end
return res
elseif #shapeA == 2 and #shapeB == 1 and shapeA[2] == shapeB[1] then
local res = { __type = "list", elements = {}, is_matrix = true, shape = shapeA }
for i = 1, shapeA[1] do
table.insert(res.elements, Matrix.elementwise_op(a.elements[i], b, op))
end
return res
elseif #shapeA == 1 and #shapeB == 2 and shapeA[1] == shapeB[2] then
local res = { __type = "list", elements = {}, is_matrix = true, shape = shapeB }
for i = 1, shapeB[1] do
table.insert(res.elements, Matrix.elementwise_op(a, b.elements[i], op))
end
return res
else
raise_err("ValueError", "operands could not be broadcast together with shapes")
end
elseif is_list_a then
local shapeA = Matrix.get_shape(a)
local res = { __type = "list", elements = {}, is_matrix = true, shape = shapeA }
for i = 1, #a.elements do table.insert(res.elements, Matrix.elementwise_op(a.elements[i], b, op)) end
return res
elseif is_list_b then
local shapeB = Matrix.get_shape(b)
local res = { __type = "list", elements = {}, is_matrix = true, shape = shapeB }
for i = 1, #b.elements do table.insert(res.elements, Matrix.elementwise_op(a, b.elements[i], op)) end
return res
else
return py_exact_binop(op, a, b)
end
end
function Matrix.unary_op(a, fn)
if type(a) == "table" and (a.__type == "list" or a.__type == "tuple") then
local shapeA = Matrix.get_shape(a)
local res = { __type = "list", elements = {}, is_matrix = true, shape = shapeA }
for i = 1, #a.elements do
table.insert(res.elements, Matrix.unary_op(a.elements[i], fn))
end
return res
end
return fn(a)
end
-- ==========================================
-- 內部超複數 2D 矩陣底層算子 (q2d: 元素皆為 {a, b, c, d})
-- ==========================================
local function q2d_eye(n)
local I = {}
for i = 1, n do
I[i] = {}
for j = 1, n do I[i][j] = { (i == j) and 1 or 0, 0, 0, 0 } end
end
return I
end
local function q2d_copy(A, m, n)
local B = {}
for i = 1, m do
B[i] = {}
for j = 1, n do
local cell = A[i][j]
B[i][j] = { cell[1], cell[2], cell[3], cell[4] }
end
end
return B
end
local function q2d_add(A, B, m, n, sign_b)
sign_b = sign_b or 1
local C = {}
for i = 1, m do
C[i] = {}
for j = 1, n do
local a, b = A[i][j], B[i][j]
C[i][j] = { a[1] + sign_b*b[1], a[2] + sign_b*b[2], a[3] + sign_b*b[3], a[4] + sign_b*b[4] }
end
end
return C
end
-- 左乘純量 s * A 與 右乘純量 A * s (在四元數中左右乘順序不同)
local function q2d_scale_left(s_q4, A, m, n)
local C = {}
for i = 1, m do
C[i] = {}
for j = 1, n do C[i][j] = q_mul(s_q4, A[i][j]) end
end
return C
end
local function q2d_scale_right(A, m, n, s_q4)
local C = {}
for i = 1, m do
C[i] = {}
for j = 1, n do C[i][j] = q_mul(A[i][j], s_q4) end
end
return C
end
local function q2d_matmul(A, B, m, k, n)
local C = {}
for i = 1, m do
C[i] = {}
for j = 1, n do
local s1, s2, s3, s4 = 0, 0, 0, 0
for p = 1, k do
local prod = q_mul(A[i][p], B[p][j])
s1 = s1 + prod[1]; s2 = s2 + prod[2]; s3 = s3 + prod[3]; s4 = s4 + prod[4]
end
C[i][j] = { s1, s2, s3, s4 }
end
end
return C
end
local function q2d_frobenius_norm(A, m, n)
local s = 0
for i = 1, m do
for j = 1, n do s = s + q_abs_sq(A[i][j]) end
end
return math.sqrt(s)
end
local function q2d_is_diagonal(A, n)
for i = 1, n do
for j = 1, n do
if i ~= j and q_abs_sq(A[i][j]) > 1e-26 then return false end
end
end
return true
end
local function q2d_has_quaternions(A, m, n)
for i = 1, m do
for j = 1, n do
if math.abs(A[i][j][3]) > 1e-15 or math.abs(A[i][j][4]) > 1e-15 then
return true
end
end
end
return false
end
-- 非交換除環左消去高斯-若爾當法 (Left-Pivoted Gauss-Jordan over Quaternions H)
-- 嚴格使用左乘列運算 E_k，求得精確雙邊反矩陣 A^-1 與行列式 (非交換矩陣回傳 Study/Dieudonne 範數模長)
local function q2d_gauss_jordan(mat, n)
local aug = {}
for i = 1, n do
aug[i] = {}
for j = 1, n do aug[i][j] = { mat[i][j][1], mat[i][j][2], mat[i][j][3], mat[i][j][4] } end
for j = 1, n do aug[i][n + j] = { (i == j) and 1 or 0, 0, 0, 0 } end
end
local det = { 1, 0, 0, 0 }
for col = 1, n do
local pivot_row = col
local max_val = q_abs(aug[col][col])
for r = col + 1, n do
local v = q_abs(aug[r][col])
if v > max_val then max_val = v; pivot_row = r end
end
if max_val <= 1e-11 then
return nil, { 0, 0, 0, 0 }
end
if pivot_row ~= col then
aug[col], aug[pivot_row] = aug[pivot_row], aug[col]
det = q_neg(det)
end
local piv = aug[col][col]
det = q_mul(det, piv)
local inv_piv = q_inv(piv)
-- 嚴格左乘 inv_piv
for j = col, 2 * n do aug[col][j] = q_mul(inv_piv, aug[col][j]) end
for r = 1, n do
if r ~= col then
local factor = aug[r][col]
if q_abs_sq(factor) > 0 then
for j = col, 2 * n do
-- 嚴格左乘 factor: Row_r = Row_r - factor * Row_col
aug[r][j] = q_sub(aug[r][j], q_mul(factor, aug[col][j]))
end
end
end
end
end
local inv = {}
for i = 1, n do
inv[i] = {}
for j = 1, n do inv[i][j] = aug[i][n + j] end
end
return inv, det
end
local function raw_transpose(A, m, n)
local T = {}
for j = 1, n do
T[j] = {}
for i = 1, m do T[j][i] = A[i][j] end
end
return T
end
local function raw_eye(n)
local I = {}
for i = 1, n do
I[i] = {}
for j = 1, n do I[i][j] = (i == j) and 1 or 0 end
end
return I
end
-- Hestenes 單邊雅可比奇異值分解 (One-Sided Jacobi SVD) 與 Moore-Penrose 偽反矩陣
local function raw_svd_and_pinv(mat, m, n)
if m < n then
local matT = raw_transpose(mat, m, n)
local pinvT, U_t, S_t, V_t, rank = raw_svd_and_pinv(matT, n, m)
return raw_transpose(pinvT, m, n), V_t, S_t, U_t, rank
end
local W = {}
for i = 1, m do
W[i] = {}
for j = 1, n do W[i][j] = mat[i][j] end
end
local V = raw_eye(n)
for _ = 1, 40 do
local converged = true
for p = 1, n - 1 do
for q = p + 1, n do
local alpha, beta, gamma = 0, 0, 0
for i = 1, m do
local wp, wq = W[i][p], W[i][q]
alpha = alpha + wp * wp
beta = beta + wq * wq
gamma = gamma + wp * wq
end
if alpha < beta then
for i = 1, m do W[i][p], W[i][q] = W[i][q], W[i][p] end
for i = 1, n do V[i][p], V[i][q] = V[i][q], V[i][p] end
alpha, beta = beta, alpha
end
if math.abs(gamma) > 1e-14 * math.sqrt(alpha * beta) and alpha > 1e-24 then
converged = false
local tau = (beta - alpha) / (2 * gamma)
local sign_tau = (tau >= 0) and 1 or -1
local t = sign_tau / (math.abs(tau) + math.sqrt(1 + tau * tau))
local c = 1 / math.sqrt(1 + t * t)
local s = c * t
for i = 1, m do
local wp, wq = W[i][p], W[i][q]
W[i][p] = c * wp - s * wq
W[i][q] = s * wp + c * wq
end
for i = 1, n do
local vp, vq = V[i][p], V[i][q]
V[i][p] = c * vp - s * vq
V[i][q] = s * vp + c * vq
end
end
end
end
if converged then break end
end
local S = {}
local max_s = 0
for j = 1, n do
local norm_sq = 0
for i = 1, m do norm_sq = norm_sq + W[i][j] * W[i][j] end
S[j] = math.sqrt(norm_sq)
if S[j] > max_s then max_s = S[j] end
end
local tol = math.max(m, n) * max_s * 1e-11
local U = {}
for i = 1, m do U[i] = {} end
local inv_S = {}
local rank = 0
for j = 1, n do
if S[j] > tol then
rank = rank + 1
inv_S[j] = 1 / S[j]
for i = 1, m do U[i][j] = W[i][j] / S[j] end
else
inv_S[j] = 0
for i = 1, m do U[i][j] = 0 end
end
end
local pinv = {}
for i = 1, n do
pinv[i] = {}
for j = 1, m do
local sum = 0
for k = 1, n do
if inv_S[k] ~= 0 then
sum = sum + V[i][k] * inv_S[k] * U[j][k]
end
end
pinv[i][j] = sum
end
end
return pinv, U, S, V, rank
end
local function q2d_inv_or_pinv(mat, m, n)
if m == n then
local inv = q2d_gauss_jordan(mat, n)
if inv then return inv end
end
if q2d_has_quaternions(mat, m, n) then
raise_err("NotImplementedError", "Moore-Penrose pseudo-inverse via SVD is not supported for non-square or singular quaternionic matrices")
end
local rmat = {}
for i = 1, m do
rmat[i] = {}
for j = 1, n do rmat[i][j] = mat[i][j][1] end
end
local rpinv = raw_svd_and_pinv(rmat, m, n)
local qpinv = {}
for i = 1, n do
qpinv[i] = {}
for j = 1, m do qpinv[i][j] = { rpinv[i][j], 0, 0, 0 } end
end
return qpinv
end
-- Denman-Beavers 縮放迭代法求方陣平方根 A^(1/2)
local function q2d_sqrtm(mat, n)
if n == 1 then return { { q_sqrt(mat[1][1]) } } end
if q2d_is_diagonal(mat, n) then
local D = q2d_eye(n)
for i = 1, n do D[i][i] = q_sqrt(mat[i][i]) end
return D
end
local Y = q2d_copy(mat, n, n)
local Z = q2d_eye(n)
for _ = 1, 35 do
local invY, detY = q2d_gauss_jordan(Y, n)
local invZ, detZ = q2d_gauss_jordan(Z, n)
if not invY or not invZ then
invY = q2d_inv_or_pinv(Y, n, n)
invZ = q2d_inv_or_pinv(Z, n, n)
detY, detZ = { 1, 0, 0, 0 }, { 1, 0, 0, 0 }
end
local prod_det = q_abs(q_mul(detY, detZ))
local gamma = 1
if prod_det > 1e-20 then
gamma = prod_det ^ (-1 / (2 * n))
if gamma < 0.1 then gamma = 0.1 elseif gamma > 10 then gamma = 10 end
end
local inv_gamma = 1 / gamma
local diff_sq, norm_sq = 0, 0
for i = 1, n do
for j = 1, n do
local yij, zij = Y[i][j], Z[i][j]
local iz, iy = invZ[i][j], invY[i][j]
local ny = {
0.5 * (gamma * yij[1] + inv_gamma * iz[1]),
0.5 * (gamma * yij[2] + inv_gamma * iz[2]),
0.5 * (gamma * yij[3] + inv_gamma * iz[3]),
0.5 * (gamma * yij[4] + inv_gamma * iz[4])
}
local nz = {
0.5 * (gamma * zij[1] + inv_gamma * iy[1]),
0.5 * (gamma * zij[2] + inv_gamma * iy[2]),
0.5 * (gamma * zij[3] + inv_gamma * iy[3]),
0.5 * (gamma * zij[4] + inv_gamma * iy[4])
}
diff_sq = diff_sq + q_abs_sq(q_sub(ny, yij))
norm_sq = norm_sq + q_abs_sq(ny)
Y[i][j] = ny
Z[i][j] = nz
end
end
if math.sqrt(diff_sq) <= 1e-12 * math.max(1, math.sqrt(norm_sq)) then break end
end
return Y
end
-- 逆縮放平方算法求矩陣主對數 logm(A)
local function q2d_logm(mat, n)
if n == 1 then return { { q_log(mat[1][1]) } } end
if q2d_is_diagonal(mat, n) then
local D = q2d_eye(n)
for i = 1, n do D[i][i] = q_log(mat[i][i]) end
return D
end
local _, det = q2d_gauss_jordan(mat, n)
if q_abs(det) <= 1e-12 then
raise_err("ValueError", "math domain error: matrix logarithm is undefined for singular matrix")
end
local A = q2d_copy(mat, n, n)
local s = 0
local function dist_from_I(M)
local d = 0
for i = 1, n do
for j = 1, n do
local cell = M[i][j]
local dr = cell[1] - ((i == j) and 1 or 0)
d = d + dr*dr + cell[2]*cell[2] + cell[3]*cell[3] + cell[4]*cell[4]
end
end
return math.sqrt(d)
end
while dist_from_I(A) > 0.18 and s < 8 do
A = q2d_sqrtm(A, n)
s = s + 1
end
local X, term, L = {}, {}, {}
for i = 1, n do
X[i], term[i], L[i] = {}, {}, {}
for j = 1, n do
local cell = A[i][j]
local xij = { cell[1] - ((i == j) and 1 or 0), cell[2], cell[3], cell[4] }
X[i][j] = xij
term[i][j] = { xij[1], xij[2], xij[3], xij[4] }
L[i][j] = { xij[1], xij[2], xij[3], xij[4] }
end
end
for k = 2, 20 do
term = q2d_matmul(term, X, n, n, n)
local coeff = ((k % 2 == 1) and 1 or -1) / k
local max_add = 0
for i = 1, n do
for j = 1, n do
local add_q = q_scale(term[i][j], coeff)
L[i][j] = q_add(L[i][j], add_q)
local mag = q_abs(add_q)
if mag > max_add then max_add = mag end
end
end
if max_add < 1e-14 then break end
end
local scale = 2 ^ s
for i = 1, n do
for j = 1, n do L[i][j] = q_scale(L[i][j], scale) end
end
return L
end
-- 縮放平方算法求矩陣指數 expm(A)
local function q2d_expm(mat, n)
if n == 1 then return { { q_exp(mat[1][1]) } } end
if q2d_is_diagonal(mat, n) then
local D = q2d_eye(n)
for i = 1, n do D[i][i] = q_exp(mat[i][i]) end
return D
end
local norm = q2d_frobenius_norm(mat, n, n)
local k = 0
while norm / (2 ^ k) > 0.25 do k = k + 1 end
local scale = 1 / (2 ^ k)
local B = {}
for i = 1, n do
B[i] = {}
for j = 1, n do B[i][j] = q_scale(mat[i][j], scale) end
end
local E = q2d_eye(n)
local term = q2d_eye(n)
for step = 1, 20 do
term = q2d_matmul(term, B, n, n, n)
local max_add = 0
for i = 1, n do
for j = 1, n do
term[i][j] = q_scale(term[i][j], 1 / step)
E[i][j] = q_add(E[i][j], term[i][j])
local mag = q_abs(term[i][j])
if mag > max_add then max_add = mag end
end
end
if max_add < 1e-14 then break end
end
for _ = 1, k do
E = q2d_matmul(E, E, n, n, n)
end
return E
end
-- ==========================================
-- 泛函微積分：三角函數、雙曲函數與反函數 (支援實數、複數與四元數矩陣)
-- ==========================================
local function ensure_square_q2d(obj, func_name, forbid_hyper)
local mA, m, n = Matrix.to_q2d(obj, forbid_hyper)
if m ~= n then
raise_err("ValueError", func_name .. "() requires a square 2D matrix")
end
return mA, n
end
-- 對含有 k, jk 的非交換四元數矩陣使用實係數冪級數計算 sinm 與 cosm (避免 1j 與 1k 不交換問題)
local function q2d_sinm_cosm_series(mA, n)
local norm = q2d_frobenius_norm(mA, n, n)
local k = 0
while norm / (2 ^ k) > 0.5 do k = k + 1 end
local X = q2d_scale_left({ 1 / (2 ^ k), 0, 0, 0 }, mA, n, n)
local X2 = q2d_matmul(X, X, n, n, n)
local S = q2d_copy(X, n, n)
local C = q2d_eye(n)
local termS = q2d_copy(X, n, n)
local termC = q2d_eye(n)
for m = 1, 12 do
termC = q2d_matmul(termC, X2, n, n, n)
termC = q2d_scale_left({ -1 / ((2*m - 1) * (2*m)), 0, 0, 0 }, termC, n, n)
C = q2d_add(C, termC, n, n, 1)
termS = q2d_matmul(termS, X2, n, n, n)
termS = q2d_scale_left({ -1 / ((2*m) * (2*m + 1)), 0, 0, 0 }, termS, n, n)
S = q2d_add(S, termS, n, n, 1)
end
for _ = 1, k do
-- sin(2X) = 2 sin(X) cos(X), cos(2X) = cos^2(X) - sin^2(X)
local nextS = q2d_scale_left({ 2, 0, 0, 0 }, q2d_matmul(S, C, n, n, n), n, n)
local nextC = q2d_add(q2d_matmul(C, C, n, n, n), q2d_matmul(S, S, n, n, n), n, n, -1)
S, C = nextS, nextC
end
return S, C
end
-- sin(M) = (e^(jM) - e^(-jM)) / 2j (對角或四元數矩陣自動派發)
local function q2d_sinm(mA, n)
if q2d_is_diagonal(mA, n) then
local D = q2d_eye(n); for i = 1, n do D[i][i] = q_sin(mA[i][i]) end; return D
end
if q2d_has_quaternions(mA, n, n) then
local S, _ = q2d_sinm_cosm_series(mA, n); return S
end
local Ep = q2d_expm(q2d_scale_left({ 0, 1, 0, 0 }, mA, n, n), n)
local Em = q2d_expm(q2d_scale_left({ 0, -1, 0, 0 }, mA, n, n), n)
return q2d_scale_right(q2d_add(Ep, Em, n, n, -1), n, n, { 0, -0.5, 0, 0 })
end
-- cos(M) = (e^(jM) + e^(-jM)) / 2
local function q2d_cosm(mA, n)
if q2d_is_diagonal(mA, n) then
local D = q2d_eye(n); for i = 1, n do D[i][i] = q_cos(mA[i][i]) end; return D
end
if q2d_has_quaternions(mA, n, n) then
local _, C = q2d_sinm_cosm_series(mA, n); return C
end
local Ep = q2d_expm(q2d_scale_left({ 0, 1, 0, 0 }, mA, n, n), n)
local Em = q2d_expm(q2d_scale_left({ 0, -1, 0, 0 }, mA, n, n), n)
return q2d_scale_right(q2d_add(Ep, Em, n, n, 1), n, n, { 0.5, 0, 0, 0 })
end
-- tan(M) = sin(M) * inv(cos(M))
local function q2d_tanm(mA, n)
if q2d_is_diagonal(mA, n) then
local D = q2d_eye(n); for i = 1, n do D[i][i] = q_tan(mA[i][i]) end; return D
end
local S = q2d_sinm(mA, n)
local C = q2d_cosm(mA, n)
return q2d_matmul(S, q2d_inv_or_pinv(C, n, n), n, n, n)
end
-- sinh(M) = (e^M - e^(-M)) / 2
local function q2d_sinhm(mA, n)
if q2d_is_diagonal(mA, n) then
local D = q2d_eye(n); for i = 1, n do D[i][i] = q_sinh(mA[i][i]) end; return D
end
local Ep = q2d_expm(mA, n)
local Em = q2d_expm(q2d_scale_left({ -1, 0, 0, 0 }, mA, n, n), n)
return q2d_scale_right(q2d_add(Ep, Em, n, n, -1), n, n, { 0.5, 0, 0, 0 })
end
-- cosh(M) = (e^M + e^(-M)) / 2
local function q2d_coshm(mA, n)
if q2d_is_diagonal(mA, n) then
local D = q2d_eye(n); for i = 1, n do D[i][i] = q_cosh(mA[i][i]) end; return D
end
local Ep = q2d_expm(mA, n)
local Em = q2d_expm(q2d_scale_left({ -1, 0, 0, 0 }, mA, n, n), n)
return q2d_scale_right(q2d_add(Ep, Em, n, n, 1), n, n, { 0.5, 0, 0, 0 })
end
-- tanh(M) = sinh(M) * inv(cosh(M))
local function q2d_tanhm(mA, n)
if q2d_is_diagonal(mA, n) then
local D = q2d_eye(n); for i = 1, n do D[i][i] = q_tanh(mA[i][i]) end; return D
end
local S = q2d_sinhm(mA, n)
local C = q2d_coshm(mA, n)
return q2d_matmul(S, q2d_inv_or_pinv(C, n, n), n, n, n)
end
-- asinh(M) = ln(M + sqrt(M^2 + I))
local function q2d_asinhm(mA, n)
if q2d_is_diagonal(mA, n) then
local D = q2d_eye(n); for i = 1, n do D[i][i] = q_asinh(mA[i][i]) end; return D
end
local I = q2d_eye(n)
local M2 = q2d_matmul(mA, mA, n, n, n)
local sq = q2d_sqrtm(q2d_add(M2, I, n, n, 1), n)
return q2d_logm(q2d_add(mA, sq, n, n, 1), n)
end
-- acosh(M) = ln(M + sqrt(M^2 - I))
local function q2d_acoshm(mA, n)
if q2d_is_diagonal(mA, n) then
local D = q2d_eye(n); for i = 1, n do D[i][i] = q_acosh(mA[i][i]) end; return D
end
local I = q2d_eye(n)
local M2 = q2d_matmul(mA, mA, n, n, n)
local sq = q2d_sqrtm(q2d_add(M2, I, n, n, -1), n)
return q2d_logm(q2d_add(mA, sq, n, n, 1), n)
end
-- atanh(M) = 0.5 * (ln(I + M) - ln(I - M))
local function q2d_atanhm(mA, n)
if q2d_is_diagonal(mA, n) then
local D = q2d_eye(n); for i = 1, n do D[i][i] = q_atanh(mA[i][i]) end; return D
end
local I = q2d_eye(n)
local num = q2d_add(I, mA, n, n, 1)
local den = q2d_add(I, mA, n, n, -1)
local q = q2d_matmul(num, q2d_inv_or_pinv(den, n, n), n, n, n)
return q2d_scale_left({ 0.5, 0, 0, 0 }, q2d_logm(q, n), n, n)
end
-- asin(M) = -j * ln(jM + sqrt(I - M^2))
local function q2d_asinm(mA, n)
if q2d_is_diagonal(mA, n) then
local D = q2d_eye(n); for i = 1, n do D[i][i] = q_asin(mA[i][i]) end; return D
end
local I = q2d_eye(n)
local M2 = q2d_matmul(mA, mA, n, n, n)
local sq = q2d_sqrtm(q2d_add(I, M2, n, n, -1), n)
local jM = q2d_scale_left({ 0, 1, 0, 0 }, mA, n, n)
local L = q2d_logm(q2d_add(jM, sq, n, n, 1), n)
return q2d_scale_left({ 0, -1, 0, 0 }, L, n, n)
end
-- acos(M) = (pi/2)*I - asin(M)
local function q2d_acosm(mA, n)
if q2d_is_diagonal(mA, n) then
local D = q2d_eye(n); for i = 1, n do D[i][i] = q_acos(mA[i][i]) end; return D
end
local pi_half_I = q2d_scale_left({ math.pi / 2, 0, 0, 0 }, q2d_eye(n), n, n)
return q2d_add(pi_half_I, q2d_asinm(mA, n), n, n, -1)
end
-- atan(M) = (j/2) * ln((I - jM) * inv(I + jM))
local function q2d_atanm(mA, n)
if q2d_is_diagonal(mA, n) then
local D = q2d_eye(n); for i = 1, n do D[i][i] = q_atan(mA[i][i]) end; return D
end
local I = q2d_eye(n)
local jM = q2d_scale_left({ 0, 1, 0, 0 }, mA, n, n)
local num = q2d_add(I, jM, n, n, -1)
local den = q2d_add(I, jM, n, n, 1)
local q = q2d_matmul(num, q2d_inv_or_pinv(den, n, n), n, n, n)
return q2d_scale_left({ 0, 0.5, 0, 0 }, q2d_logm(q, n), n, n)
end
-- ==========================================
-- Matrix 對外核心方法 (支援雙軌制：精確代數 vs 超複數泛函)
-- ==========================================
function Matrix.matmul(A, B)
-- 若 A 與 B 皆為純量超複數/四元數，轉發至四元數外積 p @ q
if not is_matrix_obj(A) and not is_matrix_obj(B) and (type(A) ~= "table" or (A.__type ~= "list" and A.__type ~= "tuple")) and (type(B) ~= "table" or (B.__type ~= "list" and B.__type ~= "tuple")) then
return Cmath.cross(A, B)
end
local shapeA = Matrix.get_shape(A)
local shapeB = Matrix.get_shape(B)
if not shapeA or not shapeB then raise_err("ValueError", "invalid matrix shape for @") end
local exact = has_exact_algebraic_elements(A) or has_exact_algebraic_elements(B)
if #shapeA == 2 and #shapeB == 2 then
if shapeA[2] ~= shapeB[1] then raise_err("ValueError", "shapes not aligned for matmul") end
local m, k, n = shapeA[1], shapeA[2], shapeB[2]
if exact then
local rows = {}
for i = 1, m do
local r = {}
local rowA = A.elements[i].elements
for j = 1, n do
local sum = py_exact_binop("*", rowA[1], B.elements[1].elements[j])
for p = 2, k do
sum = py_exact_binop("+", sum, py_exact_binop("*", rowA[p], B.elements[p].elements[j]))
end
r[j] = clean_zero(sum)
end
rows[i] = { __type = "list", elements = r, is_matrix = true, shape = {n} }
end
return { __type = "list", elements = rows, is_matrix = true, shape = {m, n} }
else
local mA = Matrix.to_q2d(A, false)
local mB = Matrix.to_q2d(B, false)
return Matrix.from_q2d(q2d_matmul(mA, mB, m, k, n), m, n)
end
elseif #shapeA == 1 and #shapeB == 2 then
if shapeA[1] ~= shapeB[1] then raise_err("ValueError", "shapes not aligned") end
local res = { __type = "list", elements = {}, is_matrix = true, shape = {shapeB[2]} }
for j = 1, shapeB[2] do
local sum = py_exact_binop("*", A.elements[1], B.elements[1].elements[j])
for k = 2, shapeA[1] do
sum = py_exact_binop("+", sum, py_exact_binop("*", A.elements[k], B.elements[k].elements[j]))
end
table.insert(res.elements, clean_zero(sum))
end
return res
elseif #shapeA == 2 and #shapeB == 1 then
if shapeA[2] ~= shapeB[1] then raise_err("ValueError", "shapes not aligned") end
local res = { __type = "list", elements = {}, is_matrix = true, shape = {shapeA[1]} }
for i = 1, shapeA[1] do
local sum = py_exact_binop("*", A.elements[i].elements[1], B.elements[1])
for k = 2, shapeA[2] do
sum = py_exact_binop("+", sum, py_exact_binop("*", A.elements[i].elements[k], B.elements[k]))
end
table.insert(res.elements, clean_zero(sum))
end
return res
elseif #shapeA == 1 and #shapeB == 1 then
if shapeA[1] ~= shapeB[1] then raise_err("ValueError", "shapes not aligned") end
local sum = py_exact_binop("*", A.elements[1], B.elements[1])
for i = 2, shapeA[1] do
sum = py_exact_binop("+", sum, py_exact_binop("*", A.elements[i], B.elements[i]))
end
return clean_zero(sum)
else
raise_err("NotImplementedError", "matmul for >2D or 0D arrays is not supported")
end
end
function Matrix.ldiv(A, B)
local is_list_a = (type(A) == "table" and (A.__type == "list" or A.__type == "tuple"))
local is_list_b = (type(B) == "table" and (B.__type == "list" or B.__type == "tuple"))
if not is_list_a and is_list_b then
return Matrix.elementwise_op(B, A, "/")
elseif is_list_a and not is_list_b then
local shapeA = Matrix.get_shape(A)
if shapeA and #shapeA == 2 then
local mA, m, n = Matrix.to_q2d(A, false)
local invA = Matrix.from_q2d(q2d_inv_or_pinv(mA, m, n), n, m)
if B == 1 then return invA end
return Matrix.elementwise_op(invA, B, "*")
end
return Matrix.elementwise_op(B, A, "/")
elseif is_list_a and is_list_b then
local shapeA = Matrix.get_shape(A)
local shapeB = Matrix.get_shape(B)
if shapeA and #shapeA == 2 then
if not shapeB or (#shapeB ~= 1 and #shapeB ~= 2) then
raise_err("TypeError", "left division (\\) requires a 1D vector or 2D matrix on the right")
end
if shapeA[1] ~= shapeB[1] then
raise_err("ValueError", "matrix dimensions not aligned for left division (\\)")
end
local mA, m, n = Matrix.to_q2d(A, false)
local pinvA = Matrix.from_q2d(q2d_inv_or_pinv(mA, m, n), n, m)
return Matrix.matmul(pinvA, B)
end
return Matrix.elementwise_op(B, A, "/")
end
return py_exact_binop("/", B, A)
end
function Matrix.rdiv(A, B)
local is_list_a = (type(A) == "table" and (A.__type == "list" or A.__type == "tuple"))
local is_list_b = (type(B) == "table" and (B.__type == "list" or B.__type == "tuple"))
if is_list_a and not is_list_b then
return Matrix.elementwise_op(A, B, "/")
elseif not is_list_a and is_list_b then
local shapeB = Matrix.get_shape(B)
if shapeB and #shapeB == 2 then
local mB, m, n = Matrix.to_q2d(B, false)
local invB = Matrix.from_q2d(q2d_inv_or_pinv(mB, m, n), n, m)
if A == 1 then return invB end
return Matrix.elementwise_op(A, invB, "*")
end
return Matrix.elementwise_op(A, B, "/")
elseif is_list_a and is_list_b then
local shapeA = Matrix.get_shape(A)
local shapeB = Matrix.get_shape(B)
if shapeB and #shapeB == 2 then
local colsA = (#shapeA == 2) and shapeA[2] or shapeA[1]
if colsA ~= shapeB[2] then
raise_err("ValueError", "matrix dimensions not aligned for right division (/)")
end
local mB, m, n = Matrix.to_q2d(B, false)
local pinvB = Matrix.from_q2d(q2d_inv_or_pinv(mB, m, n), n, m)
return Matrix.matmul(A, pinvB)
end
return Matrix.elementwise_op(A, B, "/")
end
return py_exact_binop("/", A, B)
end
function Matrix.pow(A, B)
local shapeA = Matrix.get_shape(A)
local shapeB = Matrix.get_shape(B)
local dimA = shapeA and #shapeA or -1
local dimB = shapeB and #shapeB or -1
-- [1] 2D 矩陣/陣列 ** 2D 矩陣/陣列 -> A^B = expm(B @ logm(A))
if dimA == 2 and dimB == 2 then
local mA, m1, n1 = Matrix.to_q2d(A, false)
local mB, m2, n2 = Matrix.to_q2d(B, false)
if m1 ~= n1 then
raise_err("ValueError", "base matrix in M ** N must be a square 2D matrix for logm()")
end
if n2 ~= m1 then
raise_err("ValueError", string.format("shapes (%d,%d) and (%d,%d) not aligned for M ** N = exp(N @ ln(M))", m2, n2, m1, n1))
end
if m2 ~= n1 then
raise_err("ValueError", "exponent matrix N in M ** N must produce a square matrix for expm()")
end
local logA = q2d_logm(mA, n1)
local prod = q2d_matmul(mB, logA, n1, n1, n1)
return Matrix.from_q2d(q2d_expm(prod, n1), n1, n1)
-- [2] 純量 c ** 2D 矩陣 M -> c^M = expm(M * ln(c))
elseif dimA == 0 and dimB == 2 then
local mB, m, n = Matrix.to_q2d(B, false)
if m ~= n then
raise_err("ValueError", "matrix exponent in c ** M requires a square 2D matrix")
end
local c_q4 = to_q4(A, false)
local ln_c = q_log(c_q4)
local scaled = q2d_scale_right(mB, n, n, ln_c)
return Matrix.from_q2d(q2d_expm(scaled, n), n, n)
-- [3] 2D 矩陣 M ** 純量 p
elseif dimA == 2 and dimB == 0 then
local p_q4 = to_q4(B, false)
local is_real_p = (p_q4[2] == 0 and p_q4[3] == 0 and p_q4[4] == 0)
local exp_r = p_q4[1]
local m, n = shapeA[1], shapeA[2]
if m ~= n then
if is_real_p and exp_r == -1 then
local mA = Matrix.to_q2d(A, false)
return Matrix.from_q2d(q2d_inv_or_pinv(mA, m, n), n, m)
end
raise_err("ValueError", "matrix power (other than -1 for pseudo-inverse) requires a square 2D matrix")
end
if is_real_p and exp_r == 0 then
return Matrix.from_lua_2d(raw_eye(n), n, n)
end
-- 軌道一：正整數次方且含精確型別 (Fraction / Decimal / BigInt) 時，使用精確快速冪
if is_real_p and exp_r >= 1 and math.floor(exp_r) == exp_r and has_exact_algebraic_elements(A) then
local k = math.floor(exp_r)
local res = Matrix.from_lua_2d(raw_eye(n), n, n)
local base = Matrix.as_matrix(A)
while k > 0 do
if k % 2 == 1 then res = Matrix.matmul(res, base) end
k = math.floor(k / 2)
if k > 0 then base = Matrix.matmul(base, base) end
end
return res
end
local mA = Matrix.to_q2d(A, false)
if is_real_p and math.floor(exp_r) == exp_r then
local k = math.floor(exp_r)
if k < 0 then
mA = q2d_inv_or_pinv(mA, n, n)
k = -k
end
local res = q2d_eye(n)
local base = mA
while k > 0 do
if k % 2 == 1 then res = q2d_matmul(res, base, n, n, n) end
k = math.floor(k / 2)
if k > 0 then base = q2d_matmul(base, base, n, n, n) end
end
return Matrix.from_q2d(res, n, n)
end
if is_real_p and math.abs(exp_r - 0.5) < 1e-12 then
return Matrix.from_q2d(q2d_sqrtm(mA, n), n, n)
end
-- 一般超複數/小數次方：A^p = expm(p * logm(A))
local L = q2d_logm(mA, n)
return Matrix.from_q2d(q2d_expm(q2d_scale_left(p_q4, L, n, n), n), n, n)
else
return Matrix.elementwise_op(A, B, "**")
end
end
function Matrix.transpose(A)
local shapeA = Matrix.get_shape(A)
if not shapeA or #shapeA == 0 then return A end
if #shapeA == 1 then
return { __type = "list", elements = A.elements, is_matrix = true, shape = shapeA }
elseif #shapeA == 2 then
local m, n = shapeA[1], shapeA[2]
local rows = {}
for j = 1, n do
local r = {}
for i = 1, m do r[i] = A.elements[i].elements[j] end
rows[j] = { __type = "list", elements = r, is_matrix = true, shape = {m} }
end
return { __type = "list", elements = rows, is_matrix = true, shape = {n, m} }
end
raise_err("NotImplementedError", "transpose is only supported for 1D/2D matrices")
end
function Matrix.inv(A)
local mA, m, n = Matrix.to_q2d(A, false)
if m ~= n then raise_err("ValueError", "inv() requires a square matrix; use pinv() or .I for non-square matrices") end
local inv = q2d_gauss_jordan(mA, n)
if not inv then raise_err("ZeroDivisionError", "Singular matrix") end
return Matrix.from_q2d(inv, n, n)
end
function Matrix.pinv(A)
local mA, m, n = Matrix.to_q2d(A, false)
return Matrix.from_q2d(q2d_inv_or_pinv(mA, m, n), n, m)
end
function Matrix.det(A)
local mA, m, n = Matrix.to_q2d(A, false)
if m ~= n then raise_err("ValueError", "det() requires a square 2D matrix") end
local _, d = q2d_gauss_jordan(mA, n)
return from_q4(d, false)
end
function Matrix.trace(A)
local shapeA = Matrix.get_shape(A)
if not shapeA or #shapeA ~= 2 then raise_err("ValueError", "trace() requires a 2D matrix") end
local k = math.min(shapeA[1], shapeA[2])
if k == 0 then return 0 end
local s = A.elements[1].elements[1]
for i = 2, k do
s = py_exact_binop("+", s, A.elements[i].elements[i])
end
return clean_zero(s)
end
function Matrix.rank(A)
local mA, m, n = Matrix.to_lua_2d(A)
local _, _, _, _, r = raw_svd_and_pinv(mA, m, n)
return r
end
function Matrix.svd(A)
local mA, m, n = Matrix.to_lua_2d(A)
local _, U, S, V = raw_svd_and_pinv(mA, m, n)
local k = math.min(m, n)
local s_elts = {}
for i = 1, k do s_elts[i] = clean_zero(S[i]) end
return {
__type = "tuple",
elements = {
Matrix.from_lua_2d(U, m, k),
{ __type = "list", elements = s_elts, is_matrix = true, shape = {k} },
Matrix.from_lua_2d(raw_transpose(V, n, k), k, n)
}
}
end
function Matrix.frobenius_norm(A, forbid_hyper)
local sum_sq = 0
local function accum(node)
if type(node) == "table" and (node.__type == "list" or node.__type == "tuple") then
for _, el in ipairs(node.elements) do accum(el) end
else
sum_sq = sum_sq + q_abs_sq(to_q4(node, forbid_hyper))
end
end
accum(A)
return clean_zero(math.sqrt(sum_sq))
end
function Matrix.expm(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "exp", forbid_hyper)
return Matrix.from_q2d(q2d_expm(mA, n), n, n)
end
function Matrix.expm1(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "expm1", forbid_hyper)
local E = q2d_expm(mA, n)
return Matrix.from_q2d(q2d_add(E, q2d_eye(n), n, n, -1), n, n)
end
function Matrix.logm(A, base, forbid_hyper)
local mA, n = ensure_square_q2d(A, "log", forbid_hyper)
local LA = Matrix.from_q2d(q2d_logm(mA, n), n, n)
if base == nil then return LA end
if type(base) == "table" and (base.__type == "list" or base.__type == "tuple") then
local mB, nb = ensure_square_q2d(base, "log", forbid_hyper)
local LB = Matrix.from_q2d(q2d_logm(mB, nb), nb, nb)
return Matrix.rdiv(LA, LB)
else
local bp = to_q4(base, forbid_hyper)
local lnb = from_q4(q_log(bp), false)
return Matrix.rdiv(LA, lnb)
end
end
function Matrix.log1p(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "log1p", forbid_hyper)
local IpA = q2d_add(q2d_eye(n), mA, n, n, 1)
return Matrix.from_q2d(q2d_logm(IpA, n), n, n)
end
function Matrix.sqrtm(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "sqrt", forbid_hyper)
return Matrix.from_q2d(q2d_sqrtm(mA, n), n, n)
end
function Matrix.cbrtm(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "cbrt", forbid_hyper)
local L = q2d_logm(mA, n)
return Matrix.from_q2d(q2d_expm(q2d_scale_left({ 1/3, 0, 0, 0 }, L, n, n), n), n, n)
end
function Matrix.sinm(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "sin", forbid_hyper)
return Matrix.from_q2d(q2d_sinm(mA, n), n, n)
end
function Matrix.cosm(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "cos", forbid_hyper)
return Matrix.from_q2d(q2d_cosm(mA, n), n, n)
end
function Matrix.tanm(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "tan", forbid_hyper)
return Matrix.from_q2d(q2d_tanm(mA, n), n, n)
end
function Matrix.sinhm(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "sinh", forbid_hyper)
return Matrix.from_q2d(q2d_sinhm(mA, n), n, n)
end
function Matrix.coshm(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "cosh", forbid_hyper)
return Matrix.from_q2d(q2d_coshm(mA, n), n, n)
end
function Matrix.tanhm(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "tanh", forbid_hyper)
return Matrix.from_q2d(q2d_tanhm(mA, n), n, n)
end
function Matrix.asinm(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "asin", forbid_hyper)
return Matrix.from_q2d(q2d_asinm(mA, n), n, n)
end
function Matrix.acosm(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "acos", forbid_hyper)
return Matrix.from_q2d(q2d_acosm(mA, n), n, n)
end
function Matrix.atanm(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "atan", forbid_hyper)
return Matrix.from_q2d(q2d_atanm(mA, n), n, n)
end
function Matrix.asinhm(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "asinh", forbid_hyper)
return Matrix.from_q2d(q2d_asinhm(mA, n), n, n)
end
function Matrix.acoshm(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "acosh", forbid_hyper)
return Matrix.from_q2d(q2d_acoshm(mA, n), n, n)
end
function Matrix.atanhm(A, forbid_hyper)
local mA, n = ensure_square_q2d(A, "atanh", forbid_hyper)
return Matrix.from_q2d(q2d_atanhm(mA, n), n, n)
end
MathCore.Matrix = Matrix
MathCore.linalg = {
matrix = function(args) return Matrix.as_matrix(args[1]) end,
inv = function(args) return Matrix.inv(args[1]) end,
pinv = function(args) return Matrix.pinv(args[1]) end,
det = function(args) return Matrix.det(args[1]) end,
trace = function(args) return Matrix.trace(args[1]) end,
rank = function(args) return Matrix.rank(args[1]) end,
svd = function(args) return Matrix.svd(args[1]) end,
norm = function(args) return Matrix.frobenius_norm(args[1], false) end,
transpose = function(args) return Matrix.transpose(args[1]) end,
solve = function(args) return Matrix.ldiv(args[1], args[2]) end,
mpower = function(args) return Matrix.pow(args[1], args[2]) end,
sqrtm = function(args) return Matrix.sqrtm(args[1], false) end,
expm = function(args) return Matrix.expm(args[1], false) end,
logm = function(args) return Matrix.logm(args[1], args[2], false) end,
sinm = function(args) return Matrix.sinm(args[1], false) end,
cosm = function(args) return Matrix.cosm(args[1], false) end,
tanm = function(args) return Matrix.tanm(args[1], false) end,
sinhm = function(args) return Matrix.sinhm(args[1], false) end,
coshm = function(args) return Matrix.coshm(args[1], false) end,
tanhm = function(args) return Matrix.tanhm(args[1], false) end,
asinm = function(args) return Matrix.asinm(args[1], false) end,
acosm = function(args) return Matrix.acosm(args[1], false) end,
atanm = function(args) return Matrix.atanm(args[1], false) end,
asinhm = function(args) return Matrix.asinhm(args[1], false) end,
acoshm = function(args) return Matrix.acoshm(args[1], false) end,
atanhm = function(args) return Matrix.atanhm(args[1], false) end,
dot = function(args) return Cmath.dot(args[1], args[2]) end,
cross = function(args) return Cmath.cross(args[1], args[2]) end,
eye = function(args)
local dim = require("Module:PythonScript/Runtime/Utils").unwrap_num(args[1])
return Matrix.from_lua_2d(raw_eye(dim), dim, dim)
end
}
function MathCore.get_math_module()
return require("Module:PythonScript/Librarys/Cmath").MathLib.init()
end
function MathCore.get_cmath_module()
return require("Module:PythonScript/Librarys/Cmath").CmathLib.init()
end
return MathCore
