Class: SyMath::Poly::DUP

Inherits:
SyMath::Poly show all
Defined in:
lib/symath/poly/dup.rb

Overview

Class representing a univariate polynomial. The polynomial is represented by an array in of the polynomial coefficients in ‘dense form’, i.e. zero-coefficients are included in the list. The coefficients are stored in decreasing order starting with the highest degree, and ending with the constant.

Instance Attribute Summary

Attributes inherited from SyMath::Poly

#arr, #p, #var

Instance Method Summary collapse

Methods inherited from SyMath::Poly

#%, #*, #**, #+, #-, #-@, #/, #==, #[], #degree, #lc, #sort_factors, #sort_factors_multiple, #strip!, #to_m, #zero?

Constructor Details

#initialize(args) ⇒ DUP

Returns a new instance of DUP.



12
13
14
15
16
17
18
19
20
21
22
23
24
 
# File 'lib/symath/poly/dup.rb', line 12

def initialize(args)
  if args.is_a?(SyMath::Value)
    init_from_exp(args)
    return
  end

  if args.key?(:arr)
    init_from_array(args[:arr], args[:var])
    return
  end

  raise 'Bad arguments for Poly::DUP constructor'
end

Instance Method Details

#add(g) ⇒ Object

Sum two polynomials



634
635
636
637
638
639
640
641
642
643
644
645
646
 
# File 'lib/symath/poly/dup.rb', line 634

def add(g)
  ret = @arr.clone

  if g.degree > degree
    (g.degree - degree).times { ret.unshift(0) }
  end

  (0..g.degree).each do |i|
    ret[ret.size - i - 1] += g[g.degree - i]
  end
  
  return new_dup(ret).strip!
end

#content ⇒ Object 



571
572
573
574
575
576
577
578
579
580
581
 
# File 'lib/symath/poly/dup.rb', line 571

def content()
  cont = 0

  @arr.each do |c|
    cont = cont.gcd(c)

    break if cont == 1
  end

  return cont
end

#diff ⇒ Object

Fast differentiation of polynomial



430
431
432
433
434
435
436
 
# File 'lib/symath/poly/dup.rb', line 430

def diff
  d = degree
  res = @arr.each_with_index.map { |e, i| e*(d - i) }
  res.pop

  return new_dup(res)
end

#div(g) ⇒ Object 



767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
 
# File 'lib/symath/poly/dup.rb', line 767

def div(g)
  # returns qv and r such that:
  # f = fv*qv + r
  df = degree
  dg = g.degree

  if g.zero?
    raise 'Division by zero'
  elsif df < dg
    return [zero, self.clone]  # no quotient, f is remainder
  end
  
  # Start with f as remainder, no quotient
  q = zero
  r = self
  dr = df

  lc_g = g.lc
  
  while true
    lc_r = r.lc

    if (lc_r % lc_g) != 0
      break
    end

    c = lc_r / lc_g
    j = dr - dg

    q = q + one.mul_term(c, j)
    r = r - g.mul_term(c, j)

    _dr = dr
    dr = r.degree

    if dr < dg
      break
    elsif dr >= _dr
      raise 'Polynomial division failed'
    end
  end

  return [q, r]
end

#factor ⇒ Object 



123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
 
# File 'lib/symath/poly/dup.rb', line 123

def factor
  # Transform to primitive form
  (cont, g) = primitive

  # Transform left coefficient to positive
  if g.lc < 0
    cont = -cont
    g = -g
  end

  n = g.degree

  # Handle some rivial cases
  if n <= 0
    # 0, cont
    return [cont, []]
  elsif n == 1
    # cont*(a*x + b)
    return [cont, [[g, 1]]]
  end

  # Remove square factors
  g = g.sqf_part

  # Use the zassenhaus algorithm to compute candidate factors
  h = g.zassenhaus

  # Check each of the candidate factors by dividing the original
  # polynomial with each of them until the quotient is 1.
  factors = trial_division(h)

  return [cont, factors]
end

#gcd(g) ⇒ Object

FIXME: Implement heuristic gcd?



462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
 
# File 'lib/symath/poly/dup.rb', line 462

