C++ Home

When you are done with that, try to understand this primality test, based on the same algorithm. This is a truly obscure piece of code.

Code: Select all

struct S {
  typedef unsigned long long ull;
  static ull n;
  ull a,b,c,x,y;
  S(ull a=1, ull b=0, ull c=0, ull x=1, ull y=0):a(a),b(b),c(c),x(x),y(y){}
  S operator*(S const &s){
    return S((a*s.a+b*s.c+c*s.b)%n,
        (a*s.b+b*(s.a+s.c)+c*(s.b+s.c))%n,
        (a*s.c+b*s.b+c*(s.a+s.c))%n,
        (x*s.x+y*s.y)%n,
        (x*s.y+y*(s.x+s.y))%n);
  }
  S operator^(unsigned n){
    return n?(n&1?*this*(*this^(n-1)):(*this* *this)^(n>>1)):1;
  }
};

S::ull S::n;

bool is_prime(unsigned n){
  S::n=n;
  S s=S(0,1,0,0,1)^n;
  return (3*s.a+2*s.c)%n==0 && (2*s.x+s.y)%n==1;
}



Alvaro wrote:When you are done with that, try to understand this primality test, based on the same algorithm. This is a truly obscure piece of code.

For me the first one is giving me enough trouble. Just looking at that one hurts my eyes

Alvaro wrote:It uses a faster algorithm. Can anybody figure out how it works?

I got this far:

Code: Select all

    for(int i=2;i<20;++i)
        std::cout << i << ' '
        << fib(i) << std::endl;

fib(i) will call the function with the value of 2

Code: Select all

return power(FibMatrix(0,1),n).b;

The first thing here is to call FibMatrix constructor with 0 and 1, when that returns then the power function is called and passed a object of FibMatrix and the value of n which is 2.
Now in the power function we have:

Code: Select all

return n?(n&1?x*power(x,n-1):power(x*x,n>>1)):T(1);

since n is two that is true so we go to :

Code: Select all

power(x*x,n>>1)

Here we call the operator * for x*x :

Code: Select all

return FibMatrix(a*s.a+b*s.b,a*s.b+b*(s.a+s.b));

in here it breaks down like this:
I will replace the letters with their values
(0*0 + 1*1, 0*1 + 1*(0 + 1)
so after all the math is done we end up with a call to the constructor with this
FibMatrix(1,1)
now we return this to the calling power function and end up with this:

Code: Select all

power(FibMatrix(1,1),n>>1)

now the bit shift (n>>1) leaves us with 1
so now we call the power function with a FibMatrix(1,1) object and 1.
now starts recursion and my head is spinning so I will stop for now.

^-^ wrote:now starts recursion and my head is spinning so I will stop for now.

I'll give the solution before anybody's head explodes.

Consider the pairs of consecutive Fibonacci numbers:
(0,1)
(1,1)
(1,2)
(2,3)
(3,5)
(5,8)
...
If one of the pairs is (x,y) the next is (y,x+y). This operation is linear, which means that we can encode it with a 2x2 matrix.

Code: Select all

(0 1) (x) =  (y)
(1 1) (y)   (x+y)

The magic comes when we realize that we can apply the linear transformation n times to the first pair (0,1) and that is equivalent to multiplying the vector by the matrix (on the left) n times, which is the same as computing the n-th power of the matrix and then multiplying the result by the vector. Then we notice that all the matrices involved are of the following form:

Code: Select all

(a  b )
(b a+b)

and that the product of two matrices of that type is again a matrix of that type. So I just implemented a structure that contains a and b, with the correct operator* to compute the product of matrices of that special type.

The final step of multiplying by the column (0,1) is the same as picking b.

Any head explosions so far?

Well, now to the primality test...

No my head is intact. Although I will need some time for the brain cells to regroup. That is very cleaver if you ask me. With your explanation it seems very obvious now what was going on. As for the primality test maybe with this new information I might fair better. But that remains to be seen.

I am throwing in the towel. After looking at the primality test I can see you did something similar. But now I am so confused. What I understood from the code for fibonacci is gone. For me at the moment this kind of code/thinking/math is to much for me to understand.

Yeah, it's kind of complicated. I spent a lot of time dealing with sequences defined by recursion, so I am pretty familiar with these things, but I remember it being puzzling at the beginning.

For the primality test, what my code does is check if a number is simultaneously a Perrin pseudoprime and a Lucas pseudoprime. Google is your friend. Here is one link that is easy to understand:
http://www.ai.univie.ac.at/perrin.html

In time it will make more sense to me.
Yes Google is my friend. Thanks for the link I will take a look at it.