///VIRER
/*
 * parallel ?
 * (abondance des composites dans crb)
 */

///\VIRER

/*
 * Questions or comments about this code : rao at labri dot fr
 */

#include "utils.hpp"

#define Q_INF -1

#ifndef NODOUBLECHECK
#define DOUBLECHECK
#endif

int nexteven(int a)
{
  if(a<2) return 2;
  return a+(a&1);
}

int next1mod4(int a)
{
  if(a>(1+(a/4*4))) return 5+(a/4*4);
  return 1+(a/4*4);
}

class gcdex_t { int i; };

/***** precomputation of small primes ******/

int B1=10000000;
int B2=100001000;
int B3=1000000;

unsigned short *firstprimes;
int *primes;
int *rprimes;
int nbrprimes;

void initprimes() // precompute small primes
{
  printf("B1= %Ld  B2= %Ld\n",ll_t(B1),ll_t(B2));

  int tt=1+B2/8;
  unsigned char *tmp=(unsigned char*)malloc(tt);
  for(int i=0;i<tt;i++) tmp[i]=0xAA;
  tmp[0]=0xAE;
  for(long long i=3;i*i<=B2;i+=2)
    if((tmp[i/8]>>(i%8))&1) {
      for(long long j=i*i;j<=B2;j+=2*i) 
	tmp[j/8]&=~(1<<(j%8));
    }
  nbrprimes=0;
  for(long long i=2;i<=B2;i++)
    if((tmp[i/8]>>(i%8))&1)
      nbrprimes++;
  printf("%d primes <=B2\n",nbrprimes);
  primes=new int[nbrprimes];
  int j=1;primes[0]=2;
  for(int i=3;i<B2;i+=2) 
    if((tmp[i/8]>>(i%8))&1) {
      primes[j++]=i;
    }
  printf("init firstprimes OK\n");

  firstprimes=new unsigned short[B1];
  firstprimes[1]=0;
  for(int i=2;i<B1;i++) {
    if((tmp[i/8]>>(i%8))&1) {
      firstprimes[i]=0;
    } else {
      int j;
      for(j=0;primes[j]*primes[j]<=i;j++)
	if(i%primes[j]==0) {
	  firstprimes[i]=primes[j];
	  assert(primes[j]>1 && primes[j]<i);
	  break;
	}
      assert(primes[j]*primes[j]<=i);
    }
  }
  printf("init firstprimes2 OK\n");
  
  rprimes=new int[B3];
  for(int i=0;i<B3;i++) {
    rprimes[i]=-1;
  }
  int k=0;
  for(int i=2;i<B3;i++) {
    if((tmp[i/8]>>(i%8))&1) {
      rprimes[i]=k++;
    }
  }
  printf("init rprimes OK\n");
  
  free(tmp);
  assert(j==nbrprimes);
  printf("initprimes OK\n");
}

int nbrbaseisprime=5;

bool isprimei(int a)
{
  if(a<B1) return firstprimes[a]==0;
  return mpz_probab_prime_p(int_t(a).get_mpz_t(),nbrbaseisprime);
}

bool isprime(const int_t &a)
{
  if(a<B1) return firstprimes[toint(a)]==0;
  return mpz_probab_prime_p(a.get_mpz_t(),nbrbaseisprime);
}

///VIRER
#ifdef NOFACT
///\VIRER
#include "listfact.hpp"
///VIRER
#else
#include "listfact2.hpp"
#endif
///\VIRER

/****************/

int BE=1000;
int_t LIMIT=1;
int LIMITNBRP=20;
int LIMEXPNBRDP=40;
int maxcrb=-1;

///VIRER
bool nospecialinb=false;
int gcd_error=0;

#define MOSTWANTED
map<pair<int_t,pair<int_t,int> >, int > mostwanted;
///\VIRER

int_t assumedtrialcomp=int_t("100000000"); //option 'as'. (not checked with this program)

list<path_t> forceds; //forced paths ( 'fU' option )
list<path_t> forceds_fail;

int tabB[]={127,19,7,11,331,31,97,61,13,398581,1093,3,5,307,17,23,0};
int tabP[]={3,5,7,11,0}; //,13,17,0};
int tabP3[]={3,0};
int tabP19[]={3,5,7,11,13,17,19,0};
int tabC[]={0};
int *tab=tabB;

int tabstar[]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; //tab for "special" branching
int taballexact[]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; //tab for allexact branching

int *tab_2=NULL;
int sizetab=0;
int sizetab_2=0;

int tab2(int i) // tab2(i) return the ith prime number in the our modified order
{
  if(i<sizetab_2) return tab_2[i];
  if(i+1>=nbrprimes) return 0;
  return primes[i+1];
}

bool isbefore(const int_t &p, int j) //returns true if p is before tab2(i) in our order 
{
  if(j>sizetab_2) return p<primes[j+1];
  if(p>=primes[sizetab_2+1]) return false;
  for(int i=0;i<j;i++) if(p==tab2(i)) return true;
  return false;
}

void inittab2()
{
  sizetab_2=0;
  assert(primes[0]==2);
  for(int i=0;tab[i];i++) {
    sizetab=i+1;
    while(primes[sizetab_2]<=tab[i]) sizetab_2++;
  }

  sizetab_2--;
  if(sizetab_2<=0) return;
  tab_2=new int[sizetab_2];
  int k=0;
  for(k=0;tab[k];k++) tab_2[k]=tab[k];
  for(int i=1;i<sizetab_2;i++) {
    int j;
    for(j=0;tab[j];j++) if(tab[j]==primes[i]) break;
    if(tab[j]) continue;
    tab_2[k++]=primes[i];
  }
  assert(k==sizetab_2);

  printf("sizetab_2=%d\n",sizetab_2);
///VIRER

#if 0
  for(int i=0;i<sizetab_2+10;i++) {
    printf("%d ",tab2(i));
    //if(i==10 && sizetab_2-10>10) {i= sizetab_2-10;printf("... ");}
  }
  printf("...\n");
  for(int i=0;i<sizetab_2+10;i++) {
    printf("%d ",primes[i+1]);
    //if(i==10 && sizetab_2-10>10) {i= sizetab_2-10;printf("... ");}
  }
  printf("...\n");
#endif
///\VIRER
}

#define B_BIG 1
#define B_SMALL 2
#define B_SMALL1M3 4

#define B_NOCMP 16

#define C_PROD 64
#define C_NBRP 128
#define C_NBRDP 1024

#define C_COMP 256
#define C_COMPM2 2048

#define F_DEF (C_PROD+B_BIG)
#define F_NBRP (C_NBRP+B_SMALL)
#define F_COMP (C_COMP+B_BIG)
#define F_COMPM2 (C_COMPM2+B_BIG)

#define F_COMPF (C_COMP+B_BIG+B_NOCMP)
#define F_FFF (C_PROD+B_BIG+B_NOCMP)
#define F_NBRPF (C_NBRP+B_SMALL+B_NOCMP)

#define F_NBRDP (C_NBRDP+B_SMALL)

///VIRER
#define C_NBRPM2 512
#define F_NBRPM2 (C_NBRPM2+B_SMALL1M3)
///\VIRER

int flags=F_DEF;

struct exp_t {
  exp_t(int m=0) {
    minexp=m;
    exp=0;
    prime=true;
    maxexp=-1;
    nondiv=false;
  }
  
  bool prime;
  int minexp;
  int exp;
  int maxexp;
  bool nondiv; // if true, minexp+1 is not supposed to be a divisor of exp+1
  
  int_t p;
  int q;
  
  bool special() const {return minexp%2;}
  bool br() const {return minexp!=0;}
};

struct pfact_t{
  pfact_t() {}

  pfact_t(int_t a) {
    fs[a]=exp_t(0);
  }
  
  map<int_t,exp_t> fs;
  void aff(ostream &o) const {
    for(map<int_t,exp_t>::const_iterator it=fs.begin();it!=fs.end();it++) {
      if(!it->second.prime) o<<"c";
      o<<it->first<<"^"<<it->second.exp<<" ";
    }
    o<<endl;
  }

  int_t prod() const {
    int_t r=1;
    for(map<int_t,exp_t>::const_iterator it=fs.begin();it!=fs.end();it++)
      r*=pow(it->first,it->second.exp);
    return r;
  }

