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