/* 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(y0){ 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=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 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(pb*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(x0){ "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;i0){ 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