summaryrefslogtreecommitdiff
path: root/macros/luatex/latex/luacas/tex/algebra
diff options
context:
space:
mode:
Diffstat (limited to 'macros/luatex/latex/luacas/tex/algebra')
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/_init.lua24
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/absexpression.lua80
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/equation.lua183
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/euclideandomain.lua63
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/factorialexpression.lua102
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/field.lua69
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/integer.lua952
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/integerquotientring.lua197
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/logarithm.lua186
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/polynomialring.lua860
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/polynomialring/berlekampfactoring.lua187
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/polynomialring/decomposition.lua93
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/polynomialring/zassenhausfactoring.lua220
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/rational.lua241
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/ring.lua326
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/rootexpression.lua135
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/sqrtexpression.lua196
-rw-r--r--macros/luatex/latex/luacas/tex/algebra/trigexpression.lua355
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