basis/math/polynomials/polynomials.factor

   1 ! Copyright (C) 2008 Doug Coleman.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: arrays combinators kernel make math math.bits
   4 math.vectors sequences vectors ;
   5 IN: math.polynomials
   6
   7 <PRIVATE
   8
   9 : 2pad-head ( p q n -- p q ) [ 0 pad-head ] curry bi@ ;
  10 : 2pad-tail ( p q n -- p q ) [ 0 pad-tail ] curry bi@ ;
  11 : pextend ( p q -- p q ) 2dup max-length 2pad-tail ;
  12 : pextend-left ( p q -- p q ) 2dup max-length 2pad-head ;
  13 : unempty ( seq -- seq ) [ { 0 } ] when-empty ;
  14 : 2unempty ( seq seq -- seq seq ) [ unempty ] bi@ ;
  15
  16 PRIVATE>
  17
  18 : powers ( n x -- seq )
  19     <repetition> 1 [ * ] accumulate nip ;
  20
  21 : p= ( p q -- ? ) pextend = ;
  22
  23 : ptrim ( p -- q )
  24     dup length 1 = [ [ zero? ] trim-tail ] unless ;
  25
  26 : 2ptrim ( p q -- p' q' ) [ ptrim ] bi@ ;
  27 : p+ ( p q -- r ) pextend v+ ;
  28 : p- ( p q -- r ) pextend v- ;
  29 ALIAS: n*p n*v
  30
  31 : pextend-conv ( p q -- p' q' )
  32     2dup [ length ] bi@ + 1 - 2pad-tail ;
  33
  34 : p* ( p q -- r )
  35     2unempty pextend-conv
  36     [ drop length [ <iota> ] keep ]
  37     [ nip <reversed> ]
  38     [ drop ] 2tri
  39     '[ _ _ <slice> _ vdot ] map reverse! ;
  40
  41 : p-sq ( p -- p^2 ) dup p* ; inline
  42
  43 ERROR: negative-power-polynomial p n ;
  44
  45 : (p^) ( p n  -- p^n )
  46     make-bits { 1 } [ [ over p* ] when [ p-sq ] dip ] reduce nip ;
  47
  48 : p^ ( p n -- p^n )
  49     dup 0 >= [ (p^) ] [ negative-power-polynomial ] if ;
  50
  51 <PRIVATE
  52
  53 : p/mod-setup ( p p -- p p n )
  54     2ptrim
  55     2dup [ length ] bi@ -
  56     dup 1 < [ drop 1 ] when
  57     [ over length + 0 pad-head pextend ] keep 1 + ;
  58
  59 : /-last ( seq1 seq2 -- x ) [ last ] bi@ / ;
  60
  61 : (p/mod) ( p p -- p p )
  62     2dup /-last
  63     2dup , n*p swapd
  64     p- >vector
  65     dup pop* swap rest-slice ;
  66
  67 PRIVATE>
  68
  69 : p/mod ( p q -- z w )
  70     p/mod-setup [ [ (p/mod) ] times ] V{ } make
  71     reverse nip swap 2ptrim pextend ;
  72
  73 <PRIVATE
  74
  75 : (pgcd) ( b a y x -- a d )
  76     dup V{ 0 } p= [
  77         drop nip
  78     ] [
  79         [ nip ] [ p/mod ] 2bi
  80         [ pick p* swap [ swapd p- ] dip ] dip (pgcd)
  81     ] if ;
  82
  83 PRIVATE>
  84
  85 : pgcd ( p q -- a d )
  86     [ V{ 0 } clone V{ 1 } clone ] 2dip swap (pgcd) [ >array ] bi@ ;
  87
  88 : pdiff ( p -- p' )
  89     dup length <iota> v* rest ;
  90
  91 : polyval ( x p -- p[x] )
  92     ! Horner scheme
  93     [ nip <reversed> unclip-slice swap ]
  94     [ drop ] 2bi
  95     '[ [ _ * ] dip + ] each ;
  96
  97 MACRO: polyval* ( p -- quot )
  98     reverse
  99     [ rest [ \ * swap \ + [ ] 3sequence ] map ]
 100     [ first \ drop swap [ ] 2sequence ] bi
 101     prefix \ cleave [ ] 2sequence ;