Further info. I did solve the problem (Euler Project no 592).
It needed the cumulative product of 2^23 odd numbers , 1 3 5
...33554431, modulo 2^48
This verb works, but is pretty slow on my laptop, taking 5 to 6 minutes:
cumoddprod =: 3 : 0
1 cumoddprod y
:
'lo hi' =. y
nc =. lo (>:@<.@-:@-~) hi
o =. nc#0 NB. simple (64-bit) integer output
cur=. x:x NB. extended current value, which might be initialised to
non-unity
in =. lo NB. simple integer input
for_i. i.nc do.
o =. (x:inv cur =. cur ((M&|@*) x:) in) i } o NB. multiply
extended, NB. save as integer
in =. 2 + in
end.
o
)
Small example:
4 8$64 | cumoddprod 1 63
1 3 15 41 49 27 31 17
33 51 47 57 17 11 63 33
1 35 15 9 49 59 31 49
33 19 47 25 17 43 63 1
It's prohibitive, as far as memory is concerned to store the whole
array as extended,
so I allow extended arithmetic on the two operands only.
I wouldn't need to use extended multiplication if multiply of <63-bit
integers modulo
M< ~ 2^63 was guaranteed correct! Too much to ask?
Any better ways to achieve the result of cumoddprod?
Thanks,
Mike
On 08/05/2017 15:08, 'Mike Day' via Programming wrote:
Yes, I did have a look, especially as one approach to the problem I
was examining was
to construct a set of coefficients modulo M based on polynomial
multiplication.
eg
(1 + ax)(1 + bx)(a+cx) .... = 1 + (a+b+c+...)x + (ab + bc + ca +
...)x^2 + ...
for which we should be able to use
add =: M&|@+
by =: M&|@*
ppr =: add//.@(by/) NB. adapted from the Wiki essay
with
ppr/ 1,. a,b,c,...
but it looked to me as if I'd have to invoke the fft multiplier in
scalar fashion.
I'm happy to be proved wrong.
But I'm probably on the wrong tack for the problem as a whole!
Thanks, Raul
Mike
On 08/05/2017 14:19, Raul Miller wrote:
On Mon, May 8, 2017 at 4:43 AM, 'Mike Day' via Programming
<[email protected]> wrote:
Is there a better way to do modulus multiply with largish arguments?
Hypothetically speaking, you might be able to do something based on
http://code.jsoftware.com/wiki/Essays/FFT
I hope this helps,
---
This email has been checked for viruses by Avast antivirus software.
https://www.avast.com/antivirus
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
