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
  36     [ over length pick <slice> pick [ * ] 2map sum ] map 2nip reverse ;
  37
  38 : p-sq ( p -- p^2 )
  39     dup p* ;
  40
  41 <PRIVATE
  42
  43 : p/mod-setup ( p p -- p p n )
  44     2ptrim
  45     2dup [ length ] bi@ -
  46     dup 1 < [ drop 1 ] when
  47     [ over length + 0 pad-head pextend ] keep 1 + ;
  48
  49 : /-last ( seq seq -- a )
  50     #! divide the last two numbers in the sequences
  51     [ peek ] bi@ / ;
  52
  53 : (p/mod) ( p p -- p p )
  54     2dup /-last
  55     2dup , n*p swapd
  56     p- >vector
  57     dup pop* swap rest-slice ;
  58
  59 PRIVATE>
  60
  61 : p/mod ( p q -- z w )
  62     p/mod-setup [ [ (p/mod) ] times ] V{ } make
  63     reverse nip swap 2ptrim pextend ;
  64
  65 <PRIVATE
  66
  67 : (pgcd) ( b a y x -- a d )
  68     dup V{ 0 } clone p= [
  69         drop nip
  70     ] [
  71         [ nip ] [ p/mod ] 2bi
  72         [ pick p* swap [ swapd p- ] dip ] dip (pgcd)
  73     ] if ;
  74
  75 PRIVATE>
  76
  77 : pgcd ( p q -- a d )
  78     [ V{ 0 } clone V{ 1 } clone ] 2dip swap (pgcd) [ >array ] bi@ ;
  79
  80 : pdiff ( p -- p' )
  81     dup length v* { 0 } ?head drop ;
  82
  83 : polyval ( x p -- p[x] )
  84     [ length swap powers ] [ nip ] 2bi v. ;
  85
  86 MACRO: polyval* ( p -- )
  87     reverse
  88     [ 1 tail [ \ * swap \ + [ ] 3sequence ] map ]
  89     [ first \ drop swap [ ] 2sequence ] bi
  90     prefix \ cleave [ ] 2sequence ;
  91