  void aff2(ostream &o) const {
    for(map<int_t,exp_t>::const_iterator it=fs.begin();it!=fs.end();it++) {
      if(!it->second.prime) o<<"c";
      o<<it->first<<"^("<<it->second.minexp;
      if(it->second.maxexp==it->second.minexp) o << "*";
      if(it->second.nondiv && it->second.maxexp==-1) o << "°°";
      o<<","<<it->second.exp<<") ";
      if(!it->second.prime) o<<"["<<it->second.p<<","<<it->second.q<<"] ";
    }
    o<<endl;
  }

  bool addprime(const int_t &a,int e=1) {
    if(fs[a].maxexp>=0 && (fs[a].exp+e)>fs[a].maxexp) return false;
    fs[a].exp+=e;
    return true;
  }

  void addcomp(const int_t &a,const int_t &p,int q) {
    if(fs[a].exp) {
      cerr << " FATAL double comp "<< a << endl;
      aff2(cerr);
      abort();
    }
    fs[a].exp+=1;
    fs[a].prime=false;
    fs[a].p=p;
    fs[a].q=q;
  }

  // return true if we have already supposed a special factor
  bool special() const {
    map<int_t,exp_t>::const_iterator it;
    for(it=fs.begin();it!=fs.end();it++) if(it->second.special()) return true;
    return false;
  }

  // check if gcd of every pair of factor is 1. just an assertion
  void assertgcd() const {
    map<int_t,exp_t>::const_iterator it,it2;
    for(it=fs.begin();it!=fs.end();it++)
      for(it2=fs.begin();it2!=it;++it2) {
	int_t g;
	mpz_gcd(g.get_mpz_t(),it->first.get_mpz_t(),it2->first.get_mpz_t());
	if(g!=1) {
	  cerr << "FATAL gcd = " << g << " : " << it->first << " " << it2->first << endl;
	  aff2(cerr);
	  {
	    ofstream gcdout("gcdfacts.txt",ios_base::out|ios_base::app);
#ifdef NOTH
	    if(envprof)
	      gcdout << "# " << envprof->indent() << endl;
#endif
	    if(!it->second.prime)
	      gcdout <<it->second.p<<" "<<it->second.q+1 << " " << g << endl;
	    if(!it2->second.prime)
	      gcdout <<it2->second.p<<" "<<it2->second.q+1 << " " << g << endl;
	  }
	  throw gcdex_t();
	  //abort();
	}
      }
  }

///VIRER
#ifdef NOFACT
///\VIRER

  // add a list of factors with exponents. return false if this is impossible (in "special" branchings)
  bool add(const pfact_t &o) {
    map<int_t,exp_t>::const_iterator it;
    for(it=o.fs.begin();it!=o.fs.end();it++) {
      if(fs[it->first].maxexp>=0 && fs[it->first].exp+it->second.exp>fs[it->first].maxexp) return false;
      fs[it->first].exp+=it->second.exp;
      if(!it->second.prime) {
	fs[it->first].prime=false;
	fs[it->first].p=it->second.p;
	fs[it->first].q=it->second.q;
      }
    }
#ifdef NOTH ///VIRER
    if(verb>=4) aff2(cout);
#endif ///\VIRER

    assertgcd();
    return true;
  }

///VIRER
#else 
  
  // add a list of factors with exponents. return false if this is impossible (in "special" branchings)
  bool add(const pfact_t &o) {
    bool cg=false;

    map<int_t,exp_t>::const_iterator it,it2;
    for(it=o.fs.begin();it!=o.fs.end();it++) {
      if(!cg) {
	if(!it->second.prime) {
	  for(it2=fs.begin();it2!=fs.end();it2++) {
	    if(it2->second.prime) continue;
	    int_t g;
	    mpz_gcd(g.get_mpz_t(),it->first.get_mpz_t(),it2->first.get_mpz_t());
	    if(g!=1) {
	      cg=true;break;
	    }
	  }
	} else {
	  for(it2=fs.begin();it2!=fs.end();it2++) {
	    int_t g;
	    mpz_gcd(g.get_mpz_t(),it->first.get_mpz_t(),it2->first.get_mpz_t());
	    if(g!=1) {
	      cg=true;break;
	    }
	  }
	}
      }

      if(fs[it->first].maxexp>=0) {
	if(fs[it->first].exp+it->second.exp>fs[it->first].maxexp) return false;
      }
      fs[it->first].exp+=it->second.exp;
      if(!it->second.prime) {
	fs[it->first].prime=false;
	fs[it->first].p=it->second.p;
	fs[it->first].q=it->second.q;
      }
    }
    if(verb>=4) aff2(cout);

    if(!cg) return true;
    
    int_t r;
    if(0!=(r=checkgcd())) {
      gcd_error++;
#ifdef DOUBLECHECK
      assert(isprime(r));
#endif
      if(verb>=2) aff2(cout);
      if(!split(r)) return false;
      if(verb>=2) aff2(cout);
    }
    return true;
  }

  bool split(const int_t &f) {
#ifdef NOTH
    if(verb>=2) indcout()<<"split by "<<f<<endl;
#endif
#ifdef DOUBLECHECK
    assert(f>1 && isprime(f));
#endif
    map<int_t,exp_t>::const_iterator it;
    for(it=fs.begin();it!=fs.end();++it) {
      if(it->second.prime) {continue;}
      if(isdiv(it->first,f)) {
	if(!addprime(f)) return false;
	int_t s=it->first/f;
	int_t c=it->first;
	int_t p=it->second.p;
	int q=it->second.q;
#ifdef DOUBLECHECK
	assert(it->second.exp==1);
#endif
	fs.erase(c);

	comperase(c);
	addfact(p,q,f);	
	if(isprime(s)){
	  if(!addprime(s)) return false;
	} else 
	  addcomp(s,p,q);
	return true;
      }
    }
    return true;
  }

  // check if gcd of every pair of factor is 1. just an assertion
  int_t checkgcd(bool fatal=false) const {
    map<int_t,exp_t>::const_iterator it,it2;
    for(it=fs.begin();it!=fs.end();it++)
      for(it2=fs.begin();it2!=it;++it2) {
	int_t g;
	mpz_gcd(g.get_mpz_t(),it->first.get_mpz_t(),it2->first.get_mpz_t());
	if(g!=1) {
	  if(verb>=2) {
	    cout << "gcd = " << g << " : " << it->first << " " << it2->first << endl;
	    aff2(cout);
	  }
	  if(fatal) {
	    cerr << "FATAL gcd != 1 " << endl;
	    abort();
	  }
	  return g;
	}
      }

    return 0;
  }

#endif
///\VIRER

  /******* CHECKS **********/

  int nbrdp() const {
    bool sp=special();
    int r=0,exp;
    int_t mi=0;
    
    map<int_t,exp_t>::const_iterator it;
    if(sp==false) {
      //do something ?
    }

    for(it=fs.begin();it!=fs.end();it++) {
      if(it->first==2) continue;
      if(it->second.prime) {
	r++;
      } else {
	if(isperfectpower(it->first)) {
	  cout << "FATAL: " << it->first << " is a perfect power (nbrdp)\n";
	  abort();
	}
	r+=2;
      }
    }
    return r;
  }

  bool checkNBRDP() const {
    return nbrdp()<LIMITNBRP;
  }

  int nbrp() const {
    bool sp=special();
    int r=0,exp;
    int_t mi=0;
    
    map<int_t,exp_t>::const_iterator it;
    if(sp==false) {
      for(it=fs.begin();it!=fs.end();it++) {
#ifdef DOUBLECHECK
	assert(it->second.minexp%2==0);
#endif
	if(it->second.prime==false || (it->first%4==1 && it->second.minexp==0)) { ////celui là peut etre impair
	  r--;break; 
	}
	//opt?: si composite, r-- et break aussi
      }
      //if(it==fs.end()) r++; // il faut un autre premier special 
    }

    for(it=fs.begin();it!=fs.end();it++) {
      if(it->first==2) continue;
      if(it->second.prime) {
	if((it->second.minexp%2)==0) {
	  exp=nexteven(it->second.exp);
	  if(exp<2) exp=2;
	  if(exp<it->second.minexp) exp=it->second.minexp;
	} else {
#ifdef DOUBLECHECK
	  assert((it->second.minexp%2)==1);
	  //assert(it->second.exp<2);
#endif
	  exp=next1mod4(it->second.exp);
	  if(exp<it->second.minexp) exp=it->second.minexp;
	}
      } else {
	if(isperfectpower(it->first)) {
	  //if it is a perfect power, we cannot assume exp>=4 (but this is very improbable).
	  cerr << "FATAL: " << it->first << " is a perfect power (nbrP)\n";
	  abort();
	}

	exp=4; 
      }
      
      r+=exp;
    }
    return r;
  }

