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 math.bits ;
   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     make-bits { 1 } [ [ over p* ] when [ p-sq ] dip ] reduce nip ;
  45
  46 : p^ ( p n -- p^n )
  47     dup 0 >=
  48     [ (p^) ]
  49     [ 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 ( seq seq -- a )
  60     #! divide the last two numbers in the sequences
  61     [ 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 } clone 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 v* { 0 } ?head drop ;
  92
  93 : polyval ( x p -- p[x] )
  94     [ length swap powers ] [ nip ] 2bi v. ;
  95
  96 MACRO: polyval* ( p -- )
  97     reverse
  98     [ 1 tail [ \ * swap \ + [ ] 3sequence ] map ]
  99     [ first \ drop swap [ ] 2sequence ] bi
 100     prefix \ cleave [ ] 2sequence ;
 101