View Single Post
Old 2021-09-09, 03:04   #9
bentonsar
 
"Sarah B"
Aug 2021
Washington, Seattle

2×3 Posts
Default Getting Closer

Quote:
Originally Posted by ewmayer View Post
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
bentonsar is offline   Reply With Quote