  bool checkNBR() const {
    return nbrp()<LIMITNBRP;
  }

///VIRER

  /****** f(N) -2 F(n) ******/

  int nbrP() const {
    int r=0;
    
    map<int_t,exp_t>::const_iterator it;
    for(it=fs.begin();it!=fs.end();it++) {
      if(it->first==2) continue;
      if(it->second.prime) {
	r+=1;
      } else {
	if(isperfectpower(it->first)) {
	  //if it is a perfect power, we cannot assume P>=2 (but this is very improbable).
	  cerr << "FATAL: " << it->first << "is a perfect power\n";
	  abort();
	}
	r+=2; 
      }
    }
    return r;
  }

  bool checkNBRM2() const {
    return nbrp()-2*nbrP() <LIMITNBRP;
  }
///\VIRER

  /************ COMP ****************/

  bool checkCOMP() const {   // marche pas avec 'all'
    bool sp=special();

    map<int_t,exp_t>::const_iterator it;
    
    for(it=fs.begin();it!=fs.end();it++) { 
      if(it->second.prime && it->first>2) {
	int exp=0;
	if(it->second.exp>it->second.minexp) {
	  exp=it->second.exp;
	  if(sp || it->first%4!=1) exp=nexteven(exp);
	  else {
	    if(nexteven(exp)<next1mod4(exp)) exp=nexteven(exp);
	    else exp=next1mod4(exp);
	  }
	} else exp=it->second.minexp;
	//printf("e=%d m=%d (sp=%d) (mod=%d)-> %d\n",it->second.exp,it->second.minexp,sp,(it->first%4)==1?1:3,exp);
	assert(exp>0);
	if(pow(it->first,exp)>LIMIT) return false;
      }
    }
    return true;
  }

  bool checkCOMPM2() const {   // marche pas avec 'all'
    bool sp=special();

    map<int_t,exp_t>::const_iterator it;
    
    for(it=fs.begin();it!=fs.end();it++) { 
      if(it->second.prime && it->first>2) {
	int exp=0;
	if(it->second.exp>it->second.minexp) {
	  exp=it->second.exp;
	  if(sp || it->first%4!=1) exp=nexteven(exp);
	  else {
	    if(nexteven(exp)<next1mod4(exp)) exp=nexteven(exp);
	    else exp=next1mod4(exp);
	  }
	} else exp=it->second.minexp;
	//printf("e=%d m=%d (sp=%d) (mod=%d)-> %d\n",it->second.exp,it->second.minexp,sp,(it->first%4)==1?1:3,exp);
	assert(exp>0);
	if(exp>2 && pow(it->first,exp-2)>LIMIT) return false;
      }
    }
    return true;
  }


  /************** B ***************/

  int_t B() const { // returns a lower bound for N
///VIRER
//VERIF?
///\VIRER
    int_t r=1;
    int_t spp=1;
    int_t mdfsp=1,mdfnsp=1;

#ifdef NOTH ///VIRER
    if(verb>=4)
      cout << "B: ";
#endif ///\VIRER

    bool sp=special();
    map<int_t,exp_t>::const_iterator it;

    for(it=fs.begin();it!=fs.end();it++) {
      if(it->first==2) continue;
      int exp=0;
      int expsp=0;
      int expnsp=0;
      if(it->second.prime) {
	assert(it->second.nondiv==0);
	if(it->second.minexp==0) {
	  if(it->first%4==1) {
	    expnsp=nexteven(it->second.exp);
	    expsp=next1mod4(it->second.exp);
	  } else {
	  exp=nexteven(it->second.exp);
	  }
	} else if(it->second.minexp%2==0) {
	  exp=nexteven(it->second.exp);
	  if(exp<it->second.minexp) exp=it->second.minexp;
	} else { // (it->second.minexp%2)==1
	  assert(sp);
#ifdef DOUBLECHECK
	  assert(it->second.minexp%4==1);
	  assert(it->first%4==1);
#endif
	  exp=next1mod4(it->second.exp);
	  if(exp<it->second.minexp) exp=it->second.minexp;
	}

      } else {
	expsp=expnsp=1; // le composite peut avoir des carrés...
	if(!(it->second.exp==1 && it->second.minexp==0))
	  cerr << "FATAL comp " << it->first << " minexp=" <<it->second.minexp<< " exp=" << it->second.exp<<endl;
#ifdef DOUBLECHECK
	assert(it->second.exp==1 && it->second.minexp==0);
#endif
	//if(exp<it->second.minexp) exp=it->second.minexp;
      }

      if(exp) {
#ifdef NOTH ///VIRER
	if(verb>=4)
	  cout << it->first << "^"<<exp<< " ";
#endif ///\VIRER
///VIRER
	if(nospecialinb && exp%2) {
	  assert(sp);
	} else {
///\VIRER
	  r*=pow(it->first,exp);
///VIRER
	}
///\VIRER
      } else {
	int_t tcsp=pow(it->first,expsp);
	int_t tcnsp=pow(it->first,expnsp);

#ifdef NOTH ///VIRER
	if(verb>=4)
	  cout << it->first << "^"<<expnsp << ","<<expsp;
#endif ///\VIRER
	
	if(tcnsp>spp) {
#ifdef NOTH ///VIRER
	  if(nospecialinb && verb>=4) cout<<"'";
#endif ///\VIRER
	  spp=tcnsp;
	}

	if(tcsp*mdfnsp<tcnsp*mdfsp) {
#ifdef NOTH ///VIRER
	  if(!nospecialinb && verb>=4) cout<<"*";
#endif ///\VIRER
	  mdfsp=tcsp;
	  mdfnsp=tcnsp;
	}

	r*=tcnsp;
#ifdef NOTH ///VIRER
	if(verb>=4) cout<<" ";
#endif ///\VIRER
      }
    }

#ifdef NOTH  ///VIRER
    if(verb>=4) cout<<"\n";
#endif  ///\VIRER
    
///VIRER
    if(!nospecialinb) {
///\VIRER
      if(sp) {
#ifdef NOTH  ///VIRER
	if(verb>=4)
	  cout << "ret= " << r <<endl;
#endif  ///\VIRER
	return r;
      }
#ifdef NOTH  ///VIRER
      if(verb>=4)
	cout << "ret(/)= " << r*mdfsp/mdfnsp <<endl;
#endif  ///\VIRER
      assert(mdfnsp>=mdfsp);
      return r*mdfsp/mdfnsp;
///VIRER
    }
///\VIRER

///VIRER
    if(nospecialinb) {
      if(sp) {
#ifdef NOTH
	if(verb>=4) cout << "ret= " << r <<endl;
#endif
	return r;
      }
#ifdef NOTH
      if(verb>=4) cout << "ret(/)= " << r/spp <<endl;
#endif
      return r/spp;
    }
///\VIRER

  }

  bool checkB() const {
    return (B()<LIMIT);
  }

  /****** AB ******/
  
  mpq_class  compAb() const { 
///VIRER
//VERIF ? OPT ?
///\VIRER

    int_t ma=0;
    bool sp=special();
    map<int_t,exp_t>::const_iterator it;

    mpq_class ab=1;
    
    for(it=fs.begin();it!=fs.end();it++) {
      if(it->first<=2 || !it->second.prime) continue;
      
      int exp=0;
      if(it->second.minexp) {
	if((it->second.minexp%2)==0) exp=nexteven(it->second.exp);
	else exp=next1mod4(it->second.exp);
      } else exp=nexteven(it->second.exp);
      if(exp<it->second.minexp) exp=it->second.minexp;
      
      ab*=sigma(it->first,exp);
      ab/=pow(it->first,exp);
      
      if(sp && ab>2) return ab;
    }
    if(sp) return ab;

    mpq_class abw=ab;
    for(it=fs.begin();it!=fs.end();it++) { 
      if(it->first<=2 || false==it->second.prime || it->second.minexp!=0 || 1!=(it->first%4)) continue;
      
      int oexp=nexteven(it->second.exp);
      int nexp=next1mod4(it->second.exp);
      mpq_class ab2=ab;
      
      ab2/=sigma(it->first,oexp);
      ab2*=pow(it->first,oexp);
      
      ab2*=sigma(it->first,nexp);
      ab2/=pow(it->first,nexp);
      
      if(ab2<abw) abw=ab2;
    }

    return abw;
  }

