factor: trim using lists

[factor.git] / basis / math / polynomials / polynomials.factor

diff --git a/basis/math/polynomials/polynomials.factor b/basis/math/polynomials/polynomials.factor
index ec09b366a180054b86d73501fb0252a3a47f3142..225ea075377bca8eff23beb858041c5966ee05de 100644 (file)
--- a/basis/math/polynomials/polynomials.factor
+++ b/basis/math/polynomials/polynomials.factor
@@ -1,15 +1,15 @@
 ! Copyright (C) 2008 Doug Coleman.
 ! See http://factorcode.org/license.txt for BSD license.
-USING: arrays kernel make math math.order math.vectors sequences
-    splitting vectors ;
+USING: arrays combinators kernel make math math.bits
+math.vectors sequences vectors ;
 IN: math.polynomials
 
 <PRIVATE
 
 : 2pad-head ( p q n -- p q ) [ 0 pad-head ] curry bi@ ;
 : 2pad-tail ( p q n -- p q ) [ 0 pad-tail ] curry bi@ ;
-: pextend ( p q -- p q ) 2dup [ length ] bi@ max 2pad-tail ;
-: pextend-left ( p q -- p q ) 2dup [ length ] bi@ max 2pad-head ;
+: pextend ( p q -- p q ) 2dup max-length 2pad-tail ;
+: pextend-left ( p q -- p q ) 2dup max-length 2pad-head ;
 : unempty ( seq -- seq ) [ { 0 } ] when-empty ;
 : 2unempty ( seq seq -- seq seq ) [ unempty ] bi@ ;
 
@@ -20,23 +20,33 @@ PRIVATE>
 
 : p= ( p q -- ? ) pextend = ;
 
-: ptrim ( p -- p )
+: ptrim ( p -- q )
     dup length 1 = [ [ zero? ] trim-tail ] unless ;
 
-: 2ptrim ( p q -- p q ) [ ptrim ] bi@ ;
+: 2ptrim ( p q -- p' q' ) [ ptrim ] bi@ ;
 : p+ ( p q -- r ) pextend v+ ;
 : p- ( p q -- r ) pextend v- ;
-: n*p ( n p -- n*p ) n*v ;
+ALIAS: n*p n*v
 
-: pextend-conv ( p q -- p q )
-    2dup [ length ] bi@ + 1- 2pad-tail [ >vector ] bi@ ;
+: pextend-conv ( p q -- p' q' )
+    2dup [ length ] bi@ + 1 - 2pad-tail ;
 
 : p* ( p q -- r )
-    2unempty pextend-conv <reversed> dup length
-    [ over length pick <slice> pick [ * ] 2map sum ] map 2nip reverse ;
+    2unempty pextend-conv
+    [ drop length [ <iota> ] keep ]
+    [ nip <reversed> ]
+    [ drop ] 2tri
+    '[ _ _ <slice> _ vdot ] map reverse! ;
 
-: p-sq ( p -- p^2 )
-    dup p* ;
+: p-sq ( p -- p^2 ) dup p* ; inline
+
+ERROR: negative-power-polynomial p n ;
+
+: (p^) ( p n  -- p^n )
+    make-bits { 1 } [ [ over p* ] when [ p-sq ] dip ] reduce nip ;
+
+: p^ ( p n -- p^n )
+    dup 0 >= [ (p^) ] [ negative-power-polynomial ] if ;
 
 <PRIVATE
 
@@ -44,11 +54,9 @@ PRIVATE>
     2ptrim
     2dup [ length ] bi@ -
     dup 1 < [ drop 1 ] when
-    [ over length + 0 pad-head pextend ] keep 1+ ;
+    [ over length + 0 pad-head pextend ] keep 1 + ;
 
-: /-last ( seq seq -- a )
-    #! divide the last two numbers in the sequences
-    [ peek ] bi@ / ;
+: /-last ( seq1 seq2 -- x ) [ last ] bi@ / ;
 
 : (p/mod) ( p p -- p p )
     2dup /-last
@@ -65,7 +73,7 @@ PRIVATE>
 <PRIVATE
 
 : (pgcd) ( b a y x -- a d )
-    dup V{ 0 } clone p= [
+    dup V{ 0 } p= [
         drop nip
     ] [
         [ nip ] [ p/mod ] 2bi
@@ -78,8 +86,16 @@ PRIVATE>
     [ V{ 0 } clone V{ 1 } clone ] 2dip swap (pgcd) [ >array ] bi@ ;
 
 : pdiff ( p -- p' )
-    dup length v* { 0 } ?head drop ;
-
-: polyval ( p x -- p[x] )
-    [ dup length ] dip powers v. ;
-
+    dup length <iota> v* rest ;
+
+: polyval ( x p -- p[x] )
+    ! Horner scheme
+    [ nip <reversed> unclip-slice swap ]
+    [ drop ] 2bi
+    '[ [ _ * ] dip + ] each ;
+
+MACRO: polyval* ( p -- quot )
+    reverse
+    [ rest [ \ * swap \ + [ ] 3sequence ] map ]
+    [ first \ drop swap [ ] 2sequence ] bi
+    prefix \ cleave [ ] 2sequence ;