View Single Post
2021-09-09, 03:04   #9
bentonsar

"Sarah B"
Aug 2021
Washington, Seattle

22 Posts
Getting Closer

Quote:
 Originally Posted by ewmayer Assuming your inputs are properly normalized, i.e. x,y < (2^p-1), you break each into 32/33-bit pieces according to the per-word bitlengths in your list. (Balanced-digit representation is easy to achieve here, but not necessary for learning the basics.) Then multiply each of the D words of each so-discretized input by the corresponding fractional power of 2, fwdFFT each resulting weighted input vector, dyadic-multiply the resulting length-D forward transforms, invFFT the resulting vector, and multiply the resulting D discrete convolution outputs by the D inverse weights. The resulting outputs should be close enough to pure-integers that they can be confidently rounded-to-nearest and carries propagated. Check that the carry-propagated result matches x*y (mod 2^p-1) computed using exact multiprecision integer arithmetic.
Thx for that.

My current code is below. I think I am close. I can smell victory! I know I am missing the transfer back to the balanced int (to get the actual result of the operation), but I want to get the first iteration working (4x4 mod M) before tackling that.

Code:
import numpy as np

import math

q = 521
#Q = 163

D = 16

d = np.array([np.ceil((q * i) / D) - np.ceil(q*(i-1) / D) for i in range(0, D)]) #creates the bit length array
print(d)
a = np.array([2**(np.ceil(q * j/D) - (q * j/D)) for j in range(0, D)]) #creates the weight signal array
print(a)
ainv = np.array([(1.0/a[i]/D) for i in range(0, len(a))]) #creates the inverse weight signal array
print("Inverted A: " + str(ainv))

i_hiBitArr = [int(0)] * D
dsignal = [int(0)] * D

llint_signal = [int(0)] * D
llint_signal[0] = 4

for i in range(0, q-2):
ival = llint_signal[i]
print("Ival: " + str(ival))
if (i == 0):
ival += i_hiBitArr[D-1]
else:
ival += i_hiBitArr[i]
dsignal[i] = ival * a[i]
print(str(dsignal))
C = np.fft.fftn(dsignal, norm="forward")
print("First C: " + str(C))
C = C * C
print("Squared C: " + str(C))
DD = np.fft.irfft(C, norm="backward")
print("Inversed C value: " + str(DD))
if i == 0:
sig = DD[i] * ainv[i] - 2
else:
sig = DD[i] * ainv[i]
print(sig)
llint_signal[i] = np.around(sig)
dsignal[i] = sig
#using carry steps from Crandall Book
z = llint_signal
outputarr = [None] * len(z)
carry = 0
for n in range(0, len(z)-1):
B = 2**d[n+1]
v = z[n] + carry
outputarr[n] = math.remainder(v, B)
carry = np.floor(v//B)
print("Carry " + str(n) + ": " + str(carry))
print("V " + str(n) + ": " + str(v))
print("B " + str(n) + ": " + str(B))
print("outputarr " + str(n) + ": " + str(outputarr))
print(outputarr)
exit()

Is this implemented correctly? If so, how do I convert that outputarr to the actual result of the sxs mod M operation in the lucas lehmer test?

-Sarah

Last fiddled with by paulunderwood on 2021-09-09 at 06:07