  mpq_class Ab(bool worst=false) const {
    mpq_class r;
    mpq_class rc;
    if(!worst)
      return compAb();

    map<int_t,exp_t>::const_iterator it;

    r=1;rc=1;
    for(it=fs.begin();it!=fs.end();it++) {
      if(it->first>2) {
	if(it->second.prime) {
	  mpq_class a;
	  if(it->second.maxexp==-1) a=1+mpq_class(1)/(it->first-1);
	  else a=mpq_class(sigma(it->first,it->second.maxexp))/pow(it->first,it->second.maxexp);
	  r*=a;
	} else {
///VIRER
// VERIFIER !!!!!
///\VIRER
	  int_t l=1;
	  int k=0;
	  while(l<=it->first) {k++;l*=assumedtrialcomp;}
	  mpq_class a=1+mpq_class(1)/assumedtrialcomp;
	  for(;k>0;k--) rc*=a;
	}
      }
    }
    
    return r*rc;
  }

  bool checkAb() const {
    return compAb()<=2;
  }


  /****** BESTS *********/

  void branchables(set<int_t> &s) const {
    s.clear();
    map<int_t,exp_t>::const_iterator it;
    for(it=fs.begin();it!=fs.end();it++) if(it->first!=2 && (!it->second.br()) && it->second.prime) s.insert(it->first);
  }

  void bests(list<int_t> &s) const {
    set<int_t> s2;
    branchables(s2);
    s.clear();
    for(set<int_t>::iterator it=s2.begin();it!=s2.end();++it) s.push_front(*it);
  }

  void bestssmalls(list<int_t> &s) const {
    set<int_t> s2;
    branchables(s2);
    s.clear();
    for(set<int_t>::iterator it=s2.begin();it!=s2.end();++it) s.push_back(*it);
  }

  void bestssmalls1m3(list<int_t> &s) const {
    set<int_t> s2;
    branchables(s2);
    s.clear();
    for(set<int_t>::iterator it=s2.begin();it!=s2.end();++it) if((*it%3)==1) s.push_back(*it);
    for(set<int_t>::iterator it=s2.begin();it!=s2.end();++it) if((*it%3)!=1) s.push_back(*it);
  }
};

///VIRER
string outrfactname="";
///\VIRER
string outrbcompname="";
string prefixname="";

///VIRER
#ifndef NOFACT
#include "fact4.hpp"
#else
///\VIRER


void initoracle(int &ac, char **&av)
{
  initprimes();
  initlistfactors();

  for(int i=1;i<ac;i++) {
    fprintf(stderr,"read %s\n",av[i]);
    readfactors(av[i]);
///VIRER
    facts_t::facts_stats();
///\VIRER
  }
///VIRER

  //fprintf(stderr,"read facts.txt\n");
  //readfactors("facts.txt");
///\VIRER
}

///VIRER
bool optoracle(int &ac, char **&av) {
  if(!strcmp(av[1]+1,"refact")) {
    printf("re facts to %s\n",av[2]);
    outrfactname=av[2];
    //outrfact=new ofstream(av[2],ios_base::out|ios_base::app);
    ac-=2;av+=2;return true;
  }
  return false;
}

#endif
///\VIRER


int limsize=0;
int profmax=0;

int upto=0;

/// thread encapsulation
struct th_t {

  const path_t *envprof;

  ostream &indcout() {
    if(!envprof) return tout;
    return (tout << envprof->indent());
  }

  int id; //thread id;

  // "tout" is the ostream where we write the tree.
#ifndef NOTH
  ofstream tout;
#else 
  ostream &tout;
#endif

  static string id2str(int id) {
#ifndef NOTH
    char b[10];
    snprintf(b,10,"%s%d_",prefixname.c_str(),id);
    return b;
#else
    return prefixname;
#endif
  }
  
  th_t(int i=0):id(i),
#ifndef NOTH
		tout((id2str(i)+"log.txt").c_str(),ios_base::out|ios_base::app) 
#else
		tout(cout)
#endif
  {
    stat_rec=0;
    envprof=NULL;
    outrbcomp=NULL;
    if(outrbcompname!="") {
      outrbcomp=new ofstream((id2str(i)+outrbcompname).c_str(),ios_base::out|ios_base::app);
      *outrbcomp<<"#"<<cmdline<<endl;
    }
///VIRER
    outrfact=NULL;
    if(outrfactname!="") {
      outrfact=new ofstream((id2str(i)+outrfactname).c_str(),ios_base::out|ios_base::app);
      *outrfact<<"#"<<cmdline<<endl;
    }
///\VIRER
  }

