Special code for f/.@:g
Inspired by the faster implementation of the polynomial multiplication, I developed a conjunction which does the same as f/.@:g , bus is much more efficient.
obta=: 2 : 0 NB. oblique at table
assert. 'dyad only'
:
assert. >{.(x *.&(1=#@$) y) ; 'vectors only'
if. x -: y do. (u@:(v|.)\y) , }.u@:(v|.)\.y
else.
's t'=: x ,&# y
z=. $0
if. s=t do. y=.|.y
if. x-:y do. z=.(u@:(v@:(,:|.))\y) , }.u@:(v@:(,:|.))\.y
else. NB. y=.|.y
z=. i.0 1
for_j. i.&.<: s do. z=.z, ((j{.x) u@:v (-j){.y) , ((-j){.x) u@:v j{.y end.
end.
elseif. s<t do. y=.|.y
for_j. i.&.<:s do. z=.z, (j{.x) u@:v (-j){.y end.
z=.z, |.s x&(u@:v)\y
for_j. |.i.&.<: s do. z=.z, ((-j){.x) u@:v j{.y end.
elseif. s>t do. y=.|.y
for_j. i.&.<:t do. z=.z, (j{.x) u@:v (-j){.y end.
z=.z, t (u@:v)&y\x
for_j. |.i.&.<: t do. z=.z, ((-j){.x) u@:v j{.y end.
end.
end.
) So instead of a special polynomial multiplication, one can put +/ obta * instead of +//.@:(*/) which is the normal polynomial multiplication.
Some performance figures.
'X Y'=: 3000 4000 <@(?@$"0)1000 rank' X +//.@:(*/) Y';'X +/ obta * Y' +----+-------+----+------+ |rank|tm*sz |time|size | +----+-------+----+------+ | 1 |1361.16|2.41|565.26| +----+-------+----+------+ | 0 | 1.00|1.00| 1.00| +----+-------+----+------+ (X +//.@:(*/) Y)-:X +/ obta * Y 1
It seems X +/ obta * Y is even more efficient than the dedicated X pm Y.
rank' X pm Y';'X +/ obta * Y' +----+-----+----+----+ |rank|tm*sz|time|size| +----+-----+----+----+ | 1 |1.76 |1.24|1.42| +----+-----+----+----+ | 0 |1.00 |1.00|1.00| +----+-----+----+----+ (X pm Y)-:X +/ obta * Y 1
And it definitely is. However, no difference is measured for X -: Y.
rank' pm~ Y';' (+/ obta *)~ Y' +----+-----+----+----+ |rank|tm*sz|time|size| +----+-----+----+----+ | 0 |1.00 |1.00|1.00| +----+-----+----+----+ | 1 |1.02 |1.02|1.00| +----+-----+----+----+ (pm~ Y)-: (+/ obta *)~ Y 1
Sparse matrices
Analogous to sparse_multiplication we can have special code for sparce matrices too.
Here the variables are given in the form index part,: value part, or 4&(,@:$.),:5&(,@:$.).
spobta=: 2 : 0
assert. 'dyad only'
:
(] /:"1 {.) x(,@:(+/)&{. (~.@[ ,: u/.) ,@:(v/)&{:) y
)
'X Y'=:300 500 <@(>:@?@#~@[ ,:~ ?@# + (*i.)~)"(0) 10
$&.>' Xe Ye'=:expnd &.> X;Y
+----+----+
|2992|4993|
+----+----+
rank'X spm Y';'X +/ spobta * Y'
+----+-----+----+----+
|rank|tm*sz|time|size|
+----+-----+----+----+
| 0 |1.00 |1.00|1.00|
+----+-----+----+----+
| 1 |1.02 |1.02|1.00|
+----+-----+----+----+Hardly any difference to be measured, so no reason to have a special verb for polynomial multiplication if we have special code for f/.@:g, neither for normal, nor for sparse matrices.
Extended precision
Also with extended precision the special code for f/.@:g performs equally or better than the special polynomial multiplication.
'X Y'=: 900 1200 <@(?@$"0)1000x rank' X +//.@:(*/) Y';'X pm Y';'X +/ obta * Y' +----+------+----+------+ |rank|tm*sz |time|size | +----+------+----+------+ | 2 |444.48|1.00|522.92| +----+------+----+------+ | 1 | 1.55|1.14| 1.59| +----+------+----+------+ | 0 | 1.00|1.18| 1.00| +----+------+----+------+ 'X Y'=:300 500 <@(?@#~@[ ,:~ ?@# + (*i.)~)"(0) 30x $&.>' Xe Ye'=:expnd &.> X;Y +----+-----+ |8980|14974| +----+-----+ rank'X spm Y';'X +/ spobta * Y' +----+-----+----+----+ |rank|tm*sz|time|size| +----+-----+----+----+ | 1 |1.01 |1.01|1.00| +----+-----+----+----+ | 0 |1.00 |1.00|1.00| +----+-----+----+----+