def gcd(g)
  # Trivial cases
  if zero? and g.zero?
    return [zero, zero, zero]
  end
  if zero?
    if g.lc >= 0
      return [g, zero, one]
    else
      return [-g, zero, minus_one]
    end
  end

  if g.zero?
    if lc >= 0
      return [f, one, zero]
    else
      return [-f, zero, one]
    end
  end

  (fc, ff) = primitive
  (gc, gg) = g.primitive

  c = fc.gcd(gc)
  
  h = subresultants(g)[-1]
  h = h.primitive[1]

  if (h.lc < 0)
    c = -c
  end

  h = h*c
  
  cff = self/h
  cfg = g/h

  return [h, cff, cfg]
end

#gcdext(x, y) ⇒ Object

Extended gcd for two integers



439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
 
# File 'lib/symath/poly/dup.rb', line 439

def gcdext(x, y)
  if x < 0
    g, a, b = gcdext(-x, y)
    return [g, -a, b]
  end
  if y < 0
    g, a, b = gcdext(x, -y)
    return [g, a, -b]
  end
  r0, r1 = x, y
  a0 = b1 = 1
  a1 = b0 = 0
  until r1.zero?
    q = r0 / r1
    r0, r1 = r1, r0 - q*r1
    a0, a1 = a1, a0 - q*a1
    b0, b1 = b1, b0 - q*b1
  end

  return [r0, a0, b0]
end

#hensel_lift(p, f_list, l) ⇒ Object 



211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
 
# File 'lib/symath/poly/dup.rb', line 211

def hensel_lift(p, f_list, l)
  r = f_list.size
  lcf = lc

  if r == 1
    ff = self * gcdext(lcf, p**l)[1]
    return [ ff % (p**l) ]
  end

  m = p
  k = r / 2
  d = CMath.log(l, 2).ceil

  g = new_dup([lcf]).to_galois(p)

  f_list[0..k - 1].each do |f_i|
    g = g*f_i.to_galois(p)
  end

  h = f_list[k].to_galois(p)

  f_list[k + 1..-1].each do |f_i|
    h = h*f_i.to_galois(p)
  end

  (s, t, x) = g.gcdex(h)

  g = g.to_dup
  h = h.to_dup
  s = s.to_dup
  t = t.to_dup

  d.times do
    (g, h, s, t) = hensel_step(m, g, h, s, t)
    m = m**2
  end

  return g.hensel_lift(p, f_list[0..k - 1], l) +
         h.hensel_lift(p, f_list[k..-1], l)
end

#hensel_step(m, g, h, s, t) ⇒ Object 



181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
 
# File 'lib/symath/poly/dup.rb', line 181

def hensel_step(m, g, h, s, t)
  mm = m**2

  e = (self - g*h) % mm

  (q, r) = (s*e).div(h)

  q = q % mm
  r = r % mm

  u = t*e + q*g

  gg = (g + u) % mm
  hh = (h + r) % mm

  u = s*gg + t*hh
  b = (u - one) % mm

  (c, d) = (s*b).div(hh)

  c = c % mm
  d = d % mm

  u = t*b + c*gg
  ss = (s - d) % mm
  tt = (t - u) % mm

  return [gg, hh, ss, tt]
end

#init_from_array(arr, var) ⇒ Object 



26
27
28
29
 
# File 'lib/symath/poly/dup.rb', line 26

def init_from_array(arr, var)
  @arr = arr
  @var = var
end

#init_from_exp(e) ⇒ Object 



31
32
33
34
35
36
37
38
39
40
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
95
96
97
98
 
# File 'lib/symath/poly/dup.rb', line 31

def init_from_exp(e)
  # From this point, expect e to be a SyMath::Value
  error = 'Expression ' + e.to_s + ' is not an univariate polynomial'
  max_degree = 0
  terms = {}

  # Assert that this really is an univariate polynomial, and build
  # the dup array representation.
  e.terms.each do |s|
    var = nil
    c = 1
    d = 0
    
    s.factors.each do |f|
      if f.is_a?(SyMath::Fraction)
        raise error
      end

      if f == -1
        c *= -1
      elsif f.is_number?
        c *= f.value
      else
        if !var.nil?
          raise error
        end

        var = f.base
        if !var.is_a?(SyMath::Definition::Variable)
          raise error
        end

        if !f.exponent.is_number?
          raise error
        end
        d = f.exponent.value
      end
    end

    if !var.nil?
      if @var.nil?
        @var = var
      elsif @var != var
        raise error
      end
    end

    if terms.key?(d)
      terms[d] += c
    else
      terms[d] = c
    end

    if d > max_degree
      max_degree = d
    end
  end

  @arr = (0..max_degree).to_a.reverse.map do |d|
    if terms.key?(d)
      terms[d]
    else
      0
    end
  end

  strip!