  ~th_t() {
    tout<< "thread end. stat_rec=" << stat_rec <<endl;
///VIRER
    if(outrfact) delete outrfact;
///\VIRER
    if(outrbcomp) delete outrbcomp;
  }

/* check for contracitions: abundancy and limit */
bool checkcontradiction(const pfact_t &fsp)
{
  if((flags & C_PROD) &&!fsp.checkB()) {
    if(verb>=2) indcout()<< "N>LIMIT :\n";
    return true;
  }
  
  if((flags & C_NBRP) &&!fsp.checkNBR()) {
    if(verb>=2) indcout()<< "NBRP>LIMIT :\n";
    return true;
  }

  if((flags & C_NBRDP) &&!fsp.checkNBRDP()) {
    if(verb>=2) indcout()<< "NBRDP'>LIMIT :\n";
    return true;
  }
  
///VIRER
  if((flags & C_NBRPM2) &&!fsp.checkNBRM2()) {
    if(verb>=2) indcout()<< "NBRPM2>LIMIT :\n";
    return true;
  }
///\VIRER
  
  if((flags & C_COMP) &&!fsp.checkCOMP()) {
    if(verb>=2) indcout()<< "COMP>LIMIT\n";
    return true;
  }

  if((flags & C_COMPM2) &&!fsp.checkCOMPM2()) {
    if(verb>=2) indcout()<< "COMPM2>LIMIT\n";
    return true;
  }
  
  if(!fsp.checkAb()) {
    if(verb>=2) indcout()<< "Ab>2 :\n";
    return true;
  }
  
  return false;
}

/* check for contracitions: previous + special banching and forb factors */
bool checkcontradiction2(pfact_t &fsp, const pfact_t &fs, const pfact_t &fa,ll_t borne)
{
  fsp=fs;
  
  if(!fsp.add(fa)) {
    if(verb>=2) indcout()<< "contradiction add :\n";
    return true;
  }

  if(checkcontradiction(fsp)) return true;
  
  bool skip=false;
  map<int_t,exp_t>::const_iterator it;
  for(it=fa.fs.begin();it!=fa.fs.end();it++) {
    if(it->second.prime && it->first>2 && isbefore(it->first,borne)) {
      if(fs.fs.find(it->first)!=fs.fs.end()) continue; // not a new prime factor
      if(verb>=2) indcout()<< it->first << " already eliminated:\n";
      skip=true;
      break;
    }      
  }

  if(skip) return true;

  return false;
}


path_t forced; //forced path ( 'U' option )
ofstream *outrbcomp;
///VIRER
ofstream *outrfact;
///\VIRER

ull_t stat_rec;

int th_id;


///VIRER
#ifdef NOFACT
///\VIRER
// 'find' factors of f (factor if s(p^q)), and put the known decomposition into pfact
int oraclefact(const int_t &p,int q,pfact_t &pfact,int_t &f) 
{
  path_t prof;
  if(envprof) prof=*envprof;
  
  facts_t *fa=getfacts(p,q);
  
  bool ok=0;

  facts_t::const_iterator it;
  for(it=fa->begin();it!=fa->end();++it)
    while(isdiv(f,*it)) {
///VIRER
      if(outrfact) *outrfact<<p<<" " <<q+1<<" "<<*it<<endl;
///\VIRER
      ok=1;
      pfact.addprime(*it);
      f/=*it;
    }

  if(f==1) return 1;
  
  if(isprime(f)) {
    pfact.addprime(f);
    return 1;
  }
  
  pfact.addcomp(f,p,q);
  
  if(!ok) return 0;
  return 2;
}
///VIRER
#endif
///\VIRER

#if 0
bool rec(pfact_t &fs,int_t p,const path_t &prof,ll_t borne);
#endif

bool rec_range(pfact_t &fs, int cond,const path_t &prof,ll_t borne)
{
  const path_t *oprof=envprof;
  envprof=&prof;

  if(prof.h()<limsize) {
    prof.aff(tout); 
    if(prof.h()==limsize-1) tout << " BORDERPATH";
    tout <<endl;
  }
  ll_t ostat=stat_rec;

  if(maxcrb>=0 && prof.crb()>=maxcrb) {
    if(verb>=1) indcout() << "REC_RANGE FAIL MAXCRB\n";
    envprof=oprof;
    return false;
  }

  mpq_class abw=fs.Ab(true);

  if(abw>=2) {
    if(verb>=1) indcout() << "REC_RANGE FAIL abw>=2\n";
    envprof=oprof;
    return false;
  }

  int rlnbr=0;
  if(cond &C_NBRP) { 
    rlnbr=LIMITNBRP-fs.nbrp();
    rlnbr=(rlnbr+1)/2;
  }

  if(cond &C_NBRDP) rlnbr=LIMITNBRP-fs.nbrdp();

  ll_t pom2;
  if(cond&C_NBRPM2) { 
    ll_t t=LIMITNBRP;
    ll_t F=fs.nbrP();
    /* 06/10/2012
      > stop when (p/(p-1))^((7*t+31)/4-F) < 2/A
      Dans le papier, c'est un peu plus optimisé, on a:
      stop when (p/(p-1))^((7*t+14)/4-F) < 2/A
      Tu peux modifier le programme et mettre la borne du papier ?
    */
    //pom2=((7*t+31)/4-F);
    pom2=((7*t+17)/4-F);
    if(verb>=2)  { 
      mpq_class b=abw;
      indcout() << " t=" << t << " F="<<F<< " power= " << pom2 << " abw="<<b.get_d()<<endl;
    }
  }
  
  int_t dB=0;
  if(cond&C_PROD) {
    dB=LIMIT/fs.B();
    if(verb>=2) indcout() << "dB= "<< dB <<endl;
  }

  ll_t pmax;

  if(cond) {
    int i;
    for(i=0;i<nbrprimes;i++) {
      ull_t p=primes[i];
      int_t pt=ll2mpz(p);
      if(cond&C_NBRP) {
	mpq_class a=1+mpq_class(1)/pt;
	mpq_class b=abw;
	for(int i=0;i<rlnbr;i++) b*=a;
	if(b<2) break;
      } else 
      if(cond&C_PROD) {
	int q;
	int_t r=1;
	for(q=1;;q++) {
	  r*=pt;
	  if(r>dB) break;
	}
	q=q/2;
	mpq_class b=abw;
	mpq_class a=1+mpq_class(1)/(pt-1);
	for(int i=0;i<q;i++) {
	  b*=a;
	  if(b>=2) break;
	}
	if(b<2) {
	  if(verb>=1)  indcout() << "RANGE stop p=" << pt << " q=" <<q << " ab= "<< b.get_d() << endl;
	  break;
	}
      } else 
      if(cond&C_NBRPM2) {
	mpq_class a=1+mpq_class(1)/(pt-1);
	mpq_class b=abw;
	for(int i=0;i<pom2;i++)
	  b*=a;
	if(b<2) {
	  if(verb>=1) indcout() << "RANGE stop p=" << pt << endl;
	  break;
	}
      } else 
      if(cond&C_NBRDP) {
	mpq_class a=1+mpq_class(1)/(pt-1);
	mpq_class b=abw;
	for(int i=0;i<rlnbr;i++)
	  b*=a;
	if(b<2) {
	  if(verb>=1) indcout() << "RANGE stop p=" << pt << endl;
	  break;
	}
      } else abort();
    }
    pmax=primes[i]-1;
  }

  ull_t str=stat_rec;

  if(verb>=1) indcout() << "REC_RANGE N="<< prof.crb()<< " npmax="<< pmax<< " borne="<< borne << endl;

  bool rok=true;

  ull_t p;
  ull_t pm=0;
  for(int i=0;;i++) {
    
    p=tab2(i);
    if(p>pm) pm=p;
    if(p==0) break;
    assert(p>2);

    int_t pt=ll2mpz(p);

    if(p>pmax && i>=sizetab) {
      if(verb>=4) indcout() << p << " > pmax and not in S\n";
      if(i>sizetab_2) break;
      continue;
    }

    if(i<borne) {
      if(verb>=3) indcout()<< p << " already eliminated\n";
      continue;
    }
    
    map<int_t,exp_t>::const_iterator it=fs.fs.find(pt);

    if(it!=fs.fs.end() && it->second.minexp) {
      if(verb>=2) indcout() << "already "<< p<<endl;
      continue;
    }

    {
      //tout << "RANGE: " << p <<endl;
      assert(i>=borne);
      pfact_t fs2=fs;
      path_t prof2=prof;
      prof2.add(-1);
      bool r=rec(fs2,pt,prof2,i);
      if(!r) rok=false;
    }
  }

  if(p==0) {
    if(verb>=1) indcout() << "REC_RANGE FAIL p limit (max = "<< pm<<" )"  << endl;
    rok=false;
  }

  if(verb>=1) indcout() << "\\REC_RANGE " << rok << " N="<< prof.crb() << " size=" << stat_rec-str<<endl;

  envprof=oprof;

  if(prof.h()<limsize) {
    prof.aff(tout); 
    tout << " [ size= " << stat_rec-ostat<< " ret= " << rok <<" ]\n";
  }

  return rok;
}


bool rec(pfact_t &fs,int_t p,const path_t &prof,ll_t borne)
{
  if(profmax && prof.h()>=profmax) return false;

  int forcedq=-1;
  if(forced.h()>prof.h()) {
    if(!forced.compat(prof,p)) {
      //indcout()<< "wrong forced path !!!\n";
      return true;
    }

    list<pair<int_t,int> >::const_iterator it1=forced.p.begin(),it2=prof.p.begin();
    while(true) {
      if(it1==forced.p.end() || it2==prof.p.end()) break;
      it1++;it2++;
    }
    assert(it1!=forced.p.end());
    forcedq=it1->second;
  }

  if(prof.h()<limsize) {
    prof.aff(tout); 
    if(prof.h()==limsize-1) tout << " BORDERPATH";
    tout <<endl;
  }
  ll_t ostat=stat_rec;

  stat_rec++;

  const path_t *oprof=envprof;
  envprof=&prof;

  int q=0;
  bool ok=true;
  bool last=false;
  try {
  while(1) {
    if(last) break;
    fs.fs[p].minexp=0;
    q++;
    if((q&1) && (fs.special() || (p%4)!=1))
      continue;

    if(flags & C_NBRDP) {
      fs.fs[p].nondiv=1;
      if(pow(p,q+1)>pow(int_t(10),LIMEXPNBRDP)) {
	last=true;
	fs.fs[p].minexp=q;
	fs.fs[p].maxexp=Q_INF;
      } else {
	fs.fs[p].minexp=q;
	fs.fs[p].maxexp=q;
      }
    } else {
      int allexact=false;
      if(p<sizeof(taballexact)/sizeof(taballexact[0])) allexact=taballexact[toint(p)];
      int star=0;
      if(p<sizeof(tabstar)/sizeof(tabstar[0])) star=tabstar[toint(p)];
      //if(q%2==0 && (p!=3 || q!=8) && firstprimes[q+1]) continue;
      assert(!(star!=0 && star!=2 && star!=4 ));
      assert(!(star==4 && (p%4)==1));
      if(allexact==0) {
	if(star==0) {
	  if(firstprimes[q+1]) continue;
	}
	if((star==2&&(p%4)!=1) && q!=8 && firstprimes[q+1]) continue;
	if((star==2&&(p%4)==1) && q!=8 && q!=5 && firstprimes[q+1]) continue;
	if((star==4) && q!=8 && q!=14 && q!=24 && firstprimes[q+1]) continue;
	fs.fs[p].maxexp=-1;
	fs.fs[p].minexp=q;

	if(q<=star) fs.fs[p].maxexp=q;
	if(q==1) fs.fs[p].maxexp=q;
	
      } else {
	fs.fs[p].maxexp=q;
	fs.fs[p].minexp=q;
      }
    }

    if(forcedq>=0) {
      if(q<forcedq) continue;
      if(q>forcedq) break;
    }

#if 0
    if(!forced.compat(prof,p,q)) { //for 'U' option

      /* contradictions */
      if(q>=2 && checkcontradiction(fs)) {
	tout<<"break without process full path...\n";

	break;
      }
      
      continue; 
    }
#endif

   if(verb>=2) indcout()<<"~ "<< p << " " << q <<endl;

    /* "special" branching contradiction */
    if(fs.fs[p].maxexp>=0 && fs.fs[p].exp>fs.fs[p].maxexp) {
      if(verb>=2) indcout()<<" exp>max"<<endl;
      continue;
    }

    /* contradictions */
    if(checkcontradiction(fs)) {
      if(q>=2) break;
      continue;
    }

    if(fs.fs[p].nondiv==0 || fs.fs[p].minexp==fs.fs[p].maxexp) {

    int_t bs=sigma(p,q);
    int_t b;

    if(verb>=3) {
      indcout()<<"current: ";fs.aff2(tout);
      indcout()<<"s("<<p<<"^"<<q<<")= "<<bs<<endl;
    }

///VIRER
#ifndef NOFACT
    searchfactorin(fs,p,q,bs);
#endif
///\VIRER

    pfact_t fa,fsp;
    {
      int rf;
      bool skip=false;

///VIRER
#ifndef NOFACT
      int ii;
      int tr[]={0,1,2,-2};
      if(nofact) tr[1]=-2;
      for(ii=0;tr[ii]>-2;ii++) {
	b=bs;
	fa=pfact_t();
	
	rf=tryfactorize(p,q,fa,b,tr[ii]);

#ifdef DOUBLECHECK
	if(fa.prod()!=bs) {
	  cerr << "assert prod error: "<<bs << " != "; fa.aff(cerr); tout<<endl;
	  abort();
	}
#endif
	
	if(checkcontradiction2(fsp,fs,fa,borne)) {
	  skip=1;
	  break;
	}

	if(rf) break;
      }
#else 
///\VIRER

      b=bs;
      rf=oraclefact(p,q,fa,b);

#ifdef DOUBLECHECK
      if(fa.prod()!=bs) {
	cerr << "assert prod error: "<<bs << " != "; fa.aff(cerr); tout<<endl;
	abort();
      }
#endif
      
      if(checkcontradiction2(fsp,fs,fa,borne)) {
	skip=1;
      }
      
///VIRER
#endif
///\VIRER

      if(verb>=2) {
	indcout()<<"s("<<p<<"^"<<q<<")= ";fa.aff(tout);
	indcout()<<"c: ";fsp.aff2(tout);
      }
      
      if(skip) continue;
      
      if(rf!=1) {
	if(verb>=3) indcout()<<"cannot fact "<<p << " " <<q+1<<" "<< b <<endl;
///VIRER
#ifndef NOFACT
	//if(thcomp) (*thcomp)<<p<<" " << q << " " << b<<endl;
#endif
///\VIRER
	
	if(flags&B_NOCMP) { //Pascal's favourite option
	  ok=false;
	  continue;
	}
      }
    }

    list<int_t> lb;
    if(flags&B_BIG) fsp.bests(lb);
    else if(flags&B_SMALL) fsp.bestssmalls(lb);
    else if(flags&B_SMALL1M3) fsp.bestssmalls1m3(lb);
    else abort();
///VIRER
#ifndef NOFACT

    /* if no factor, 'lastchance' factoring... */
    
    if(!lb.size()) {
      if(verb>=1) {
	indcout()<< "cannot branch : "; // << endl;
	fsp.aff2(tout);
	indcout()<< "try fact"<< endl;
      }
      if(tryfactorize(fsp)) {
	if(verb>=1) indcout()<< " OK no RB\n";
	if(flags&B_BIG) fsp.bests(lb);
	else if(flags&B_SMALL) fsp.bestssmalls(lb);
	else if(flags&B_SMALL1M3) fsp.bestssmalls1m3(lb);
	else abort();
      }
    }

#endif

///\VIRER

    if(outrbcomp && !lb.size())
      for(map<int_t,exp_t>::const_iterator it=fsp.fs.begin();it!=fsp.fs.end();it++) {
	if(it->second.prime) continue;
	(*outrbcomp)<<it->second.p<<" "<<it->second.q<<" "<<it->first<<endl;
      }

    /* if no factor, try CRB */
    
    if(!lb.size()) {

///VIRER

      int mwosize=stat_rec;
      
      if(1) {
	char fn[100];
	snprintf(fn,100,"rb_lev_%d.txt",prof.crb());
	ofstream rb(fn,ios_base::out|ios_base::app);
	path_t p2=prof;
	p2.add(p,q);
	p2.aff(rb);
	rb<<endl;
      }

///\VIRER

      bool r=false;
      if(verb>=1) {
	indcout()<< "no branch ! : "; // << endl;
	fsp.aff2(tout);
      }

      { //aff
	mpq_class ab=fsp.Ab(),abw=fsp.Ab(true);
	double fl=ab.get_d(),flw=abw.get_d();
	if(verb>=1) indcout() << "Ab= " << fl << " Abw=" << flw << endl;
      }

      if(flags & C_NBRP) {
	{ //aff 
	  mpq_class ab=fsp.Ab(),abw=fsp.Ab(true);
	  double fl=ab.get_d(),flw=abw.get_d();
	  int re=LIMITNBRP-fsp.nbrp();
	  int np=(re+1)/2;
	  assert(np>0);
	  double lala=2/fl;
	  double lim=pow(lala,1./np);
	  assert(lim>1.);
	  double lp=1./(lim-1.);
	  if(verb>=1) indcout() << "P max: " << lp << endl; 
	  if(flw<2) {
	    double lala=2/flw;
	    double lim=pow(lala,1./np);
	    assert(lim>1.);
	    double lp=1./(lim-1.);
	    if(verb>=1) indcout() << "P max w: " << lp << endl; 
	  }
	}

	// CRB:
	path_t p2=prof;
	p2.add(p,q);
	r=rec_range(fsp,C_NBRP,p2,borne);

	if(r==false) ok=false;
	if(forced.h()>= p2.h()) break;
	
      } else 
      if(flags & C_NBRPM2) {
	// CRB:
	path_t p2=prof;
	p2.add(p,q);
	r=rec_range(fsp,C_NBRPM2,p2,borne);

	if(r==false) ok=false;
	if(forced.h()>= p2.h()) break;
	
      } else 
      if(flags & C_PROD) {
	// CRB:
	path_t p2=prof;
	p2.add(p,q);
	r=rec_range(fsp,C_PROD,p2,borne);

	if(r==false) ok=false;
	if(forced.h()>= p2.h()) break;
      } else if(flags & C_NBRDP) {
	path_t p2=prof;
	p2.add(p,q);
	r=rec_range(fsp,C_NBRDP,p2,borne);
	
	if(r==false) ok=false;
	if(forced.h()>= p2.h()) break;
      } else { //CBR not handled
	if(verb>=0) indcout()<< " RB !!!\n";

	abort();
	path_t p2=prof;
	p2.add(p,q);

	ok=false;

	if(forced.h()>= p2.h()) break;
      }

      ///VIRER

      int mwsize=stat_rec-mwosize;

#ifdef MOSTWANTED
      for(map<int_t,exp_t>::const_iterator it=fsp.fs.begin();it!=fsp.fs.end();it++) {
	if(it->second.prime) continue;
	pair<int_t,pair<int_t,int> > x=make_pair(it->first,make_pair(it->second.p,it->second.q));
	mostwanted[x]+=mwsize;
      }
#endif
      ///\VIRER

      continue;
    }

    /* if factor(s), branch on the "best" prime */

    list<int_t>::const_iterator it=lb.begin();
    path_t p2=prof;
    p2.add(p,q);
    if(!rec(fsp,*it,p2,borne)) ok=false;

    if(forced.h()>= p2.h()) break;

    } else { //NBRDP 

      pfact_t fsp;
      fsp=fs;
      
      if(verb>=3) {
	indcout()<<"current: ";fs.aff2(tout);
	indcout()<<p<<"^(>="<<q<<") "<<endl;
      }
      
      list<int_t> lb;
      if(flags&B_BIG) fsp.bests(lb);
      else if(flags&B_SMALL) fsp.bestssmalls(lb);
      else abort();
      
      if(outrbcomp && !lb.size())
	for(map<int_t,exp_t>::const_iterator it=fsp.fs.begin();it!=fsp.fs.end();it++) {
	  if(it->second.prime) continue;
	  (*outrbcomp)<<it->second.p<<" "<<it->second.q<<" "<<it->first<<endl;
	}
      
      /* if no factor, try CRB */
      
      if(!lb.size()) {
	bool r=false;
	if(verb>=1) {
	  indcout()<< "no branch ! : "; // << endl;
	  fsp.aff2(tout);
	}
	
	{ //aff
	  mpq_class ab=fsp.Ab(),abw=fsp.Ab(true);
	  double fl=ab.get_d(),flw=abw.get_d();
	  if(verb>=1) indcout() << "Ab= " << fl << " Abw=" << flw << endl;
	}
	
	if(flags & C_NBRDP) {
	  path_t p2=prof;
	  p2.add(p,q);
	  r=rec_range(fsp,C_NBRDP,p2,borne);

	  if(r==false) ok=false;
	  if(forced.h()>= p2.h()) break;
	  
	} else { //CBR not handled
	  if(verb>=0) indcout()<< " RB !!!\n";

	  abort();
	  path_t p2=prof;
	  p2.add(p,q);
	  
	  ok=false;
	  
	  if(forced.h()>= p2.h()) break;
	}
	
	continue;
      }
      
      /* if factor(s), branch on the "best" prime */
      
      list<int_t>::const_iterator it=lb.begin();
      path_t p2=prof;
      p2.add(p,q);
      if(!rec(fsp,*it,p2,borne)) ok=false;
      
    }

  }
  } catch(gcdex_t e) {
    indcout()<< "GCD error..."<<endl;
    ok=false;
  }

  if(verb>=2) indcout()<<"end " << p << " " << q << endl;

  envprof=oprof;

  if(prof.h()<limsize) {
    prof.aff(tout); 
    tout << " [ size= " << stat_rec-ostat<< " ret= " << ok <<" ]\n";
  }
  
  return ok;
}

bool main_proc(path_t forced_,int uniq=0, int after=0)
{
  bool ret=1;
  forced=forced_;
  if(forced.h()) {
    tout << "*** "; forced.aff(tout); tout << " ***\n";
    tout.flush();
    uniq=toint(forced.p.begin()->first);
    after=0;
  }
  for(int mit=0;upto || tab[mit] ;mit++) {
    int i=tab2(mit);
    ///VIRER
#if 0
    if(!i) {
      static int io=tab2(nbrprimes-2);
      io++;
      while(!isprime(io)) io++;
      i=io;
    }
#endif
    ///\VIRER
    assert(i>2);
    if(upto && i>upto) break;
    if(uniq && i!=uniq) continue;
    if(after) uniq=0;
#ifdef NOTH
    if(mit<100 || mit%5000==0) 
      fprintf(stderr,"i=%d        \r",tab2(mit));
#endif
    path_t p;
    pfact_t pf = pfact_t(int_t(i));
    ll_t osr=stat_rec;
    bool rr=rec(pf,i,p,mit);
    if(!rr) ret=0;
    tout << "\n--- RETURN FOR " << i;
    if(forced.h()) {tout<<" / "; forced.aff(tout);}
    tout<< " : " << rr<< " (ar : " << ret << " size= " <<stat_rec-osr<<" ) "<<endl;
    tout.flush();
  }
  
  tout<< "stat_rec=" << stat_rec <<endl;
  return ret;
}

bool main_proc_uniq(int ii)
{
  forced=path_t();

  tout << "*** "<< tab2(ii) <<  " ***\n";
  tout.flush();
  
  path_t p;
  pfact_t pf = pfact_t(int_t(tab2(ii)));
  bool rr=rec(pf,tab2(ii),p,ii);
  tout << "\n--- RETURN FOR " << tab2(ii) << " : " << rr<<endl;
  tout.flush();
  return rr;
}

};

