#include <stdio.h>
#define LIMIT 1000000 //Numerically verify the conjecture up to this limit
//Integer square root code by Wilco Dijkstra
typedef unsigned uint32;
#define iter1(N) \
    try = root + (1 << (N)); \
    if (n >= try << (N))   \
    {   n -= try << (N);   \
        root |= 2 << (N); \
    }
uint32 isqrt (uint32 n)
{
    uint32 root = 0, try;
    iter1 (15);    iter1 (14);    iter1 (13);    iter1 (12);
    iter1 (11);    iter1 (10);    iter1 ( 9);    iter1 ( 8);
    iter1 ( 7);    iter1 ( 6);    iter1 ( 5);    iter1 ( 4);
    iter1 ( 3);    iter1 ( 2);    iter1 ( 1);    iter1 ( 0);
    return root >> 1;
}
//generates a "boolean" list of primes
char* sieve (int n)
{
  char* sieve;
  int i, j, m;
  m = isqrt(n);
  /* calloc initializes to zero */
  sieve = (char*)calloc (n+1, sizeof(char));
  sieve[0] = 1;
  sieve[1] = 1;
  for (i = 2; i <= m; i++)
    if (sieve[i] == 0)
      for (j = i*i; j <= n; j += i)
        if (sieve[j] == 0)
          sieve[j] = 1; 
  return sieve;
}
//calculates the nth lucas number mod m
int lucasmod(int n, int m){
	int x=2,y=1,temp;
	while(n){
		temp=(x+y)%m;
		x=y;
		y=temp;
		--n;
	}
	return x;
}
//calculates 2^e mod m
int mod2pow(int e, int m){
	unsigned long long int b=2;
	unsigned long long int r=1;
    while (e){
		if (e&1)
			r=(r*b)%m;
		e >>= 1;
		b=(b*b)%m;
	}
	return r;
}
//calculates (-1)^pow
inline char powm1(int pow){
	if(pow&1) return -1;
	else return 1;
}
int main(){
	char* composites=sieve(LIMIT);
	int p=3,p2,incrementTick=0;
	while(p<LIMIT){
		if(!composites[p]){
			p2=p<<1;
			int xpow2=1,x;
			for(x=1;xpow2<=p2;++x){
				xpow2=x*x;
				if((p2-xpow2)%5 == 0){
					int ypow2=(p2-xpow2)/5, y=isqrt(ypow2);
					if(ypow2==y*y){
						//we have found our integers: 2*p=x^2+5y^2
						int ind=(p+1)>>2, l=lucasmod(ind,p), val = mod2pow(ind,p);
						//y%8 is functionally equivalent to y&7
						//given p odd, (p-1)/2 is also functionally equivalent to p/2=p>>1
						//Also, note that ((-x) mod m) is congruent to -(x mod m)+m 
						if((y&7)==((p>>1)&7) || (y&7)==-((p>>1)&7)+8)
							val= val*powm1(ind);
						else
							val=-val*powm1(ind);
						if(val<0)
							val+=p;
						if(val!=l)
							printf("2*%d=%d^2+5*%d^2:%d,%d\n",p,x,y,l,val);
						break;
					}
				}
			}
		}
		//compute only for numbers = 3 or 7 mod 20 
		if(incrementTick=!incrementTick) p+=4;
		else p+=16;
	}
	return 0;
}