end

#l1_norm ⇒ Object 



627
628
629
630
631
 
# File 'lib/symath/poly/dup.rb', line 627

def l1_norm()
  return 0 if zero?

  return @arr.map { |e| e.abs }.inject(:+)
end

#lshift(n) ⇒ Object 



583
584
585
586
587
588
589
 
# File 'lib/symath/poly/dup.rb', line 583

def lshift(n)
  if zero?
    return zero
  else
    return new_dup(@arr + [0]*n)
  end
end

#max_norm ⇒ Object 



823
824
825
826
827
828
829
 
# File 'lib/symath/poly/dup.rb', line 823

def max_norm()
  if zero?
    return 0
  else
    return @arr.map { |e| e.abs }.max
  end
end

#minus_one ⇒ Object 



119
120
121
 
# File 'lib/symath/poly/dup.rb', line 119

def minus_one()
  return new_dup([-1])
end

#mul(g) ⇒ Object 



668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
 
# File 'lib/symath/poly/dup.rb', line 668

def mul(g)
  if self == g
    return self**2
  end

  if zero? and g.zero?
    return zero
  end

  df = degree
  dg = g.degree

  n = [df, dg].max + 1

  if n < 100
    h = []

    (0..df + dg).each do |i|
      coeff = 0

      a = [0, i - dg].max
      b = [df, i].min
      (a..b).each do |j|
        coeff += self[j]*g[i - j]
      end
      h << coeff
    end

    return new_dup(h).strip!
  else
    # Use Karatsuba's algorithm for large polygons.
    n2 = n/2

    fl = slice(0, n2)
    gl = g.slice(0, n2)

    fh = slice(n2, n).rshift(n2)
    gh = g.slice(n2, n).rshift(n2)

    lo = fl*gl
    hi = fh*gh

    mid = (fl + fh)*(gl + gh)
    mid -= lo + hi

    return lo + mid.lshift(n2) + hi.lshift(2*n2)
  end
end

#mul_ground(c) ⇒ Object 



717
718
719
720
 
# File 'lib/symath/poly/dup.rb', line 717

def mul_ground(c)
  ret = @arr.map { |t| t*c }
  return new_dup(ret).strip!
end

#mul_term(c, j) ⇒ Object

Multiply a polynomial with a single term.



723
724
725
726
727
728
 
# File 'lib/symath/poly/dup.rb', line 723

def mul_term(c, j)
  ret = @arr.map { |t| t*c }
  j.times { ret.push(0) }

  return new_dup(ret).strip!
end

#neg ⇒ Object

Return the negative of the polynomial



664
665
666
 
# File 'lib/symath/poly/dup.rb', line 664

def neg()
  return new_dup(@arr.map { |t| -t })
end

#new_dup(arr) ⇒ Object

Convenience method. Returns a a new instance from array, on the same variable as this instance.



106
107
108
 
# File 'lib/symath/poly/dup.rb', line 106

def new_dup(arr)
  return self.class.new({ :arr => arr, :var => @var })
end

#one ⇒ Object 



115
116
117
 
# File 'lib/symath/poly/dup.rb', line 115

def one()
  return new_dup([1])
end

#primitive ⇒ Object

Compute content and primitive form



557
558
559
560
561
562
563
564
565
566
567
568
569
 
# File 'lib/symath/poly/dup.rb', line 557

def primitive()
  if zero?
    return [0, self.clone]
  end

  cont = content

  if cont == 1
    return [cont, self.clone]
  else
    return [cont, self/cont]
  end
end

#pseudo_rem(g) ⇒ Object

Compute pseudo remainder of self / g



731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
 
# File 'lib/symath/poly/dup.rb', line 731

def pseudo_rem(g)
  df = self.degree
  dg = g.degree

  r = self
  dr = df

  if g.zero?
    raise 'Division by zero'
  elsif df < dg
    return self.clone # self is remainder
  end

  n = df - dg + 1

  while true
    j = dr - dg
    n -= 1

    rr = r*g.lc
    gg = g.mul_term(r.lc, j)
    r = rr - gg

    _dr = dr
    dr = r.degree

    if dr < dg
      break
    elsif !(dr < _dr)
      raise 'Polynomial division failed'
    end
  end  

  return r*g.lc**n