#ifndef NOTH

pthread_t thtab[128];
pthread_mutex_t mutex;
int thmode=0;
int thmit=0;
int thupto=0;

void *thread_main(void *pt)
{
  int id=*(int*)pt;
  th_t th(id);
  
  printf("START THREAD %d\n",id);fflush(stdout);

  bool ret=1;
  while(1) {
    path_t f;
    bool end=false;

    pthread_mutex_lock(&mutex);

    int ip=-1;
    if(!thmode) {
      if(forceds.empty()) end=true;
      else {f=forceds.front();forceds.pop_front();}
    } else {
      if(thupto || tab[thmit]!=0) {
	ip=thmit;
	int p=tab2(thmit);
	assert(p>2);
	if(thupto && p>thupto) ip=-1;
      }
      if(ip>=0) thmit++;
      else end=true;
    }
    
    pthread_mutex_unlock(&mutex);
    
    if(end) break;
    
    int rr=(thmode ? th.main_proc_uniq(ip) : th.main_proc(f));
    if(!rr) {
      ret=0;
      pthread_mutex_lock(&mutex);
      forceds_fail.push_back(f);
      pthread_mutex_unlock(&mutex);
    }
  }

  printf("END THREAD %d RET=%d\n",id,ret);fflush(stdout);

  return NULL;
}
#endif

