Gene, Thanks a lot for all the information you provided on the LL test. With your clear explanations I was able to get my program mostly working. You explained it better than anything else I could find on the web!
Now I am having a problem with my IFFT implementation. My problem is with the x values getting negative values. After a couple iterations, the complex part of the fft output is negative and stays negative after the squaring. When the IFFT is performed, the output contains negative real parts for some elements. I have included the output of my program for M(37). Please take a look at iteration 4. I think this is a bug in my IFFT routine. I used Matlab to perform iteration 4, and the output of my IFFT is the same for the even elements, but different for the odd elements. I am using a recursive in-place radix-2 DIF for the forward transform and a DIT for the inverse transform. My algorithm is based on the one here: http://www.eptools.com/tn/T0001/PT04.HTM#Head283. In the DIT, I am multiplying the odd elements by [cos(2*PI*n/N) + sin(2*PI*n/N) i]. Since half of the elements are correct, I think the twiddle factor I'm using might be incorrect. Could anyone provide some insight on where I might be going wrong? Thanks, Chris Clark Output of my program below: a[0] = 1.000000 a[1] = 1.681793 a[2] = 1.414214 a[3] = 1.189207 b[0] = 10 b[1] = 9 b[2] = 9 b[3] = 9 sum[0] = 0 sum[1] = 10 sum[2] = 19 sum[3] = 28 convert: x[3] = 0.000000 x[2] = 0.000000 x[1] = 0.000000 x[0] = 4.000000 --------------- ITERATION 1 --------------- multiply input signal by weight signal y[0] = 4.000000 y[1] = 0.000000 y[2] = 0.000000 y[3] = 0.000000 fft squared y[0] 4.0000 + 0.0000i 16.0000 + 0.0000i y[1] 4.0000 + 0.0000i 16.0000 + 0.0000i y[2] 4.0000 + 0.0000i 16.0000 + 0.0000i y[3] 4.0000 + 0.0000i 16.0000 + 0.0000i ifft divide round y[0] 16.0000 + 0.0000i 16.0000 + 0.0000i 16.0000 + 0.0000i y[1] 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i y[2] 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i y[3] 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i propagate carries y[0] = 14.000000 carry = 0 y[1] = 0.000000 carry = 0 y[2] = 0.000000 carry = 0 y[3] = 0.000000 carry = 0 --------------- ITERATION 2 --------------- multiply input signal by weight signal y[0] = 14.000000 y[1] = 0.000000 y[2] = 0.000000 y[3] = 0.000000 fft squared y[0] 14.0000 + 0.0000i 196.0000 + 0.0000i y[1] 14.0000 + 0.0000i 196.0000 + 0.0000i y[2] 14.0000 + 0.0000i 196.0000 + 0.0000i y[3] 14.0000 + 0.0000i 196.0000 + 0.0000i ifft divide round y[0] 196.0000 + 0.0000i 196.0000 + 0.0000i 196.0000 + 0.0000i y[1] 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i y[2] 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i y[3] 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i propagate carries y[0] = 194.000000 carry = 0 y[1] = 0.000000 carry = 0 y[2] = 0.000000 carry = 0 y[3] = 0.000000 carry = 0 --------------- ITERATION 3 --------------- multiply input signal by weight signal y[0] = 194.000000 y[1] = 0.000000 y[2] = 0.000000 y[3] = 0.000000 fft squared y[0] 194.0000 + 0.0000i 37636.0000 + 0.0000i y[1] 194.0000 + 0.0000i 37636.0000 + 0.0000i y[2] 194.0000 + 0.0000i 37636.0000 + 0.0000i y[3] 194.0000 + 0.0000i 37636.0000 + 0.0000i ifft divide round y[0] 37636.0000 + 0.0000i 37636.0000 + 0.0000i 37636.0000 + 0.0000i y[1] 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i y[2] 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i y[3] 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i propagate carries y[0] = 770.000000 carry = 36 y[1] = 36.000000 carry = 0 y[2] = 0.000000 carry = 0 y[3] = 0.000000 carry = 0 --------------- ITERATION 4 --------------- multiply input signal by weight signal y[0] = 770.000000 y[1] = 60.544542 y[2] = 0.000000 y[3] = 0.000000 fft squared y[0] 830.5445 + 0.0000i 689804.2361 + 0.0000i y[1] 709.4555 + 0.0000i 503327.0470 + 0.0000i y[2] 770.0000 + -60.5445i 589234.3584 + -93238.5945i y[3] 770.0000 + 60.5445i 589234.3584 + 93238.5945i ifft divide round y[0] 592900.0000 + 0.0000i 592900.0000 + 0.0000i 592900.0000 + 0.0000i y[1] 139857.8918 + -0.0000i 83160.0000 + -0.0000i 83160.0000 + 0.0000i y[2] 3665.6416 + 0.0000i 2592.0000 + 0.0000i 2592.0000 + 0.0000i y[3] -46619.2973 + 0.0000i -39201.9999 + 0.0000i -39202.0000 + 0.0000i propagate carries y[0] = 2.000000 carry = 0 y[1] = 283.000000 carry = 0 y[2] = 195.000000 carry = 0 y[3] = -285.000000 carry = 0 This is what Matlab produces for the output of the IFFT: ifft(fft(y .* a).^2) 592900.0000 93238.5945 3665.6416 -0.0000 As you can see, my elements 1 and 3 are not correct. -----Original Message----- From: Ravvin, Gennadiy (Marlboro) [mailto:[EMAIL PROTECTED]] Sent: Sunday, April 21, 2002 1:35 PM To: 'Chris Clark ' Subject: RE: Mersenne: DWT help Hi Chris i recently had done the same thing you did and got a working version of the program (lots of man - hours!!!!). i won't go into details of that example but will try to explain what is the principle behind Lucas Lehmer test in all programs using DWT. Lucas Lehmer is very simple: let S(0) = 4. Let S(n) = S(n-1)^2 - 2. let q = 2^p-1 then if S(p-2) = 0 Mod 2^p-1 => q is prime. To implement this test we need p-2 squarings mod 2^P-1. for that we use DWT squaring because if we used ordinary FFT squaring we would need to pad the digits of array x with lots of zeroes and double run-length (but it would still work!!!). DWT squaring allows you to decrease FFT length by 2, plus get Mod 2^p-1 step almost for free!!! So your array x will contain the digits of squared numbers S(n). choose the dimension of x (i'll call it N)such that none of the elements of x would exceed 16 bits (otherwise loose precision) this is done by taking p and dividing by 16 then increase the result to the next power of 2 (since FFT is used on N which are powers of 2). now we do the squaring mod 2^p-1 to do the DWT squaring 0.Choose run length N so that none of the elements of b(i) 1.first multiply vector x by vector a (pointwise) 2.perform FFT on your new vector. the result is a complex vector x'. 3.do a pointwise-complex multiplication of vector x' with itself (i'll call the result y' - i.e. Re(y') = a^2 - b^2; Im(y') = 2*a*b where a = Re(x') and b = Im(x')) 4.do an inverse FFT on vector y' to obtain vector y. 5.divide all elements of y by vector a pointwise. your result is a vector which contains your digits of x^2 mod 2^p-1, however - in different bases - those bases are contained in vector b. so as a very last step make sure that elements of x(i) are smaller then elements in b(i) by doing a carry propagation step by in a loop: (this representation is what Dr. crandall and Dr. fagin call standard base) x(i) = x(i) mod b(i) carry = floor(x(i) / b(i)) x(i+1) = x(i+1) + carry and so on.... very important to note that if carry > 0 after the whole loop, then we need to add carry back to the element x[0] otherwise the results are going to be wrong and (well, at least i did it in my program and it works) possibly do a propagate carry step again in case x[0] + carry > b[0] now subtract 2 and repeat .... very sorry if i'm going into too much or too little detail for you - i simply don't know how familiar you are with this subject. will be glad to help again if needed Regards, \ Gene -----Original Message----- From: Chris Clark To: [EMAIL PROTECTED] Sent: 4/21/2002 4:00 AM Subject: Mersenne: DWT help I am trying to understand the DWT algorithm in the Crandall/Fagin paper so that I can implement it in C. I am aware of the many implementations available, and I have downloaded and looked at the source code for several. However, the optimizations performed by the code and the sparse comments make it difficult to follow the code without understanding the theory. I wasn't able to fine much by searching the list archives either. Can anyone point me to a reference with a detailed explanation of the DWT and LL algorithms (either on the web or in print), preferably with a worked example? Specifically, I am trying to reproduce the example from the paper for M(37). I understand how to compute the x, b, and a vectors. The next step is to perform the DWT. The paper says to calculate X = DWT(N,a)x = DFT(N)ax, but I don't completely understand the notation here. What is the input to the DFT? As I understand it, I should perform a standard FFT and then do a point-by-point multiply on the FFT result, the a signal, and the x signal. Is this correct? I appreciate any help you can provide. Thanks, Chris Clark ________________________________________________________________________ _ Unsubscribe & list info -- http://www.ndatech.com/mersenne/signup.htm Mersenne Prime FAQ -- http://www.tasam.com/~lrwiman/FAQ-mers _________________________________________________________________________ Unsubscribe & list info -- http://www.ndatech.com/mersenne/signup.htm Mersenne Prime FAQ -- http://www.tasam.com/~lrwiman/FAQ-mers
