More timings, not dumping the file just yet, but showing interface samples. I
made a direct pointer locale that still manages memory. Its much faster than
managing J objects, but will segfault if you think you are passing a pointer
when you are not. The J object version gives an innocent locale error. So
there is room for both.
Miller Rabin prime test in J
MRabin =: 4 : 0 NB. millerrabin test. x is random witnesses < y
e =. <. @ -:^:(0=2&|)^:a: y-1
for_a. x do. if. (+./c=y-1) +: 1={:c=. a y&|@^"0 e do. 0 return. end. end.
1
)
In oop BN mode,
MRabinBN=: 4 : 0 NB. millerrabin text. x is random witnesses < y
Y =. newM y
E =. newM e =. <. @ -:^:(0=2&|)^:a: y-1
for_a. 'I' inl newM x do.
BNmod_exp"1 EEi ;"0 1 a ;"1 EEi;"0 1 I__Y ; (CTX)[EEi =. 'I' inl EE =. dupM
E
'toX a:' inl EE
if. (+./ (y-1) = 'toX a:' inl EE) +: 1= BNget_word {: EEi do. out=.0 return.
end.
end. 1
)
In new raw BN mode (still managed memory)
MRabinBN=: 4 : 0 NB. millerrabin text. x is random witnesses < y
Y =. newM y
E =. newM e =. <. @ -:^:(0=2&|)^:a: y-1
for_a. newM x do.
BNmod_exp"1 EEi ;"0 1 a ;"1 EEi;"0 1 Y ; (CTX)[EEi =. dupM E
NB.pD toX EEi
if. (+./ (y-1) = toX EEi) +: 1= BNget_word {: EEi do. out=.0 return. end.
end. 1
)
testing a prime with 3 witnesses
C is OOP locale, D is direct locale
Pure J,
timespacex '123 129830012 12382902983209832098338901x MRabin__C
16513016020288569463710715785129617461891649734274611646334861564198564772449744004786622524842236224099696226154815067523x'
0.22579 3.07098e6
Direct BN: 100+x faster
timespacex '123 129830012 12382902983209832098338901x MRabinBN__D
16513016020288569463710715785129617461891649734274611646334861564198564772449744004786622524842236224099696226154815067523x'
0.00196416 18816
Object BN: 40x faster,
timespacex '123 129830012 12382902983209832098338901x MRabinBN__C
16513016020288569463710715785129617461891649734274611646334861564198564772449744004786622524842236224099696226154815067523x'
0.00597568 68096
J is faster at pure multiplication (of many smaller numbers) though. But its
easy, as in the miller rabin code above to mix and match J operations with BN
operations modular exponentiation is the main thing it does fast.
ProductW =: 3 : 0
acc =. newM {: y
(}: y) 4 : 'y[ BNfree a [ BNmul_word y;a =. newI x' reduce acc
)
timespacex 'ProductW__D p: i.300'
0.013215 79488
timespacex '*/ p: i.300x'
0.00380256 109312
There's actually not much gain to erasing things as you go along, instead of
all at once when you are done (if you have the memory). This is my new sum
code, XC displays the Jxint values and releases the BNs
sum =: (][ [: BNadd_ctx_ ] ; ;)/@:newM
XC__D sum__D i.10
45
product just needs to add a context (only variable in the manager)
product =: (][ [: BNmul_ctx_ ] ;wCTX ;)/@:newM
XC__D product__D p: i.10
6469693230
speed gains over J as numbers get larger (timings include memory free)
timespacex'( ] * 32 ^ ~ [ ) / p: i. 50x'
0.0053535982868485486 339840
timespacex'''XC (][[:BNmul ];wCTX ];[ [ [:BNexp [;wCTX a ;~[)/@:newM p:i.50 [ a
=: newM 32'' inl D'
0.0057142381714437853 63232
timespacex'''XC (][[:BNmul ];wCTX ];[ [ [:BNexp [;wCTX a ;~[)/@:newM p:i.200 [
a =: newM 47'' inl D'
0.0611437 182528
timespacex'( ] * 47 ^ ~ [ ) / p: i. 200x'
0.135313 259840
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