end

#quo(g) ⇒ Object 



812
813
814
 
# File 'lib/symath/poly/dup.rb', line 812

def quo(g)
  return div(g)[0]
end

#quo_ground(c) ⇒ Object

Quotient by constant for each coefficient



817
818
819
820
821
 
# File 'lib/symath/poly/dup.rb', line 817

def quo_ground(c)
  ret = @arr.map { |t| t.to_i/c.to_i }

  return new_dup(ret).strip!
end

#rshift(n) ⇒ Object 



591
592
593
 
# File 'lib/symath/poly/dup.rb', line 591

def rshift(n)
  return new_dup(@arr[0..-n - 1])
end

#slice(a, b) ⇒ Object

Slice of polynomial between two degrees, >= a and < b



596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
 
# File 'lib/symath/poly/dup.rb', line 596

def slice(a, b)
  s = @arr.size
  aa = [0, s - a].max
  bb = [0, s - b].max

  if aa <= 0
    return zero
  end
  
  ret = @arr[bb..aa-1]

  if ret == []
    return zero
  else
    return new_dup(ret + [0]*a)
  end
end

#sqf_list ⇒ Object

Decompose a polynomial into square free components



391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
 
# File 'lib/symath/poly/dup.rb', line 391

def sqf_list()
  (coeff, f) = primitive

  if f.lc < 0
    f = -f
    coeff = -coeff
  end

  # Trivial case, constant polynomial
  if f.degree <= 0
    return coeff, []
  end

  res = []
  i = 1

  (g, p, q) = f.gcd(f.diff)

  while true
    h = q - p.diff

    if h.zero?
      res << [p, i]
      break
    end

    (g, p, q) = p.gcd(h)

    if g.degree > 0
      res << [g, i]
    end

    i += 1
  end

  return [coeff, res]
end

#sqf_part ⇒ Object

Return square free part of f



373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
 
# File 'lib/symath/poly/dup.rb', line 373

def sqf_part()
  # Trivial case
  if zero?
    return self.clone
  end

  if lc < 0
    f = -self
  else
    f = self
  end

  sqf = f/f.gcd(f.diff)[0]

  return sqf.primitive[1]
end

#sub(g) ⇒ Object

Subtract a polynomial from this one



649
650
651
652
653
654
655
656
657
658
659
660
661
 
# File 'lib/symath/poly/dup.rb', line 649

def sub(g)
  ret = @arr.clone

  if g.degree > degree
    (g.degree - degree).times { ret.unshift(0) }
  end

  (0..g.degree).each do |i|
    ret[ret.size - i - 1] -= g[g.degree - i]
  end
  
  return new_dup(ret).strip!
end

#subresultants(g) ⇒ Object

Calculate subresultants of polynomials self and g



504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
 
# File 'lib/symath/poly/dup.rb', line 504

def subresultants(g)
  f = self
  
  n = f.degree
  m = g.degree

  if n < m
    f, g = g, f
    n, m = m, n
  end
  
  if f.zero?
    return []
  end

  if g.zero?
    return [f.clone]
  end

  r = [f.clone, g.clone]
  d = n - m

  b = (-1)**(d + 1)

  h = f.pseudo_rem(g)*b

  lc = g.lc
  c = -(lc**d)

  while !h.zero?
    k = h.degree
    r << h

    f, g, m, d = g, h, k, m - k

    b = -lc * c**d

    h = f.pseudo_rem(g)/b

    lc = g.lc

    if d > 1
      q = c**(d - 1)
      c = ((-lc)**d).to_i/q.to_i
    else
      c = -lc
    end
  end
  
  return r          
end

#to_galois(p) ⇒ Object 



100
101
102
 
# File 'lib/symath/poly/dup.rb', line 100

def to_galois(p)
  return SyMath::Poly::Galois.new({ :dup => self, :p => p })      
end

#to_s ⇒ Object 



831
832
833
 
# File 'lib/symath/poly/dup.rb', line 831

def to_s()
  return @arr.to_s
end

#trial_division(factors) ⇒ Object 



157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
 
