module Algebra::MPolynomialFactorization

Public Instance Methods

_factorize(na = 0) click to toggle source
# File lib/algebra/m-polynomial-factor.rb, line 41
def _factorize(na = 0)
  pring = Algebra.Polynomial(ground, 'x')

  f = self
  f_one = nil
  pb0 = nil
  # f     : m-poly. moved by pb0 reduced onto f_one
  # f_one : sqfree poly.
  # pb0   : parallel moving vector

  if ground == Integer
    each_point_int(vars.size, na) do |pb0_dash|
      pb0 = pb0_dash
      f = move_const(pb0) if pb0.find { |x| !x.zero? }
      f_one = f.reduce2onevar(pring, na)
      break if f_one.psqfree?
    end
    fact_one, char = f_one._factorize
    hensel = :hensel_lift_int
  else # ground == Zp
    char = ground.char
    each_point_zp(vars.size, na, char) do |pb0_dash|
      pb0 = pb0_dash
      f = move_const(pb0) if pb0.find { |x| !x.zero? }
      f_one = f.reduce2onevar(pring, na)
      break if f_one.sqfree?
    end
    fact_one = f_one.factorize_modp
    hensel = :hensel_lift_zp
  end

  factors = Factors.new
  if fact_one.fact_size <= 1
    factors *= self
    return factors
  end
  height = totdeg_away_from(0)

  f1 = f
  fact_one.each_product do |g0|
    f1_one = f1.reduce2onevar(pring, na)
    next if g0.deg > f1_one.deg
    #      fact_q = f1.send(hensel, g0, f1_one, char, height)
    fact_q = f1.send(hensel, g0, f1_one, char, height, na)
    next unless fact_q
    #        if fact_q = f1.hensel_lift(g0, f1_one, char, height)
    f0, f1 = fact_q
    factors.push [f0, 1]
    break if f1.totdeg <= 0
    true
  end

  factors.collect { |f3, i| [f3.move_const(pb0, -1), i] }
end
_hensel_lift(g0, f0, _char, height, where) { |ha, ary, cofacts| ... } click to toggle source

def _hensel_lift(g0, f0, char, height)

# File lib/algebra/m-polynomial-factor.rb, line 97
def _hensel_lift(g0, f0, _char, height, where)
  # self in MPolynomial/ZorZp
  # g0 in Polyomial/ZorZp, candidate of factor of f0
  # f0 in Polyomial/ZoZzp, one variable reduction of self

  ring = self.class
  ring_one = g0.class

  h0, r0 = f0.divmod g0
  raise 'each_product does not work well' unless r0.zero?

  #    where = 0
  ary = [g0, h0]
  cofacts = mk_cofacts(ary)
  fk = ary.collect { |fx| fx.value_on_one(ring, where) } # MPolynomial

  height.times do |k|
    c = self - fk[0] * fk[1]
    h = c.annihilate_away_from(where, k + 1)
    h_alpha0 = h.collect_away_from(where, ring_one) # Hash: ind0=>Polynomial
    h_alpha = {}

    h_alpha0.each do |ind, ha|
      h_alpha[ind] = yield(ha, ary, cofacts)
    end

    hias = ary.collect { {} }
    h_alpha.each do |ind, ha_i|
      ha_i.each_with_index do |fx, i|
        next if fx.zero?
        hias[i][ind] = fx
      end
    end

    hi = hias.collect do |hia|
      e = ring.zero
      hia.each do |ind, fx|
        e += ring.monomial(ind) * fx.value_on_one(ring, where)
      end
      e
    end
    fk_new = []
    hi.each_with_index do |fx, i|
      fk_new.push fk[i] + fx
    end
    fk = fk_new
  end
  fk
end
_sqfree() click to toggle source
# File lib/algebra/m-polynomial-factor.rb, line 147
def _sqfree
  dvs = []
  vars.each_with_index do |_v, i|
    dv = derivate_at(i)
    next if dv.zero?
    dvs.push dv
  end
  g = gcd_all(*dvs)
  self / g
end
centorize(f, n) click to toggle source
# File lib/algebra/m-polynomial-factor-int.rb, line 92
def centorize(f, n)
  half = n / 2
  #    f.class.new(*f.collect{|c| c > half ? c - n : c })
  f.project(f.class) do |c, _ind|
    c = c % n
    c > half ? c - n : c
  end
end
each_point_int(n0, na) { |pb0| ... } click to toggle source
# File lib/algebra/m-polynomial-factor.rb, line 17
def each_point_int(n0, na)
  # obscure algorithm
  pblock = 7
  m = 0
  loop do
    each_point_zp(n0, na, m + pblock) do |pb0|
      next if (1...n0 - na).find { |k| pb0[na + k] < m }
      yield pb0
    end
    m += pblock
    $stderr.puts "warning (each_point_int): it's hard to find adequate zero point"
  end
end
each_point_zp(n0, na, char) { |pb0| ... } click to toggle source
# File lib/algebra/m-polynomial-factor.rb, line 31
def each_point_zp(n0, na, char)
  avoid_indice = indice_of_constant
  Combinatorial.power_cubic(char, n0 - na - 1) do |a|
    pb0 = (0..na).collect { 0 } + a
    next if avoid_indice.find { |i| !pb0[i].zero? }
    yield(pb0)
  end
  raise 'each_point(limit over)'
end
factorize_int(an = 0) click to toggle source

