/* sqf.c ssw 4-4-91 */
#include <stdio.h>
#include <math.h>

#define MPDIM		50		/* size of MP arrays		*/
#define	SIZE(x)		x[0]		/* Size of x in base RADIX	*/
#define LAST(x)		x[1]		/* Last digit of x		*/
#define SINGLE		2		/* single precision size of MP	*/

int qqueue[20],qpoint;
char s3465[3465];

main()
{
int N[MPDIM];
long fac1,fac2;

printf("Give N: ");
get_file(N,stdin);			/* read the number		*/
sqfprep();
(void)squfof(N,&fac1,&fac2);
}

squfof(n,fac1,fac2)
int n[];
long *fac1, *fac2;
{	/* begin SQUFOF: factor n = fac1*fac2 */
int temp[MPDIM], temp1[MPDIM];
register int iq,l,l2,p,pnext,q,qlast,r,s,t,i;
long jter,iter;

qlast = 1;
mpsqrt(n,temp); if (SIZE(temp) > SINGLE) return(0);
s = LAST(temp);
p = s;
mult(temp,temp,temp1);
subt(n,temp1,temp);	/* temp = n - floor(sqrt(n))^2 */
if (SIZE(temp) <= SINGLE && LAST(temp) == 0)
	{	/* Here n is a square */
	*fac1 = s;
	*fac2 = s;
	return(1);
	}
q = LAST(temp);	/* q = excess of n over next smaller square */
l = 1 + 2*(int)sqrt((double)(p+p));
l2 = l/2;
qpoint = 0;
for (jter=0; jter < 20000; jter++)
	{	/* begin main loop */
	iq = (s + p)/q;
	pnext = iq*q - p;
	if (q <= l)
		{
		if ((q & 1) == 0) enqu(q/2,&jter);
		else if (q <= l2) enqu(q,&jter);
		if (jter < 0) return(0); /* queue overflow: give up */
		}
	t = qlast + iq*(p - pnext);
	qlast = q;
	q = t;
	p = pnext;
	if (jter & 1 == 1) continue;	/* jter is odd: omit square test */
	if (t = (q & 3) > 1) continue;	/* skip 2,3 mod 4 */
	if (q & 7 == 5) continue;	/* skip 5 mod 8 */
	if (t == 0 && (q & 15) > 5) continue;	/* skip 8, 12 mod 16 */
	if (s3465[q%3465] == 'n') continue;	/* skip QNR's mod 3465 */
	r = (int)sqrt((double)q);	/* r = floor(sqrt(q)) */
	if (q != r*r) continue;
	if (qpoint == 0) goto gotit;
	for (i=0; i<qpoint; i++)	/* treat queue as list for simplicity */
		if (r == qqueue[i]) goto contin;
	goto gotit;
contin: ; }	/* end of main loop */
return(0);	/* takes too long: give up */
gotit:	;
qlast = r;
p = p + r*((s - p)/r);
SIZE(temp1) = SINGLE; LAST(temp1) = p;
mult(temp1,temp1,temp);
subt(n,temp,temp);
div_single(temp,qlast,temp1);
q = LAST(temp1);	/* q = (n - p*p)/qlast (division is exact) */
for (iter=0; iter<20000; iter++)
	{	/* begin second main loop */
	iq = (s + p)/q;
	pnext = iq*q - p;
	if (p == pnext) goto gotfac;
	t = qlast + iq*(p - pnext);
	qlast = q;
	q = t;
	p = pnext;
	}
return(0);	/* this shouldn't happen */
gotfac:	; if ((q & 1) == 0) q/=2;	/* q was factor or 2*factor */
*fac1 = q;
div_single(n,q,temp);
if (SIZE(temp) > SINGLE) return (0);
*fac2 = LAST(temp);
printf("sqf success: iter1,iter2,facts = %d %d %d %d\n",jter,iter,*fac1,*fac2);
return(1);
}

sqfprep()
{
int i;
for (i=0; i<3465; i++) s3465[i] = 'n';
for (i=0; i<1733; i++) s3465[(i*i)%3465] = 'r';
}

enqu(q,iter)
int q;
long *iter;
{
qqueue[qpoint] = q;
printf("sqf: enqu %d at %d\n",q,*iter);
(void)fflush(stdout);
if (++qpoint < 20) return;
printf("sqf: qqueue overflow\n");
*iter = -1;
}

get_file(number,ptr)
int number[];
FILE *ptr;
{
char string[100];
(void)fscanf(ptr,"%s",string);
mpctoi(string,number);
}