extra/math/polynomials/polynomials.factor

   1 USING: arrays kernel sequences vectors math math.vectors namespaces
   2 shuffle splitting sequences.lib math.order ;
   3 IN: math.polynomials
   4
   5 ! Polynomials are vectors with the highest powers on the right:
   6 ! { 1 1 0 1 } -> 1 + x + x^3
   7 ! { } -> 0
   8
   9 : powers ( n x -- seq )
  10     #! Output sequence has n elements, { 1 x x^2 x^3 ... }
  11     <array> 1 [ * ] accumulate nip ;
  12
  13 <PRIVATE
  14 : 2pad-left ( p p n -- p p ) 0 [ pad-left swap ] 2keep pad-left swap ;
  15 : 2pad-right ( p p n -- p p ) 0 [ pad-right swap ] 2keep pad-right swap ;
  16 : pextend ( p p -- p p ) 2dup [ length ] bi@ max 2pad-right ;
  17 : pextend-left ( p p -- p p ) 2dup [ length ] bi@ max 2pad-left ;
  18 : unempty ( seq -- seq ) [ { 0 } ] when-empty ;
  19 : 2unempty ( seq seq -- seq seq ) [ unempty ] bi@ ;
  20
  21 PRIVATE>
  22 : p= ( p p -- ? ) pextend = ;
  23
  24 : ptrim ( p -- p )
  25     dup length 1 = [ [ zero? ] trim-right ] unless ;
  26
  27 : 2ptrim ( p p -- p p ) [ ptrim ] bi@ ;
  28 : p+ ( p p -- p ) pextend v+ ;
  29 : p- ( p p -- p ) pextend v- ;
  30 : n*p ( n p -- n*p ) n*v ;
  31
  32 ! convolution
  33 : pextend-conv ( p p -- p p )
  34     #! extend to: p_m + p_n - 1
  35     2dup [ length ] bi@ + 1- 2pad-right [ >vector ] bi@ ;
  36
  37 : p* ( p p -- p )
  38     #! Multiply two polynomials.
  39     2unempty pextend-conv <reversed> dup length
  40     [ over length pick <slice> pick [ * ] 2map sum ] map 2nip reverse ;
  41
  42 : p-sq ( p -- p-sq )
  43     dup p* ;
  44
  45 <PRIVATE
  46
  47 : p/mod-setup ( p p -- p p n )
  48     2ptrim
  49     2dup [ length ] bi@ -
  50     dup 1 < [ drop 1 ] when
  51     [ over length + 0 pad-left pextend ] keep 1+ ;
  52
  53 : /-last ( seq seq -- a )
  54     #! divide the last two numbers in the sequences
  55     [ peek ] bi@ / ;
  56
  57 : (p/mod) ( p p -- p p )
  58     2dup /-last
  59     2dup , n*p swapd
  60     p- >vector
  61     dup pop* swap rest-slice ;
  62
  63 PRIVATE>
  64
  65 : p/mod ( a b -- / mod )
  66     p/mod-setup [ [ (p/mod) ] times ] V{ } make
  67     reverse nip swap 2ptrim pextend ;
  68
  69 : (pgcd) ( b a y x -- a d )
  70     dup V{ 0 } clone p= [
  71         drop nip
  72     ] [
  73         tuck p/mod >r pick p* swap >r swapd p- r> r> (pgcd)
  74     ] if ;
  75
  76 : pgcd ( p p -- p q )
  77     swap V{ 0 } clone V{ 1 } clone 2swap (pgcd) [ >array ] bi@ ;
  78
  79 : pdiff ( p -- p' )
  80     #! Polynomial derivative.
  81     dup length v* { 0 } ?head drop ;
  82
  83 : polyval ( p x -- p[x] )
  84     #! Evaluate a polynomial.
  85     >r dup length r> powers v. ;
  86