include PolynomialFactorization::Q

# File lib/algebra/m-polynomial-factor-int.rb, line 21
def factorize_int(an = 0)
  return Factors.new([[self, 1]]) if constant?
  f = self
  f = f.sqfree if an == 0
  an += 1 while f.deg_at(an) <= 0
  f = f.pp(an)
  lc_of_f = f.lc_at(an)
  f = f.monic_int(an)
  facts = f._factorize(an)
  facts = facts.map! do |f2|
    f3 = f2.project(self.class) { |c, j| c * lc_of_f**j[an] }
    f3.pp(an)
  end
  facts = facts.fact_all_u(self)
  facts.mcorrect_lc!(self, an) do |fac|
    fac.factorize_int(an + 1)
  end
  facts
end
factorize_modp(an = 0) click to toggle source

include PolynomialFactorization::Z

# File lib/algebra/m-polynomial-factor-zp.rb, line 19
def factorize_modp(an = 0)
  f = self
  f = f.sqfree_zp if an == 0
  an += 1 while f.deg_at(an) <= 0
  f = f.pp(an)
  lc_of_f = f.lc_at(an)
  f = f.monic_int(an)
  facts = f._factorize(an)
  facts = facts.map! do |f2|
    f3 = f2.project(self.class) { |c, j| c * lc_of_f**j[an] }
    f3.pp(an)
  end
  facts = facts.fact_all_u(self)
  facts.mcorrect_lc!(self, an) do |fac|
    fac.factorize_modp(an + 1)
  end
  facts
end
gelfond_bound(char) click to toggle source
# File lib/algebra/m-polynomial-factor-int.rb, line 45
def gelfond_bound(char)
  n = vars.size
  ds = (0...n).collect { 0 }
  c0 = 0
  each do |ind, c|
    next if c.zero?
    ind.each_with_index do |d0, i|
      ds[i] = d0 if d0 > ds[i]
    end
    c0 = c.abs if c.abs > c0
  end
  tot = 0
  ds.each do |x|
    tot += x
  end
  b = Math::E**tot * c0
  l = 0
  s = 1
  while s <= 2 * b
    s *= char
    l += 1
  end
  l
end
hensel_lift_int(g0, f0, char, height, where) click to toggle source

def hensel_lift_int(g0, f0, char, height)

# File lib/algebra/m-polynomial-factor-int.rb, line 71
def hensel_lift_int(g0, f0, char, height, where)
  # self in MPolynomial/int
  # g0 in Polyomial/Z, candidate of factor of f0
  # f0 in Polyomial/Z, one variable reduction of self

  pheight = gelfond_bound(char)

  fk = _hensel_lift(g0, f0, char, height, where) do |ha, ary, cofacts|
    #    fk =  _hensel_lift(g0, f0, char, height) {|ha, ary, cofacts|
    decompose_on_cofactors_p_adic(ha, ary, cofacts, char, pheight)
  end
  mod = char**pheight
  g = centorize(fk[0], mod)
  h = centorize(fk[1], mod)

  r = self - g * h
  return [g, h] if r.zero?

  nil
end
hensel_lift_zp(g0, f0, char, height, where) click to toggle source

def hensel_lift_zp(g0, f0, char, height)

# File lib/algebra/m-polynomial-factor-zp.rb, line 59
def hensel_lift_zp(g0, f0, char, height, where)
  # self in MPolynomial/int
  # g0 in Polyomial/Zp, candidate of factor of f0
  # f0 in Polyomial/Zp, one variable reduction of self

  #    fk = _hensel_lift(g0, f0, char, height) {|ha, ary, cofacts|
  fk = _hensel_lift(g0, f0, char, height, where) do |ha, ary, cofacts|
    decompose_on_cofactors_modp(ha, ary, cofacts)
  end

  r = self - fk[0] * fk[1]
  return fk if r.zero?

  nil
end
monic_int(na = 0) click to toggle source
# File lib/algebra/m-polynomial-factor.rb, line 164
def monic_int(na = 0)
  d = deg_at(na)
  a = lc_at(na)
  ary = self.class.vars
  m = project0(self.class) do |c, ind|
    n = ind[na]
    if n == d
      zero
    elsif n < d
      c * a**(d - n - 1) * index_eval(self.class, ary, ind)
    elsif c.zero?
      zero
    else
      raise 'fatal contradiction!'
    end
  end

  index_eval(self.class, ary, MIndex.monomial(na, d)) + m
end
pp(an) click to toggle source
# File lib/algebra/m-polynomial-factor.rb, line 158
def pp(an)
  x, *a = coeffs_at(an)
  content = x.gcd_all(*a)
  self / content
end
sqfree() click to toggle source
# File lib/algebra/m-polynomial-factor-int.rb, line 41
def sqfree
  _sqfree
end
sqfree_zp() click to toggle source
# File lib/algebra/m-polynomial-factor-zp.rb, line 47
def sqfree_zp
  f = self
  s = unity
  until f.constant?
    g = f._sqfree
    s *= g
    f = (f / g).verschiebung
  end
  s / s.lc
end
verschiebung() click to toggle source
# File lib/algebra/m-polynomial-factor-zp.rb, line 38
def verschiebung
  ary = vars
  ch = ground.char
  e = zero
  project0 do |c, ind|
    c * index_eval(self.class, ary, ind.collect { |i| i / ch })
  end
end