# Multivariate polynomial ring class over arbitary ring
#
# by Shin-ichiro Hara
#
# Version 1.6 (2002.02.05)
require "algebra/import-module-single-thread"
require "algebra/m-index"
require "algebra/algebraic-system"
module Algebra
autoload :Groebner, "algebra/groebner-basis"
def MPolynomial(ground, *vars)
Algebra::MPolynomial.create(ground, *vars)
end
module_function :MPolynomial
class MPolynomial
auto_req_init
auto_req(:value_on, "algebra/polynomial-converter"){
include MPolynomialConverter
}
auto_req_s_init
auto_req_s(:convert_to, "algebra/polynomial-converter"){
extend MPolynomialConvertTo
}
=begin
class MonomialOrder
def initialize(ord = :lex, t = nil)
# ord ||= :lex
case ord
when Symbol
@ord_str = ord.id2name
@ord_sym = ord
when String
@ord_sym = ord.intern
@ord_str = ord
else
raise "want Symbol or String"
end
if t
@t = t
instance_eval("alias cmp #{@ord_str}_t")
else
instance_eval("alias cmp #@ord_str")
end
end
def to_s; @ord_str; end
def to_sym; @ord_sym; end
def v_ord; @t; end
#standerd ordering
def lex(me, other)
me.to_a <=> other.to_a
end
def grlex(me, other)
s = (me.totdeg <=> other.totdeg)
return s unless s.zero?
lex(me, other)
end
def grevlex(me, other)
s = (me.totdeg <=> other.totdeg)
return s unless s.zero?
n = [me.size, other.size].max
(n-1).downto 0 do |i|
x = other[i] - me[i]
return x unless x.zero?
end
0
end
#ordering by transition array
def lex_t(me, other)
n = [me.size, other.size].max
0.upto (n-1) do |i|
x = me[@t[i]] - other[@t[i]]
return x unless x.zero?
end
0
end
def grlex_t(me, other)
s = (me.totdeg <=> other.totdeg)
return s unless s.zero?
lex_t(me, other)
end
def grevlex_t(me, other)
s = (me.totdeg <=> other.totdeg)
return s unless s.zero?
n = [me.size, other.size].max
(n-1).downto 0 do |i|
x = other[@t[i]] - me[@t[i]]
return x unless x.zero?
end
0
end
end
=end
include Enumerable
include Comparable
extend AlgebraCreator
include AlgebraBase
auto_req :factorize, "algebra/m-polynomial-factor"
auto_req :divmod_s, "algebra/groebner-basis"
auto_req :S_pair, "algebra/groebner-basis"
auto_req :S_pair_coeff, "algebra/groebner-basis"
def initialize(ha = {})
@bone = ha
method_cash_clear
end
def self.create(ground, *vars)
klass = super(ground)
klass.sysvar(:order, nil)
klass.sysvar(:v_order, nil)
klass.set_ord(:lex)
klass.sysvar(:variables, vars)
klass.sysvar(:display_type, :norm)
klass
end
def variables; self.class.variables; end
def vars; self.class.vars; end
def self.vars(*vs)
if vs.size > 0
vs = vs[0].scan(/[a-zA-Z]_?\d*/) if vs.size == 1 && vs[0].is_a?(String)
vs.collect{|var| indeterminate(var)}
else
variables.collect{|var| indeterminate(var) }
end
end
def self.to_ary
[self, *vars]
end
def self.indeterminate(obj)
ind = if i = variables.index(obj)
[0] * i + [1]
# Algebra::MIndex.monomial(i)
else
variables.push obj
# [0] * (variables.size-1) + [1]
Algebra::MIndex.monomial(variables.size-1)
end
monomial(ind)
end
def self.var(obj)
indeterminate(obj)
end
def self.mvar(*array)
if array.empty?
variables.collect{|var| indeterminate(var) }
else
array.collect{|var| indeterminate(var)}
end
end
def self.monomial(ind = [1], c = ground.unity)
new({Algebra::MIndex[*ind] => c})
end
def monomial(ind = [1], c = ground.unity)
self.class.monomial(ind, c)
end
def self.const(x)
new({Algebra::MIndex::Unity => x})
end
def self.regulate(x)
if x.is_a? self
x
elsif y = ground.regulate(x)
const(y)
else
nil
end
end
def self.[](*x)
new(Hash[*x])
end
def self.euclidian?
ground.euclidian? && variables.size <= 1
end
def self.ufd?
euclidian? || ground.field?
end
def each(&b); @bone.each(&b); end
def keys; @bone.keys; end
def values; @bone.values; end
def [](ind)
@bone[ind] || 0
end
def []=(ind, c)
@bone[ind] = c
end
def monomial?
flag = false
each do |m, c|
unless c.zero?
if flag
return false
else
flag = true
end
end
end
flag
end
def constant?
each do |m, c|
if m.totdeg > 0 && !c.zero?
return false
end
end
true
end
def unit?
constant? && lc.unit?
end
def constant_coeff
# we assume {[]=>2, [0]=>3} is not appeared.
# it must be reduced to {[]=>6}.
@bone[Algebra::MIndex::Unity]
# k = ground.unity
# @bone.each do |ind, c|
# unless c.zero?
# k *= c if ind.totdeg == 0
# end
# end
# k
end
def include?(k)
@bone.include?(k)
end
def self.zero; new; end
def self.unity; const(ground.unity); end
def zero?; !find{|m, c| !c.zero?}; end
def self.set_ord(ord = nil, v_ord = nil)
# @@order = MonomialOrder.new(ord, v_ord)
if ord.is_a? String
ord = ord.intern
elsif ord.is_a?(Symbol) || ord.nil?
else
raise "#{ord} must be String or Symbol"
end
if v_ord.is_a? String
v_ord = v_ord.intern
elsif v_ord.is_a?(Symbol) || v_ord.nil?
else
raise "#{v_ord} must be String or Symbol"
end
if ord.nil?
ord = get_ord
end
# self.order = MonomialOrder.new(ord, v_ord) #o
self.order = ord if ord #n
self.v_order = v_ord if v_ord #n
MIndex.set_ord(ord) if ord #n
MIndex.set_v_ord(v_ord) if v_ord #n
end
def self.ord=(ord)
self.set_ord(ord, nil)
end
def self.get_ord
order
end
def self.ord
get_ord
end
def self.get_v_ord
v_order
end
def self.with_ord_old(ord = nil, v_ord = nil, methods = [])
method_cash_clear(*methods)
o_ord, o_v_ord = get_ord, get_v_ord
result = nil
begin
set_ord(ord, v_ord)
result = yield
ensure
set_ord(o_ord, o_v_ord)
end
result
end
def self.with_ord(ord = nil, v_ord = nil, methods = [])
method_cash_clear(*methods)
o_ord, o_v_ord = get_ord, get_v_ord
ord = o_ord unless ord
result = nil
if v_ord
MIndex.set_v_ord(ord, v_ord)
MIndex.import_module(MIndex.get_module(ord || self.order, v_ord)) do
result = yield
end
MIndex.set_v_ord(o_ord, o_v_ord)
else
MIndex.import_module(MIndex.get_module(ord || self.order)) do
result = yield
end
end
result
end
def self.method_cash_clear(*m)
# $stderr.puts "method_cash_clear is given no parameter" if m.empty?
m.each do |x|
x.method_cash_clear
end
end
def method_cash_clear
@lc = @lm = @lt = @rt = @multideg = nil
end
def compact!
a = {}
each do |m, c|
a[m.compact] = c unless c.zero?
end
@bone = a
self
end
def multideg
if @multideg
@multideg
else
d = Algebra::MIndex::Unity
each do |m, c|
# d = m if !c.zero? && m > d
# d = m if !c.zero? && self.class.order.cmp(m, d) > 0 #o
d = m if !c.zero? && (m <=> d) > 0 #n
end
@multideg = d
end
end
def totdeg
if zero?
-1
else
max = 0
each do |ind, c|
unless c.zero?
t = ind.totdeg
max = t if t > max
end
end
max
end
end
def deg; multideg; end
def ==(other)
super{ |other|
each do |m, c|
return false unless c == other[m]
end
other.each do |m, c|
if !include? m
return false unless c.zero?
end
end
true
}
end
def <=>(other)
if o = regulate(other)
# THIS IS RATHER VAGUE
# lm.ind <=> o.lm.ind
# self.class.order.cmp(lm.ind, o.lm.ind) #o
lm.ind <=> o.lm.ind #n
else
x , y = other.coerce(self)
x <=> y
end
end
def +(other)
super{ |other|
a = zero
each do |m, c|
a[m] = c + other[m]
end
other.each do |m, c|
if !include? m
a[m] = self[m] + c
# a[m] = c + self[m]
end
end
a
}
end
def -(other)
super{ |other|
a = zero
each do |m, c|
a[m] = c - other[m]
end
other.each do |m, c|
if !include? m
a[m] = self[m] - c
end
end
a
}
end
def *(other)
super{ |other|
a = zero
each do |m0, c0|
other.each do |m1, c1|
a[m0 + m1] += c0 * c1
end
end
# a.compact!
a
}
end
def /(other) # other asume to be scalar
# if (o = regulate(other))
super{ |o|
if o.constant?
# self * (ground.unity / o.constant_coeff)
project{|c, ind|
c/o.constant_coeff
}
else
a, r = divmod(other)
if r.zero?
a.first
else
raise "#{other} not divides #{self}"
end
end
}
# else
# raise "other(#{o}, #{o.class}) must be monomial"
# end
end
def divmod(*f)
g = clone # not dup
f_lts = f.collect{|x| x.lt}
a = (0...f.size).collect{zero}
r = zero
until g.zero?
flag = false
f_lts.each_with_index do |f_lt, i|
g_lt = g.lt
if f_lt.divide? g_lt
d = g_lt / f_lt
if d == 0 then exit; end
a[i] += d
g = g.rt - (d * f[i]).rt
g.compact!
flag = true
break
end
end
unless flag
r += g.lt
g = g.rt #g.rt!
g.compact!
end
end
[a, r]
end
def divmod_variant(*f)
g = clone # not dup
f_lts = f.collect{|x| x.lt}
a = (0...f.size).collect{zero}
r = zero
until g.zero?
flag = false
f_lts.each_with_index do |f_lt, i|
g_lt = g.lt
if f_lt.divide? g_lt
d = g_lt / f_lt
a[i] += d
g = g.rt - (d * f[i]).rt
g.compact!
flag = true
g.zero? ? break : redo
end
end
unless flag
r += g.lt
g = g.rt #g.rt!
end
end
[a, r]
end
def %(others)
divmod(*others).last
end
def gcd(other)
fact = factorize
fact.gcd(other).pi
end
def lc
@lc ||= self[multideg]
end
def lm
#@lm ||= monomial(multideg).extend(Monomial)
monomial(multideg).extend(Monomial)
end
def lt
md = multideg
@lt ||= monomial(md, self[md]).extend(Monomial)
end
def rt
@rt ||= self - lt
end
def rt! # deepy destructive, do NOT use
@bone.delete(multideg)
method_cash_clear
self
end
def index_eval(ring, ary, ind)
e = ring.unity
ind.each_with_index do |n, i|
e *= ary[i]**n if n > 0
end
e
end
def project0(ring = self.class, ary = ring.vars)
e = ring.zero
each do |ind, c|
e += yield(c, ind)
end
e
end
alias map_to project0
def project(ring = self.class, ary = ring.vars)
e = ring.zero
each do |ind, c|
e += index_eval(ring, ary, ind) * yield(c, ind)
end
e
end
def evaluate(*ary)
project(ground, ary){|c, ind| c}
# project(self.class, ary){|c, ind| c}
end
def sub(var, value)
vs = self.class.vars
if i = vs.index(var)
vs = vs.dup
vs[i] = value
evaluate(*vs)
else
raise "#{var} is not a variable"
end
end
alias call evaluate
def derivate(v)
if an = self.class.variables.index(v) || vars.index(v)
derivate_at(an)
else
raise "#{var} is not a variable"
end
end
def derivate_at(an)
ary = vars
project0 {|c, ind|
if ind[an] < 1
zero
else
ind0 = ind - Algebra::MIndex.monomial(an)
index_eval(self.class, ary, ind0) * ind[an] * c
end
}
end
def move_const(pb, coef = ground.unity)
vs = self.class.vars.dup
a = []
vs.each_with_index do |x, i|
a << x + pb[i]*coef
end
evaluate(*a)
end
def convert_to(ring)
project(ring){|c, ind| c}
end
def need_paren_in_coeff?
if constant?
c = constant_coeff
if c.respond_to?(:need_paren_in_coeff?)
c.need_paren_in_coeff?
elsif c.is_a?(Numeric)
false
else
true
end
elsif !monomial?
true
else
false
end
end
def to_s
case self.class.display_type
when :code
pr, po = "*", "**"
when :tex
pr, po = "", "}^"
else
pr, po = "", "^"
end
a = []
# (k = keys).sort!
# k = keys.sort{|x, y| self.class.order.cmp(x, y)} #o
k = keys.sort{|x, y| x <=> y} #n
k.each do |m| c = self[m]
next if c.zero?
s = if m.unity?
c.to_s
else
case c
when Numeric
c = c == 1 ? "" : c == -1 ? "-" : c.to_s
when self.class
c = "(#{c})"
else
c = if c == 1
""
elsif c == -1
"-"
elsif c.respond_to?(:need_paren_in_coeff?)
if c.need_paren_in_coeff?
"(#{c})"
else
c.to_s
end
elsif c.is_a?(Numeric)
c.to_s
else
"(#{c})"
end
end
ms = m.to_s(variables, pr, po)
c.empty? ? ms : c + pr + ms
end
a.unshift s
end
a.unshift "0" if a.empty?
a.join(" + ").gsub(/\+ -/, "- ")
end
alias inspect to_s unless $DEBUG
end
class MPolynomial
module Monomial# < MPolynomial
attr_reader :ind, :coeff
def self.extend_object(obj)
# p [obj, obj.class, obj.monomial?, obj == 0]
# exit unless obj.monomial?
raise "'#{obj}' must be monomial" unless obj.monomial?
super
obj.ind
obj.coeff
self
end
def ind
return @ind if @ind
each do |@ind, @coeff|
return @ind
end
raise "internal error"
end
def coeff
return @coeff if @coeff
each do |@ind, @coeff|
return @coeff
end
raise "internal error"
end
def <=>(other)
#@ind <=> other.ind
# self.class.order.cmp(@ind, other.ind) #o
@ind <=> other.ind #n
end
def divide?(other)
return true if other.zero?
case other
when Monomial
ind.devide?(other.ind) && (ground.field? ||
ground.euclidian? &&
coeff.devide?(other.coeff))
else
raise "unkown self.class #{other}(#{other.class})"
end
end
def divide_or?(other0, other1)
return true if other0.zero? && other1.zero?
if other0.is_a? Monomial and other1.is_a? Monomial
ind.devide_or?(other0.ind, other1.ind)
else
raise "unkown self.class #{other0.class}, #{other1.class}"
end
end
def prime_to?(other)
return false if other.zero?
case other
when Monomial
ind.prime_to?(other.ind)
else
raise "unkown self.class #{other}(#{other.class})"
end
end
def /(other)
case other
when self.class
self.class[(ind - other.ind), coeff / other.coeff]
else
super
end
end
def lcm(other)
monomial(ind.lcm(other.ind)).extend(Monomial)
end
end
end
end
if $0 == __FILE__
include Algebra
Foo = MPolynomial(Integer)
f = Foo[Algebra::MIndex[1,1]=>1, Algebra::MIndex[1,3]=>4, Algebra::MIndex[]=>9]
g = Foo[Algebra::MIndex[1,1]=>1, Algebra::MIndex[1,3]=>4, Algebra::MIndex[]=>9,
Algebra::MIndex[0,2,3]=>6, Algebra::MIndex[2,3]=>-1, Algebra::MIndex[2]=>-5,Algebra::MIndex[0,1]=>7]
$DEBUG = true
puts f, g, f*g
$DEBUG = false
puts
require "algebra/polynomial-converter"
require "algebra/mathn"
P = MPolynomial(Rational)
x, y, z = P.vars("xyz")
f = x**2*y + x*y**2 + y*2 + z**3
g = x*y-z**3
h = y*2-6*z
P.set_ord(:lex)
puts "LEX:"
puts "(#{f}).divmod([#{g}, #{h}]) =>"
fgh = f.divmod(g, h)
puts fgh
puts
P.method_cash_clear(f, g, h)
P.set_ord(:grlex)
puts "GRLEX:"
puts "(#{f}).divmod([#{g}, #{h}]) =>"
fgh = f.divmod(g, h)
puts fgh
puts
P.method_cash_clear(f, g, h)
P.set_ord(:grevlex)
puts "GREVLEX:"
puts "(#{f}).divmod([#{g}, #{h}]) =>"
fgh = f.divmod(g, h)
puts fgh
puts
P.variables.clear
z, y, x = P.vars("zyx")
f = x**2*y + x*y**2 + y*2 + z**3
g = x*y-z**3
h = y*2-6*z
puts "GREVLEX: order z, y, x"
puts "(#{f}).divmod([#{g}, #{h}]) =>"
fgh = f.divmod(g, h)
puts fgh
puts
end
syntax highlighted by Code2HTML, v. 0.9.1