#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
#include <string.h>
#include <gmp.h>


#include "../../prmsieve/atkins/primegen.h"
#include "math52.h"

#if !defined(__alpha__)
//PFGW modular reduction
extern int32 _Imod2(uint32 * const a,const int32 n1,const int32 n2,const int32 len, int32 *p1,int32 *p2);
#endif


#if defined(__alpha__)
uint64 getrpc(uint64 from)
{
    register uint64 result, scratch;
    asm ( "rpcc %0" : "=r" (result) );
    scratch=result>>32;
    result+=scratch-from;
    return result&0xffffffff;
}
#endif


#define FMIN 30000
#define FMAX 100000
uint32 plus[(FMAX-FMIN)/32];
uint32 minus[(FMAX-FMIN)/32];

void addbackbit(uint32* p, int i)
{
    i-=FMIN;
    uint32 om=p[i>>5];
    uint32 b=1u<<(i&31);
    p[i>>5]=om&~b;
}
int removebit(uint32* p, int i, uint64 pr)
{
    i-=FMIN;
    uint32 om=p[i>>5];
    uint32 b=1u<<(i&31);
    if(om&b) 
	return 0;
    if(pr!=0)
    { 
	printf("%llu divides %i!%c1\n", pr, i+FMIN, (p==plus)?'+':'-'); 
	fflush(stdout); 
    }
    p[i>>5]=om|b;
    return 1;
}
void removeallbits()
{
    memset(plus, 0xff, sizeof(plus));
    memset(minus, 0xff, sizeof(minus));
}
int main(int argc, char**argv)
{
    primegen pg;
    mpz_t start;

    uint64 lastlog=0;
    int contents=(FMAX-FMIN)*2;

    // sensible defaults:
    char*filename="factsv.factors";
    uint64 startp=FMIN;
    int abcmode=0;

    if(argv[1])
    {
	filename=argv[1];
    }
    if(!access(filename, R_OK))
    {
	FILE* fp=fopen(filename, "r");
	if(fp)
	{
	    uint64 tmp64;
	    int tmp1, tmp2;
	    char buff[50];
	    while(fgets(buff, 50, fp))
	    {
		if(1==sscanf(buff, "ABC $a!+$b // %llu", &tmp64))
		{
		    abcmode=1;
		    startp=tmp64;
		    fprintf(stderr, "ABC mode\n");
		    contents=0;
		    removeallbits();
		}
		else if(3==sscanf(buff, "%llu divides %i!%i", &tmp64, &tmp1, &tmp2))
		{
		    if(startp<tmp64) startp=tmp64;
		    if(tmp2==-1) removebit(minus, tmp1, 0);
		    else if(tmp2==1) removebit(plus, tmp1, 0);
		    else { fprintf(stderr, "Bad: %s", buff); exit(1); }
		    --contents;
		}
		else if(abcmode && (2==sscanf(buff, "%i %i", &tmp1, &tmp2)))
		{
		    if(tmp2==-1) addbackbit(minus, tmp1);
		    else if(tmp2==1) addbackbit(plus, tmp1);
		    else { fprintf(stderr, "Bad: %s", buff); exit(1); }
		    ++contents;
		}
		else
		{
		    fprintf(stderr, "What's line : %s", buff);
		}
	    }
	}
	fclose(fp);
    }

    mpz_init(start);
    mpz_fac_ui(start, FMIN);

    primegen_init(&pg);
    primegen_skipto(&pg, startp-1);
    fprintf(stderr, "Starting a %llu, contents=%i\n", startp, contents);
    while(1)
    {
	int64 p1=primegen_next(&pg); 
	int64 p2=primegen_next(&pg); 

	if((lastlog^p2)>0xffff)
	{
	    fprintf(stderr, "p=%llu c=%i\n", p2, contents);
	    lastlog=p2;
	}

#if !(defined(__alpha__))
	if(p2<(1u<<31))
	{
	    int32 ip1=p1, ip2=p2;
	    int32 res1, res2;
	    _Imod2(start[0]._mp_d, ip1, ip2, start[0]._mp_size, &res1, &res2);
	    if(res1<0) res1+=ip1;
	    if(res2<0) res2+=ip2;
	    //printf("%i->%i, %i->%i\n", ip1, res1, ip2, res2);
	    int32 f=FMIN;
	    int32 pm1=ip1-1, pm2=ip2-1;
	    while(f<FMAX)
	    {
		if(res1==1)   { contents-=removebit(minus, f, p1); }
		else if(res1==pm1) { contents-=removebit(plus, f, p1); }
		if(res2==1)   { contents-=removebit(minus, f, p2); }
		else if(res2==pm2) { contents-=removebit(plus, f, p2); }
		++f;
		res1 = ((int64)res1*f)%ip1;
		res2 = ((int64)res2*f)%ip2;
	    }
	}	    
	else
	{
	    break;
	}
#else
	int64 res1=mpz_tdiv_ui(start, p1);
	int64 res2=mpz_tdiv_ui(start, p2);
	//fprintf(stderr, "%lli->%lli, %lli->%lli\n", p1, res1, p2, res2);
	int f=FMIN;
	
	//uint64 t1=getrpc(0);

#if 1
	double ff=(double)f;
	double fres1=res1;
	double fres2=res2;
	const double fone=1.0;
	const double fminusone=-fone;
	const double frounder=6755399441055744.0;
#else
	int64 pm1=p1-1;
	int64 pm2=p2-1;
#endif
	const double fp1=p1;
	const double fp2=p2;	
	const double fpr1=1.0/fp1;
	const double fpr2=1.0/fp2;
	
#define ALLFLOAT 1

	//printf("p=%llu, 30000!%%p=%llu\n", p, res);
	while(f<FMAX)
	{
#if 1
	    double ft1,ft2;
	    if(fres1==fone)   { contents-=removebit(minus, f, p1); }
	    else if(fres1==fminusone) { contents-=removebit(plus, f, p1); }
	    if(fres2==fone)   { contents-=removebit(minus, f, p2); }
	    else if(fres2==fminusone) { contents-=removebit(plus, f, p2); }
#else
	    if(res1==1)   { contents-=removebit(minus, f, p1); }
	    else if(res1==pm1) { contents-=removebit(plus, f, p1); }
	    if(res2==1)   { contents-=removebit(minus, f, p2); }
	    else if(res2==pm2) { contents-=removebit(plus, f, p2); }
#endif
	    ++f;
#if ALLFLOAT
	    ff+=fone;
	    fres1*=ff;      
	    fres2*=ff;
	    ft1=fres1*fpr1; 
	    ft2=fres2*fpr2;
	    ft1+=frounder;  
	    ft2+=frounder;
	    ft1-=frounder;  
	    ft2-=frounder;
	    ft1*=fp1;       
	    ft2*=fp2;
	    fres1-=ft1;     
	    fres2-=ft2;
#else
	    res1 = mulmod52_x_pr(res1, f, p1, fpr1);
	    res2 = mulmod52_x_pr(res2, f, p2, fpr2);
#endif
#if 0
	    if(((fres1>=0) && res1!=(int64)fres1) || ((fres1<0) && res1!=p1+(int64)fres1))
	    { fprintf(stderr, "fres1=%f res1=%llu\n", fres1, res1); }
	    if(((fres2>=0) && res2!=(int64)fres2) || ((fres2<0) && res2!=p2+(int64)fres2))
	    { fprintf(stderr, "fres2=%f res2=%llu\n", fres2, res2); }
#endif		
	}
	//t1=getrpc(t1);
	//fprintf(stderr, "t=%llu ticks\n", t1);
#endif
	
    }
    return 0;
}