/* gnubc program: factors */

"
typing factor(n) attempts to factors n using Brent-Pollard.
Consequently In general not very effective on numbers of more
than say 20 digits.
"
pglobal[0]=2
pglobal[1]=3
pglobal[2]=5
pglobal[3]=7
pglobal[4]=11
pglobal[5]=13
pglobal[6]=17
pglobal[7]=19
pglobal[8]=23
pglobal[9]=29
pglobal[10]=31
pglobal[11]=37
pglobal[12]=41
pglobal[13]=43
pglobal[14]=47
pglobal[15]=53
pglobal[16]=59
pglobal[17]=61
pglobal[18]=67
pglobal[19]=71
pglobal[20]=73
pglobal[21]=79
pglobal[22]=83
pglobal[23]=89
pglobal[24]=97
pglobal[25]=101
pglobal[26]=103
pglobal[27]=107
pglobal[28]=109
pglobal[29]=113
pglobal[30]=127
pglobal[31]=131
pglobal[32]=137
pglobal[33]=139
pglobal[34]=149
pglobal[35]=151
pglobal[36]=157
pglobal[37]=163
pglobal[38]=167
pglobal[39]=173
pglobal[40]=179
pglobal[41]=181
pglobal[42]=191
pglobal[43]=193
pglobal[44]=197
pglobal[45]=199
pglobal[46]=211
pglobal[47]=223
pglobal[48]=227
pglobal[49]=229
pglobal[50]=233
pglobal[51]=239
pglobal[52]=241
pglobal[53]=251
pglobal[54]=257
pglobal[55]=263
pglobal[56]=269
pglobal[57]=271
pglobal[58]=277
pglobal[59]=281
pglobal[60]=283
pglobal[61]=293
pglobal[62]=307
pglobal[63]=311
pglobal[64]=313
pglobal[65]=317
pglobal[66]=331
pglobal[67]=337
pglobal[68]=347
pglobal[69]=349
pglobal[70]=353
pglobal[71]=359
pglobal[72]=367
pglobal[73]=373
pglobal[74]=379
pglobal[75]=383
pglobal[76]=389
pglobal[77]=397
pglobal[78]=401
pglobal[79]=409
pglobal[80]=419
pglobal[81]=421
pglobal[82]=431
pglobal[83]=433
pglobal[84]=439
pglobal[85]=443
pglobal[86]=449
pglobal[87]=457
pglobal[88]=461
pglobal[89]=463
pglobal[90]=467
pglobal[91]=479
pglobal[92]=487
pglobal[93]=491
pglobal[94]=499
pglobal[95]=503
pglobal[96]=509
pglobal[97]=521
pglobal[98]=523
pglobal[99]=541
pglobal[100]=547
pglobal[101]=557
pglobal[102]=563
pglobal[103]=569
pglobal[104]=571
pglobal[105]=577
pglobal[106]=587
pglobal[107]=593
pglobal[108]=599
pglobal[109]=601
pglobal[110]=607
pglobal[111]=613
pglobal[112]=617
pglobal[113]=619
pglobal[114]=631
pglobal[115]=641
pglobal[116]=643
pglobal[117]=647
pglobal[118]=653
pglobal[119]=659
pglobal[120]=661
pglobal[121]=673
pglobal[122]=677
pglobal[123]=683
pglobal[124]=691
pglobal[125]=701
pglobal[126]=709
pglobal[127]=719
pglobal[128]=727
pglobal[129]=733
pglobal[130]=739
pglobal[131]=743
pglobal[132]=751
pglobal[133]=757
pglobal[134]=761
pglobal[135]=769
pglobal[136]=773
pglobal[137]=787
pglobal[138]=797
pglobal[139]=809
pglobal[140]=811
pglobal[141]=821
pglobal[142]=823
pglobal[143]=827
pglobal[144]=829
pglobal[145]=839
pglobal[146]=853
pglobal[147]=857
pglobal[148]=859
pglobal[149]=863
pglobal[150]=877
pglobal[151]=881
pglobal[152]=883
pglobal[153]=887
pglobal[154]=907
pglobal[155]=911
pglobal[156]=919
pglobal[157]=929
pglobal[158]=937
pglobal[159]=941
pglobal[160]=947
pglobal[161]=953
pglobal[162]=967
pglobal[163]=971
pglobal[164]=977
pglobal[165]=983
pglobal[166]=991
pglobal[167]=997

