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|
+----+-----+----+----+

RE Boss/J-blog/Special code for f/.@:g (last edited 2010-12-26 17:58:32 by RE Boss)