20210819, 16:23  #1 
"Sarah B"
Aug 2021
Washington, Seattle
2^{2} 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 
20210819, 17:37  #2 
Sep 2002
Database er0rr
2^{3}·5·97 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 20210819 at 18:59 
20210826, 21:52  #3 
"Sarah B"
Aug 2021
Washington, Seattle
4_{10} 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 
20210828, 10:33  #4 
Jun 2003
2^{5}×7×23 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 20210828 at 10:33 
20210828, 18:04  #5 
"Sarah B"
Aug 2021
Washington, Seattle
4_{16} 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 
20210907, 00:38  #6  
∂^{2}ω=0
Sep 2002
República de California
2·7^{3}·17 Posts 
Assuming your inputs are properly normalized, i.e. x,y < (2^p1), you break each into 32/33bit pieces according to the perword bitlengths in your list. (Balanceddigit representation is easy to achieve here, but not necessary for learning the basics.) Then multiply each of the D words of each sodiscretized input by the corresponding fractional power of 2, fwdFFT each resulting weighted input vector, dyadicmultiply the resulting lengthD 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 pureintegers that they can be confidently roundedtonearest and carries propagated. Check that the carrypropagated result matches x*y (mod 2^p1) computed using exact multiprecision integer arithmetic.
Quote:


20210907, 02:51  #7 
If I May
"Chris Halsall"
Sep 2002
Barbados
13^{2}×59 Posts 

20210907, 20:50  #8 
∂^{2}ω=0
Sep 2002
República de California
2×7^{3}×17 Posts 
BTW, a useful identity for computing inverse weights for modM(p) transform is as follows  for lengthN 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 primalityornot of p remains to be established ... I instead use M(p) = 2^p1 below):
w[j] = 2^[ceiling(j*p/N)  j*p/N], j = 0,...,N1 .[*] 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,...,N1, where sw denotes the number of "smallwords" in our variablebase 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[Nj] = 2 for j = 0,...,N . Thus w[j] = 2/w[Nj], and the jth inverse weight can be computed as 1/w[j] = w[Nj]/2 . The 1/N multiplier needed for inversetransform outputs can be lumped into the denominator of the RHS, thus one can define "effective inverse weights" by winv[j] = w[Nj]/(2*N); this needs just a single reciprocal 1/2N to be computed. Last fiddled with by ewmayer on 20210907 at 21:00 
20210909, 03:04  #9  
"Sarah B"
Aug 2021
Washington, Seattle
2^{2} Posts 
Getting Closer
Quote:
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*(i1) / 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, q2): ival = llint_signal[i] print("Ival: " + str(ival)) if (i == 0): ival += i_hiBitArr[D1] 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 20210909 at 06:07 

20210909, 05:27  #10 
∂^{2}ω=0
Sep 2002
República de California
2·7^{3}·17 Posts 
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 modM(p) multiply rather than a squaring), the resulting N mixedbase inputs and forwardweights, and the resulting outputs after inverseweighting but before carry propagation, and then after you've done the carries.

Thread Tools  
Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Is there utility in implementing an LLTlike sequence for GIMPS?  kriesel  Math  11  20190121 18:16 
Implementing Factoring Algorithms  Raman  Hobbies  45  20090511 05:11 
Implementing MPQS: SOS!  smoking81  Factoring  10  20071002 12:30 
Implementing algorithms, did I do this right?  ShiningArcanine  Programming  18  20051229 21:47 
improvement needed in SSE2 for Pentium 4 for nonIBDWT case for LLR/PRP  wfgarnett3  Hardware  3  20040629 03:03 