\\
\\---------------------------------------------------------------------------------
\\ Testing of prime numbers
\\---------------------------------------------------------------------------------
\\	input: NO  -  pr1()  N=4k+1 Only by using a new method
\\	input: NO  -  pr1ok() N=4k+1 new method and tested 2^(N-1)=1
\\	input: NO  -  pr3()  N=4k-1 Only by using a new method
\\	input: NO  -  pr3ok() N=4k-1 new method and tested 2^(N-1)=1
\\---------------------------------------------------------------------------------
\\
\\	allocatemem(1024000000);
\\
\\---------------------------------------------------------------------------------
pr1()={
\\
	local(pr,nmin,nnum,npr,fmoc,flOK);
\\
	print("first number  -  q=4k+1 =?");
	nmin=input();
	if(nmin%2==0,nmin++);
	if(nmin%4==3,nmin=nmin+2);
	if(nmin%4==1, ,print("ERROR - nmin%4!=1");return());
	print("number of primes =?");
	nnum=input();
	print(" ");
	print("nmin=",nmin);
	print("nnum=",nnum);
	print(" ");
\\
	flOK=1;
	pr=nmin;
	npr=0;
	while(1,
		if(isprime(pr),
			npr++;
			fmoc=moc(pr,1,2,(pr-1)/2);
			if(fmoc==0 || fmoc==-pr, ,flOK=0; printf("ERROR1  -  pr=%20d  fmoc=%20d\n",pr,fmoc));
			if(npr==nnum,print("\nlast prime= ",pr);break());
		,
			fmoc=moc(pr,1,2,(pr-1)/2);
			if(fmoc==0 || fmoc==-pr,flOK=0; printf("Pseudoprime  -  pr=%20d  fmoc=%20d\n",pr,fmoc));
		);
		pr=pr+4;
	);
\\
	if(flOK,print("\nTEST - OK\n\n"));
	return();
\\
}
\\---------------------------------------------------------------------------------
pr3()={
\\
	local(pr,nmin,nnum,npr,fmoc,flOK);
\\
	print("first number  -  q=4k-1 =?");
	nmin=input();
	if(nmin%2==0,nmin++);
	if(nmin%4==1,nmin=nmin+2);
	if(nmin%4==3, ,print("ERROR - nmin%4!=3");return());
	print("number of primes =?");
	nnum=input();
	print(" ");
	print("nmin=",nmin);
	print("nnum=",nnum);
	print(" ");
\\
	flOK=1;
	pr=nmin;
	npr=0;
	while(1,
		if(isprime(pr),
			npr++;
			fmoc=moc(pr,1,2,(pr+1)/2);
			if(fmoc==0 || fmoc==-pr, ,flOK=0; printf("ERROR1  -  pr=%20d  fmoc=%20d\n",pr,fmoc));
			if(npr==nnum,print("\nlast prime= ",pr);break());
		,
			fmoc=moc(pr,1,2,(pr+1)/2);
			if(fmoc==0 || fmoc==-pr,flOK=0; printf("Pseudoprime  -  pr=%20d  fmoc=%20d\n",pr,fmoc));
		);
		pr=pr+4;
	);
\\
	if(flOK,print("\nTEST - OK\n\n"));
	return();
\\
}
\\---------------------------------------------------------------------------------
pr1ok()={
\\
	local(pr,nmin,nnum,npr,fmoc,flOK);
\\
	print("first number  -  q=4k+1 =?");
	nmin=input();
	if(nmin%2==0,nmin++);
	if(nmin%4==3,nmin=nmin+2);
	if(nmin%4==1, ,print("ERROR - nmin%4!=3");return());
	print("number of primes  =?");
	nnum=input();
	print(" ");
	print("nmin=",nmin);
	print("nnum=",nnum);
	print(" ");
\\
	flOK=1;
	pr=nmin;
	npr=0;
	while(1,
		if(isprime(pr),
			npr++;
			fmoc=moc(pr,1,2,(pr-1)/2);
			if(fmoc==0 || fmoc==-pr, ,flOK=0; printf("ERROR1  -  pr=%20d  fmoc=%20d\n",pr,fmoc));
			if(npr==nnum,print("\nlast prime= ",pr);break());
		,
			fmoc=moc(pr,1,2,(pr-1)/2);
			if(fmoc==0 || fmoc==-pr,
				m2=Mod(2,pr)^(pr-1);
				if(component(m2,2)==1,
					flOK=0;
					printf("Pseudoprime  -  pr=%20d  fmoc=%20d\n",pr,fmoc);
				,
				);
			);
		);
		pr=pr+4;
	);
\\
	if(flOK,print("\n\nTEST - OK\n\n"));
	return();
\\
}
\\---------------------------------------------------------------------------------
pr3ok()={
\\
	local(pr,nmin,nnum,npr,fmoc,flOK);
\\
	print("first number  -  q=4k-1 =?");
	nmin=input();
	if(nmin%2==0,nmin++);
	if(nmin%4==1,nmin=nmin+2);
	if(nmin%4==3, ,print("ERROR - nmin%4!=3");return());
	print("number of primes =?");
	nnum=input();
	print(" ");
	print("nmin=",nmin);
	print("nnum=",nnum);
	print(" ");
\\
	flOK=1;
	pr=nmin;
	npr=0;
	while(1,
		if(isprime(pr),
			npr++;
			fmoc=moc(pr,1,2,(pr+1)/2);
			if(fmoc==0 || fmoc==-pr, ,flOK=0; printf("ERROR1  -  pr=%20d  fmoc=%20d\n",pr,fmoc));
			if(npr==nnum,print("\nlast prime= ",pr);break());
		,
			fmoc=moc(pr,1,2,(pr+1)/2);
			if(fmoc==0 || fmoc==-pr,
				m2=Mod(2,pr)^(pr-1);
				if(component(m2,2)==1,
					flOK=0;
					printf("Pseudoprime  -  pr=%20d  fmoc=%20d\n",pr,fmoc);
				,
				);
			);
		);
		pr=pr+4;
	);
\\
	if(flOK,print("\n\nTEST - OK\n\n"));
	return();
\\
}
\\---------------------------------------------------------------------------------
\\---------------------------------------------------------------------------------
qcomp(x1,x2,N,cc)={
	local(xsuma,gc);
	xsuma=x1+x2;
	gc=gcd(xsuma,N);
	if(gc==1, ,return(-gc));
	return(((x1*x2-cc^2)/xsuma)%N);
}
\\---------------------------------------------------------------------------------
moc(N,cc,y,ex)={
	local(i,bex,y0);
\\
	bex=binary(ex);
	y0=y;
	for(i=2,matsize(bex)[2],
		if(bex[i],
			y=qcomp(y,y,N,cc);
			if(y<0,return(y));
\\			if(y<0,printf("i=%5d  factorizacia-y0  %35d\n\n",i,-y);return());
			y=qcomp(y,y0,N,cc);
			if(y<0,return(y));
\\			if(y<0,printf("i=%5d  factorizacia-y1  %35d\n\n",i,-y);return());
		,
			y=qcomp(y,y,N,cc);
			if(y<0,return(y));
\\			if(y<0,printf("i=%5d  factorizacia-y2  %35d\n\n",i,-y);return());
		);
	);
\\
	return(y);
}
\\---------------------------------------------------------------------------------