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 ;
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@ ;
18 : powers ( n x -- seq )
19 <repetition> 1 [ * ] accumulate nip ;
21 : p= ( p q -- ? ) pextend = ;
24 dup length 1 = [ [ zero? ] trim-tail ] unless ;
26 : 2ptrim ( p q -- p' q' ) [ ptrim ] bi@ ;
27 : p+ ( p q -- r ) pextend v+ ;
28 : p- ( p q -- r ) pextend v- ;
31 : pextend-conv ( p q -- p' q' )
32 2dup [ length ] bi@ + 1 - 2pad-tail ;
36 [ drop length [ iota ] keep ]
39 '[ _ _ <slice> _ v* sum ] map reverse ;
41 : p-sq ( p -- p^2 ) dup p* ; inline
43 ERROR: negative-power-polynomial p n ;
46 make-bits { 1 } [ [ over p* ] when [ p-sq ] dip ] reduce nip ;
51 [ negative-power-polynomial ] if ;
55 : p/mod-setup ( p p -- p p n )
58 dup 1 < [ drop 1 ] when
59 [ over length + 0 pad-head pextend ] keep 1 + ;
61 : /-last ( seq1 seq2 -- x ) [ last ] bi@ / ;
63 : (p/mod) ( p p -- p p )
67 dup pop* swap rest-slice ;
71 : p/mod ( p q -- z w )
72 p/mod-setup [ [ (p/mod) ] times ] V{ } make
73 reverse nip swap 2ptrim pextend ;
77 : (pgcd) ( b a y x -- a d )
82 [ pick p* swap [ swapd p- ] dip ] dip (pgcd)
88 [ V{ 0 } clone V{ 1 } clone ] 2dip swap (pgcd) [ >array ] bi@ ;
91 dup length iota v* rest ;
93 : polyval ( x p -- p[x] )
95 [ nip <reversed> unclip-slice swap ]
97 '[ [ _ * ] dip + ] each ;
99 MACRO: polyval* ( p -- )
101 [ rest [ \ * swap \ + [ ] 3sequence ] map ]
102 [ first \ drop swap [ ] 2sequence ] bi
103 prefix \ cleave [ ] 2sequence ;