Module: Algebra::PolynomialFactorization::Z

Included in:
MPolynomialFactorization, Algebra::PolynomialFactorization
Defined in:
lib/algebra/polynomial-factor-int.rb

Instance Method Summary collapse

Instance Method Details

#_factorize ⇒ Object 



135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
 
# 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) ⇒ Object 



25
26
27
28
 
# 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) ⇒ Object

divmod in the ring over Euclidian Domain

Raises:

  • (ZeroDivisionError) 


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
 
# 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 ⇒ Object 



165
166
167
 
# File 'lib/algebra/polynomial-factor-int.rb', line 165

def factorize_int
  _factorize.first
end

#hensel_lift(g0) ⇒ Object 



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
 
# 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) ⇒ Object 



97
98
99
 
# 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) ⇒ Object 



84
85
86
87
88
89
90
91
92
93
94
95
 
# 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) ⇒ Object

find [a, b] s.t. a * g + b * h == c



13
14
15
16
17
18
19
20
21
22
23
 
# 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) ⇒ Object

ground ring must be a field



80
81
82
 
# 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 ⇒ Object 



44
45
46
47
48
49
 
# 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) ⇒ Object 



30
31
32
33
34
35
36
37
38
39
40
41
42
 
# 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