mersenneforum.org Implementing the IBDWT in Python?
 Register FAQ Search Today's Posts Mark Forums Read

 2021-08-19, 16:23 #1 bentonsar   "Sarah B" Aug 2021 Washington, Seattle 68 Posts Implementing the IBDWT in Python? Hi fellow crunchers, I am trying to implement the IBDWT in Python (for clarity not for speed) and I am having some trouble. I read through the relevant code in gwnum, gpuowl, gpulucas, and Prime Numbers: a computational perspective and I still can't seem to get it running correctly. I am using Numpy right now for the FFTs but I am fine changing that if need be. Has somebody already implemented one in python? Any help would be greatly appreciated. -Sarah
 2021-08-19, 17:37 #2 paulunderwood     Sep 2002 Database er0rr 102016 Posts See the posts here. It should be relatively easy to take the expressive Pari/GP code and convert into Python, The code there might be for Wagstaff numbers. I am sure someone here will point out the differences for it an that for Mersenne numbers. Last fiddled with by paulunderwood on 2021-08-19 at 18:59
 2021-08-26, 21:52 #3 bentonsar   "Sarah B" Aug 2021 Washington, Seattle 2·3 Posts Thanks Paul. That helped a ton. The part I am stuck on now is the digit to balanced and then back again (after the FFT) part. Is there any pseudocode anywhere for that? -Sarah
 2021-08-28, 10:33 #4 axn     Jun 2003 14F016 Posts This thread contains some description for balanced digit representation Also wikipedia Finally, have a look at cudaLucas source code for a GPU based implementation of LL test. The core FFT is based on Nvidia cuFFT library, but all the weightings, balancing, carry etc are done there. Ask if you have any specific questions/difficulties Last fiddled with by axn on 2021-08-28 at 10:33
 2021-08-28, 18:04 #5 bentonsar   "Sarah B" Aug 2021 Washington, Seattle 2·3 Posts Hi Axn, Thx for the links. Take the example from page 498 in Prime Numbers: A Computational Exploration where p = 2^521 - 1, q = 521, D=16, d = (33, 33, 32, 33, 32, 33, 32, 33, 33, 32, 33, 32, 33, 32, 33, 32), a = ... Right now I can successfully generate d and a based on those input values and my FFTs seem to be working correctly. My question is, let's say we want to have the x and y integers be 10 and 10 (or some other integer, it doesn't really matter). How do you apply those weights that were created to the integers so that they are ready for the FFT (of size D)? I read through the x = formula that is listed in 9.5.18 but I couldn't quite get it working correctly. -Sarah
2021-09-07, 00:38   #6
ewmayer
2ω=0

Sep 2002
República de California

1172710 Posts

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.

Quote:
 Originally Posted by bentonsar Hi Axn, Thx for the links. Take the example from page 498 in Prime Numbers: A Computational Exploration where p = 2^521 - 1, q = 521, D=16, d = (33, 33, 32, 33, 32, 33, 32, 33, 33, 32, 33, 32, 33, 32, 33, 32), a = ... Right now I can successfully generate d and a based on those input values and my FFTs seem to be working correctly. My question is, let's say we want to have the x and y integers be 10 and 10 (or some other integer, it doesn't really matter). How do you apply those weights that were created to the integers so that they are ready for the FFT (of size D)? I read through the x = formula that is listed in 9.5.18 but I couldn't quite get it working correctly. -Sarah

2021-09-07, 02:51   #7
chalsall
If I May

"Chris Halsall"
Sep 2002

243238 Posts

Quote:
 Originally Posted by ewmayer Assuming your inputs are properly normalized, i.e. x,y < (2^p-1).
What if that isn't a safe assumption?

 2021-09-07, 20:50 #8 ewmayer ∂2ω=0     Sep 2002 República de California 32×1,303 Posts BTW, a useful identity for computing inverse weights for mod-M(p) transform is as follows - for length-N transform, start with the weights as per the original Crandall/Fagin paper (I dislike their p = 2^q - 1 notation, since "p" in this sort of context generally implies prime, but their q is prime and primality-or-not of p remains to be established ... I instead use M(p) = 2^p-1 below): w[j] = 2^[ceiling(j*p/N) - j*p/N], j = 0,...,N-1 .[*] This gives w[0] = 1, followed by a sequence of nonrepeating fractional powers of 2. E.g. for p = 263 and N = 16 we have w[j] = 2^[0,9,2,11,4,13,6,15,8,1,10,3,12,5,14,7]/16 = 2^[j*sw (mod N)] for j = 0,...,N-1, where sw denotes the number of "smallwords" in our variable-base representation; if bw = p%N is the number of "bigwords", then sw = N - bw. If we extend the formula[*] to j = N, we get w[N] = w[0] = 1. Instead define w[N] := 2, then observe that w[j] * w[N-j] = 2 for j = 0,...,N . Thus w[j] = 2/w[N-j], and the jth inverse weight can be computed as 1/w[j] = w[N-j]/2 . The 1/N multiplier needed for inverse-transform outputs can be lumped into the denominator of the RHS, thus one can define "effective inverse weights" by winv[j] = w[N-j]/(2*N); this needs just a single reciprocal 1/2N to be computed. Last fiddled with by ewmayer on 2021-09-07 at 21:00
2021-09-09, 03:04   #9
bentonsar

"Sarah B"
Aug 2021
Washington, Seattle

68 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

2021-09-09, 05:27   #10
ewmayer
2ω=0

Sep 2002
República de California

32·1,303 Posts

Quote:
 Originally Posted by bentonsar 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
It'll be easier for the rest of us to tell if it's implemented correctly if you would provide a sample input, say a random number of the desired bitlength p you are working (or pair of such if you're doing a general mod-M(p) multiply rather than a squaring), the resulting N mixed-base inputs and forward-weights, and the resulting outputs after inverse-weighting but before carry propagation, and then after you've done the carries.

 Similar Threads Thread Thread Starter Forum Replies Last Post kriesel Math 11 2019-01-21 18:16 Raman Hobbies 45 2009-05-11 05:11 smoking81 Factoring 10 2007-10-02 12:30 ShiningArcanine Programming 18 2005-12-29 21:47 wfgarnett3 Hardware 3 2004-06-29 03:03

All times are UTC. The time now is 11:32.

Sun May 22 11:32:27 UTC 2022 up 38 days, 9:33, 0 users, load averages: 0.84, 1.07, 1.12