contrib/math/polynomials.factor

   1 IN: polynomials-internals
   2 USING: arrays kernel sequences vectors math math-internals namespaces arrays
   3     sequences-contrib ;
   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 : 2pad-left ( p p n -- p p ) 0 [ pad-left swap ] 2keep pad-left swap ;
  10 : 2pad-right ( p p n -- p p ) 0 [ pad-right swap ] 2keep pad-right swap ;
  11 : pextend ( p p -- p p ) 2dup [ length ] 2apply max 2pad-right ;
  12 : pextend-left ( p p -- p p ) 2dup [ length ] 2apply max 2pad-left ;
  13 : unempty ( seq -- seq ) dup empty? [ drop { 0 } ] when ;
  14 : 2unempty ( seq seq -- seq seq ) [ unempty ] 2apply ;
  15
  16 IN: math-contrib
  17 : p= ( p p -- ? ) pextend = ;
  18
  19 : ptrim ( p -- p ) [ zero? ] rtrim* ;
  20
  21 : 2ptrim ( p p -- p p ) [ ptrim ] 2apply ;
  22 : p+ ( p p -- p ) pextend v+ ;
  23 : p- ( p p -- p ) pextend v- ;
  24 : n*p ( n p -- n*p ) n*v ;
  25
  26 ! convolution
  27 : pextend-conv ( p p -- p p )
  28     #! extend to: p_m + p_n - 1
  29     2dup [ length ] 2apply + 1- 2pad-right [ >vector ] 2apply ;
  30
  31 : p* ( p p -- p )
  32     #! Multiply two polynomials.
  33     2unempty pextend-conv <reversed> dup length
  34     [ over length pick <slice> pick [ * ] 2map sum ] map 2nip reverse ;
  35
  36 : p-sq ( p -- p-sq )
  37     dup p* ;
  38
  39 IN: polynomials-internals
  40
  41 : pop-front ( seq -- seq )
  42     1 tail-slice ;
  43
  44 : /-last ( seq seq -- a )
  45     #! divide the last two numbers in the sequences
  46     [ peek ] 2apply / ;
  47
  48 : p/mod-setup ( p p -- p p n )
  49     2ptrim 2dup [ length ] 2apply - dup 1 < [ drop 1 ] when
  50     dup >r over length + 0 pad-left pextend r> 1+ ;
  51
  52 : (p/mod)
  53     2dup /-last 2dup , n*p swapd p- >vector dup pop drop swap pop-front ;
  54
  55 IN: math-contrib
  56 : p/mod
  57     p/mod-setup [ [ (p/mod) ] times ] V{ } make
  58     reverse nip swap 2ptrim pextend ;
  59
  60 : (pgcd) ( b a y x -- a d )
  61     dup V{ 0 } clone p= [
  62         drop nip
  63     ] [
  64         tuck p/mod >r pick p* swap >r swapd p- r> r> (pgcd)
  65     ] if ;
  66
  67 : pgcd ( p p -- p q )
  68     swap V{ 0 } clone V{ 1 } clone 2swap (pgcd) [ >array ] 2apply ;
  69
  70 : pdiff ( p -- p' )
  71     #! Polynomial derivative.
  72     dup length v* { 0 } ?head drop ;
  73
  74 : polyval ( p x -- p[x] )
  75     #! Evaluate a polynomial.
  76     >r dup length r> powers v. ;