module Algebra::PolynomialFactorization::Z
Public Instance Methods
_factorize()
click to toggle source
# File lib/algebra/polynomial-factor-int.rb, line 135 def _factorize if constant? return Algebra::Factors.new([[self, 1]]), 0 end f = sqfree_over_integral flc = f.lc f = f.monic_int(flc) f0 = reduce_over_prime_field(f) fa = f0.factorize_modp f1 = f factors = Algebra::Factors.new # fa.each_product0 do |g0| fa.each_product do |g0| next if g0.deg > f1.deg # if fact_q = f1.hensel_lift(g0, f0.class) if fact_q = f1.hensel_lift(g0) fact, f1 = fact_q factors.push [fact, 1] break if f1.deg <= 0 true end end factors.map!{|f2| f2.project(self.class){|c, j| c*flc**j}.pp } facts = factors.fact_all(self) [facts.correct_lc!(self), f0.ground.char] end
centorize(f, n)
click to toggle source
# File lib/algebra/polynomial-factor-int.rb, line 25 def centorize(f, n) half = n / 2 f.class.new(*f.collect{|c| c > half ? c - n : c }) end
divmod_ED(o)
click to toggle source
# File lib/algebra/polynomial-factor-int.rb, line 51 def divmod_ED(o) # divmod in the ring over Euclidian Domain raise ZeroDivisionError, "divisor is zero" if o.zero? d0, d1 = deg, o.deg if d1 <= 0 a = collect{|c| q, r = (ground.field? ? [ground.zero, c/o.lc] : c.divmond(o.lc)) if r.zero? [q, zero] else return [zero, self] end } [self.class.new(*a), zero] elsif d0 < d1 [zero, self] else q, r = (ground.field? ? [lc/o.lc, ground.zero] : lc.divmod(o.lc)) if r.zero? tk = monomial(d0 - d1) * q e = rt - o.rt * tk q, r = e.divmod_ED(o) [q + tk, r] else [zero, self] end end end
factorize_int()
click to toggle source
# File lib/algebra/polynomial-factor-int.rb, line 165 def factorize_int _factorize.first end
hensel_lift(g0)
click to toggle source
# File lib/algebra/polynomial-factor-int.rb, line 101 def hensel_lift(g0) ring, pf, char = self.class, g0.class, g0.class.ground.char f0 = pf.reduce(self) h0, r0 = f0.divmod g0 unless r0.zero? raise "each_product does not work well" #return nil end g = lifting(g0, ring) h = lifting(h0, ring) a, b = 0, 0 0.upto try_e(g, char) do |e| mod = char ** e mod1 = mod * char g = centorize(g + b * mod, mod1) h = centorize(h + a * mod, mod1) q, r = divmod(g) if r.zero? return [g, q] end c = (self - g * h) / mod1 g0, h0, c0 = pf.reduce(g), pf.reduce(h), pf.reduce(c) a0, b0 = resolve_c(g0, h0, c0) a, b = lifting(a0, ring), lifting(b0, ring) end return nil end
lifting(f, ring)
click to toggle source
# File lib/algebra/polynomial-factor-int.rb, line 97 def lifting(f, ring) f.project(ring){|c, i| c.lift} end
reduce_over_prime_field(f)
click to toggle source
# File lib/algebra/polynomial-factor-int.rb, line 84 def reduce_over_prime_field(f) Primes.new.each do |prime| if prime > 10 pf = Algebra::Polynomial.reduction(Integer, prime, "x") f0 = pf.reduce(f) r = f0.resultant_by_elimination(f0.derivate) unless r.zero? return f0 end end end end
resolve_c(g, h, c)
click to toggle source
# File lib/algebra/polynomial-factor-int.rb, line 13 def resolve_c(g, h, c) # find [a, b] s.t. a * g + b * h == c d, a, b = g.gcd_coeff(h) q, r = c.divmod(d) raise "#{c} is not devided by (#{g}, #{h})" unless r.zero? a = a * q b = b * q q, r = a.divmod(h) a = r #= a - q * h b = b + q * g [a, b] end
resultant_by_elimination(o)
click to toggle source
# File lib/algebra/polynomial-factor-int.rb, line 80 def resultant_by_elimination(o) #ground ring must be a field sylvester_matrix(o).determinant_by_elimination end
sqfree_over_integral()
click to toggle source
# File lib/algebra/polynomial-factor-int.rb, line 44 def sqfree_over_integral f = pp g = f.derivate.pp h = f.pgcd(g) f / h.pp end
try_e(g, u)
click to toggle source
# File lib/algebra/polynomial-factor-int.rb, line 30 def try_e(g, u) m = g.deg c = 1.0 n = (m/2).to_i 0.upto n - 1 do |i| c *= (m - i) / (n - i) end s = 0 each do |c| s += c**2 end (Math.log(c * Math.sqrt(s))/Math.log(u)).to_i + 1 end