diff options
Diffstat (limited to 'macros/luatex/latex/luacas/tex/algebra')
18 files changed, 4469 insertions, 0 deletions
diff --git a/macros/luatex/latex/luacas/tex/algebra/_init.lua b/macros/luatex/latex/luacas/tex/algebra/_init.lua new file mode 100644 index 0000000000..0a9b114cb9 --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/_init.lua @@ -0,0 +1,24 @@ +-- Loads algebra files in the correct order. +require("_lib.table") + +require("core._init") + +require("algebra.ring") +require("algebra.euclideandomain") +require("algebra.field") +require("algebra.polynomialring") +require("algebra.integer") +require("algebra.rational") +require("algebra.integerquotientring") +require("algebra.sqrtexpression") + +require("algebra.absexpression") +require("algebra.equation") +require("algebra.factorialexpression") +require("algebra.logarithm") +require("algebra.rootexpression") +require("algebra.trigexpression") + +require("algebra.polynomialring.berlekampfactoring") +require("algebra.polynomialring.zassenhausfactoring") +require("algebra.polynomialring.decomposition")
\ No newline at end of file diff --git a/macros/luatex/latex/luacas/tex/algebra/absexpression.lua b/macros/luatex/latex/luacas/tex/algebra/absexpression.lua new file mode 100644 index 0000000000..2bf0b6a7ef --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/absexpression.lua @@ -0,0 +1,80 @@ +--- @class AbsExpression +--- The absolute value of an expression. +--- @field expression Expression +AbsExpression = {} +__AbsExpression = {} + +---------------------------- +-- Instance functionality -- +---------------------------- + +--- Creates a new absolute value expression with the given expression. +--- @param expression Expression +--- @return AbsExpression +function AbsExpression:new(expression) + local o = {} + local __o = Copy(__ExpressionOperations) + + o.expression = expression + + __o.__index = AbsExpression + __o.__tostring = function(a) + return '|' .. tostring(a.expression) .. '|' + end + + o = setmetatable(o, __o) + return o +end + +--- @return Expression +function AbsExpression:evaluate() + if self.expression:isconstant() then + if self.expression >= Integer.zero() then + return self.expression + end + return -self.expression + end + return self +end + +--- @return Expression +function AbsExpression:autosimplify() + return AbsExpression(self.expression:autosimplify()):evaluate() +end + +--- @return table<number, Expression> +function AbsExpression:subexpressions() + return {self.expression} +end + +--- @param subexpressions table<number, Expression> +--- @return AbsExpression +function AbsExpression:setsubexpressions(subexpressions) + return AbsExpression(subexpressions[1]) +end + +--- @param other Expression +--- @return boolean +function AbsExpression:order(other) + return FunctionExpression("abs", self.expression):order(other) +end + +--- @return string +function AbsExpression:tolatex() + return "\\left|" .. self.expression:tolatex() .. "\\right|" +end + +----------------- +-- Inheritance -- +----------------- + +__AbsExpression.__index = CompoundExpression +__AbsExpression.__call = AbsExpression.new +AbsExpression = setmetatable(AbsExpression, __AbsExpression) + +---------------------- +-- Static constants -- +---------------------- +ABS = function (a) + return AbsExpression(a) +end
\ No newline at end of file diff --git a/macros/luatex/latex/luacas/tex/algebra/equation.lua b/macros/luatex/latex/luacas/tex/algebra/equation.lua new file mode 100644 index 0000000000..4e81d5ca95 --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/equation.lua @@ -0,0 +1,183 @@ +--- @class Equation +--- An expression that represents an equation of the form lhs = rhs. +--- @field lhs Expression +--- @field rhs Expression +Equation = {} +__Equation = {} + +-------------------------- +-- Static functionality -- +-------------------------- + +--- Attempts to isolate the variable var in lhs by moving expressions to rhs. Ony performs a single step. +--- @param lhs Expression +--- @param rhs Expression +--- @param var SymbolExpression +--- @return Expression, Expression +function Equation.isolatelhs(lhs, rhs, var) + if lhs:type() == BinaryOperation then + local stay = Integer.zero() + local switch = Integer.zero() + if lhs.operation == BinaryOperation.ADD then + for _, exp in ipairs(lhs:subexpressions()) do + if exp:freeof(var) then + switch = switch + exp + else + stay = stay + exp + end + end + if switch == Integer.zero() then + lhs = lhs:factor() -- TODO: Replace with collect for efficiency reasons + else + return stay:autosimplify(), (rhs - switch):autosimplify() + end + end + if lhs.operation == BinaryOperation.MUL then + stay = Integer.one() + switch = Integer.one() + for _, exp in ipairs(lhs:subexpressions()) do + if exp:freeof(var) then + switch = switch * exp + else + stay = stay * exp + end + end + return stay:autosimplify(), (rhs / switch):autosimplify() + end + if lhs.operation == BinaryOperation.POW then + if lhs:subexpressions()[1]:freeof(var) then + return lhs:subexpressions()[2]:autosimplify(), Logarithm(lhs:subexpressions()[1], rhs):autosimplify() + elseif lhs:subexpressions()[2]:freeof(var) then + return lhs:subexpressions()[1]:autosimplify(), (rhs ^ (Integer.one()/lhs:subexpressions()[2])):autosimplify() + end + end + elseif lhs:type() == Logarithm then + if lhs.base:freeof(var) then + return lhs.expression:autosimplify(), (lhs.base ^ rhs):autosimplify() + elseif lhs.expression:freeof(var) then + return lhs.base:autosimplify(), (lhs.expression ^ (Integer.one()/rhs)):autosimplify() + end + elseif lhs:type() == TrigExpression then + return lhs.expression:autosimplify(), TrigExpression(TrigExpression.INVERSES[lhs.name], rhs):autosimplify() + end + + return lhs, rhs +end + +---------------------------- +-- Instance functionality -- +---------------------------- + +--- Creates a new equation with the given expressions. +--- @param lhs Expression +--- @param rhs Expression +--- @return Equation +function Equation:new(lhs, rhs) + + if lhs:type() == Equation or rhs:type() == Equation then + error("Sent parameter of wrong type: cannot nest equations or inequalities") + end + + local o = {} + local __o = Copy(__ExpressionOperations) -- TODO: Ensure only one metatable for each instance of a class + + o.lhs = lhs + o.rhs = rhs + + __o.__index = Equation + __o.__tostring = function(a) + return tostring(a.lhs) .. ' = ' .. tostring(a.rhs) + end + __o.__eq = function(a, b) + -- This shouldn't be needed, since __eq should only fire if both metamethods have the same function, but for some reason Lua always runs this anyway + if not b:type() == Equation then + return false + end + return a.lhs == b.lhs and a.rhs == b.rhs + end + o = setmetatable(o, __o) + + return o +end + +--- Evaluation in this case just checks for structural equality, or guarenteed inequality in the case of constants +--- @return Equation|boolean +function Equation:evaluate() + if self.lhs == self.rhs then + return true -- TODO: Add Boolean Expressions + end + if self.lhs:isconstant() and self.rhs:isconstant() and self.lhs ~= self.rhs then + return false + end + return self +end + +--- @return Equation|boolean +function Equation:autosimplify() + local lhs = self.lhs:autosimplify() + local rhs = self.rhs:autosimplify() + + return Equation(lhs, rhs):evaluate() +end + +--- @return table<number, Expression> +function Equation:subexpressions() + return {self.lhs, self.rhs} +end + +--- Attempts to solve the equation for a particular variable. +--- @param var SymbolExpression +--- @return Equation +function Equation:solvefor(var) + local lhs = self.lhs + local rhs = self.rhs + + if lhs:freeof(var) and rhs:freeof(var) then + return self + end + + -- Check for monovariate polynomial expressions + local root = (lhs - rhs):autosimplify() + local poly, status = root:expand():topolynomial() + if status then + -- TODO: Add Set expressions + return Equation(var, poly:roots()[1]) + end + + local newlhs, newrhs = root, Integer(0) + local oldlhs + while newlhs ~= var and oldlhs ~= newlhs do + oldlhs = newlhs + newlhs, newrhs = Equation.isolatelhs(newlhs, newrhs, var) + end + + return Equation(newlhs, newrhs) +end + +--- @param subexpressions table<number, Expression> +--- @return Equation +function Equation:setsubexpressions(subexpressions) + return Equation(subexpressions[1], subexpressions[2]) +end + +--- @param other Expression +--- @return boolean +function Equation:order(other) + if other:isatomic() then + return false + end + + return self.lhs:order(other) +end + +--- @return string +function Equation:tolatex() + return self.lhs:tolatex() .. '=' .. self.rhs:tolatex() +end + +----------------- +-- Inheritance -- +----------------- +__Equation.__index = CompoundExpression +__Equation.__call = Equation.new +Equation = setmetatable(Equation, __Equation)
\ No newline at end of file diff --git a/macros/luatex/latex/luacas/tex/algebra/euclideandomain.lua b/macros/luatex/latex/luacas/tex/algebra/euclideandomain.lua new file mode 100644 index 0000000000..fab2c5c7c2 --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/euclideandomain.lua @@ -0,0 +1,63 @@ +--- @class EuclideanDomain +--- Interface for an element of a euclidean domain. +EuclideanDomain = {} +__EuclideanDomain = {} + +---------------------- +-- Required methods -- +---------------------- + +--- @param b EuclideanDomain +--- @return EuclideanDomain, EuclideanDomain +function EuclideanDomain:divremainder(b) + error("Called unimplemented method : divremainder()") +end + +---------------------------- +-- Instance functionality -- +---------------------------- + +--- @return boolean +function EuclideanDomain:iscommutative() + return true +end + +-------------------------- +-- Instance metamethods -- +-------------------------- + +__EuclideanOperations = Copy(__RingOperations) + +-- Division with remainder +-- Unfortunately, this can only return 1 result, so it returns the quotient - for the remainder use a % b, or a:divremainder(b) +__EuclideanOperations.__idiv = function(a, b) + if(b == b:zero()) then + error("Cannot divide by zero.") + end + local aring, bring = a:getring(), b:getring() + local oring = Ring.resultantring(aring, bring) + if not oring then + error("Attempted to divide two elements of incompatable rings") + end + return a:inring(oring):divremainder(b:inring(oring)) +end + +__EuclideanOperations.__mod = function(a, b) + if(b == b:zero()) then + error("Cannot divide by zero.") + end + local aring, bring = a:getring(), b:getring() + local oring = Ring.resultantring(aring, bring) + if not oring then + error("Attempted to divide two elements of incompatable rings") + end + local _,q = a:inring(oring):divremainder(b:inring(oring)) + return q +end + +----------------- +-- Inheritance -- +----------------- + +__EuclideanDomain.__index = Ring +EuclideanDomain = setmetatable(EuclideanDomain, __EuclideanDomain)
\ No newline at end of file diff --git a/macros/luatex/latex/luacas/tex/algebra/factorialexpression.lua b/macros/luatex/latex/luacas/tex/algebra/factorialexpression.lua new file mode 100644 index 0000000000..1c5b824ff7 --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/factorialexpression.lua @@ -0,0 +1,102 @@ +--- @class FactorialExpression +--- The factorial of an expression. +--- @field expression Expression +FactorialExpression = {} +__FactorialExpression = {} + +---------------------------- +-- Instance functionality -- +---------------------------- + +--- Creates a new factorial expression with the given expression. +--- @param expression Expression +--- @return FactorialExpression +function FactorialExpression:new(expression) + local o = {} + local __o = Copy(__ExpressionOperations) + + o.expression = expression + + __o.__index = FactorialExpression + __o.__tostring = function(a) + return '(' .. tostring(a.expression) .. ')!' + end + + o = setmetatable(o, __o) + return o +end + +--- @return Expression +function FactorialExpression:evaluate() + if self.expression:type() == Integer then + if self.expression < Integer.zero() then + error("Aritmetic Error: Factorials of negative integers are not defined.") + end + + if not FactorialExpression.LIMIT then + FactorialExpression.LIMIT = Integer(5000) + end + + if self.expression > FactorialExpression.LIMIT then + return self + end + -- TODO: More efficient factorial computations. + local out = Integer.one() + local i = Integer.zero() + while i < self.expression do + i = i + Integer.one() + out = out * i + end + return out + end + return self +end + +--- @return Expression +function FactorialExpression:autosimplify() + return FactorialExpression(self.expression:autosimplify()):evaluate() +end + +--- @return table<number, Expression> +function FactorialExpression:subexpressions() + return {self.expression} +end + +--- @param subexpressions table<number, Expression> +--- @return AbsExpression +function FactorialExpression:setsubexpressions(subexpressions) + return FactorialExpression(subexpressions[1]) +end + +--- @param other Expression +--- @return boolean +function FactorialExpression:order(other) + return FunctionExpression("fact", self.expression):order(other) +end + +--- @return string +function FactorialExpression:tolatex() + if self.expression:isatomic() then + return self.expression:tolatex() .. "!" + end + return "(" .. self.expression:tolatex() .. ")!" +end + +----------------- +-- Inheritance -- +----------------- + +__FactorialExpression.__index = CompoundExpression +__FactorialExpression.__call = FactorialExpression.new +FactorialExpression = setmetatable(FactorialExpression, __FactorialExpression) + +---------------------- +-- Static constants -- +---------------------- + +-- Do not attempt to compute factorials larger than this. +FactorialExpression.LIMIT = Integer(5000) + +FACT = function (a) + return FactorialExpression(a) +end
\ No newline at end of file diff --git a/macros/luatex/latex/luacas/tex/algebra/field.lua b/macros/luatex/latex/luacas/tex/algebra/field.lua new file mode 100644 index 0000000000..c845e94c55 --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/field.lua @@ -0,0 +1,69 @@ +--- @class Field +--- Interface for an element of a field. +Field = {} +__Field = {} + +---------------------- +-- Required methods -- +---------------------- + +--- @return Field +function Field:div(b) + return self:mul(b:inv()) +end + +--- @return Field +function Field:inv() + error("Called unimplemented method: inv()") +end + +---------------------------- +-- Instance functionality -- +---------------------------- + +--- Field exponentiation based on the definition. Specific rings may implement more efficient methods. +--- @return Field +function Field:pow(n) + local base = self + if(n < Integer.zero()) then + n = -n + base = base:inv() + end + local k = Integer.zero() + local b = self.getring().one() + while k < n do + b = b.mul(base) + k = k + Integer.one() + end + return b +end + +-------------------------- +-- Instance metamethods -- +-------------------------- + +__FieldOperations = Copy(__EuclideanOperations) + +__FieldOperations.__div = function(a, b) + if not b.getring and not b:isconstant() then + return BinaryOperation.DIVEXP({a, b}) + end + + if(b == b:zero()) then + error("Arithmetic Error: Cannot divide by zero.") + end + + local aring, bring = a:getring(), b:getring() + local oring = Ring.resultantring(aring, bring) + if not oring then + error("Attempted to divide two elements of incompatable rings") + end + return a:inring(oring):div(b:inring(oring)) +end + +----------------- +-- Inheritance -- +----------------- + +__Field.__index = EuclideanDomain +Field = setmetatable(Field, __Field)
\ No newline at end of file diff --git a/macros/luatex/latex/luacas/tex/algebra/integer.lua b/macros/luatex/latex/luacas/tex/algebra/integer.lua new file mode 100644 index 0000000000..a3c54f498c --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/integer.lua @@ -0,0 +1,952 @@ +--- @class Integer +--- Represents an element of the ring of integers. +--- @field self table<number, number> +--- @field sign number +Integer = {} +__Integer = {} + +-------------------------- +-- Static functionality -- +-------------------------- + +-- The length of each digit in base 10. 10^15 < 2^53 < 10^16, so 15 is the highest value that will work with double-percision numbers. +-- For multiplication to work properly, however, this also must be even so we can take the square root of the digit size exactly. +-- 10^14 is still larger than 2^26, so it is still efficient to do multiplication this way. +Integer.DIGITLENGTH = 14 +-- The maximum size for a digit. While this doesn't need to be a power of 10, it makes implementing converting to and from strings much easier. +Integer.DIGITSIZE = 10 ^ Integer.DIGITLENGTH +-- Partition size for multiplying integers so we can get both the upper and lower bits of each digits +Integer.PARTITIONSIZE = math.floor(math.sqrt(Integer.DIGITSIZE)) + +--- Method for computing the gcd of two integers using Euclid's algorithm. +--- @param a Integer +--- @param b Integer +--- @return Integer +function Integer.gcd(a, b) + while b ~= Integer.zero() do + a, b = b, a%b + end + return a +end + +--- Method for computing the gcd of two integers using Euclid's algorithm. +--- Also returns Bezout's coefficients via extended gcd. +--- @param a Integer +--- @param b Integer +--- @return Integer, Integer, Integer +function Integer.extendedgcd(a, b) + local oldr, r = a, b + local olds, s = Integer.one(), Integer.zero() + local oldt, t = Integer.zero(), Integer.one() + while r ~= Integer.zero() do + local q = oldr // r + oldr, r = r, oldr - q*r + olds, s = s, olds - q*s + oldt, t = t, oldt - q*t + end + return oldr, olds, oldt +end + +--- Method for computing the larger of two integers. +--- Also returns the other integer for sorting purposes. +--- @param a Integer +--- @param b Integer +--- @return Integer, Integer +function Integer.max(a, b) + if a > b then + return a, b + end + return b, a +end + +--- Method for computing the smaller of two integers. +--- Also returns the other integer for sorting purposes. +--- @param a Integer +--- @param b Integer +--- @return Integer, Integer +function Integer.min(a, b) + if a < b then + return a, b + end + return b, a +end + +--- Methods for computing the larger magnitude of two integers. +--- Also returns the other integer for sorting purposes, and the number -1 if the two values were swapped, 1 if not. +--- @param a Integer +--- @param b Integer +--- @return Integer, Integer, number +function Integer.absmax(a, b) + if b:ltabs(a) then + return a, b, 1 + end + return b, a, -1 +end + +-- Returns the ceiling of the log base (defaults to 10) of a. +-- In other words, returns the least n such that base^n > a. +--- @param a Integer +--- @param base Integer +--- @return Integer +function Integer.ceillog(a, base) + base = base or Integer(10) + local k = Integer.zero() + + while (base ^ k) < a do + k = k + Integer.one() + end + + return k +end + +--- Returns a ^ b (mod n). This should be used when a ^ b is potentially large. +--- @param a Integer +--- @param b Integer +--- @param n Integer +--- @return Integer +function Integer.powmod(a, b, n) + if n == Integer.one() then + return Integer.zero() + else + local r = Integer.one() + a = a % n + while b > Integer.zero() do + if b % Integer(2) == Integer.one() then + r = (r * a) % n + end + a = (a ^ Integer(2)) % n + b = b // Integer(2) + end + return r + end +end + +--- @return RingIdentifier +local t = {ring=Integer} +t = setmetatable(t, {__index = Integer, __eq = function(a, b) + return a["ring"] == b["ring"] +end, __tostring = function(a) + return "ZZ" +end}) +function Integer.makering() + return t +end + + +---------------------------- +-- Instance functionality -- +---------------------------- + +-- So we don't have to copy the Euclidean operations each time we create an integer. +local __o = Copy(__EuclideanOperations) +__o.__index = Integer +__o.__tostring = function(a) -- Only works if the digit size is a power of 10 + local out = "" + for i, digit in ipairs(a) do + local pre = tostring(math.floor(digit)) + if i ~= #a then + while #pre ~= Integer.DIGITLENGTH do + pre = "0" .. pre + end + end + out = pre .. out + end + if a.sign == -1 then + out = "-" .. out + end + return out +end +__o.__div = function(a, b) -- Constructor for a rational number disguised as division + if not b.getring then + return BinaryOperation.DIVEXP({a, b}) + end + if(a:getring() == Integer:getring() and b:getring() == Integer:getring()) then + return Rational(a, b) + end + return __FieldOperations.__div(a, b) +end +__o.__concat = function(a, b) -- Like a decimal, but fancier. Used mainly for the parser with decimal numbers. + return a + b / (Integer(10) ^ Integer.ceillog(b)) +end + +--- Creates a new integer given a string or number representation of the integer. +--- @param n number|string|Integer +--- @return Integer +function Integer:new(n) + local o = {} + o = setmetatable(o, __o) + + if not n then + o[1] = 0 + o.sign = 0 + return o + end + + -- Can convert any floating-point number into an integer, though we generally only want to pass whole numbers into this. + -- This will only approximate very large floating point numbers to a small proportion of the total significant digits + -- After that the result will just be nonsense - strings should probably be used for big numbers + if type(n) == "number" then + n = math.floor(n) + if n == 0 then + o[1] = 0 + o.sign = 0 + else + if n < 0 then + n = -n + o.sign = -1 + else + o.sign = 1 + end + local i = 1 + while n >= Integer.DIGITSIZE do + o[i] = n % Integer.DIGITSIZE + n = n // Integer.DIGITSIZE + i = i + 1 + end + o[i] = n + end + -- Only works on strings that are exact (signed) integers + elseif type(n) == "string" then + if not tonumber(n) then + error("Sent parameter of wrong type: " .. n .. " is not an integer.") + end + if n == "0" then + o[1] = 0 + o.sign = 0 + else + local s = 1 + if string.sub(n, 1, 1) == "-" then + s = s + 1 + o.sign = -1 + else + o.sign = 1 + end + + while string.sub(n, s, s) == "0" do + s = s + 1 + end + + local e = #n + local i = 1 + while e > s + Integer.DIGITLENGTH - 1 do + o[i] = tonumber(string.sub(n, e - Integer.DIGITLENGTH + 1, e)) + e = e - Integer.DIGITLENGTH + i = i + 1 + end + o[i] = tonumber(string.sub(n, s, e)) or 0 + end + -- Copying is expensive in Lua, so this constructor probably should only sparsely be called with an Integer argument. + elseif type(n) == "table" then + o = Copy(n) + else + error("Sent parameter of wrong type: Integer does not accept " .. type(n) .. ".") + end + + return o +end + +--- Returns the ring this object is an element of. +--- @return RingIdentifier +function Integer:getring() + return t +end + +--- @param ring RingIdentifier +--- @return Ring +function Integer:inring(ring) + if ring == self:getring() then + return self + end + + if ring == PolynomialRing:getring() then + return PolynomialRing({self:inring(ring.child)}, ring.symbol) + end + + if ring == Rational:getring() then + if ring.child then + return Rational(self:inring(ring.child), self:inring(ring.child):one(), true) + end + return Rational(self, Integer.one(), true):inring(ring) + end + + if ring == IntegerModN:getring() then + return IntegerModN(self, ring.modulus) + end + + error("Unable to convert element to proper ring.") +end + +--- @param b Integer +--- @return Integer +function Integer:add(b) + if self.sign == 1 and b.sign == -1 then + return self:usub(b, 1) + end + if self.sign == -1 and b.sign == 1 then + return self:usub(b, -1) + end + + local sign = self.sign + if sign == 0 then + sign = b.sign + end + return self:uadd(b, sign) +end + +--- Addition without sign so we don't have to create an entire new integer when switching signs. +--- @param b Integer +--- @param sign number +--- @return Integer +function Integer:uadd(b, sign) + local o = Integer() + o.sign = sign + + local c = 0 + local n = math.max(#self, #b) + for i = 1, n do + local s = (self[i] or 0) + (b[i] or 0) + c + if s >= Integer.DIGITSIZE then + o[i] = s - Integer.DIGITSIZE + c = 1 + else + o[i] = s + c = 0 + end + end + if c == 1 then + o[n + 1] = c + end + return o +end + +--- @param b Integer +--- @return Integer +function Integer:sub(b) + if self.sign == 1 and b.sign == -1 then + return self:uadd(b, 1) + end + if self.sign == -1 and b.sign == 1 then + return self:uadd(b, -1) + end + + local sign = self.sign + if sign == 0 then + sign = b.sign + end + return self:usub(b, sign) +end + +-- Subtraction without sign so we don't have to create an entire new integer when switching signs. +-- Uses subtraction by compliments. +--- @param b Integer +--- @param sign number +--- @return Integer +function Integer:usub(b, sign) + local a, b, swap = Integer.absmax(self, b) + local o = Integer() + o.sign = sign * swap + + local c = 0 + local n = #a + for i = 1, n do + local s = (a[i] or 0) + Integer.DIGITSIZE - 1 - (b[i] or 0) + c + if i == 1 then + s = s + 1 + end + if s >= Integer.DIGITSIZE then + o[i] = s - Integer.DIGITSIZE + c = 1 + else + o[i] = s + c = 0 + end + end + + -- Remove leading zero digits, since we want integer representations to be unique. + while o[n] == 0 do + o[n] = nil + n = n - 1 + end + + if not o[1] then + o[1] = 0 + o.sign = 0 + end + + return o +end + +--- @return Integer +function Integer:neg() + local o = Integer() + o.sign = -self.sign + for i, digit in ipairs(self) do + o[i] = digit + end + return o +end + +--- @param b Integer +--- @return Integer +function Integer:mul(b) + local o = Integer() + o.sign = self.sign * b.sign + if o.sign == 0 then + o[1] = 0 + return o + end + + -- Fast single-digit multiplication in the most common case + if #self == 1 and #b == 1 then + o[2], o[1] = self:mulone(self[1], b[1]) + + if o[2] == 0 then + o[2] = nil + end + + return o + end + + -- "Grade school" multiplication algorithm for numbers with small numbers of digits works faster than Karatsuba + local n = #self + local m = #b + o[1] = 0 + o[2] = 0 + for i = 2, n+m do + o[i + 1] = 0 + for j = math.max(1, i-m), math.min(n, i-1) do + local u, l = self:mulone(self[j], b[i - j]) + o[i - 1] = o[i - 1] + l + o[i] = o[i] + u + if o[i - 1] >= Integer.DIGITSIZE then + o[i - 1] = o[i - 1] - Integer.DIGITSIZE + o[i] = o[i] + 1 + end + if o[i] >= Integer.DIGITSIZE then + o[i] = o[i] - Integer.DIGITSIZE + o[i + 1] = o[i + 1] + 1 + end + end + end + + -- Remove leading zero digits, since we want integer representations to be unique. + if o[n+m+1] == 0 then + o[n+m+1] = nil + end + + if o[n+m] == 0 then + o[n+m] = nil + end + + return o +end + +--- Multiplies two single-digit numbers and returns two digits. +--- @param a number +--- @param b number +--- @return number, number +function Integer:mulone(a, b) + local P = Integer.PARTITIONSIZE + + local a1 = a // P + local a2 = a % P + local b1 = b // P + local b2 = b % P + + local u = a1 * b1 + local l = a2 * b2 + + local m = ((a1 * b2) + (b1 * a2)) + local mu = m // P + local ml = m % P + + u = u + mu + l = l + ml * P + + if l >= Integer.DIGITSIZE then + l = l - Integer.DIGITSIZE + u = u + 1 + end + + return u, l +end + +--- Naive exponentiation is slow even for small exponents, so this uses binary exponentiation. +--- @param b Integer +--- @return Integer +function Integer:pow(b) + if b < Integer.zero() then + return Integer.one() / (self ^ -b) + end + + if b == Integer.zero() then + return Integer.one() + end + + -- Fast single-digit exponentiation + if #self == 1 and #b == 1 then + local test = (self.sign * self[1]) ^ b[1] + if test < Integer.DIGITSIZE and test > -Integer.DIGITSIZE then + return Integer(test) + end + end + + local x = self + local y = Integer.one() + while b > Integer.one() do + if b[1] % 2 == 0 then + x = x * x + b = b:divbytwo() + else + y = x * y + x = x * x + b = b:divbytwo() + end + end + + return x * y +end + +-- Fast integer division by two for binary exponentiation. +--- @return Integer +function Integer:divbytwo() + local o = Integer() + o.sign = self.sign + for i = #self, 1, -1 do + if self[i] % 2 == 0 then + o[i] = self[i] // 2 + else + o[i] = self[i] // 2 + if i ~= 1 then + o[i - 1] = self[i - 1] * 2 + end + end + end + return o +end + +--- Division with remainder over the integers. Uses the standard base 10 long division algorithm. +--- @param b Integer +--- @return Integer, Integer +function Integer:divremainder(b) + if self >= Integer.zero() and b > self or self <= Integer.zero() and b < self then + return Integer.zero(), Integer(self) + end + + if #self == 1 and #b == 1 then + return Integer((self[1]*self.sign) // (b[1]*b.sign)), Integer((self[1]*self.sign) % (b[1]*b.sign)) + end + + local Q = Integer() + local R = Integer() + + Q.sign = self.sign * b.sign + R.sign = 1 + local negativemod = false + if b.sign == -1 then + b.sign = -b.sign + negativemod = true + end + + for i = #self, 1, -1 do + local s = tostring(math.floor(self[i])) + while i ~= #self and #s ~= Integer.DIGITLENGTH do + s = "0" .. s + end + Q[i] = 0 + for j = 1, #s do + R = R:mulbyten() + R[1] = R[1] + tonumber(string.sub(s, j, j)) + if R[1] > 0 then + R.sign = 1 + end + while R >= b do + R = R - b + Q[i] = Q[i] + 10^(#s - j) + end + end + end + + -- Remove leading zero digits, since we want integer representations to be unique. + while Q[#Q] == 0 do + Q[#Q] = nil + end + + if negativemod then + R = -R + b.sign = -b.sign + elseif self.sign == -1 then + R = b - R + end + + return Q, R +end + +--- Fast in-place multiplication by ten for the division algorithm. This means the number IS MODIFIED by this method unlike the rest of the library. +--- @return Integer +function Integer:mulbyten() + local DIGITSIZE = Integer.DIGITSIZE + for i, _ in ipairs(self) do + self[i] = self[i] * 10 + end + for i, _ in ipairs(self) do + if self[i] > DIGITSIZE then + local msd = self[i] // DIGITSIZE + if self[i+1] then + self[i+1] = self[i+1] + msd + else + self[i+1] = msd + end + self[i] = self[i] - DIGITSIZE*msd + end + end + return self +end + +--- @param b Integer +--- @return boolean +function Integer:eq(b) + for i, digit in ipairs(self) do + if not b[i] or not (b[i] == digit) then + return false + end + end + return #self == #b and self.sign == b.sign +end + +--- @param b Integer +--- @return boolean +function Integer:lt(b) + local selfsize = #self + local bsize = #b + if selfsize < bsize then + return b.sign == 1 + end + if selfsize > bsize then + return self.sign == -1 + end + local n = selfsize + while n > 0 do + if self[n]*self.sign < b[n]*b.sign then + return true + end + if self[n]*self.sign > b[n]*b.sign then + return false + end + n = n - 1 + end + return false +end + +--- Same as less than, but ignores signs. +--- @param b Integer +--- @return boolean +function Integer:ltabs(b) + if #self < #b then + return true + end + if #self > #b then + return false + end + local n = #self + while n > 0 do + if self[n] < b[n] then + return true + end + if self[n] > b[n] then + return false + end + n = n - 1 + end + return false +end + +--- @param b Integer +--- @return boolean +function Integer:le(b) + local selfsize = #self + local bsize = #b + if selfsize < bsize then + return b.sign == 1 + end + if selfsize > bsize then + return self.sign == -1 + end + local n = selfsize + while n > 0 do + if self[n]*self.sign < b[n]*b.sign then + return true + end + if self[n]*self.sign > b[n]*b.sign then + return false + end + n = n - 1 + end + return true +end + +local zero = Integer:new(0) +--- @return Integer +function Integer:zero() + return zero +end + +local one = Integer:new(1) +--- @return Integer +function Integer:one() + return one +end + +--- Returns this integer as a floating point number. Can only approximate the value of large integers. +--- @return number +function Integer:asnumber() + local n = 0 + for i, digit in ipairs(self) do + n = n + digit * Integer.DIGITSIZE ^ (i - 1) + end + return self.sign*math.floor(n) +end + +--- Returns all positive divisors of the integer. Not guarenteed to be in any order. +--- @return table<number, Integer> +function Integer:divisors() + local primefactors = self:primefactorizationrec() + local divisors = {} + + local terms = {} + for prime in pairs(primefactors) do + if prime == Integer(-1) then + primefactors[prime] = nil + end + terms[prime] = Integer.zero() + end + + local divisor = Integer.one() + + while true do + divisors[#divisors+1] = divisor + for prime, power in pairs(primefactors) do + if terms[prime] < power then + terms[prime] = terms[prime] + Integer.one() + divisor = divisor * prime + break + else + terms[prime] = Integer.zero() + divisor = divisor / (prime ^ power) + end + end + if divisor == Integer.one() then + break + end + end + + return divisors +end + +--- Returns whether this integer is a prime power, of the form p^a for prime p and positive integer a. +--- If it is a prime power, also returns the prime and the power. +--- @return boolean, Expression|nil, Expression|nil +function Integer:isprimepower() + if self <= Integer.one() then + return false + end + local factorization = self:primefactorization() + if factorization:type() == BinaryOperation and #factorization:subexpressions() == 1 then + return true, factorization.expressions[1].expressions[2], factorization.expressions[1].expressions[1] + end + return false +end + +--- Returns whether this integer is a perfect power, of the form a^b for positive integers a and b. +--- If it is a prime power, also returns the prime and the power. +--- @return boolean, Expression|nil, Expression|nil +function Integer:isperfectpower() + if self <= Integer.one() then + return false + end + local factorization = self:primefactorization() + if factorization:type() ~= BinaryOperation then + return false + end + local power = Integer.zero() + for _, term in ipairs(factorization:subexpressions()) do + power = Integer.gcd(power, term.expressions[2]) + if power == Integer.one() then + return false + end + end + local base = Integer.one() + for _, term in ipairs(factorization:subexpressions()) do + base = base * term.expressions[1] ^ (term.expressions[2] / power) + end + return true, base, power +end + +--- Returns the prime factorization of this integer as a expression. +--- @return Expression +function Integer:primefactorization() + if not Integer.FACTORIZATIONLIMIT then + Integer.FACTORIZATIONLIMIT = Integer(Integer.DIGITSIZE) + end + if self > Integer.FACTORIZATIONLIMIT then + return self + end + local result = self:primefactorizationrec() + local mul = {} + local i = 1 + for factor, degree in pairs(result) do + mul[i] = BinaryOperation.POWEXP({factor, degree}) + i = i + 1 + end + return BinaryOperation.MULEXP(mul):lock(Expression.NIL) +end + +--- Recursive part of prime factorization using Pollard Rho. +function Integer:primefactorizationrec() + if self < Integer.zero() then + return Integer.mergefactors({[Integer(-1)]=Integer.one()}, (-self):primefactorizationrec()) + end + if self == Integer.one() then + return {[Integer.one()]=Integer.one()} + end + local result = self:findafactor() + if result == self then + return {[result]=Integer.one()} + end + local remaining = self / result + return Integer.mergefactors(result:primefactorizationrec(), remaining:primefactorizationrec()) +end + + +function Integer.mergefactors(a, b) + local result = Copy(a) + + for factor, degree in pairs(b) do + for ofactor, odegree in pairs(result) do + if factor == ofactor then + result[ofactor] = degree + odegree + goto continue + end + end + result[factor] = degree + ::continue:: + end + return result +end + +-- Return a non-trivial factor of n via Pollard Rho, or returns n if n is prime. +function Integer:findafactor() + if self:isprime() then + return self + end + + if self % Integer(2) == Integer.zero() then + return Integer(2) + end + + if self % Integer(3) == Integer.zero() then + return Integer(3) + end + + if self % Integer(5) == Integer.zero() then + return Integer(5) + end + + local g = function(x) + local temp = Integer.powmod(x, Integer(2), self) + return temp + end + + local xstart = Integer(2) + while xstart < self do + local x = xstart + local y = xstart + local d = Integer.one() + while d == Integer.one() do + x = g(x) + y = g(g(y)) + d = Integer.gcd((x - y):abs(), self) + end + + if d < self then + return d + end + + xstart = xstart + Integer.one() + end +end + +--- Uses Miller-Rabin to determine whether a number is prime up to a very large number. +local smallprimes = {Integer:new(2), Integer:new(3), Integer:new(5), Integer:new(7), Integer:new(11), Integer:new(13), Integer:new(17), +Integer:new(19), Integer:new(23), Integer:new(29), Integer:new(31), Integer:new(37), Integer:new(41), Integer:new(43), Integer:new(47)} + +function Integer:isprime() + if self % Integer(2) == Integer.zero() then + if self == Integer(2) then + return true + end + return false + end + + if self == Integer.one() then + return false + end + + for _, value in pairs(smallprimes) do + if value == self then + return true + end + end + + local r = Integer.zero() + local d = self - Integer.one() + while d % Integer(2) == Integer.zero() do + r = r + Integer.one() + d = d / Integer(2) + end + + for _, a in ipairs(smallprimes) do + local s = r + local x = Integer.powmod(a, d, self) + if x == Integer.one() or x == self - Integer.one() then + goto continue + end + + while s > Integer.zero() do + x = Integer.powmod(x, Integer(2), self) + if x == self - Integer.one() then + goto continue + end + s = s - Integer.one() + end + do + return false + end + ::continue:: + end + + return true +end + +--- Returns the absolute value of an integer. +--- @return Integer +function Integer:abs() + if self.sign >= 0 then + return Integer(self) + end + return -self +end + +----------------- +-- Inheritance -- +----------------- + +__Integer.__index = EuclideanDomain +__Integer.__call = Integer.new +Integer = setmetatable(Integer, __Integer) + +---------------------- +-- Static constants -- +---------------------- + +Integer.FACTORIZATIONLIMIT = Integer(Integer.DIGITSIZE)
\ No newline at end of file diff --git a/macros/luatex/latex/luacas/tex/algebra/integerquotientring.lua b/macros/luatex/latex/luacas/tex/algebra/integerquotientring.lua new file mode 100644 index 0000000000..4c2188c9e5 --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/integerquotientring.lua @@ -0,0 +1,197 @@ +--- @class IntegerModN +--- Represents an element of the ring of integers mod n (this is also a field iff n is prime). +--- @field element Integer +--- @field modulus Integer + +IntegerModN = {} +__IntegerModN = {} + +-- Metatable for ring objects. +local __obj = {__index = IntegerModN, __eq = function(a, b) + return a["ring"] == b["ring"] and (a["modulus"] == b["modulus"] or a["modulus"] == nil or b["modulus"] == nil) +end, __tostring = function(a) + if a.modulus then return "Z/Z" .. tostring(a.modulus) else return "(Generic Integer Mod Ring)" end +end} + +-------------------------- +-- Static functionality -- +-------------------------- + +--- Creates a new ring with the given modulus. +--- @param modulus Integer +--- @return RingIdentifier +function IntegerModN.makering(modulus) + local t = {ring = IntegerModN} + t.modulus = modulus + t = setmetatable(t, __obj) + return t +end + +-- Shorthand constructor for a ring with a particular modulus. +function IntegerModN.R(modulus) + return IntegerModN.makering(modulus) +end + +---------------------------- +-- Instance functionality -- +---------------------------- + +-- So we don't have to copy the field operations each time +local __o +__o = Copy(__FieldOperations) + +__o.__index = IntegerModN +__o.__tostring = function(a) + return tostring(a.element) +end + +--- Creates a new integer i in Z/nZ. +--- @param i Integer +--- @param n Integer +--- @return IntegerModN +function IntegerModN:new(i, n) + local o = {} + + if n:getring() ~= Integer:getring() or n < Integer.one() then + error("Argument error: modulus must be an integer greater than 0.") + end + + o = setmetatable(o, __o) + + if i < Integer.zero() or i >= n then + i = i % n + end + + o.element = i + o.modulus = n + + return o +end + +--- @return RingIdentifier +function IntegerModN:getring() + local t = {ring = IntegerModN} + if self then + t.modulus = self.modulus + end + t = setmetatable(t, __obj) + return t +end + +--- @param ring RingIdentifier +--- @return Ring +function IntegerModN:inring(ring) + if ring == IntegerModN:getring() then + if ring.modulus then + return IntegerModN(self.element, ring.modulus) + end + return self + end + + if ring == PolynomialRing:getring() then + return PolynomialRing({self:inring(ring.child)}, ring.symbol) + end + + if ring == Rational:getring() and ring.symbol then + return Rational(self:inring(ring.child), self:inring(ring.child):one(), true) + end + + if ring == Integer:getring() then + return self.element:inring(ring) + end + + error("Unable to convert element to proper ring.") +end + +--- @param b IntegerModN +--- @return IntegerModN +function IntegerModN:add(b) + return IntegerModN(self.element + b.element, self.modulus) +end + +--- @return IntegerModN +function IntegerModN:neg() + return IntegerModN(-self.element, self.modulus) +end + +--- @param b IntegerModN +--- @return IntegerModN +function IntegerModN:mul(b) + return IntegerModN(self.element * b.element, self.modulus) +end + +-- Overrides the generic power method with powmod. +--- @param b IntegerModN +--- @return IntegerModN +function IntegerModN:pow(b) + return IntegerModN(Integer.powmod(self.element, b.element, self.modulus), self.modulus) +end + +-- Returns the multiplicative inverse of this number if it exists. +--- @return IntegerModN +function IntegerModN:inv() + local r, t, _ = Integer.extendedgcd(self.element, self.modulus) + + if r > Integer.one() then + error("Element does not have an inverse in this ring") + end + + return IntegerModN(t, self.modulus) +end + +--- @param b IntegerModN +--- @return IntegerModN +function IntegerModN:div(b) + return self:mul(b:inv()) +end + +--- @param b IntegerModN +--- @return boolean +function IntegerModN:eq(b) + return self.element == b.element +end + +--- @param b IntegerModN +--- @return boolean +function IntegerModN:lt(b) + return self.element < b.element +end + +--- @param b IntegerModN +--- @return boolean +function IntegerModN:le(b) + return self.element <= b.element +end + +--- @return IntegerModN +function IntegerModN:zero() + if not self or not self.modulus then + return Integer.zero() + end + return IntegerModN(Integer.zero(), self.modulus) +end + +--- @return IntegerModN +function IntegerModN:one() + if not self or not self.modulus then + return Integer.one() + end + return IntegerModN(Integer.one(), self.modulus) +end + +--- @return string +function IntegerModN:tolatex(mod) + mod = mod or false + if mod then + return self.element:tolatex() .. "\\bmod{" .. self.modulus:tolatex() .. "}" + else + return self.element:tolatex() + end +end +----------------- +-- Inheritance -- +----------------- + +__IntegerModN.__index = Field +__IntegerModN.__call = IntegerModN.new +IntegerModN = setmetatable(IntegerModN, __IntegerModN)
\ No newline at end of file diff --git a/macros/luatex/latex/luacas/tex/algebra/logarithm.lua b/macros/luatex/latex/luacas/tex/algebra/logarithm.lua new file mode 100644 index 0000000000..07a9db84d1 --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/logarithm.lua @@ -0,0 +1,186 @@ +--- @class Logarithm +--- An expression for the logarithm of an expression with respect to another. +--- Currently, logarithms are not being evaluated since we are just doing symbolic computation. +--- @field base Expression +--- @field expression Expression +Logarithm = {} +__Logarithm = {} + +---------------------------- +-- Instance functionality -- +---------------------------- + +--- Creates a new logarithm expression with the given symbol and expression. +--- @param base Expression +--- @param expression Expression +--- @ +function Logarithm:new(base, expression) + local o = {} + local __o = Copy(__ExpressionOperations) + + o.base = Copy(base) + o.expression = Copy(expression) + + __o.__index = Logarithm + __o.__tostring = function(a) + return 'log(' .. tostring(base) .. ', ' .. tostring(expression) .. ')' + end + __o.__eq = function(a, b) + -- This shouldn't be needed, since __eq should only fire if both metamethods have the same function, but for some reason Lua always runs this anyway + if not b:type() == Logarithm then + return false + end + return a.base == b.base and a.expression == b.expression + end + o = setmetatable(o, __o) + + return o +end + +--- @return Expression +function Logarithm:evaluate() + + if (self.base:isconstant() and (self.base <= Integer.zero() or self.base == Integer.one())) or + (self.expression:isconstant() and self.expression <= Integer.zero()) then + error("Arithmetic error: division by zero") + end + + if not self.base:isconstant() or not self.expression:isconstant() then + return self + end + + local power = Integer.one() + local base = self.base + if base:type() == Integer then + local pp, b, p = base:isperfectpower() + if pp then + base = b + power = p + end + end + + if base:type() == Rational then + local ppn, bn, pn = base.numerator:isperfectpower() + local ppd, bd, pd = base.denominator:isperfectpower() + if base.numerator == Integer.one() then + ppn = true + bn = Integer.one() + pn = pd + end + if ppn and ppd and pn == pd then + base = bn / bd + power = pn + end + end + + local result = Integer.one() + local expression = self.expression + local sign = Integer.one() + if base < Integer.one() then + base = Integer.one() / base + sign = -sign + end + + + local current = base + while current < expression do + current = current * base + result = result + Integer.one() + end + if current == expression then + return sign * result / power + else + while current > expression do + current = current / base + result = result - Integer.one() + end + if current == expression then + return sign * result / power + end + end + + return self +end + +--- @return Expression +function Logarithm:autosimplify() + + local base = self.base:autosimplify() + local expression = self.expression:autosimplify() + + local evaluated = Logarithm(base, expression):evaluate() + if evaluated:type() ~= Logarithm then + return evaluated + end + + -- Uses the property that log(b, 1) = 0 + if expression == Integer.one() then + return Integer.zero() + end + + -- Uses the property that log(b, b) = 1 + if expression == base then + return Integer.one() + end + + -- Uses the propery that log(b, x^y) = y * log(b, x) + if expression.operation == BinaryOperation.POW then + return BinaryOperation.MULEXP({expression.expressions[2], Logarithm(base, expression.expressions[1])}):autosimplify() + end + + if expression:type() == Rational and expression.numerator == Integer.one() then + return (-Logarithm(base, expression.denominator)):autosimplify() + end + + -- Our expression cannot be simplified + return Logarithm(base, expression) +end + +--- @return Expression +function Logarithm:expand() + return Logarithm(self.base:expand(), self.expression:expand()):autosimplify() +end + +--- @return table<number, Expression> +function Logarithm:subexpressions() + return {self.base, self.expression} +end + +--- @param subexpressions table<number, Expression> +--- @return Logarithm +function Logarithm:setsubexpressions(subexpressions) + return Logarithm(subexpressions[1], subexpressions[2]) +end + +--- @param other Expression +--- @return boolean +function Logarithm:order(other) + return FunctionExpression("log", {self.base, self.expression}):order(other) +end + +--- @return string +function Logarithm:tolatex() + if self.base == E then + return '\\ln\\mathopen{}\\left(' .. self.expression:tolatex() .. '\\right)' + end + return '\\log_' .. self.base:tolatex() .. '\\mathopen{}\\left(' .. self.expression:tolatex() .. '\\right)' +end + +----------------- +-- Inheritance -- +----------------- +__Logarithm.__index = CompoundExpression +__Logarithm.__call = Logarithm.new +Logarithm = setmetatable(Logarithm, __Logarithm) + +---------------------- +-- Static constants -- +---------------------- + +LOG = function(base, expression) + return Logarithm(base, expression) +end + +LN = function(expression) + return Logarithm(E, expression) +end
\ No newline at end of file diff --git a/macros/luatex/latex/luacas/tex/algebra/polynomialring.lua b/macros/luatex/latex/luacas/tex/algebra/polynomialring.lua new file mode 100644 index 0000000000..568c21c921 --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/polynomialring.lua @@ -0,0 +1,860 @@ +--- @class PolynomialRing +--- Represents an element of a polynomial ring. +--- @field coefficients table<number, Ring> +--- @field symbol SymbolExpression +--- @field ring RingIdentifier +PolynomialRing = {} +__PolynomialRing = {} + +-- Metatable for ring objects. +local __obj = {__index = PolynomialRing, __eq = function(a, b) + return a["ring"] == b["ring"] and + (a["child"] == b["child"] or a["child"] == nil or b["child"] == nil) and + (a["symbol"] == b["symbol"] or a["child"] == nil or b["child"] == nil) +end, __tostring = function(a) + if a.child and a.symbol then return tostring(a.child) .. "[" .. a.symbol .. "]" else return "(Generic Polynomial Ring)" end +end} + +-------------------------- +-- Static functionality -- +-------------------------- + +--- Creates a new ring with the given symbol and child ring. +--- @param symbol SymbolExpression +--- @param child RingIdentifier +--- @return RingIdentifier +function PolynomialRing.makering(symbol, child) + local t = {ring = PolynomialRing} + t.symbol = symbol + t.child = child + t = setmetatable(t, __obj) + return t +end + +-- Shorthand constructor for a polynomial ring with integer or integer mod ring coefficients. +function PolynomialRing.R(symbol, modulus) + if modulus then + return PolynomialRing.makering(symbol, IntegerModN.makering(modulus)) + end + return PolynomialRing.makering(symbol, Integer.getring()) +end + +--- Returns the GCD of two polynomials in a ring, assuming both rings are euclidean domains. +--- @param a PolynomialRing +--- @param b PolynomialRing +--- @return PolynomialRing +function PolynomialRing.gcd(a, b) + if a.symbol ~= b.symbol then + error("Cannot take the gcd of two polynomials with different symbols") + end + while b ~= Integer.zero() do + a, b = b, a % b + end + return a // a:lc() +end + +-- Returns the GCD of two polynomials in a ring, assuming both rings are euclidean domains. +-- Also returns bezouts coefficients via extended gcd. +--- @param a PolynomialRing +--- @param b PolynomialRing +--- @return PolynomialRing, PolynomialRing, PolynomialRing +function PolynomialRing.extendedgcd(a, b) + local oldr, r = a, b + local olds, s = Integer.one(), Integer.zero() + local oldt, t = Integer.zero(), Integer.one() + while r ~= Integer.zero() do + local q = oldr // r + oldr, r = r, oldr - q*r + olds, s = s, olds - q*s + oldt, t = t, oldt - q*t + end + return oldr // oldr:lc(), olds // oldr:lc(), oldt // oldr:lc() +end + +-- Returns the resultant of two polynomials in the same ring, whose coefficients are all part of a field. +--- @param a PolynomialRing +--- @param b PolynomialRing +--- @return Field +function PolynomialRing.resultant(a, b) + + if a.ring == PolynomialRing.getring() or b.ring == PolynomialRing.getring() then + return PolynomialRing.resultantmulti(a, b) + end + + local m, n = a.degree, b.degree + if n == Integer.zero() then + return b.coefficients[0]^m + end + + local r = a % b + if r == Integer.zero() then + return r.coefficients[0] + end + + local s = r.degree + local l = b:lc() + + return Integer(-1)^(m*n) * l^(m-s) * PolynomialRing.resultant(b, r) +end + +-- Returns the resultant of two polynomials in the same ring, whose coefficients are not part of a field. +--- @param a PolynomialRing +--- @param b PolynomialRing +--- @return Ring +function PolynomialRing.resultantmulti(a, b) + local m, n = a.degree, b.degree + + if m < n then + return Integer(-1) ^ (m * n) * PolynomialRing.resultantmulti(b, a) + end + if n == Integer.zero() then + return b.coefficients[0]^m + end + + local delta = m - n + Integer(1) + local _ , r = PolynomialRing.pseudodivide(a, b) + if r == Integer.zero() then + return r.coefficients[0] + end + + local s = r.degree + local w = Integer(-1)^(m*n) * PolynomialRing.resultant(b, r) + local l = b:lc() + local k = delta * n - m + s + local f = l ^ k + return w // f +end + +-- Given two polynomials a and b, returns a list of the remainders generated by the monic Euclidean algorithm. +--- @param a PolynomialRing +--- @param b PolynomialRing +--- @return table<number, Ring> +function PolynomialRing.monicgcdremainders(a, b) + if a.symbol ~= b.symbol then + error("Cannot take the gcd of two polynomials with different symbols") + end + + local remainders = {a / a:lc(), b / b:lc()} + while true do + local q = remainders[#remainders - 1] // remainders[#remainders] + local c = remainders[#remainders - 1] - q*remainders[#remainders] + if c ~= Integer.zero() then + remainders[#remainders+1] = c/c:lc() + else + break + end + end + + return remainders +end + +-- Returns the partial fraction decomposition of the rational function g/f +-- given g, f, and some (not nessecarily irreducible) factorization of f. +-- If the factorization is omitted, the irreducible factorization is used. +-- The degree of g must be less than the degree of f. +--- @param g PolynomialRing +--- @param f PolynomialRing +--- @param ffactors Expression +--- @return Expression +function PolynomialRing.partialfractions(g, f, ffactors) + + if g.degree >= f.degree then + error("Argument Error: The degree of g must be less than the degree of f.") + end + + -- Converts f to a monic polynomial. + g = g * f:lc() + f = f / f:lc() + + ffactors = ffactors or f:factor() + + local expansionterms = {} + + for _, factor in ipairs(ffactors.expressions) do + local k + local m + if factor.getring and factor:getring() == PolynomialRing:getring() then + m = factor + k = Integer.one() + elseif not factor:isconstant() then + m = factor.expressions[1] + k = factor.expressions[2] + end + + if not factor:isconstant() then + -- Uses Chinese Remainder Theorem for each factor to determine the numerator of the term in the decomposition + local mk = m^k + local v = g % mk + local _, minv, _ = PolynomialRing.extendedgcd(f // mk, mk) + local c = v*minv % mk + + + if k == Integer.one() then + expansionterms[#expansionterms+1] = BinaryOperation.ADDEXP({BinaryOperation.DIVEXP({c, BinaryOperation.POWEXP({m, Integer.one()})})}) + else + -- Uses the p-adic expansion of c to split terms with repeated roots. + local q = c + local r + local innerterms = {} + for i = k:asnumber(), 1, -1 do + q, r = q:divremainder(m) + innerterms[#innerterms+1] = BinaryOperation.DIVEXP({r, BinaryOperation.POWEXP({m, Integer(i)})}) + end + expansionterms[#expansionterms+1] = BinaryOperation.ADDEXP(innerterms) + end + end + end + + return BinaryOperation.ADDEXP(expansionterms) + +end + +---------------------------- +-- Instance functionality -- +---------------------------- + +-- So we don't have to copy the Euclidean operations each time +local __o = Copy(__EuclideanOperations) +__o.__index = PolynomialRing +__o.__tostring = function(a) + local out = "" + local loc = a.degree:asnumber() + while loc >= 0 do + if a.ring == PolynomialRing.getring() or (a.ring == Rational.getring() and a.ring.symbol) then + out = out .. "(" .. tostring(a.coefficients[loc]) .. ")" .. a.symbol .. "^" .. tostring(math.floor(loc)) .. "+" + else + out = out .. tostring(a.coefficients[loc]) .. a.symbol .. "^" .. tostring(math.floor(loc)) .. "+" + end + loc = loc - 1 + end + return string.sub(out, 1, string.len(out) - 1) +end +__o.__div = function(a, b) + if not b.getring then + return BinaryOperation.DIVEXP({a, b}) + end + if Ring.resultantring(a.ring, b:getring()) ~= Ring.resultantring(a:getring(), b:getring()) then + return a:div(b:inring(Ring.resultantring(a:getring(), b:getring()))) + end + if b.ring and b:getring() == Rational:getring() and a.symbol == b.ring.symbol then + return a:inring(Ring.resultantring(a:getring(), b:getring())):div(b) + end + if a:getring() == b:getring() then + return Rational(a, b, true) + end + -- TODO: Fix this for arbitrary depth + if a:getring() == PolynomialRing:getring() and b:getring() == PolynomialRing:getring() and a.symbol == b.symbol then + local oring = Ring.resultantring(a:getring(), b:getring()) + return Rational(a:inring(oring), b:inring(oring), true) + end + return BinaryOperation.DIVEXP({a, b}) +end + +function PolynomialRing:tolatex() + local out = '' + local loc = self.degree:asnumber() + if loc == 0 then + return self.coefficients[loc]:tolatex() + end + if self.ring == Rational.getring() or self.ring == Integer.getring() or self.ring == IntegerModN.getring() then + if self.coefficients[loc] ~= Integer.one() then + out = out .. self.coefficients[loc]:tolatex() .. self.symbol + else + out = out .. self.symbol + end + if loc ~=1 then + out = out .. "^{" .. loc .. "}" + end + loc = loc -1 + while loc >=0 do + local coeff = self.coefficients[loc] + if coeff == Integer.one() then + if loc == 0 then + out = out .. "+" .. coeff:tolatex() + goto skip + else + out = out .. "+" + goto continue + end + end + if coeff == Integer(-1) then + if loc == 0 then + out = out .. "-" .. coeff:neg():tolatex() + goto skip + else + out = out .. "-" + goto continue + end + end + if coeff < Integer.zero() then + out = out .. "-" .. coeff:neg():tolatex() + end + if coeff == Integer.zero() then + goto skip + end + if coeff > Integer.zero() then + out = out .. "+" .. coeff:tolatex() + end + ::continue:: + if loc > 1 then + out = out .. self.symbol .. "^{" .. loc .. "}" + end + if loc == 1 then + out = out .. self.symbol + end + ::skip:: + loc = loc-1 + end + else + while loc >=0 do + if loc >=1 then + out = out .. self.coefficients[loc]:tolatex() .. self.symbol .. "^{" .. loc .. "} + " + else + out = out .. self.coefficients[loc]:tolatex() .. self.symbol .. "^{" .. loc .. "}" + end + loc = loc-1 + end + end + return out +end + +function PolynomialRing:isatomic() + --if self.degree >= Integer.one() then + -- return false + --else + return false + --end +end +--test + +-- Creates a new polynomial ring given an array of coefficients and a symbol +function PolynomialRing:new(coefficients, symbol, degree) + local o = {} + o = setmetatable(o, __o) + + if type(coefficients) ~= "table" then + error("Sent parameter of wrong type: Coefficients must be in an array") + end + o.coefficients = {} + o.degree = degree or Integer(-1) + + if type(symbol) ~= "string" and not symbol.symbol then + error("Symbol must be a string") + end + o.symbol = symbol.symbol or symbol + + -- Determines what ring the polynomial ring should have as its child + for index, coefficient in pairs(coefficients) do + if type(index) ~= "number" then + error("Sent parameter of wrong type: Coefficients must be in an array") + end + if not coefficient.getring then + error("Sent parameter of wrong type: Coefficients must be elements of a ring") + end + if not o.ring then + o.ring = coefficient:getring() + else + local newring = coefficient:getring() + local combinedring = Ring.resultantring(o.ring, newring) + if combinedring == newring then + o.ring = newring + elseif not o.ring == combinedring then + error("Sent parameter of wrong type: Coefficients must all be part of the same ring") + end + end + end + + if not coefficients[0] then + -- Constructs the coefficients when a new polynomial is instantiated as an array + for index, coefficient in ipairs(coefficients) do + o.coefficients[index - 1] = coefficient + o.degree = o.degree + Integer.one() + end + else + -- Constructs the coefficients from an existing polynomial of coefficients + local loc = o.degree:asnumber() + while loc > 0 do + if not coefficients[loc] or coefficients[loc] == coefficients[loc]:zero() then + o.degree = o.degree - Integer.one() + else + break + end + loc = loc - 1 + end + + while loc >= 0 do + o.coefficients[loc] = coefficients[loc] + loc = loc - 1 + end + end + + -- Each value of the polynomial greater than its degree is implicitly zero + o.coefficients = setmetatable(o.coefficients, {__index = function (table, key) + return o:zeroc() + end}) + return o +end + +-- Returns the ring this object is an element of +function PolynomialRing:getring() + local t = {ring = PolynomialRing} + if self then + t.child = self.ring + t.symbol = self.symbol + end + t = setmetatable(t, __obj) + return t +end + +-- Explicitly converts this element to an element of another ring +function PolynomialRing:inring(ring) + + -- Faster equality check + if ring == self:getring() then + return self + end + + if ring == Rational:getring() and ring.symbol then + return Rational(self:inring(ring.child), self:inring(ring.child):one(), true) + end + + if ring.symbol == self.symbol then + local out = {} + for i = 0, self.degree:asnumber() do + out[i + 1] = self.coefficients[i]:inring(ring.child) + end + return PolynomialRing(out, self.symbol) + end + + -- TODO: Allow re-ordering of polynomial rings, so from R[x][y] -> R[y][x] for instance + if ring == PolynomialRing:getring() then + return PolynomialRing({self:inring(ring.child)}, ring.symbol) + end + + error("Unable to convert element to proper ring.") +end + + +-- Returns whether the ring is commutative +function PolynomialRing:iscommutative() + return true +end + +function PolynomialRing:add(b) + local larger + + if self.degree > b.degree then + larger = self + else + larger = b + end + + local new = {} + local loc = 0 + while loc <= larger.degree:asnumber() do + new[loc] = self.coefficients[loc] + b.coefficients[loc] + loc = loc + 1 + end + + return PolynomialRing(new, self.symbol, larger.degree) +end + +function PolynomialRing:neg() + local new = {} + local loc = 0 + while loc <= self.degree:asnumber() do + new[loc] = -self.coefficients[loc] + loc = loc + 1 + end + return PolynomialRing(new, self.symbol, self.degree) +end + +function PolynomialRing:mul(b) + -- Grade-school multiplication is actually faster up to a very large polynomial size due to Lua's overhead. + local new = {} + + local sd = self.degree:asnumber() + local bd = b.degree:asnumber() + + for i = 0, sd+bd do + new[i] = self:zeroc() + for j = math.max(0, i-bd), math.min(sd, i) do + new[i] = new[i] + self.coefficients[j]*b.coefficients[i-j] + end + end + return PolynomialRing(new, self.symbol, self.degree + b.degree) + -- return PolynomialRing(PolynomialRing.mul_rec(self.coefficients, b.coefficients), self.symbol, self.degree + b.degree) +end + +-- Performs Karatsuba multiplication without constructing new polynomials recursively +function PolynomialRing.mul_rec(a, b) + if #a==0 and #b==0 then + return {[0]=a[0] * b[0], [1]=Integer.zero()} + end + + local k = Integer.ceillog(Integer.max(Integer(#a), Integer(#b)) + Integer.one(), Integer(2)) + local n = Integer(2) ^ k + local m = n / Integer(2) + local nn = n:asnumber() + local mn = m:asnumber() + + local a0, a1, b0, b1 = {}, {}, {}, {} + + for e = 0, mn - 1 do + a0[e] = a[e] or Integer.zero() + a1[e] = a[e + mn] or Integer.zero() + b0[e] = b[e] or Integer.zero() + b1[e] = b[e + mn] or Integer.zero() + end + + local p1 = PolynomialRing.mul_rec(a1, b1) + local p2a = Copy(a0) + local p2b = Copy(b0) + for e = 0, mn - 1 do + p2a[e] = p2a[e] + a1[e] + p2b[e] = p2b[e] + b1[e] + end + local p2 = PolynomialRing.mul_rec(p2a, p2b) + local p3 = PolynomialRing.mul_rec(a0, b0) + local r = {} + for e = 0, mn - 1 do + p2[e] = p2[e] - p1[e] - p3[e] + r[e] = p3[e] + r[e + mn] = p2[e] + r[e + nn] = p1[e] + end + for e = mn, nn - 1 do + p2[e] = p2[e] - p1[e] - p3[e] + r[e] = r[e] + p3[e] + r[e + mn] = r[e + mn] + p2[e] + r[e + nn] = p1[e] + end + + return r +end + +-- Uses synthetic division. +function PolynomialRing:divremainder(b) + local n, m = self.degree:asnumber(), b.degree:asnumber() + + if m > n then + return self:zero(), self + end + + local o = Copy(self.coefficients) + local lc = b:lc() + for i = n, m, -1 do + o[i] = o[i] / lc + + if o[i] ~= self:zeroc() then + for j = 1, m do + o[i-j] = o[i-j] - b.coefficients[m - j] * o[i] + end + end + end + + local q = {} + local r = {} + for i = 0, m-1 do + r[i] = o[i] + end + + r[0] = r[0] or self:zeroc() + + for i = m, #o do + q[i - m] = o[i] + end + + return PolynomialRing(q, self.symbol, self.degree), PolynomialRing(r, self.symbol, Integer.max(Integer.zero(), b.degree-Integer.one())) +end + +-- Performs polynomial pseudodivision of this polynomial by another in the same ring, +-- and returns both the pseudoquotient and pseudoremainder. +-- In the case where both coefficients are fields, this is equivalent to division with remainder. +function PolynomialRing:pseudodivide(b) + + local p = self:zero() + local s = self + local m = s.degree + local n = b.degree + local delta = Integer.max(m - n + Integer.one(), Integer.zero()) + + local lcb = b:lc() + local sigma = Integer.zero() + + while m >= n and s ~= Integer.zero() do + local lcs = s:lc() + p = p * lcb + self:one():multiplyDegree((m-n):asnumber()) * lcs + s = s * lcb - b * self:one():multiplyDegree((m-n):asnumber()) * lcs + sigma = sigma + Integer.one() + m = s.degree + end + + if delta - sigma == Integer.zero() then + return p,s + else + return lcb^(delta - sigma) * p, lcb^(delta - sigma) * s + end +end + +-- Polynomial rings are never fields, but when dividing by a polynomial by a constant we may want to use / instead of // +function PolynomialRing:div(b) + return self:divremainder(b) +end + +function PolynomialRing:zero() + return self.coefficients[0]:zero():inring(self:getring()) +end + +function PolynomialRing:zeroc() + return self.coefficients[0]:zero() +end + +function PolynomialRing:one() + return self.coefficients[0]:one():inring(self:getring()) +end + +function PolynomialRing:onec() + return self.coefficients[0]:one() +end + +function PolynomialRing:eq(b) + for i=0,math.max(self.degree:asnumber(), b.degree:asnumber()) do + if self.coefficients[i] ~= b.coefficients[i] then + return false + end + end + return true +end + +-- Returns the leading coefficient of this polynomial +function PolynomialRing:lc() + return self.coefficients[self.degree:asnumber()] +end + +--- @return boolean +function PolynomialRing:isconstant() + return false +end + +-- This expression is free of a symbol if and only if the symbol is not the symbol used to create the ring. +function PolynomialRing:freeof(symbol) + return symbol.symbol ~= self.symbol +end + +-- Replaces each expression in the map with its value. +function PolynomialRing:substitute(map) + return self:tocompoundexpression():substitute(map) +end + +-- Expands a polynomial expression. Polynomials are already in expanded form, so we just need to autosimplify. +function PolynomialRing:expand() + return self:tocompoundexpression():autosimplify() +end + +function PolynomialRing:autosimplify() + return self:tocompoundexpression():autosimplify() +end + +-- Transforms from array format to an expression format. +function PolynomialRing:tocompoundexpression() + local terms = {} + for exponent, coefficient in pairs(self.coefficients) do + terms[exponent + 1] = BinaryOperation(BinaryOperation.MUL, {coefficient:tocompoundexpression(), + BinaryOperation(BinaryOperation.POW, {SymbolExpression(self.symbol), Integer(exponent)})}) + end + return BinaryOperation(BinaryOperation.ADD, terms) +end + +-- Uses Horner's rule to evaluate a polynomial at a point +function PolynomialRing:evaluateat(x) + local out = self:zeroc() + for i = self.degree:asnumber(), 1, -1 do + out = out + self.coefficients[i] + out = out * x + end + return out + self.coefficients[0] +end + +-- Multiplies this polynomial by x^n +function PolynomialRing:multiplyDegree(n) + local new = {} + for e = 0, n-1 do + new[e] = self:zeroc() + end + local loc = n + while loc <= self.degree:asnumber() + n do + new[loc] = self.coefficients[loc - n] + loc = loc + 1 + end + return PolynomialRing(new, self.symbol, self.degree + Integer(n)) +end + +-- Returns the formal derivative of this polynomial +function PolynomialRing:derivative() + if self.degree == Integer.zero() then + return PolynomialRing({self:zeroc()}, self.symbol, Integer(-1)) + end + local new = {} + for e = 1, self.degree:asnumber() do + new[e - 1] = Integer(e) * self.coefficients[e] + end + return PolynomialRing(new, self.symbol, self.degree - Integer.one()) +end + +-- Returns the square-free factorization of a polynomial +function PolynomialRing:squarefreefactorization() + local terms + if self.ring == Rational.getring() or self.ring == Integer.getring() then + terms = self:rationalsquarefreefactorization() + elseif self.ring == IntegerModN.getring() then + if not self.ring.modulus:isprime() then + error("Cannot compute a square-free factorization of a polynomial ring contructed from a ring that is not a field.") + end + terms = self:modularsquarefreefactorization() + end + + local expressions = {self:lc()} + local j = 1 + for index, term in ipairs(terms) do + if term.degree ~= Integer.zero() or term.coefficients[0] ~= Integer.one() then + j = j + 1 + expressions[j] = BinaryOperation.POWEXP({term, Integer(index)}) + end + end + + return BinaryOperation.MULEXP(expressions) +end + +-- Factors a polynomial into irreducible terms +function PolynomialRing:factor() + -- Square-free factorization over an integral domain (so a polynomial ring constructed from a field) + local squarefree = self:squarefreefactorization() + local squarefreeterms = {} + local result = {squarefree.expressions[1]} + for i, expression in ipairs(squarefree.expressions) do + if i > 1 then + -- Converts square-free polynomials with rational coefficients to integer coefficients so Rational Roots / Zassenhaus can factor them + if expression.expressions[1].ring == Rational.getring() then + local factor, integerpoly = expression.expressions[1]:rationaltointeger() + result[1] = result[1] * factor ^ expression.expressions[2] + squarefreeterms[i - 1] = integerpoly + else + squarefreeterms[i - 1] = expression.expressions[1] + end + end + end + + for i, expression in ipairs(squarefreeterms) do + local terms + if expression.ring == Integer.getring() then + -- Factoring over the integers first uses the rational roots test to factor out monomials (for efficiency purposes) + local remaining, factors = expression:rationalroots() + terms = factors + -- Then applies the Zassenhaus algorithm if there entire polynomial has not been factored into monomials + if remaining ~= Integer.one() then + remaining = remaining:zassenhausfactor() + for _, exp in ipairs(remaining) do + terms[#terms+1] = exp + end + end + end + if expression.ring == IntegerModN.getring() then + -- Berlekamp factorization is used for rings with integers mod a prime as coefficients + terms = expression:berlekampfactor() + end + for _, factor in ipairs(terms) do + result[#result+1] = BinaryOperation.POWEXP({factor, squarefree.expressions[i + 1].expressions[2]}) + end + end + return BinaryOperation.MULEXP(result) +end + +-- Uses the Rational Root test to factor out monomials of a square-free polynomial. +function PolynomialRing:rationalroots() + local remaining = self + local roots = {} + if self.coefficients[0] == Integer.zero() then + roots[1] = PolynomialRing({Integer.zero(), Integer.one()}, self.symbol) + remaining = remaining // roots[1] + end + -- This can be slower than Zassenhaus if the digits are large enough, since factoring integers is slow + -- if self.coefficients[0] > Integer(Integer.DIGITSIZE - 1) or self:lc() > Integer(Integer.DIGITSIZE - 1) then + -- return remaining, roots + -- end + while remaining ~= Integer.one() do + :: nextfactor :: + local a = remaining.coefficients[0] + local b = remaining:lc() + local afactors = a:divisors() + local bfactors = b:divisors() + for _, af in ipairs(afactors) do + for _, bf in ipairs(bfactors) do + local testroot = Rational(af, bf, true) + if remaining:evaluateat(testroot) == Integer.zero() then + roots[#roots+1] = PolynomialRing({-testroot.numerator, testroot.denominator}, self.symbol) + remaining = remaining // roots[#roots] + goto nextfactor + end + if remaining:evaluateat(-testroot) == Integer.zero() then + roots[#roots+1] = PolynomialRing({testroot.numerator, testroot.denominator}, self.symbol) + remaining = remaining // roots[#roots] + goto nextfactor + end + end + end + break + end + + return remaining, roots +end + +-- Returns a list of roots of the polynomial, simplified up to cubics. +function PolynomialRing:roots() + local roots = {} + local factorization = self:factor() + + for i, factor in ipairs(factorization.expressions) do + if i > 1 then + local decomp = factor.expressions[1]:decompose() + for _, poly in ipairs(decomp) do + if poly.degree > Integer(3) then + table.insert(roots,RootExpression(factor.expressions[1])) + goto nextfactor + end + end + local factorroots = RootExpression(decomp[#decomp]):autosimplify() + if factorroots == true then + return true + end + if factorroots == false then + goto nextfactor + end + local replaceroots = {} + for j = #decomp - 1,1,-1 do + for _, root in ipairs(factorroots) do + local temp = RootExpression(decomp[j]):autosimplify(root) + if temp == true then + return true + end + if factorroots == false then + goto nextfactor + end + replaceroots = JoinArrays(replaceroots, temp) + end + factorroots = replaceroots + end + roots = JoinArrays(roots, factorroots) + end + end + ::nextfactor:: + return roots +end + +----------------- +-- Inheritance -- +----------------- + +__PolynomialRing.__index = Ring +__PolynomialRing.__call = PolynomialRing.new +PolynomialRing = setmetatable(PolynomialRing, __PolynomialRing)
\ No newline at end of file diff --git a/macros/luatex/latex/luacas/tex/algebra/polynomialring/berlekampfactoring.lua b/macros/luatex/latex/luacas/tex/algebra/polynomialring/berlekampfactoring.lua new file mode 100644 index 0000000000..feabb61a3f --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/polynomialring/berlekampfactoring.lua @@ -0,0 +1,187 @@ +-- Methods related to the Berlekamp factoring algorithm. + +-- Square-free factorization in the modular field Zp. +function PolynomialRing:modularsquarefreefactorization() + local monic = self / self:lc() + local terms = {} + terms[0] = PolynomialRing.gcd(monic, monic:derivative()) + local b = monic // terms[0] + local c = monic:derivative() // terms[0] + local d = c - b:derivative() + local i = 1 + while b ~= Integer.one() do + terms[i] = PolynomialRing.gcd(b, d) + b, c = b // terms[i], d // terms[i] + i = i + 1 + d = c - b:derivative() + end + + if not (terms[i-1]:derivative().degree == Integer.zero() and terms[i-1]:derivative().coefficients[0] == Integer.zero()) then + return terms + end + + local recursiveterms = terms[i-1]:collapseterms(self.ring.modulus):modularsquarefreefactorization() + for k, poly in ipairs(recursiveterms) do + recursiveterms[k] = poly:expandterms(self.ring.modulus) + end + return JoinArrays(terms, recursiveterms) +end + +-- Returns a new polnomial consisting of every nth term of the old one - helper method for square-free factorization +function PolynomialRing:collapseterms(n) + local new = {} + local loc = 0 + local i = 0 + local nn = n:asnumber() + while loc <= self.degree:asnumber() do + new[i] = self.coefficients[loc] + loc = loc + nn + i = i + 1 + end + + return PolynomialRing(new, self.symbol, self.degree // n) +end + +-- Returns a new polnomial consisting of every nth term of the old one - helper method for square-free factorization +function PolynomialRing:expandterms(n) + local new = {} + local loc = 0 + local i = 0 + local nn = n:asnumber() + while i <= self.degree:asnumber() do + new[loc] = self.coefficients[i] + for j = 1, nn do + new[loc + j] = IntegerModN(Integer.zero(), n) + end + loc = loc + nn + i = i + 1 + end + + return PolynomialRing(new, self.symbol, self.degree * n) +end + +-- Uses Berlekamp's Algorithm to factor polynomials in mod p +function PolynomialRing:berlekampfactor() + if self.degree == 0 or self.degree == 1 then + return {self} + end + + local R = self:RMatrix() + local S = self:auxillarybasis(R) + if #S == 1 then + return {self} + end + return self:findfactors(S) +end + +-- Gets the R Matrix for Berlekamp factorization +function PolynomialRing:RMatrix() + local R = {} + for i = 1, self.degree:asnumber() do + R[i] = {} + end + for i = 0, self.degree:asnumber()-1 do + local remainder = PolynomialRing({IntegerModN(Integer.one(), self.ring.modulus)}, self.symbol):multiplyDegree(self.ring.modulus:asnumber()*i) % self + for j = 0, self.degree:asnumber()-1 do + R[j + 1][i + 1] = remainder.coefficients[j] + if j == i then + R[j + 1][i + 1] = R[j + 1][i + 1] - IntegerModN(Integer.one(), self.ring.modulus) + end + end + end + return R +end + +-- Creates an auxillary basis using the R matrix +function PolynomialRing:auxillarybasis(R) + local P = {} + local n = self.degree:asnumber() + for i = 1, n do + P[i] = 0 + end + S = {} + local q = 1 + for j = 1, n do + local i = 1 + local pivotfound = false + while not pivotfound and i <= n do + if R[i][j] ~= self:zeroc() and P[i] == 0 then + pivotfound = true + else + i = i + 1 + end + end + if pivotfound then + P[i] = j + local a = R[i][j]:inv() + for l = 1, n do + R[i][l] = a * R[i][l] + end + for k = 1, n do + if k ~= i then + local f = R[k][j] + for l = 1, n do + R[k][l] = R[k][l] - f*R[i][l] + end + end + end + else + local s = {} + s[j] = self:onec() + for l = 1, j - 1 do + local e = 0 + i = 1 + while e == 0 and i <= n do + if l == P[i] then + e = i + else + i = i + 1 + end + end + if e > 0 then + local c = -R[e][j] + s[l] = c + else + s[l] = self:zeroc() + end + end + S[#S+1] = PolynomialRing(s, self.symbol) + end + end + return S +end + +-- Uses the auxilary basis to find the irreirrducible factors of the polynomial. +function PolynomialRing:findfactors(S) + local r = #S + local p = self.ring.modulus + local factors = {self} + for k = 2,r do + local b = S[k] + local old_factors = Copy(factors) + for i = 1,#old_factors do + local w = old_factors[i] + local j = 0 + while j <= p:asnumber() - 1 do + local g = PolynomialRing.gcd(b-IntegerModN(Integer(j), p), w) + if g == Integer.one() then + j = j + 1 + elseif g == w then + j = p:asnumber() + else + factors = Remove(factors, w) + local q = w // g + factors[#factors+1] = g + factors[#factors+1] = q + if #factors == r then + return factors + else + j = j + 1 + w = q + end + end + + end + end + end +end
\ No newline at end of file diff --git a/macros/luatex/latex/luacas/tex/algebra/polynomialring/decomposition.lua b/macros/luatex/latex/luacas/tex/algebra/polynomialring/decomposition.lua new file mode 100644 index 0000000000..84f2af1efc --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/polynomialring/decomposition.lua @@ -0,0 +1,93 @@ +-- Methods related to polynomial decomposition. + +-- Returns a list of polynomials that form a complete decomposition of a polynomial. +function PolynomialRing:decompose() + local U = self - self.coefficients[0] + local S = U:divisors() + local decomposition = {} + local C = PolynomialRing({Integer.zero(), Integer.one()}, self.symbol) + local finalcomponent + + while S[1] do + local w = S[1] + for _, poly in ipairs(S) do + if poly.degree < w.degree then + w = poly + end + end + S = Remove(S, w) + if C.degree < w.degree and w.degree < self.degree and self.degree % w.degree == Integer.zero() then + local g = w:polyexpand(C, self.symbol) + local R = self:polyexpand(w, self.symbol) + if g.degree == Integer.zero() and R.degree == Integer.zero() then + g.symbol = self.symbol + decomposition[#decomposition+1] = g.coefficients[0] + decomposition[#decomposition].symbol = self.symbol + C = w + finalcomponent = R.coefficients[0] + end + end + end + + if not decomposition[1] then + return {self} + end + + finalcomponent.symbol = self.symbol + decomposition[#decomposition+1] = finalcomponent + return decomposition +end + +-- Returns a list of all monic divisors of positive degree of the polynomial, assuming the polynomial ring is a Euclidean Domain. +function PolynomialRing:divisors() + local factors = self:factor() + -- Converts each factor to a monic factor (we don't need to worry updating the constant term) + for i, factor in ipairs(factors.expressions) do + if i > 1 then + factor.expressions[1] = factor.expressions[1] / factor.expressions[1]:lc() + end + end + + local terms = {} + for i, _ in ipairs(factors.expressions) do + if i > 1 then + terms[i] = Integer.zero() + end + end + + local divisors = {} + local divisor = PolynomialRing({self:onec()}, self.symbol) + while true do + for i, factor in ipairs(factors.expressions) do + if i > 1 then + local base = factor.expressions[1] + local power = factor.expressions[2] + if terms[i] < power then + terms[i] = terms[i] + Integer.one() + divisor = divisor * base + break + else + terms[i] = Integer.zero() + divisor = divisor // (base ^ power) + end + end + end + if divisor == Integer.one() then + break + end + divisors[#divisors+1] = divisor + end + + return divisors + +end + +-- Polynomial expansion as a subroutine of decomposition. +function PolynomialRing:polyexpand(v, x) + local u = self + if u == Integer.zero() then + return Integer.zero() + end + local q,r = u:divremainder(v) + return PolynomialRing({PolynomialRing({Integer.zero(), Integer.one()}, "_")}, x) * q:polyexpand(v, x) + r +end
\ No newline at end of file diff --git a/macros/luatex/latex/luacas/tex/algebra/polynomialring/zassenhausfactoring.lua b/macros/luatex/latex/luacas/tex/algebra/polynomialring/zassenhausfactoring.lua new file mode 100644 index 0000000000..be19298e64 --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/polynomialring/zassenhausfactoring.lua @@ -0,0 +1,220 @@ +-- Methods related to the Zassenhaus factorization algorithm. + + +-- Square-free factorization in the rational field. +function PolynomialRing:rationalsquarefreefactorization(keeplc) + local monic = self / self:lc() + local terms = {} + terms[0] = PolynomialRing.gcd(monic, monic:derivative()) + local b = monic // terms[0] + local c = monic:derivative() // terms[0] + local d = c - b:derivative() + local i = 1 + while b.degree ~= Integer.zero() or b.coefficients[0] ~= Integer.one() do + terms[i] = PolynomialRing.gcd(b, d) + b, c = b // terms[i], d // terms[i] + i = i + 1 + d = c - b:derivative() + end + if keeplc and terms[1] then + terms[1] = terms[1] * self:lc() + end + return terms +end + +-- Factors the largest possible constant out of a polynomial whos underlying ring is a Euclidean domain but not a field +function PolynomialRing:factorconstant() + local gcd = Integer.zero() + for i = 0, self.degree:asnumber() do + gcd = self.ring.gcd(gcd, self.coefficients[i]) + end + if gcd == Integer.zero() then + return Integer.one(), self + end + return gcd, self / gcd +end + +-- Converts a polynomial in the rational polynomial ring to the integer polynomial ring +function PolynomialRing:rationaltointeger() + local lcm = Integer.one() + for i = 0, self.degree:asnumber() do + if self.coefficients[i]:getring() == Rational:getring() then + lcm = lcm * self.coefficients[i].denominator / Integer.gcd(lcm, self.coefficients[i].denominator) + end + end + return Integer.one() / lcm, self * lcm +end + +-- Uses Zassenhaus's Algorithm to factor sqaure-free polynomials over the intergers +function PolynomialRing:zassenhausfactor() + + -- Creates a monic polynomial V with related roots + local V = {} + local n = self.degree:asnumber() + local l = self:lc() + for i = 0, n - 1 do + V[i] = l ^ Integer(n - 1 - i) * self.coefficients[i] + end + V[n] = Integer.one() + V = PolynomialRing(V, "y", self.degree) + + -- Performs Berlekamp Factorization in a sutable prime base + local p = V:findprime() + local S = V:inring(PolynomialRing.R("y", p)):berlekampfactor() + + -- If a polynomial is irreducible with coefficients in mod p, it is also irreducible over the integers + if #S == 1 then + return {self} + end + + -- Performs Hensel lifting on the factors mod p + local k = V:findmaxlifts(p) + local W = V:henselift(S, k) + local M = {} + + -- Returns the solutions back to the original from the monic transformation + for i, factor in ipairs(W) do + local w = {} + for j = 0, factor.degree:asnumber() do + w[j] = factor.coefficients[j]:inring(Integer.getring()) * l ^ Integer(j) + end + _, M[i] = PolynomialRing(w, self.symbol, factor.degree):factorconstant() + end + + return M + +end + +-- Finds the smallest prime such that this polynomial with coefficients in mod p is square-free +function PolynomialRing:findprime() + + local smallprimes = {Integer(2), Integer(3), Integer(5), Integer(7), Integer(11), Integer(13), Integer(17), Integer(19), Integer(23), + Integer(29), Integer(31), Integer(37), Integer(41), Integer(43), Integer(47), Integer(53), Integer(59)} + + for _, p in pairs(smallprimes) do + local P = PolynomialRing({IntegerModN(Integer.one(), p)}, self.symbol) + local s = self:inring(P:getring()) + if PolynomialRing.gcd(s, s:derivative()) == P then + return p + end + end + + error("Execution error: No suitable prime found for factoring.") +end + +-- Finds the maximum number of times Hensel Lifting will be applied to raise solutions to the appropriate power +function PolynomialRing:findmaxlifts(p) + local n = self.degree:asnumber() + local h = self.coefficients[0] + for i=0 , n do + if self.coefficients[i] > h then + h = self.coefficients[i] + end + end + + local B = 2^n * math.sqrt(n) * h:asnumber() + return Integer(math.ceil(math.log(2*B, p:asnumber()))) +end + +-- Uses Hensel lifting on the factors of a polynomial S mod p to find them in the integers +function PolynomialRing:henselift(S, k) + local p = S[1].ring.modulus + if k == Integer.one() then + return self:truefactors(S, p, k) + end + G = self:genextendsigma(S) + local V = S + for j = 2, k:asnumber() do + local Vp = V[1]:inring(PolynomialRing.R("y")) + for i = 2, #V do + Vp = Vp * V[i]:inring(PolynomialRing.R("y")) + end + local E = self - Vp:inring(PolynomialRing.R("y")) + if E == Integer.zero() then + return V + end + E = E:inring(PolynomialRing.R("y", p ^ Integer(j))):inring(PolynomialRing.R("y")) + F = E / p ^ (Integer(j) - Integer.one()) + R = self:genextendR(V, G, F) + local Vnew = {} + for i, v in ipairs(V) do + local vnew = v:inring(PolynomialRing.R("y", p ^ Integer(j))) + local rnew = R[i]:inring(PolynomialRing.R("y", p ^ Integer(j))) + Vnew[i] = vnew + (p) ^ (Integer(j) - Integer.one()) * rnew + end + V = Vnew + end + return self:truefactors(V, p, k) +end + +-- Gets a list of sigma polynomials for use in hensel lifting +function PolynomialRing:genextendsigma(S) + local v = S[1] * S[2] + local _, A, B = PolynomialRing.extendedgcd(S[2], S[1]) + local SIGMA = {A, B} + for i, _ in ipairs(S) do + if i >= 3 then + v = v * S[i] + local sum = SIGMA[1] * (v // S[1]) + for j = 2, i-1 do + sum = sum + SIGMA[j] * (v // S[j]) + end + _, A, B = PolynomialRing.extendedgcd(sum, v // S[i]) + for j = 1, i-1 do + SIGMA[j] = SIGMA[j] * A + end + SIGMA[i] = B + end + end + + return SIGMA +end + +-- Gets a list of r polynomials for use in hensel lifting +function PolynomialRing:genextendR(V, G, F) + R = {} + for i, v in ipairs(V) do + local pring = G[1]:getring() + R[i] = F:inring(pring) * G[i] % v:inring(pring) + end + return R +end + +-- Updates factors of the polynomial to the correct ones in the integer ring +function PolynomialRing:truefactors(l, p, k) + local U = self + local L = l + local factors = {} + local m = 1 + while m <= #L / 2 do + local C = Subarrays(L, m) + while #C > 0 do + local t = C[1] + local prod = t[1] + for i = 2, #t do + prod = prod * t[i] + end + local T = prod:inring(PolynomialRing.R("y", p ^ k)):inring(PolynomialRing.R("y")) + -- Convert to symmetric representation - this is the only place it actually matters + for i = 0, T.degree:asnumber() do + if T.coefficients[i] > p ^ k / Integer(2) then + T.coefficients[i] = T.coefficients[i] - p^k + end + end + local Q, R = U:divremainder(T) + if R == Integer.zero() then + factors[#factors+1] = T + U = Q + L = RemoveAll(L, t) + C = RemoveAny(C, t) + else + C = Remove(C, t) + end + end + m = m + 1 + end + if U ~= Integer.one() then + factors[#factors+1] = U + end + return factors +end
\ No newline at end of file diff --git a/macros/luatex/latex/luacas/tex/algebra/rational.lua b/macros/luatex/latex/luacas/tex/algebra/rational.lua new file mode 100644 index 0000000000..811909a948 --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/rational.lua @@ -0,0 +1,241 @@ +--- @class Rational +--- Represents an element of the field of rational numbers or rational functions. +--- @field numerator Ring +--- @field denominator Ring +--- @field ring RingIdentifier +Rational = {} +__Rational = {} + + +-------------------------- +-- Static functionality -- +-------------------------- + +-- Metatable for ring objects. +local __obj = {__index = Rational, __eq = function(a, b) + return a["ring"] == b["ring"] and + (a["child"] == b["child"] or a["child"] == nil or b["child"] == nil) and + (a["symbol"] == b["symbol"] or a["child"] == nil or b["child"] == nil) +end, __tostring = function(a) + if a.symbol then + return tostring(a.child.child) .. "(" .. a.symbol .. ")" + end + if a.child then + return "QQ" + end + return "(Generic Fraction Field)" + end} + +--- @param symbol SymbolExpression +--- @param child RingIdentifier +--- @return RingIdentifier +function Rational.makering(symbol, child) + local t = {ring = Rational} + t.symbol = symbol + t.child = child + t = setmetatable(t, __obj) + return t +end + +--- Converts a string of the form -?[0-9]+ or -?[0-9]+\/[0-9]+ to a rational number. +--- @param str string +--- @return Rational|Integer +function Rational.fromstring(str) + local divloc = string.find(str, "/"); + if not divloc then + return Integer(str) + end + return Rational(Integer(string.sub(str, 1, divloc - 1)), Integer(string.sub(str, divloc + 1, #str))) +end + + +---------------------------- +-- Instance functionality -- +---------------------------- + +-- So we don't have to copy the field operations each time. +local __o = Copy(__FieldOperations) +__o.__index = Rational +__o.__tostring = function(a) + if a.ring.symbol then + return "(" .. tostring(a.numerator)..")/("..tostring(a.denominator) .. ")" + end + return tostring(a.numerator).."/"..tostring(a.denominator) +end + +--- Creates a new rational given a numerator and denominator that are part of the same ring. +--- Rational numbers are represented uniquely. +--- @param n Ring +--- @param d Ring +--- @param keep boolean +function Rational:new(n, d, keep) + local o = {} + o = setmetatable(o, __o) + + if n:getring() == PolynomialRing.getring() then + o.symbol = n.symbol + end + + if d:getring() == PolynomialRing.getring() then + o.symbol = d.symbol + end + + if d == Integer(0) then + error("Arithmetic error: division by zero") + end + + n = n or Integer.zero() + d = d or Integer.one() + o.numerator = n + o.denominator = d + o:reduce() + + if (not keep) and o.denominator == Integer.one() or (not keep) and o.numerator == Integer.zero() then + return o.numerator + end + + return o +end + +--- Reduces a rational expression to standard form. This method mutates its object. +function Rational:reduce() + if self.numerator:getring() == Integer.getring() then + if self.denominator < Integer.zero() then + self.denominator = -self.denominator + self.numerator = -self.numerator + end + local gcd = Integer.gcd(self.numerator, self.denominator) + self.numerator = self.numerator//gcd + self.denominator = self.denominator//gcd + self.ring = Integer.getring() + elseif self.numerator:getring() == PolynomialRing.getring() then + local lc = self.denominator:lc() + self.denominator = self.denominator/lc + self.numerator = self.numerator/lc + local gcd = PolynomialRing.gcd(self.numerator, self.denominator) + self.numerator = self.numerator//gcd + self.denominator = self.denominator//gcd + self.ring = Ring.resultantring(self.numerator:getring(), self.denominator:getring()) + end +end + + +--- @return RingIdentifier +function Rational:getring() + local t = {ring=Rational} + if self then + t.child = self.ring + t.symbol = self.symbol + end + t = setmetatable(t, __obj) + return t +end + +--- @param ring RingIdentifier +--- @return Ring +function Rational:inring(ring) + if ring == self:getring() then + return self + end + + if ring == Rational:getring() and ring.symbol then + if not self:getring().symbol then + return Rational(self:inring(ring.child), self:inring(ring.child):one(), true) + end + return Rational(self.numerator:inring(ring.child), self.denominator:inring(ring.child), true) + end + + if ring == PolynomialRing:getring() then + return PolynomialRing({self:inring(ring.child)}, ring.symbol) + end + + error("Unable to convert element to proper ring.") +end + +--- @return boolean +function Rational:isconstant() + if self.symbol then + return false + end + return true +end + +--- @return Expression +function Rational:tocompoundexpression() + return BinaryOperation(BinaryOperation.DIV, {self.numerator:tocompoundexpression(), self.denominator:tocompoundexpression()}) +end + +--- Returns this rational as a floating point number. Can only approximate the value of most rationals. +--- @return number +function Rational:asnumber() + return self.numerator:asnumber() / self.denominator:asnumber() +end + +function Rational:add(b) + return Rational(self.numerator * b.denominator + self.denominator * b.numerator, self.denominator * b.denominator) +end + +function Rational:neg() + return Rational(-self.numerator, self.denominator, true) +end + +function Rational:mul(b) + return Rational(self.numerator * b.numerator, self.denominator * b.denominator) +end + +-- function Rational:inv(b) +-- return Rational(self.numerator * b.numerator, self.denominator * b.denominator) +-- end + +function Rational:pow(b) + return (self.numerator ^ b) / (self.denominator ^ b) +end + +function Rational:div(b) + return Rational(self.numerator * b.denominator, self.denominator * b.numerator) +end + +function Rational:eq(b) + return self.numerator == b.numerator and self.denominator == b.denominator +end + +function Rational:lt(b) + if self.numerator < Integer.zero() and b.numerator > Integer.zero() then + return true + end + if self.numerator > Integer.zero() and b.numerator < Integer.zero() then + return false + end + + if (self.numerator >= Integer.zero() and b.numerator >= Integer.zero()) or (self.numerator <= Integer.zero() and b.numerator <= Integer.zero()) then + return self.numerator * b.denominator < self.denominator * b.numerator + end + return self.numerator * b.denominator > self.denominator * b.numerator +end + +function Rational:le(b) + return self:eq(b) or self:lt(b) +end + +function Rational:zero() + return Integer.zero() +end + +function Rational:one() + return Integer.one() +end + +function Rational:tolatex() + if string.sub(self.numerator:tolatex(),1,1) == '-' then + return "- \\frac{" .. string.sub(self.numerator:tolatex(),2,-1) .. "}{" .. self.denominator:tolatex() .. "}" + end + return "\\frac{" .. self.numerator:tolatex() .."}{".. self.denominator:tolatex().. "}" +end + +----------------- +-- Inheritance -- +----------------- + +__Rational.__index = Field +__Rational.__call = Rational.new +Rational = setmetatable(Rational, __Rational)
\ No newline at end of file diff --git a/macros/luatex/latex/luacas/tex/algebra/ring.lua b/macros/luatex/latex/luacas/tex/algebra/ring.lua new file mode 100644 index 0000000000..4903b9c912 --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/ring.lua @@ -0,0 +1,326 @@ +--- @class Ring +--- Interface for an element of a ring with unity. +Ring = {} +__Ring = {} + +-------------------------- +-- Static functionality -- +-------------------------- + +--- Determines which ring the output of a binary operation with inputs in ring1 and ring2 should be, if such a ring exists. +--- If one of the rings is a subring of another ring, the result should be one of the two rings. +--- @param ring1 RingIdentifier +--- @param ring2 RingIdentifier +--- @return RingIdentifier +function Ring.resultantring(ring1, ring2) + if ring1 == ring2 then + return ring1 + end + + if ((ring1 == PolynomialRing.getring() and ring2 == Rational.getring()) or + (ring2 == PolynomialRing.getring() and ring1 == Rational.getring())) + and ring1.symbol == ring2.symbol then + return Rational.makering(ring1.symbol, Ring.resultantring(ring1.child, ring2.child)) + end + + if ring1 == PolynomialRing.getring() or ring2 == PolynomialRing.getring() then + if ring1 == ring2.child then + return ring2 + end + if ring2 == ring1.child then + return ring1 + end + + if ring1 == PolynomialRing.getring() and ring2 == PolynomialRing.getring() and ring1.symbol == ring2.symbol then + return PolynomialRing.makering(ring1.symbol, Ring.resultantring(ring1.child, ring2.child)) + end + + -- If none of the above conditions are satisfied, recusion is a pain, so we just strip all of the variables off of both rings. + -- TODO: Make this properly recursive, or just use a multivariable polynomial ring class + local symbols = {} + while ring1 == PolynomialRing.getring() do + symbols[#symbols+1] = ring1.symbol + ring1 = ring1.child + end + while ring2 == PolynomialRing.getring() do + if not Contains(symbols, ring2.symbol) then + symbols[#symbols+1] = ring2.symbol + end + ring2 = ring2.child + end + local ring = Ring.resultantring(ring1, ring2) + + if ring == Rational.getring() and Contains(symbols, ring.symbol) then + symbols = Remove(symbols, ring.symbol) + end + for i = #symbols, 1, -1 do + ring = PolynomialRing.makering(symbols[i], ring) + end + return ring + end + + if ring1 == Integer.getring() then + if ring2 == Integer.getring() then + return ring2 + end + + if ring2 == Rational.getring() then + return ring2 + end + + if ring2 == IntegerModN.getring() then + return ring2 + end + end + + if ring1 == Rational.getring() then + if ring2 == Integer.getring() then + return ring1 + end + + if ring2 == Rational.getring() then + if not ring1.symbol then + return Rational.makering(ring2.symbol, Ring.resultantring(ring1, ring2.child)) + end + if not ring2.symbol then + return Rational.makering(ring1.symbol, Ring.resultantring(ring1.child, ring2)) + end + if ring1.symbol and ring2.symbol and ring1.symbol == ring2.symbol then + return Rational.makering(ring1.symbol, Ring.resultantring(ring1.child, ring2.child)) + end + return ring2 + end + + if ring2 == IntegerModN.getring() then + return nil + end + end + + if ring1 == IntegerModN.getring() then + if ring2 == Integer.getring() then + return ring1 + end + + if ring2 == Rational.getring() then + return nil + end + + if ring2 == IntegerModN.getring() then + return IntegerModN.makering(Integer.gcd(ring1.modulus, ring2.modulus)) + end + end + + return nil +end + +--- Returns a particular instantiation of a ring. +--- Does the same thing as getring() if there is only one possible ring for a class, i.e., the integers and rationals. +--- @return RingIdentifier +function Ring.makering() + error("Called unimplemented method : makering()") +end + +---------------------- +-- Required methods -- +---------------------- + +--- Returns the ring this element is part of. +--- @return RingIdentifier +function Ring:getring() + error("Called unimplemented method : getring()") +end + +--- Explicitly converts this element to an element of another ring. +--- @param ring RingIdentifier +--- @return Ring +function Ring:inring(ring) + error("Called unimplemented method : in()") +end + +--- Returns whether the ring is commutative. +--- @return boolean +function Ring:iscommutative() + error("Called unimplemented method : iscommutative()") +end + +--- @return Ring +function Ring:add(b) + error("Called unimplemented method : add()") +end + +--- @return Ring +function Ring:sub(b) + return(self:add(b:neg())) +end + +--- @return Ring +function Ring:neg() + error("Called unimplemented method : neg()") +end + +--- @return Ring +function Ring:mul(b) + error("Called unimplemented method : mul()") +end + +--- Ring exponentiation by definition. Specific rings may implement more efficient methods. +--- @return Ring +function Ring:pow(n) + if(n < Integer.zero()) then + error("Execution error: Negative exponentiation is undefined over general rings") + end + local k = Integer.zero() + local b = self:one() + while k < n do + b = b * self + k = k + Integer.one() + end + return b +end + +--- @return boolean +function Ring:eq(b) + error("Execution error: Ring does not have a total order") +end + +--- @return boolean +function Ring:lt(b) + error("Execution error: Ring does not have a total order") +end + +--- @return boolean +function Ring:le(b) + error("Execution error: Ring does not have a total order") +end + +--- The additive identitity of the ring. +--- @return Ring +function Ring:zero() + error("Called unimplemented method : zero()") +end + +--- The multiplicative identitity of the ring. +--- @return Ring +function Ring:one() + error("Called unimplemented method : one()") +end + +-------------------------- +-- Instance metamethods -- +-------------------------- +__RingOperations = {} + +-- Each of these methods just handles coverting each element in the ring to an instance of the proper ring, if possible, +-- then passing the arguments to the function in a specific ring. + +__RingOperations.__unm = function(a) + return a:neg() +end + +__RingOperations.__add = function(a, b) + if not b.getring then + return BinaryOperation.ADDEXP({a, b}) + end + + local aring, bring = a:getring(), b:getring() + local oring = Ring.resultantring(aring, bring) + if not oring then + error("Attempted to add two elements of incompatable rings") + end + return a:inring(oring):add(b:inring(oring)) +end + +__RingOperations.__sub = function(a, b) + if not b.getring then + return BinaryOperation.SUBEXP({a, b}) + end + + local aring, bring = a:getring(), b:getring() + local oring = Ring.resultantring(aring, bring) + if not oring then + error("Attempted to subtract two elements of incompatable rings") + end + return a:inring(oring):sub(b:inring(oring)) +end + +-- Allows for multiplication by writing two expressions next to each other. +__RingOperations.__call = function (a, b) + return a * b +end + +__RingOperations.__mul = function(a, b) + if not b.getring then + return BinaryOperation.MULEXP({a, b}) + end + + local aring, bring = a:getring(), b:getring() + local oring = Ring.resultantring(aring, bring) + if not oring then + error("Attempted to muliply two elements of incompatable rings") + end + return a:inring(oring):mul(b:inring(oring)) +end + +__RingOperations.__pow = function(a, n) + if (not n.getring) or (n.getring and n:getring().ring ~= Integer) then + return BinaryOperation.POWEXP({a, n}) + end + + -- if a == a:zero() and n == Integer.zero() then + -- error("Cannot raise 0 to the power of 0") + -- end + + return a:pow(n) +end + +-- Comparison operations assume, of course, that the ring operation is equipped with a total order +-- All elements of all rings need these metamethods, since in Lua comparisons on tables only fire if both objects have the table +__RingOperations.__eq = function(a, b) + -- This shouldn't be needed, since __eq should only fire if both metamethods have the same function, but for some reason Lua always runs this anyway + if not a.getring or not b.getring then + return false + end + local aring, bring = a:getring(), b:getring() + if aring == bring then + return a:eq(b) + end + local oring = Ring.resultantring(aring, bring) + if not oring then + error("Attempted to compare two elements of incompatable rings") + end + return a:inring(oring):eq(b:inring(oring)) +end + +__RingOperations.__lt = function(a, b) + local aring, bring = a:getring(), b:getring() + if aring == bring then + return a:lt(b) + end + local oring = Ring.resultantring(aring, bring) + if not oring then + error("Attempted to compare two elements of incompatable rings") + end + return a:inring(oring):lt(b:inring(oring)) +end + +__RingOperations.__le = function(a, b) + local aring, bring = a:getring(), b:getring() + if aring == bring then + return a:le(b) + end + local oring = Ring.resultantring(aring, bring) + if not oring then + error("Attempted to compare two elements of incompatable rings") + end + return a:inring(oring):le(b:inring(oring)) +end + +----------------- +-- Inheritance -- +----------------- + +__Ring.__index = ConstantExpression +Ring = setmetatable(Ring, __Ring) + +--- Used for comparing and converting between rings. +--- @class RingIdentifier diff --git a/macros/luatex/latex/luacas/tex/algebra/rootexpression.lua b/macros/luatex/latex/luacas/tex/algebra/rootexpression.lua new file mode 100644 index 0000000000..a66182579a --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/rootexpression.lua @@ -0,0 +1,135 @@ +--- @class RootExpression +--- An expression that represents the solutions to expression = 0. +--- @field expression Expression +RootExpression = {} +__RootExpression = {} + +---------------------------- +-- Instance functionality -- +---------------------------- + +--- Creates a new root expression with the given expression. +--- @param expression Expression +--- @return RootExpression +function RootExpression:new(expression) + local o = {} + local __o = Copy(__ExpressionOperations) + + o.expression = Copy(expression) + + __o.__index = RootExpression + __o.__tostring = function(a) + return 'Root Of: (' .. tostring(a.expression) .. ')' + end + __o.__eq = function(a, b) + -- This shouldn't be needed, since __eq should only fire if both metamethods have the same function, but for some reason Lua always rungs this anyway + if not b:type() == RootExpression then + return false + end + return a.expression == b.expression + end + o = setmetatable(o, __o) + + return o +end + +--- @return Expression +function RootExpression:autosimplify(subpart) + local simplified = self.expression:autosimplify() + local simplified, ispoly = simplified:topolynomial() + + if simplified:isconstant() then + -- 0 = 0 is always true (obviously). + return simplified == simplified:zero() + end + + if ispoly then + if simplified.degree == Integer.zero() then + return simplified == simplified:zero() + end + if simplified.degree == Integer.one() then + return {-simplified.coefficients[0] / simplified.coefficients[1]} + end + if simplified.degree == Integer(2) then + local a = simplified.coefficients[2] + local b = simplified.coefficients[1] + local c = simplified.coefficients[0] + -- This is a hack until we can get more expression manipulation working, but that's okay. + if subpart then + c = (c - subpart):autosimplify() + end + return {((-b + sqrt(b^Integer(2) - Integer(4) * a * c)) / (Integer(2) * a)):autosimplify(), + ((-b - sqrt(b^Integer(2) - Integer(4) * a * c)) / (Integer(2) * a)):autosimplify()} + end + if simplified.degree == Integer(3) then + local a = simplified.coefficients[3] + local b = simplified.coefficients[2] + local c = simplified.coefficients[1] + local d = simplified.coefficients[0] + -- This is a hack until we can get more expression manipulation working, but that's okay. + if subpart then + d = (d - subpart):autosimplify() + end + + local delta0 = (b^Integer(2) - Integer(3)*a*c):autosimplify() + local delta1 = (Integer(2) * b^Integer(3) - Integer(9)*a*b*c+Integer(27)*a^Integer(2)*d):autosimplify() + + local C = sqrt((delta1 + sqrt(delta1 ^ Integer(2) - Integer(4) * delta0 ^ Integer(3))) / Integer(2), Integer(3)):autosimplify() + + if C == Integer.zero() then + C = sqrt((delta1 - sqrt(delta1 ^ Integer(2) - Integer(4) * delta0 ^ Integer(3))) / Integer(2), Integer(3)):autosimplify() + end + + if C == Integer.zero() then + C = (-b/(Integer(3)*a)):autosimplify() + end + + local eta = ((Integer(-1) + sqrt(Integer(-3))) / Integer(2)):autosimplify() + + return {((-Integer.one() / (Integer(3) * a)) * (b + C + delta0 / C)):autosimplify(), + ((-Integer.one() / (Integer(3) * a)) * (b + C*eta + delta0 / (C*eta))):autosimplify(), + ((-Integer.one() / (Integer(3) * a)) * (b + C*eta^Integer(2) + delta0 / (C*eta^Integer(2)))):autosimplify()} + end + end + if ispoly then + simplified = simplified:autosimplify() + end + if subpart then + simplified = (simplified - subpart):autosimplify() + end + return {RootExpression(simplified)} +end + +--- @return table<number, Expression> +function RootExpression:subexpressions() + return {self.expression} +end + +--- @param subexpressions table<number, Expression> +--- @return RootExpression +function RootExpression:setsubexpressions(subexpressions) + return RootExpression(subexpressions[1]) +end + +--- @param other Expression +--- @return boolean +function RootExpression:order(other) + --- TODO: Fix ordering on new expression types + if other:type() ~= RootExpression then + return false + end + + return self.expression:order(other.expression) +end + +--- @return string +function RootExpression:tolatex() + return '\\operatorname{RootOf}\\left(' .. self.expression:tolatex() .. '\\right)' +end + +----------------- +-- Inheritance -- +----------------- +__RootExpression.__index = CompoundExpression +__RootExpression.__call = RootExpression.new +RootExpression = setmetatable(RootExpression, __RootExpression)
\ No newline at end of file diff --git a/macros/luatex/latex/luacas/tex/algebra/sqrtexpression.lua b/macros/luatex/latex/luacas/tex/algebra/sqrtexpression.lua new file mode 100644 index 0000000000..bc96f4e1d8 --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/sqrtexpression.lua @@ -0,0 +1,196 @@ +--- @class SqrtExpression +--- An expression that represents the positive real solution to x^n = a where n is a positive integer and a is constant. +--- @field expression Expression +SqrtExpression = {} +__SqrtExpression = {} + +---------------------------- +-- Instance functionality -- +---------------------------- + +--- Creates a new sqrt expression with the given expression. +--- @param expression Expression +--- @param root Integer +--- @return SqrtExpression +function SqrtExpression:new(expression, root) + root = root or Integer(2) + local o = {} + local __o = Copy(__ExpressionOperations) + + o.expression = Copy(expression) + o.root = root + + __o.__index = SqrtExpression + __o.__tostring = function(a) + return tostring(a.expression) .. ' ^ (1/' .. tostring(a.root) .. ')' + end + __o.__eq = function(a, b) + -- This shouldn't be needed, since __eq should only fire if both metamethods have the same function, but for some reason Lua always rungs this anyway + if not b:type() == SqrtExpression then + return false + end + return a.expression == b.expression and a.root == b.root + end + o = setmetatable(o, __o) + + return o +end + + +--- @return table<number, Expression> +function SqrtExpression:subexpressions() + return {self.expression} +end + +--- @param subexpressions table<number, Expression> +--- @return SqrtExpression +function SqrtExpression:setsubexpressions(subexpressions) + return SqrtExpression(subexpressions[1], self.root) +end + +--- @param other Expression +--- @return boolean +function SqrtExpression:order(other) + return self:topower():order(other) +end + +function SqrtExpression:topower() + local exponent = BinaryOperation(BinaryOperation.DIV,{Integer.one(),self.root}):autosimplify() + local base = self.expression + return BinaryOperation(BinaryOperation.POW,{base,exponent}):autosimplify() +end + +function SqrtExpression:autosimplify() + local expression = self.expression:autosimplify() + local root = self.root:autosimplify() + + if root == Integer.one() then + return expression + end + + if root:type() == Rational then + return SqrtExpression(BinaryOperation(BinaryOperation.POW,{expression,root.denominator}):autosimplify(), root.numerator):autosimplify() + end + + if not root:isconstant() then + return BinaryOperation(BinaryOperation.POW,{expression,Integer.one() / root}):autosimplify() + end + + if not expression:isconstant() then + if expression.operation == BinaryOperation.MUL and expression.expressions[1]:isconstant() then + local coeff = SqrtExpression(expression.expressions[1],root):autosimplify() + expression.expressions[1] = BinaryOperation(BinaryOperation.MUL,{Integer.one()}) + expression = expression:autosimplify() + local sqrtpart = SqrtExpression(expression,root):autosimplify() + local result = coeff*sqrtpart + return result:autosimplify() + end + return BinaryOperation(BinaryOperation.POW,{expression,Integer.one() / root}):autosimplify() + end + + if expression:type() == Rational then + local result = BinaryOperation(BinaryOperation.MUL, {SqrtExpression(expression.numerator,root):autosimplify(),BinaryOperation(BinaryOperation.POW,{SqrtExpression(expression.denominator,root):autosimplify(),Integer(-1)})}) + return result:autosimplify() + end + + if expression:type() == Integer then + if expression == Integer.zero() then + return Integer.zero() + end + if expression == Integer.one() then + return Integer.one() + end + if expression < Integer.zero() then + if root == Integer(2) then + local result = SqrtExpression(expression:neg(),root):autosimplify() + result = I*result + return result:autosimplify() + end + if root % Integer(2) == Integer.one() then + local result = SqrtExpression(expression:neg(),root):autosimplify() + result = -result + return result:autosimplify() + end + end + local primes = expression:primefactorization() + local coeffresult = {} + local exprresult = {} + local reduction = root + for _, term in ipairs(primes.expressions) do + local primepower = term.expressions[2] + reduction = Integer.gcd(primepower,reduction) + if reduction == Integer.one() then + goto skip + end + end + ::skip:: + local newroot = root / reduction + for index, term in ipairs(primes.expressions) do + local prime = term.expressions[1] + local primepower = term.expressions[2] / reduction + local coeffpower = primepower // newroot + coeffresult[index] = prime ^ coeffpower + local exprpower = primepower - coeffpower*newroot + exprresult[index] = prime ^ exprpower + end + local newexpression = BinaryOperation(BinaryOperation.MUL,exprresult):autosimplify() + local coeff = BinaryOperation(BinaryOperation.MUL,coeffresult):autosimplify() + if coeff == Integer.one() then + if reduction == Integer.one() then + goto stop + end + return SqrtExpression(newexpression,newroot) + end + if newroot == Integer.one() then + return coeff + end + return BinaryOperation(BinaryOperation.MUL,{coeff,SqrtExpression(newexpression,newroot)}):autosimplify() + end + ::stop:: + + if expression.operation == BinaryOperation.POW and expression.expressions[2]:type() == Integer then + local exponent = expression.expressions[2] + local power = exponent // root + local newexponent = (exponent / root) - power + local coeff = expression.expressions[1] ^ power + coeff = coeff:evaluate() + if newexponent == Integer.zero() then + return coeff + else + local num = newexponent.numerator + local den = newexponent.denominator + local newexpression = expression ^ num + newexpression = newexpression:autosimplify() + local result = coeff * SqrtExpression(newexpression,den) + return result + end + end + + return SqrtExpression(expression,root) +end + +function SqrtExpression:tolatex() + local printout = '\\sqrt' + if self.root == Integer(2) then + printout = printout .. '{' .. self.expression:tolatex() .. '}' + else + printout = printout .. '[' .. self.root:tolatex() .. ']' .. '{' .. self.expression:tolatex() .. '}' + end + return printout +end + + +----------------- +-- Inheritance -- +----------------- +__SqrtExpression.__index = CompoundExpression +__SqrtExpression.__call = SqrtExpression.new +SqrtExpression = setmetatable(SqrtExpression, __SqrtExpression) + +---------------------- +-- Static constants -- +---------------------- + +sqrt = function(expression, root) + return SqrtExpression(expression, root) +end
\ No newline at end of file diff --git a/macros/luatex/latex/luacas/tex/algebra/trigexpression.lua b/macros/luatex/latex/luacas/tex/algebra/trigexpression.lua new file mode 100644 index 0000000000..307bfe1737 --- /dev/null +++ b/macros/luatex/latex/luacas/tex/algebra/trigexpression.lua @@ -0,0 +1,355 @@ +--- @class TrigExpression +--- Represents a trigonometric function from one expression to another. +--- @field name string +--- @field expression Expression +TrigExpression = {} +__TrigExpression = {} + +---------------------------- +-- Instance functionality -- +---------------------------- + +--- Creates a new trig expression with the given name and expression. +--- @param name string|SymbolExpression +--- @param expression Expression +--- @return TrigExpression +function TrigExpression:new(name, expression) + local o = {} + local __o = Copy(__ExpressionOperations) + + if not TrigExpression.NAMES[name] then + error("Argument error: " .. name .. " is not the name of a trigonometric function.") + end + + o.name = name + o.expression = expression + o.expressions = {expression} + if expression:isatomic() then + o.variables = {expression} + else + o.variables = {SymbolExpression('x')} + end + o.derivatives = {Integer.zero()} + + __o.__index = TrigExpression + __o.__tostring = function(a) + return tostring(a.name) .. '(' .. tostring(a.expression) .. ')' + end + __o.__eq = function(a, b) + -- if b:type() == FunctionExpression then + -- return a:tofunction() == b + -- end + -- This shouldn't be needed, since __eq should only fire if both metamethods have the same function, but for some reason Lua always runs this anyway + if not b:type() == TrigExpression then + return false + end + return a.name == b.name and a.expression == b.expression + end + + o = setmetatable(o, __o) + return o +end + +--- @return TrigExpression +function TrigExpression:evaluate() + local expression = self.expression:autosimplify() + + if expression == Integer.zero() then + if self.name == "cos" or self.name == "sec" then + return Integer.one() + end + if self.name == "sin" or self.name == "tan" then + return Integer.zero() + end + if self.name == "arctan" or self.name == "arcsin" then + return Integer.zero() + end + if self.name == "arccos" or self.name == "arccot" then + return PI / Integer(2) + end + end + + if expression == PI then + if self.name == "cos" or self.name == "sec" then + return Integer(-1) + end + if self.name == "sin" or self.name == "tan" then + return Integer.zero() + end + end + + if expression:ismulratlPI() then + local coeff = expression.expressions[1] + if TrigExpression.COSVALUES[tostring(coeff)] ~= nil then + if self.name == "cos" then + return TrigExpression.COSVALUES[tostring(coeff)]:autosimplify() + end + if self.name == "sin" then + local sign = Integer.one() + if coeff > Integer.one() then + sign = Integer(-1) + end + return (sign*sqrt(Integer.one()-cos(expression)^Integer(2))):autosimplify() + end + if self.name == "tan" then + return (sin(expression) / cos(expression)):autosimplify() + end + if self.name == "sec" then + return (Integer.one() / cos(expression)):autosimplify() + end + if self.name == "csc" then + return (Integer.one() / sin(expression)):autosimplify() + end + if self.name == "cot" then + return (cos(expression) / sin(expression)):autosimplify() + end + end + end + + if TrigExpression.ACOSVALUES[tostring(expression)] ~= nil then + if self.name == "arccos" then + return TrigExpression.ACOSVALUES[tostring(expression)]:autosimplify() + end + if self.name == "arcsin" then + if expression == Integer(-1) then + return TrigExpression.ACOSVALUES["-1"]:autosimplify() + elseif expression.expressions and expression.expressions[1] == Integer(-1) then + local expr = (Integer(-1)*sqrt(Integer.one() - expression ^ Integer(2))):autosimplify() + return TrigExpression.ACOSVALUES[tostring(expr)]:autosimplify() + else + local expr = (sqrt(Integer.one() - expression ^ Integer(2))):autosimplify() + return TrigExpression.ACOSVALUES[tostring(expr)]:autosimplify() + end + end + end + + if self.name == "arctan" and TrigExpression.ATANVALUES[tostring(expression)] ~= nil then + return TrigExpression.ATANVALUES[tostring(expression)]:autosimplify() + end + + return self +end + +--- checks if expression is a rational multiple of pi +--- @return boolean +function Expression:ismulratlPI() + if self.operation == BinaryOperation.MUL and #self.expressions == 2 and (self.expressions[1]:type() == Integer or self.expressions[1]:type() == Rational) and self.expressions[2] == PI then + return true + end + + return false +end + +--- @return TrigExpression +function TrigExpression:autosimplify() + local expression = self.expression:autosimplify() + + -- even and odd properties of trig functions + if (self.name == "sin" or self.name == "tan" or self.name == "csc" or self.name == "cot") and + expression.operation == BinaryOperation.MUL and expression.expressions[1]:isconstant() and expression.expressions[1] < Integer(0) then + return (-Integer.one() * TrigExpression(self.name, -expression)):autosimplify() + end + + if (self.name == "cos" or self.name == "sec") and + expression.operation == BinaryOperation.MUL and expression.expressions[1]:isconstant() and expression.expressions[1] < Integer(0) then + expression = (-expression):autosimplify() + end + + -- uses periodicity of sin and cos and friends + if self.name == "sin" or self.name == "cos" or self.name == "csc" or self.name == "sec" then + if expression == Integer.zero() or expression == PI then + goto skip + end + if expression.operation ~= BinaryOperation.ADD then + expression = BinaryOperation(BinaryOperation.ADD,{expression}) + end + for index,component in ipairs(expression.expressions) do + if component:ismulratlPI() then + local coeff = component.expressions[1] + if coeff:type() == Integer then + coeff = coeff % Integer(2) + coeff = coeff:autosimplify() + end + if coeff:type() == Rational then + local n = coeff.numerator + local d = coeff.denominator + local m = {n:divremainder(d)} + coeff = (m[1] % Integer(2)) + m[2]/d + coeff = coeff:autosimplify() + end + expression.expressions[index].expressions[1] = coeff + end + expression = expression:autosimplify() + end + ::skip:: + end + + -- uses periodicity of tan and cot + if self.name == "tan" or self.name == "cot" then + if expression == Integer.zero() or expression == PI then + goto skip + end + if expression.operation ~= BinaryOperation.ADD then + expression = BinaryOperation(BinaryOperation.ADD,{expression}) + end + for index,component in ipairs(expression.expressions) do + if component:ismulratlPI() then + local coeff = component.expressions[1] + if coeff:type() == Integer then + coeff = Integer.zero() + end + if coeff:type() == Rational then + local n = coeff.numerator + local d = coeff.denominator + local m = {n:divremainder(d)} + coeff = m[2]/d + coeff = coeff:autosimplify() + end + expression.expressions[index].expressions[1] = coeff + end + if component == PI then + expression.expressions[index] = Integer.zero() + end + end + expression = expression:autosimplify() + ::skip:: + end + + return TrigExpression(self.name, expression):evaluate() +end + +--- @return table<number, Expression> +function TrigExpression:subexpressions() + return {self.expression} +end + +--- @param subexpressions table<number, Expression> +--- @return TrigExpression +function TrigExpression:setsubexpressions(subexpressions) + return TrigExpression(self.name, subexpressions[1]) +end + +-- function TrigExpression:freeof(symbol) +-- return self.expression:freeof(symbol) +-- end + +-- function TrigExpression:substitute(map) +-- for expression, replacement in pairs(map) do +-- if self == expression then +-- return replacement +-- end +-- end +-- return TrigExpression(self.name, self.expression:substitute(map)) +-- end + +-- function TrigExpression:order(other) +-- return self:tofunction():order(other) +-- end + +-- function TrigExpression:tofunction() +-- return FunctionExpression(self.name, {self.expression}, true) +-- end + +----------------- +-- Inheritance -- +----------------- + +__TrigExpression.__index = FunctionExpression +__TrigExpression.__call = TrigExpression.new +TrigExpression = setmetatable(TrigExpression, __TrigExpression) + +---------------------- +-- Static constants -- +---------------------- +TrigExpression.NAMES = {sin=1, cos=2, tan=3, csc=4, sec=5, cot=6, + arcsin=7, arccos=8, arctan=9, arccsc=10, arcsec=11, arccot=12} + +TrigExpression.INVERSES = {sin="arcsin", cos="arccos", tan="arctan", csc="arccsc", sec="arcsec", cot="arccot", + arcsin="sin", arccos="cos", arctan="tan", arccsc="csc", arcsec="sec", arccot="cot"} + +TrigExpression.COSVALUES = { + ["0"] = Integer.one(), + ["1/6"] = sqrt(Integer(3))/Integer(2), + ["1/4"] = sqrt(Integer(2))/Integer(2), + ["1/3"] = Integer.one()/Integer(2), + ["1/2"] = Integer.zero(), + ["2/3"] = -Integer.one()/Integer(2), + ["3/4"] = -sqrt(Integer(2))/Integer(2), + ["5/6"] = -sqrt(Integer(3))/Integer(2), + ["1"] = -Integer.one(), + ["7/6"] = -sqrt(Integer(3))/Integer(2), + ["5/4"] = -sqrt(Integer(2))/Integer(2), + ["4/3"] = -Integer.one()/Integer(2), + ["3/2"] = Integer.zero(), + ["5/3"] = Integer.one()/Integer(2), + ["7/4"] = sqrt(Integer(2))/Integer(2), + ["11/6"] = sqrt(Integer(3))/Integer(2), +} +TrigExpression.ACOSVALUES = { + ["1"] = Integer.zero(), + ["(1/2 * sqrt(3,2))"] = PI * Integer(6) ^ Integer(-1), + ["(1/2 * sqrt(2,2))"] = PI * Integer(4) ^ Integer(-1), + ["1/2"] = PI * Integer(3) ^ Integer(-1), + ["0"] = PI * Integer(2) ^ Integer(-1), + ["-1/2"] = PI * Integer(2) * Integer(3) ^ Integer(-1), + ["(-1/2 * sqrt(2,2))"]= PI * Integer(3) * Integer(4) ^ Integer(-1), + ["(-1/2 * sqrt(3,2))"]= PI * Integer(5) * Integer(6) ^ Integer(-1), + ["-1"] = Integer(-1)*PI, +} +TrigExpression.ATANVALUES = { + ["(-1 * sqrt(3,2))"] = Integer(-1) * PI * Integer(3) ^ Integer(-1), + ["-1"] = Integer(-1) * PI * Integer(4) ^ Integer(-1), + ["(-1/3 * sqrt(3,2))"] = Integer(-1) * Integer(6) ^ Integer(-1), + ["0"] = Integer.zero(), + ["(1/3 * sqrt(3,2))"] = PI * Integer(6) ^ Integer(-1), + ["1"] = PI * Integer(4) ^ Integer(-1), + ["sqrt(3,2)"] = PI * Integer(3) ^ Integer(-1) +} + +SIN = function (a) + return TrigExpression("sin", a) +end + +COS = function (a) + return TrigExpression("cos", a) +end + +TAN = function (a) + return TrigExpression("tan", a) +end + +CSC = function (a) + return TrigExpression("csc", a) +end + +SEC = function (a) + return TrigExpression("sec", a) +end + +COT = function (a) + return TrigExpression("cot", a) +end + +ARCSIN = function (a) + return TrigExpression("arcsin", a) +end + +ARCCOS = function (a) + return TrigExpression("arccos", a) +end + +ARCTAN = function (a) + return TrigExpression("arctan", a) +end + +ARCCSC = function (a) + return TrigExpression("arccsc", a) +end + +ARCSEC = function (a) + return TrigExpression("arcsec", a) +end + +ARCCOT = function (a) + return TrigExpression("arccot", a) +end
\ No newline at end of file |