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 fry kernel macros make math math.bits
   4 math.order math.vectors sequences splitting 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 [ length ] bi@ max 2pad-tail ;
  12 : pextend-left ( p q -- p q ) 2dup [ length ] bi@ max 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> _ v* sum ] 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 >=
  50     [ (p^) ]
  51     [ negative-power-polynomial ] if ;
  52
  53 <PRIVATE
  54
  55 : p/mod-setup ( p p -- p p n )
  56     2ptrim
  57     2dup [ length ] bi@ -
  58     dup 1 < [ drop 1 ] when
  59     [ over length + 0 pad-head pextend ] keep 1 + ;
  60
  61 : /-last ( seq1 seq2 -- x ) [ last ] bi@ / ;
  62
  63 : (p/mod) ( p p -- p p )
  64     2dup /-last
  65     2dup , n*p swapd
  66     p- >vector
  67     dup pop* swap rest-slice ;
  68
  69 PRIVATE>
  70
  71 : p/mod ( p q -- z w )
  72     p/mod-setup [ [ (p/mod) ] times ] V{ } make
  73     reverse nip swap 2ptrim pextend ;
  74
  75 <PRIVATE
  76
  77 : (pgcd) ( b a y x -- a d )
  78     dup V{ 0 } p= [
  79         drop nip
  80     ] [
  81         [ nip ] [ p/mod ] 2bi
  82         [ pick p* swap [ swapd p- ] dip ] dip (pgcd)
  83     ] if ;
  84
  85 PRIVATE>
  86
  87 : pgcd ( p q -- a d )
  88     [ V{ 0 } clone V{ 1 } clone ] 2dip swap (pgcd) [ >array ] bi@ ;
  89
  90 : pdiff ( p -- p' )
  91     dup length iota v* rest ;
  92
  93 : polyval ( x p -- p[x] )
  94     ! Horner scheme
  95     [ nip <reversed> unclip-slice swap ]
  96     [ drop ] 2bi
  97     '[ [ _ * ] dip + ] each ;
  98
  99 MACRO: polyval* ( p -- )
 100     reverse
 101     [ rest [ \ * swap \ + [ ] 3sequence ] map ]
 102     [ first \ drop swap [ ] 2sequence ] bi
 103     prefix \ cleave [ ] 2sequence ;
 104