/* absolute value of an integer n */

define a(n){
	if(n>=0)return(n)
	return(-n)
}

/*  the least non-negative remainder when an integer a is divided by a positive
integer b */
/* a%b=m(a,b) if a<=0 or a<0 and b divides a */
/* a%b=m(a,b)-b if a<0, b>0, a not divisible by b */

define m(a,b){
	auto c
	c=a%b
	if(a>=0)return(c)
	if(c==0)return(0)
	return(c+b)	
}

/* min(x,y) */

define n(x,y){
	if(y<x)return(y)
	return(x)
}

/* returns a raised to the power b, where a is an integer
and b is a non-negative integer. */

define z(a,b){
	auto x,y,z
	x=a
	y=b
	z=1
	while(y>0){
		while(y%2==0){
			y=y/2
			x=x*x
		}
		y=y-1
		z=z*x
	}
	return(z)
}

/*  gcd(m,n) for any integers m and n */
/* Euclid's division algorithm is used. */

define g(m,n){
	auto a,b,c
	a=a(m)
	if(n==0)return(a)
	b=a(n)
	c=a%b
	while(c>0){
		a=b
		b=c
		c=a%b
	}
	return(b)
}   

/* the Brent-Pollard method for returning a proper factor of a composite n 
see R. Brent, Bit 20, 176-184 */

define b(n){
	auto a,y,r,g,x,i,k,f,z
	"BRENT-POLLARD: SEARCHING FOR A PROPER FACTOR OF ";n
	a=1;y=0;r=1;g=1;prod=1
	while(g==1){
		x=y
		for(i=1;i<=r;i++)y=(y*y+a)%n
		k=0
		f=0
		while(f==0){
			y=(y*y+a)%n
			k=k+1
			z=a(x-y)
			prod=prod*z
			if(k%10==0){
				g=g(prod,n)
				prod=1
				if (g>1)f=1
			}
			if(k>=r)f=1
		}
		r=2*r;print "r=",r,"\n"
		if(g==n || r == 16384){
			g=1
			r=1
			y=0
			a=a+1
			print "increasing a\n"
			if(a==5){
				print "Brent-Pollard failed to find a factor\n"	
				return(0)
			}
		}
	}
	"--
"
	"FINISHED ";g
	"is a proper factor of ";n
	"--
"
	return(g)
}

define pollard(n){
	auto i,p,t,b
	b=t=2
	p=1
	for(i=2;i<=10^4;i++){
			if(i%10==0){"i=";i}
			t=e(t,i,n)/* now t=b ^(i!) */
			p=g(n,t-1)
			if(p>1){
				if(p<n){
					print "With b = ",b,", i = ",i,", n = ",n,",\n"
					print "gcd(b^i!-1,n)=",p, " is a proper factor of ",n,"\n"
					return(p)
				}
			}
			if(p==n){
			i=2;t=b=b+1
			"b= ";b
			if(b==5) {
			print "Pollard p-1 failed to deliver a factor of ",n,"\n"
				return(0)}
			}
	}
			print "Pollard p-1 failed to deliver a factor of ",n,"\n"
                         return(0)
}
/* returns n/m, where m is the part of n composed of primes 
 qglobal[i],...,qglobal[jglobal-1] < 1000, where i is the value of the global
 variable jglobal before d(n) is called.  qglobal[] and exponent array
 kglobal[] are global variables. */