# File 'lib/symath/poly/dup.rb', line 157

def trial_division(factors)
  result = []
  f = self

  factors.each do |factor|
    k = 0

    while true
      (q, r) = f.div(factor)

      if r.zero?
        f = q
        k += 1
      else
        break
      end
    end

    result << [factor, k]
  end

  return sort_factors_multiple(result)
end

#trunc(p) ⇒ Object 



614
615
616
617
618
619
620
621
622
623
624
625
 
# File 'lib/symath/poly/dup.rb', line 614

def trunc(p)
  ret = @arr.map do |e|
    ep = e % p
    if ep > p / 2
      ep - p
    else
      ep
    end
  end
                
  return new_dup(ret).strip!
end

#zassenhaus ⇒ Object

Zassenhaus algorithm for factorizing square free polynomial



253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
 
# File 'lib/symath/poly/dup.rb', line 253

def zassenhaus()
  n = degree

  # Trivial case, a*x + b
  if n == 1
    return [self.clone]
  end

  # Calculate bound of px:
  # n = deg(f)
  # B = sqrt(n + 1)*2^n*max_norm(f)*lc(f)
  # C = (n + 1)^2n*A^(2n - 1)
  # gm = 2log(cc,2)
  # bound = 2*gm*ln(gm)
  fc = self[-1]
  aa = max_norm
  b = lc
  bb = (CMath.sqrt(n + 1).floor*2**n*aa*b).abs.to_i # Integer square root??
  cc = ((n + 1)**(2*n)*aa**(2*n - 1)).to_i
  gamma = (2*CMath.log(cc, 2)).ceil
  bound = (2*gamma*CMath.log(gamma)).to_i

  a = []

  # Choose a prime number p such that f be square free in Z_p
  # if there are many factors in Z_p, choose among a few different p
  # the one with fewer factors
  (3..bound).each do |px|
    # Skip non prime px and px which do not divide lc(f)
    if !Prime.prime?(px) or (b % px) == 0
      next
    end

    # px = convert(px) ???

    # Convert f to a galois field of order px
    ff = self.to_galois(px)

    # Skip if ff has square factors
    if !ff.sqf_p
      next
    end

    # Factorize ff and store all factors together with its order px
    fsqfx = ff.factor_sqf[1]
    a << [px, fsqfx]

    if fsqfx.size < 15 or a.size > 4
      break
    end
  end

  # Select the factor list with the fewest factors.
  (p, fsqf) = a.min { |x| x[1].size }
  l = CMath.log(2*bb + 1, p).ceil

  # Convert the factors back to integer polynomials
  modular = fsqf.map { |ff| ff.to_dup }

  # Hensel lift of modular -> g
  g = hensel_lift(p, modular, l)

  # Start with T as the set of factors in array g.
  tt = (0..g.size - 1).to_a
  factors = []
  s = 1
  pl = p**l  # pl =~ 2*bb + 1

  f = self

  while 2*s <= tt.size
    inc_s = true

    tt.combination(s).each do |ss|
      # Calculate G as the product of the subset S of factors. Lift
      # the constant coefficient of G.
      gg = new_dup([b])
      ss.each { |i| gg = gg*g[i] }
      gg = (gg % pl).primitive[1]
      q = gg[-1]

      # If it does not divide the input polynomial constant (fc), G
      # does not divide the input polynomial.
      if q != 0 and fc % q != 0
        next
      end
      
      tt_new = tt - ss

      # Calculate H as the product of the remaining factors in T.
      hh = new_dup([b])
      tt_new.each { |i| hh = hh*g[i] }
      hh = hh % pl

      # If the norm of the candidate G and the remaining H are bigger than
      # the bound B, we have a valid candidate.
      # - Store it in the factors list
      # - Remove its corresponding selection from T
      # - Continue with H as the remaining polynomial
      if gg.l1_norm*hh.l1_norm <= bb
        tt = tt_new

        gg = gg.primitive[1]
        f = hh.primitive[1]

        factors << gg
        b = f.lc

        inc_s = false
        break
      end
    end

    s += 1 if inc_s
  end

  return factors + [f]
end

#zero ⇒ Object

Convenience methods for creating some commonly used constant polynomials



111
112
113
 
# File 'lib/symath/poly/dup.rb', line 111

def zero()
  return new_dup([0])
end