Module: Algebra::MPolynomialFactorization

Includes:
ChineseRemainderTheorem, Q, PolynomialFactorization::Z
Included in:
MPolynomial
Defined in:
lib/algebra/m-polynomial-factor.rb,
lib/algebra/m-polynomial-factor-zp.rb,
lib/algebra/m-polynomial-factor-int.rb

Defined Under Namespace

Modules: Q

Instance Method Summary collapse

Methods included from Q

#contQ, #factorize_rational, #ppQ

Methods included from PolynomialFactorization::Z

#divmod_ED, #hensel_lift, #lifting, #reduce_over_prime_field, #resolve_c, #resultant_by_elimination, #sqfree_over_integral, #try_e

Methods included from ChineseRemainderTheorem

#decompose_on_cofactors_modp, #decompose_on_cofactors_p_adic, #decompose_on_factors, #mk_cofacts

Instance Method Details

#_factorize(na = 0) ⇒ Object 



41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
 
# 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) ⇒ Object

def _hensel_lift(g0, f0, char, height)



97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
 
# 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 ⇒ Object 



147
148
149
150
151
152
153
154
155
156
 
# 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) ⇒ Object 



92
93
94
95
96
97
98
99
 
# 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) ⇒ Object 



17
18
19
20
21
22
23
24
25
26
27
28
29
 
# 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) ⇒ Object 



31
32
33
34
35
36
37
38
39
 
# 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) ⇒ Object

include PolynomialFactorization::Q



21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
 
# 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) ⇒ Object

include PolynomialFactorization::Z



19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
 
# 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) ⇒ Object 



45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
 
# 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) ⇒ Object

def hensel_lift_int(g0, f0, char, height)



71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
 
# 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) ⇒ Object

def hensel_lift_zp(g0, f0, char, height)



59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
 
# 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) ⇒ Object 



164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
 
# 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) ⇒ Object 



158
159
160
161
162
 
# 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 ⇒ Object 



41
42
43
 
# File 'lib/algebra/m-polynomial-factor-int.rb', line 41

def sqfree
  _sqfree
end

#sqfree_zp ⇒ Object 



47
48
49
50
51
52
53
54
55
56
 
# 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 ⇒ Object 



38
39
40
41
42
43
44
45
 
# 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