define d(n){
auto a,u,x,k,i
	a=167
	x=n
	i=0
	while(i<= a){
		k=0
		if(x<pglobal[i])break
		if(x%pglobal[i]==0){
			while(x%pglobal[i]==0){
				k=k+1
				x=x/pglobal[i]
			}
			pglobal[i]
			"is a prime factor of ";n
			qglobal[jglobal]=pglobal[i]
			"exponent ";k
			kglobal[jglobal]=k
			jglobal=jglobal+1
			"--
	"
		}
		i=i+1
	}
return(x)
}

/* a^b (mod c) */
/* a,b,c integers, a,b>=0,c>=1 */

define e(a,b,c){
	auto x,y,z
	x=a
	y=b
	z=1
	while(y>0){
		while(y%2==0){
			y=y/2
			x=(x*x)%c
		}
		y=y-1
		z=(z*x)%c
	}
	return(z)
}

/*
 n>1 and odd and n does not divide b, b > 0, n-1=2^s*t, t odd.
 If r(n,b)=1, then n passes Miller's test for base b and n is either
 a prime or a strong pseudoprime to base b.
 If r(n,b)=0, then n is composite.
 */

define r(n,b){
	auto a,q,i,j
	if(n%2==0){
	"n is even
";return
}
	if(b%n==0){
	"b divides n
";return
	}
	q=(n-1)/2
	a=q
	i=0
	while(a%2==0){a=a/2;i=i+1}
	/* a = t, i=s-1 */
	b=e(b,a,n) /* "b" = b^t(mod n) */
	if(b==1)return(1)
	/* b^t != 1 (mod n) */
	j=0
	while(1){
		/* "b=b^{(2^j)*t}(mod n) */
		if(b==n-1)return(1)
		j=j+1
		if(i < j)return(0)/* j = s */
		b=(b*b)%n
	}
}

/* n>1 is distinct from pglobal[0],...,pglobal[15]. */
/* if q(n)=1, then n passes Miller's test for bases pglobal[0],...,pglobal[4]
   and is likely to be prime. If q(n)=0, then n is composite. */

define q(n){
	auto i
	"MILLER'S TEST IN PROGRESS FOR ";n
	"--
"
	for(i=0;i<=15;i++){
		if(r(n,pglobal[i])==0){
		"FINISHED: ";n
		"is composite
"
		"--
"
		return(0)
		}
	}
	"FINISHED: ";n
	"passes Miller's test 
"
	"--
"
	return(1)
}

/* n>1 is not divisible by pglobal[0],...,pglobal[167]. */
/* f(n) returns a factor of n which is < 1,000,000 (and hence prime)
 or which passes Miller's test for bases pglobal[0],...,pglobal[15] 
 and is therefore likely to be prime. */

define f(n){
	auto f,x,y,b
	b=1000
	if(n<b*b){
		return(n)
	}
	if(q(n)==1)return(n)
	f=1
	x=n
	while(f!=0){
		y=b(x)
		if(y==0){
			y=pollard(x)
			if(y==0){
				print "Exiting factoring program\n"
				halt
			}
		}
		x=n(x/y,y)
		if(x<b*b)return(x)
		if(q(x)==1)f=0
	}
	return(x)
}

/* A quasi-prime (q-prime) factor of n is > 1,000,000, passes Miller's test 
and is not divisible by pglobal[0],...,pglobal[167]. It is likely to be prime.
h(n) returns the number lglobal-t of q-prime factors of n, (rglobal[] and
lglobal are global variables.
The prime and q-prime factors of n are qglobal[i],...,qglobal[jglobal-1]
with exponents kglobal[i],...,kglobal[jglobal-1], where i is the value of the
global variable jglobal before h(n) is called.
 */

define h(n){
	auto b,p,x,k,t
	b=1000 
	"FACTORIZING ";n
	"--
"
	x=d(n)
	t=lglobal
	while(x!=1){
		k=0
		p=f(x)
		while(x%p==0){
			k=k+1
			x=x/p
		}
		if(p<b*b){
			p
			"is a prime factor of ";n
		}
		if(p>b*b){
			p
			"is a q-prime factor of ";n
			rglobal[lglobal]=p
			lglobal=lglobal+1
		}
		qglobal[jglobal]=p
		"exponent ";k
		kglobal[jglobal]=k
		jglobal=jglobal+1
		"--
"
	}
	"FINISHED FACTORIZATION INTO PRIMES AND Q-PRIMES FOR ";n
	"--
"
	return(lglobal-t)
}

