find_a_ORDOWSKI(n=48, a=1, B_START=48, LIM=10^11,DEBUG=1)={
  my(i,j,B,FF,LL,CYCLO,p_intrinsic);
  my(fn="_THOMAS_ORDOWSKI_b_a_n.txt");
  
  if(isprime(n)==1 && a==1, return([n,2]));
  [FF,LL,p_intrinsic,CYCLO]=ALGEBRAIC_PART_DECOMPOSITION('b,a,n);

  print(FF);
  print(LL);
  print(p_intrinsic);
  print(CYCLO);

  if(Mod(n,2)==1,  \\ n-ODD
      B=FIND_BASE(n,BSTART=B_START,LIM,STEP=1,FF); 
  ,                \\ n-EVEN
      B=FIND_BASE(n,BSTART=n,LIM,STEP=n,FF);  	 
  );
  if(B>0,
	print_and_write(fn,n,B,a,CYCLO,p_intrinsic,LL,FF,DEBUG);
	return([n,B,[subst(FF,'b,B)]]);    
  );
  return(0);
}
\\----------------------------------------
 ALGEBRAIC_PART_DECOMPOSITION(base='b,a=1,n)={
  my(LL,CYCLO,p_intrinsic,F);
  
  LL=R(base,a,n);           \\ R(b,a,n)=(b^n-a^n)/(b-a) 
  CYCLO=polcyclo(n,base);   \\primitive part, p|CYCLO/intrinsic*->rank(p)=n
  p_intrinsic=gcd(n,polcyclo(n,base)); \\intrinsic divisor p*
  LL=LL/CYCLO;              \\algebraic part   
  LL=LL/p_intrinsic;
  F=factor(LL);
  
 return([F,LL,p_intrinsic,CYCLO]); 
}  
\\-------------------------------
FIND_BASE(n,BSTART,LIM,STEP,F)={
my(i,ll=0,sss=0,OK=0,f,fs,ROWS,b=0,B=0);
      
ROWS=matsize(F)[1];

forstep(b=BSTART, LIM, STEP,
	for(i=1,ROWS,
	  ll++;
	  f=F[i,1];
	  sss=subst(f,'b, b);
	  fs=factor(sss);  \\???
	  if(Mod(ll,100000)==0,
	    print(sss);
	    print(f);
		print(fs);
	  );
	  OK=check_factors(fs,b,n,ll);
      if(OK==0,
	    break; 
	  );	  
	); 
	if(OK==1, B=b; break);
 ); 
 return(B);
}
\\---------------------
check_factors(f,b,n,ll)={
  my(r,j,rows,p,OK=1,s="");	

  rows=matsize(f)[1]; 	  
  for(j=1,rows,
	p=f[j,1];
	if(Mod(ll,10000)==0,
  	  s=Strprintf("? b=%d p=%d",b,p); \\ processing && debug
	  print(s);
    );  	
    r=Mod(p-1,n);
    if(!(r==0),
	   OK=0; break;   		  
	);	  
  ); 
  return(OK);
}
\\---------------------
print_and_write(fn,n,b,a,CYCLO,p_intrinsic,LL,F,DEBUG=1)={
	my(s="",v=[]);
	
	if(debug==1,
	  s=Strprintf("b=%d",b);
      print(s);
      write(fn,s);
	  	  
      s=Strprintf("cyclo(%d,b)=%s",n,CYCLO);
      print(s);
      write(fn,s);
	  s=Strprintf("intrinsic=%s",p_intrinsic);
      print(s);
	  write(fn,s);
      s=Strprintf("R(b,%d,%d)/intrinsic/cyclo=%s",a,n,LL);
      print(s);
	  write(fn,s);
      s=Strprintf("F=%s",F);
      print(s);
	  write(fn,s);
	);
    v=[subst(F,'b,b)];
   	write(fn,[n,b,v]);
    if(DEBUG==1,	  
	  write(fn,"");
	);
};