basis/math/polynomials/polynomials.factor

   1 ! Copyright (C) 2008 Doug Coleman.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: arrays kernel make math math.order math.vectors sequences
   4     splitting vectors macros combinators ;
   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 -- p )
  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 : n*p ( n p -- n*p ) n*v ;
  30
  31 : pextend-conv ( p q -- p q )
  32     2dup [ length ] bi@ + 1 - 2pad-tail [ >vector ] bi@ ;
  33
  34 : p* ( p q -- r )
  35     2unempty pextend-conv <reversed> dup length iota
  36     [ over length pick <slice> pick [ * ] 2map sum ] map 2nip reverse ;
  37
  38 : p-sq ( p -- p^2 )
  39     dup p* ;
  40
  41 ERROR: negative-power-polynomial p n ;
  42
  43 : p^ ( p n -- p^n )
  44     {
  45         { [ dup 0 > ] [ 1 - dupd [ p* ] with times ] }
  46         { [ dup 0 = ] [ 2drop { 1 } ] }
  47         { [ dup 0 < ] [ negative-power-polynomial ] }
  48     } cond ;
  49
  50 <PRIVATE
  51
  52 : p/mod-setup ( p p -- p p n )
  53     2ptrim
  54     2dup [ length ] bi@ -
  55     dup 1 < [ drop 1 ] when
  56     [ over length + 0 pad-head pextend ] keep 1 + ;
  57
  58 : /-last ( seq seq -- a )
  59     #! divide the last two numbers in the sequences
  60     [ last ] bi@ / ;
  61
  62 : (p/mod) ( p p -- p p )
  63     2dup /-last
  64     2dup , n*p swapd
  65     p- >vector
  66     dup pop* swap rest-slice ;
  67
  68 PRIVATE>
  69
  70 : p/mod ( p q -- z w )
  71     p/mod-setup [ [ (p/mod) ] times ] V{ } make
  72     reverse nip swap 2ptrim pextend ;
  73
  74 <PRIVATE
  75
  76 : (pgcd) ( b a y x -- a d )
  77     dup V{ 0 } clone p= [
  78         drop nip
  79     ] [
  80         [ nip ] [ p/mod ] 2bi
  81         [ pick p* swap [ swapd p- ] dip ] dip (pgcd)
  82     ] if ;
  83
  84 PRIVATE>
  85
  86 : pgcd ( p q -- a d )
  87     [ V{ 0 } clone V{ 1 } clone ] 2dip swap (pgcd) [ >array ] bi@ ;
  88
  89 : pdiff ( p -- p' )
  90     dup length v* { 0 } ?head drop ;
  91
  92 : polyval ( x p -- p[x] )
  93     [ length swap powers ] [ nip ] 2bi v. ;
  94
  95 MACRO: polyval* ( p -- )
  96     reverse
  97     [ 1 tail [ \ * swap \ + [ ] 3sequence ] map ]
  98     [ first \ drop swap [ ] 2sequence ] bi
  99     prefix \ cleave [ ] 2sequence ;
 100