int main(int ac, char **av)
{
  path_t forced;
  int uniq=0,after=0;
  int nbrthreads=1;

  ///VIRER
  int retry=0;
  ///\VIRER

  {
    strcpy(cmdline,"");
    for(int i=0;i<ac;i++) {
      if(i) strcat(cmdline," ");
      strcat(cmdline,av[i]);
    }
    printf("nohup %s > log.log &\n",cmdline);
  }

///VIRER
  cout<< "define: ";
#ifdef DOUBLECHECK
  cout << "DOUBLECHECK ";
#endif

  cout<<endl;
///\VIRER

  int addflags=0;
  while(ac>1 && av[1][0]=='-') {
    if(!strcmp(av[1]+1,"v")) {verb++;ac-=1;av+=1;continue;}
    if(!strcmp(av[1]+1,"q")) {verb--;ac-=1;av+=1;continue;}

    if(!strcmp(av[1]+1,"SIZE")) {BE=atoi(av[2]);ac-=2;av+=2;continue;}
    if(!strcmp(av[1]+1,"NBRP")) {
      tab=tabP;      
      flags=F_NBRP;
      LIMITNBRP=atoi(av[2]); 
      assert(LIMITNBRP%2);
      ac-=2;av+=2;
      continue;
    }

    if(!strcmp(av[1]+1,"NBRDP")) {
      tab=tabP3;
      flags=F_NBRDP;
      LIMITNBRP=atoi(av[2]); 
      if(av[1][2]=='F') {assert(0); /*flags=F_NBRPF;*/}
      ac-=2;av+=2;continue;
    }
    
    if(!strcmp(av[1]+1,"ED")) {LIMEXPNBRDP=atoi(av[2]);ac-=2;av+=2;continue;}

///VIRER
    if(!strcmp(av[1]+1,"NBRPM2")) {
      tab=tabP19;      
      flags=F_NBRPM2;
      LIMITNBRP=atoi(av[2]); 
      assert(LIMITNBRP%2);
      ac-=2;av+=2;
      continue;
    }
///\VIRER

    if(!strcmp(av[1]+1,"COMP")) {
      tab=tabC;
      flags=F_COMP;
      BE=atoi(av[2]);
      ac-=2;av+=2;continue;
    }

    if(!strcmp(av[1]+1,"COMPM2")) {
      tab=tabC;
      flags=F_COMPM2;
      BE=atoi(av[2]);
      ac-=2;av+=2;continue;
    }

    if(!strcmp(av[1]+1,"u")) {uniq=atoi(av[2]); ac-=2;av+=2;continue;}
    if(!strcmp(av[1]+1,"a")) {uniq=atoi(av[2]); after=1;ac-=2;av+=2;continue;}
    if(!strcmp(av[1]+1,"U")) {
      int_t p=int_t(av[2]);
      int q=atoi(av[3]);
      forced.add(p,q);
      av+=3;ac-=3;
      continue;
    }

    if(!strcmp(av[1]+1,"pm")) {profmax=atoi(av[2]);ac-=2;av+=2;continue;}
    if(!strcmp(av[1]+1,"lm")) {limsize=atoi(av[2]);ac-=2;av+=2;continue;}

    if(!strcmp(av[1]+1,"t")) {nbrthreads=atoi(av[2]);ac-=2;av+=2;continue;}

    if(!strcmp(av[1]+1,"as")) {assumedtrialcomp=int_t(av[2]);ac-=2;av+=2;continue;}

    if(!strcmp(av[1]+1,"fileprefix")) {prefixname=(av[2]);ac-=2;av+=2;continue;}

    if(!strcmp(av[1]+1,"fU")) {
      FILE *in=fopen(av[2],"rt");
      assert(in);
      char bf[100000];
      while(fgets(bf,100000,in)) {
	assert(strlen(bf)<99900);
	if(bf[0]=='#' || bf[0]==0 || bf[0]=='\n') continue;
	forceds.push_back(path_t(bf));
      }
      fclose(in);
      av+=2;ac-=2;
      continue;
    }

    if(!strcmp(av[1]+1,"F")) {addflags|=B_NOCMP;ac-=1;av+=1;continue;}

    ///VIRER
    if(!strcmp(av[1]+1,"retry")) {retry++;ac-=1;av+=1;continue;}
    ///\VIRER

    if(!strcmp(av[1]+1,"crb")) {maxcrb=atoi(av[2]);ac-=2;av+=2;continue;}

    if(!strcmp(av[1]+1,"ut")) {
      double lala=atoi(av[2]);
      upto=pow(10.,lala);
      printf("upto %d\n",upto);
      ac-=2;av+=2;continue;
    }

    if(!strcmp(av[1]+1,"UT")) {
      upto=atoll(av[2]);
      printf("upto %d\n",upto);
      ac-=2;av+=2;continue;
    }
    
    if(!strncmp(av[1]+1,"X",1)) {
      int k=atoi(av[1]+2);
      int i=atoi(av[2]); 
      assert(k==2 || k==4);
      assert(i<sizeof(tabstar)/sizeof(tabstar[0]));
      tabstar[i]=k;
      ac-=2;av+=2;
      continue;
    }

    if(!strncmp(av[1]+1,"Y",1)) {
      int i=atoi(av[2]); 
      assert(i<sizeof(taballexact)/sizeof(taballexact[0]));
      taballexact[i]=1;
      ac-=2;av+=2;
      continue;
    }

    if(!strcmp(av[1]+1,"rbcomp")) {
      printf("rbcomp to %s\n",av[2]);
      outrbcompname=av[2];
      //outrbcomp=new ofstream(av[2],ios_base::out|ios_base::app);
      ac-=2;av+=2;
      continue;
    }

///VIRER
#ifndef NOFACT
    if(!strcmp(av[1]+1,"writecomp")) {
      printf("write comp to %s\n",av[2]);
      outcomp=new ofstream(av[2],ios_base::out|ios_base::app);
      ac-=2;av+=2;
      continue;
    }

    /*if(!strcmp(av[1]+1,"writethcomp")) {
      printf("write thcomp to %s\n",av[2]);
      thcomp=new ofstream(av[2],ios_base::out|ios_base::app);
      ac-=2;av+=2;
      continue;
    }//*/
#endif

    if(!strcmp(av[1]+1,"ns")) {nospecialinb=true;ac-=1;av+=1;continue;}
    if(optoracle(ac,av)) continue;
    ///\VIRER
       
    printf("bad option: %s\n",av[1]);
    exit(1);
  }
  flags|=addflags;

  printf("list of factors: ");
  for(int i=0;tab[i];i++) {printf("%d ",tab[i]);}
  printf("\n");

  LIMIT=1;
  for(int i=0;i<BE;i++) LIMIT*=10;

  nbrbaseisprime=0; // only 2 bases for a quick check when read facts 

  initoracle(ac,av);

  makesmallfacts(3,8); makesmallfacts(3,14); makesmallfacts(3,24); makesmallfacts(3,48);
  makesmallfacts(5,5); makesmallfacts(5,8);
  makesmallfacts(7,8);
  makesmallfacts(11,8);
  makesmallfacts(13,5); makesmallfacts(13,8);
  makesmallfacts(17,5); makesmallfacts(17,8);
  makesmallfacts(19,8);
  makesmallfacts(23,8);    

  inittab2();

  nbrbaseisprime=5;

///VIRER
  facts_t::facts_stats();
  
///\VIRER
  if(forced.h()) forceds.push_front(forced);

  ///VIRER

  while(1) {

    gcd_error=0;

    list<path_t> forceds_copy;
    if(retry==2) forceds_copy=forceds;
  ///\VIRER
    
    printf("forceds.size= %d\n",int(forceds.size()));

#ifdef NOTH

  if(forceds.empty()) forceds.push_back(path_t());

  printf("START\n");
  fprintf(stderr,"START      \r");

  {
    th_t th;
    
    for(list<path_t>::iterator it=forceds.begin();it!=forceds.end();++it) th.main_proc(*it,uniq,after);
  }
#else

  if(forceds.empty()) {
    thmode=1;
    assert(!uniq || after);
    if(uniq) {
      for(thmit=0;tab2(thmit) && tab2(thmit)!=uniq;thmit++);
      assert(tab2(thmit));
    }
    if(upto) thupto=upto;
  }

  printf("START\n");fprintf(stderr,"START      \r");

  if(nbrthreads<=1) {
    int zero=0;
    thread_main(&zero);
  } else {
    int tid[nbrthreads];
    for(int i=0;i<nbrthreads;i++) {
      tid[i]=i;
      pthread_create(&thtab[i],NULL,thread_main,tid+i);
    }
    
    printf("WAIT...\n");fflush(stdout);

    for(int i=0;i<nbrthreads;i++) {
      pthread_join(thtab[i],NULL);
      printf("JOIN THREAD %d\n",i);fflush(stdout);
    }
  }

#endif

#if 1 ///VIRER

  printf("gcd_error=%d\n",gcd_error);

  if(retry) {
    if(retry ==2) {
      forceds=forceds_copy;
      while(1) {
	printf("do what ?\n");
	char bf[10000];
	bf[0]=0;

	if(fgets(bf,10000,stdin)) for(int i=0;bf[i];i++) {
	  if(bf[i]=='\n') {
	    bf[i]=0;
	    break;
	  }
	}
	if(bf[0]==0) continue;

#ifndef NOFACT
	if(!strncmp(bf,"C ",2)) {
	  ecmmaxrb=atoi(bf+2);
	  printf("C: %d\n",ecmmaxrb);
	  continue;
	}

	if(!strncmp(bf,"c ",2)) {
	  ecmmax=atoi(bf+2);
	  printf("c: %d\n",ecmmax);
	  continue;
	}

	if(!strcmp(bf,"thclear")) {
	  fstatus.clear();
	  printf("th={}\n");
	  continue;
	}

#endif
	if(!strncmp(bf,"crb ",4)) {
	  maxcrb=atoi(bf+4);printf("crb=%d\n",maxcrb);continue;
	}

	if(!strcmp(bf,"Uclear")) {
	  forceds.clear();
	  printf("forceds={}\n");
	  continue;
	}

	if(!strcmp(bf,"Ulist")) {
	  for(list<path_t>::iterator it=forceds.begin();it!=forceds.end();++it) {
	    it->aff(cout);cout<<endl;
	  }
	  continue;
	}

	if(!strncmp(bf,"U ",2)) {
	  forceds.push_back(path_t(bf+2));
	  continue;
	}

	if(!strncmp(bf,"fU ",3)) {
	  forceds.clear();
	  FILE *in=fopen(bf+3,"rt");
	  if(in==NULL) printf("file not found\n");
	  else {
	    char bf[100000];
	    while(fgets(bf,100000,in)) {
	      assert(strlen(bf)<99900);
	      if(bf[0]=='#' || bf[0]==0 || bf[0]=='\n') continue;
	      forceds.push_back(path_t(bf));
	    }
	    fclose(in);
	  }
	  printf("forceds.size= %d\n",int(forceds.size()));
	  continue;
	}

	if(!strcmp(bf,"cont")) break;
	if(!strcmp(bf,"continue")) break;
	
	if(!strncmp(bf,"read ",5)) {
	  readfactors(bf+5);
	  continue;
	}
	
	printf("unknown cmd '%s'\n",bf);

	continue;
      }
      forceds_fail.clear();
      continue;
    }
    if(forceds_fail.size()) {
      readfactors("gcdfacts.txt");
      forceds.swap(forceds_fail);
      continue;
    }
  }
#endif 

  break;
  }
///\VIRER

  printf("END\n");fflush(stdout);
  fprintf(stderr,"END\n");fflush(stderr);

#ifdef MOSTWANTED ///VIRER

  if(!(flags&B_NOCMP)){
    ofstream mwout("mostwantedrb.txt",ios_base::out|ios_base::app);
    mwout<<"#"<<cmdline<<endl;
    map<pair<int_t,pair<int_t,int> >,int >::iterator it;
    map<int,list<pair<int_t,pair<int_t,int> > > > mwback;
    map<int,list<pair<int_t,pair<int_t,int> > > >::iterator it2;
    for(it=mostwanted.begin();it!=mostwanted.end();it++) mwback[it->second].push_back(it->first);
    for(it2=mwback.begin();it2!=mwback.end();++it2)
      for(list<pair<int_t,pair<int_t,int> > >::iterator it3=it2->second.begin();it3!=it2->second.end();it3++)
        mwout<<it3->second.first<<" " <<it3->second.second<<" "<<it3->first<<" " <<it2->first<<endl;
  }
#endif ///\VIRER

  desinitoracle (ac, av);

  return 0;
}