On Sat, Mar 29, 2008 at 8:11 PM, Zsbán Ambrus <[EMAIL PROTECTED]> wrote:
> For multivariate polynomials, I've written a library in perl once that
> calculates with a sparse representation of them.
Actually, it's quite easy to do the basic operations in J too.
Let us represent a polynomial as a matrix with boxed coefficients in
the second columns and boxed list of variables in the first column.
(The following code might have bugs, I threw it together hastily.)
num2pol =: 3 :',:(,a:),(,<y)'"0 : [: NB. constant polynomial from number
]one =: num2pol 1
++-+
||1|
++-+
NB. identity in single variable from name of variable
name2pol =: 3 :',:(,<,<,>y),(,<1)' : [:
name2pol 'a'
+---+-+
|+-+|1|
||a|| |
|+-+| |
+---+-+
NB. normalize a polynomial
normpol1 =: ({:,~[:/:~&.:>{.)"1
normpol2 =: {."1 ({.@:{.,[:+/&.:>{:"1)/. ]
normpol3 =: #~ (0~:[:>{:"1)
normpol4 =: /:~
normpol =: normpol4@:normpol3@:normpol2@:normpol1 : [:
normpol a,(num2pol 1),a,(num2pol _1)
NB. scale a polynomial with a number
scalpol =: [: : ([: normpol3 ({.@:] , [*&.:>{:@:])"1)
2 scalpol a
+---+-+
|+-+|2|
||a|| |
|+-+| |
+---+-+
addpol =: [: : (normpol@:,) NB. add two polynomials
a addpol one
+---+-+
| |1|
+---+-+
|+-+|1|
||a|| |
|+-+| |
+---+-+
a addpol one addpol a addpol num2pol _1
+---+-+
|+-+|2|
||a|| |
|+-+| |
+---+-+
negpol =: _1 scalpol ] NB. negate
subpol =: [ addpol negpol@:] NB. subtract two polynomials
a subpol one
+---+--+
| |_1|
+---+--+
|+-+|1 |
||a|| |
|+-+| |
+---+--+
NB. multiply polynomials
mulpol1 =: ,&:(>@:{.) ,&:< *&:(>@:{:)
mulpol =: [: : ([: normpol [: ,/ mulpol1"1/)
]l =: (a addpol one) mulpol (a subpol one) NB. (a+1)*(a-1)
+-----+--+
| |_1|
+-----+--+
|+-+-+|1 |
||a|a|| |
|+-+-+| |
+-----+--+
]r =: (mulpol~ a) subpol one NB. a^2 - 1
+-----+--+
| |_1|
+-----+--+
|+-+-+|1 |
||a|a|| |
|+-+-+| |
+-----+--+
eqpol =: [: : (-:&:normpol) NB. polynomial equality
l eqpol r
1
You'd also need power, substitution, and a pretty-printer. I leave
those as homework for the reader.
Ambrus
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
