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 max-length 2pad-tail ;
12 : pextend-left ( p q -- p q ) 2dup max-length 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 ;
49 dup 0 >= [ (p^) ] [ negative-power-polynomial ] if ;
53 : p/mod-setup ( p p -- p p n )
56 dup 1 < [ drop 1 ] when
57 [ over length + 0 pad-head pextend ] keep 1 + ;
59 : /-last ( seq1 seq2 -- x ) [ last ] bi@ / ;
61 : (p/mod) ( p p -- p p )
65 dup pop* swap rest-slice ;
69 : p/mod ( p q -- z w )
70 p/mod-setup [ [ (p/mod) ] times ] V{ } make
71 reverse nip swap 2ptrim pextend ;
75 : (pgcd) ( b a y x -- a d )
80 [ pick p* swap [ swapd p- ] dip ] dip (pgcd)
86 [ V{ 0 } clone V{ 1 } clone ] 2dip swap (pgcd) [ >array ] bi@ ;
89 dup length iota v* rest ;
91 : polyval ( x p -- p[x] )
93 [ nip <reversed> unclip-slice swap ]
95 '[ [ _ * ] dip + ] each ;
97 MACRO: polyval* ( p -- )
99 [ rest [ \ * swap \ + [ ] 3sequence ] map ]
100 [ first \ drop swap [ ] 2sequence ] bi
101 prefix \ cleave [ ] 2sequence ;