- **Programming**
(*https://www.mersenneforum.org/forumdisplay.php?f=29*)

- - **Implementing the IBDWT in Python?**
(*https://www.mersenneforum.org/showthread.php?t=27089*)

Implementing the IBDWT in Python?Hi fellow crunchers,
I am trying to implement the [URL="http://www.faginfamily.net/barry/Papers/Discrete%20Weighted%20Transforms.pdf"]IBDWT[/URL] 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 |

See the [URL="https://mersenneforum.org/showthread.php?p=531780#post531780"]posts here[/URL]. It should be relatively easy to take the expressive Pari/GP code and convert into Python, :smile:
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. |

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 |

[URL="https://mersenneforum.org/showthread.php?t=27012"]This[/URL] thread contains some description for balanced digit representation
Also [URL="https://en.wikipedia.org/wiki/Signed-digit_representation"]wikipedia[/URL] Finally, have a look at [URL="https://sourceforge.net/p/cudalucas/code/HEAD/tree/trunk/"]cudaLucas source code[/URL] 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 |

Hi Axn,
Thx for the links. Take the example from page 498 in [URL="http://thales.doa.fmph.uniba.sk/macaj/skola/teoriapoli/primes.pdf"]Prime Numbers: A Computational Exploration[/URL] 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 |

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=bentonsar;586739]Hi Axn, Thx for the links. Take the example from page 498 in [URL="http://thales.doa.fmph.uniba.sk/macaj/skola/teoriapoli/primes.pdf"]Prime Numbers: A Computational Exploration[/URL] 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[/QUOTE] |

[QUOTE=ewmayer;587417]Assuming your inputs are properly normalized, i.e. x,y < (2^p-1).[/QUOTE]
What if that isn't a safe assumption? |

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. |

Getting Closer[QUOTE=ewmayer;587417]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]
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()[/Code] 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 |

[QUOTE=bentonsar;587541]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[/QUOTE] 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. |

All times are UTC. The time now is 08:40. |

Powered by vBulletin® Version 3.8.11

Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.