/* Selfridge's test for primality - see "Prime Numbers and Computer
Methods for Factorization" by H. Riesel, Theorem 4.4, p.106. */

/* input: n (q-prime) */
/* first finds the prime and q-prime factors of n-1. If no q-prime factor is
present and 1 is returned, then n is prime. However if at least one q-prime is
present and 1 is returned, then n retains "q-prime" status.
If 0 is returned, then either n or one of the q-prime factors of n-1 is 
composite. */

define s(n){
	auto i,x,s,t,u
	"SELFRIDGE'S TEST FOR PRIMALITY OF QUASI-PRIME ";n
	"--
"
	i=jglobal
	u=h(n-1)
	cglobal=u+cglobal
	/* cglobal,jglobal,lglobal are global variables. */
/* h(n) returns jglobal-i primes and q-primes 
qglobal[i],...,qglobal[jglobal-1] */
/* and q-primes rglobal[lglobal-u],...,rglobal[lglobal-1], where u>=0. */
	"SELFRIDGE'S EXPONENT TEST IN PROGRESS FOR Q-PRIME ";n
	"--
"
	while(i<=jglobal-1){
		x=2
		while(x<n){
			s=e(x,n-1,n)
			if(s!=1){
				"SELFRIDGE'S TEST FINISHED: q-prime ";n
				"is composite
"
				return(0)
			}
			t=e(x,(n-1)/qglobal[i],n)
			if(t!=1){
				i=i+1
				break
			}
		x=x+1
		}
		if(x==n){ /* unlikely to be tested */
			"SELFRIDGE'S TEST INCONCLUSIVE
"
				return(0)
		}
	}
	if(u==0){
		"SELFRIDGE' TEST FINISHED: q-prime ";n
		"is prime
"
		"--
"
	return(1)
	}
	if(u>0){
		"FINISHED: q-prime n = ";n
		"retains its q-prime status.
"
		"--
"
	return(1)
	}
}

/* a factorization of n into prime and q-prime factors is first obtained.
Selfridge's primality test is then applied to any q-prime factors; the test 
is applied repeatedly until either a q-prime is found to be composite
(in which case the initial factorization is doubtful and an extra base should be
used in Miller's test) or we run out of q-primes, in which case every q-prime 
factor of n is certified as prime.*/
/* the number of distinct prime factors of n is returned. */

define factor(n){
	auto i,s,t
	"--
"
	jglobal=lglobal=0
	cglobal=h(n)
	/* cglobal gives a progress count of the number of q-prime factors */
	/* yet to be tested by Selfridges' method */
	s=jglobal
	/* s is the number of prime and q-prime factors of n */
	if(cglobal==0){
		"NO Q-PRIMES, FACTORIZATION VALID:
"
		"
"
		"prime factors: 
"
		for(i=0;i<s;i++) qglobal[i]
		"
"
		"exponents:
"
		for(i=0;i<s;i++) kglobal[i]
		"
"
		print "number of prime factors of ",n, " is ";return(s)
	}
	"TESTING Q-PRIMES 
"
	"--
"
	while(cglobal>0){
		t=s(rglobal[i])
		if(t==0){
			"do factor(n) again  with an extra base in q(n) 
"
			return(0)
		}
		i=i+1
		cglobal=cglobal-1
	}
		"ALL Q-PRIMES ARE PRIMES: FACTORIZATION VALID: 
"
		"
"
		"prime factors: 
"
		for(i=0;i<s;i++) qglobal[i]
		"
"
		"exponents:
"
		for(i=0;i<s;i++) kglobal[i]
		"
"
		print "number of prime factors of ",n," is ";return(s)
}