[net.sources] alternative factor

colonel@sunybcs.UUCP (George Sicherman) (12/08/83)

/*
 *	factor - factor an integer into primes.  VAX version.
 *	the Colonel. 82/09/09.
 *
 *	This program is an imitation of factor(1) in
 *	UNIX V. 7.
 */

#include <stdio.h>
#include <signal.h>
#include <ctype.h>

int pflag = 0;
int grind = 0;
long int a[2], b[2];
int blog; /* NUMBER OF BITS IN DIVISOR */
long int barrier; /* BLOG CHANGE BOUNDARY */
char number[19];
#define MASK 01777777777L
#define BUMP(i) {if (barrier<(i+=y)){barrier<<=1;blog++;} y=z=6-z;}

main(argc,argv)
int argc;
char **argv;
{
	int fudge;
	int factor(), primes();
	int (*func)() = factor;
	void bailout();
	signal(SIGINT,bailout);
	while (--argc) {
		if ('-'==**++argv) switch(*++*argv) {
		case 'p':
			pflag++;
			break;
		default:
			bomb();
		}
		else break;
	}
	if (pflag) {
		func=primes;
		if (argc<2) grind++;
	}
	if (argc) while (argc--) (*func)(*argv++);
	else while (EOF!=(fudge=scanf("%18s",number))) {
		if (fudge) (*func)(number);
		else printf("invalid number\n");
	}
}

factor(s)
char *s;
{
	int y, z;
	long int i, j, divide();
	double sqrt();
	if (getbign(s)) {
		i=3;
		z=4;
		y=2;
		blog=2;
		barrier=4;
		j=((long int)sqrt((double)(1+a[0])))<<14;
		while (a[0]) {
			if (i>j) {
				printbig();
				printf(" \n");
				return;
			}
			if (!(a[1]&1))
			do {
				a[1]>>=1;
				a[1]|=(a[0]&1)<<27;
				a[0]>>=1;
				printf("2 ");
			} 
			while (!(a[1]&1));
			else if (divide(i)) BUMP(i)
			else printf("%ld ",i);
		}
		if (a[1]==1) return;
		if (a[1]<4) {
			printf("%d\n",a[1]);
			return;
		} 
		else while (!(a[1]&1)) {
			a[1]>>=1;
			printf("2 ");
		}
		j=(long int)sqrt((double)a[1]);
		while (a[1]>1) {
			if (i>j) i=a[1];
			if (!(a[1]%i)) {
				a[1]/=i;
				printf("%ld ",i);
			}
			else BUMP(i)
		}
		printf("\n");
	}
	else printf("invalid number\n");
}

getbign(s)
char *s;
{
	a[0]=a[1]=0;
	while (isdigit(*s)) {
		a[0]=(a[0]<<1)+(a[0]<<3);
		if (a[0]&~MASK) return(0);
		a[1]=(a[1]<<1)+(a[1]<<3);
		a[1]+= *s++-'0';
		a[0]+=(017&((a[1]&~MASK)>>28));
		a[1]&=MASK;
	}
	return(1);
}

primes(s)
char *s;
{
	if (getbign(s)) for (;;) {
		if (!a[0] && a[1]<3) {
			printf("2\n");
			a[1]=3;
			continue;
		} else {
			a[1]|=1;
			if (isprime()) {
				printbig();
				putchar('\n');
				if (!grind) break;
			}
			a[1]+=2;
			if (a[1]&~MASK) {
				a[1]&=MASK;
				if (++a[0]&~MASK) break;
			}
		}
	}
	else fprintf(stderr,"primes: bad number %s\n",s);
}

isprime()
/*
 *	contents of a[] should be odd.
 */
{
	int y, z;
	long int i, j, divide();
	double sqrt();
	z=4;
	y=2;
	blog=2;
	barrier=4;
	if (a[0]) {
		j=((long int)sqrt((double)(1+a[0])))<<14;
		for (i=3;;) {
			if (i>j) return 1;
			if (!divide(i)) return 0;
			BUMP(i)
		}
	} 
	else {
		if (a[1]<2) return 0;
		if (a[1]<4) return 1;
		if (~a[1]&1) return 0;
		j=(long int)sqrt((double)a[1]);
		for (i=3;;) {
			if (i>j) return 1;
			if (!(a[1]%i)) return 0;
			BUMP(i)
		}
	}
}

long int divide(dvsr)
long int dvsr;
{
	long int c, d;
	int q;
	b[0]=a[0]/dvsr;
	c=a[0]-b[0]*dvsr;
	b[1]=0;
	d=a[1];
	for (q=blog; q; q--) {
		b[1]<<=1;
		c<<=1;
		c+=d>>27;
		d<<=1;
		d&=MASK;
		if (c>=dvsr) {
			b[1]++;
			c-=dvsr;
		}
	}
	d>>=blog;
	d|=c<<(28-blog);
	c=d/dvsr;
	d=d-c*dvsr;
	b[1]<<=(28-blog);
	b[1]+=c;
	if (!(pflag||d)) for (q=0; q<2; q++) a[q]=b[q];
	return d;
}

printbig()
{
	char *sp;
	long int u[2];
	int q;
	for (q=0; q<2; q++) u[q]=a[q];
	sp = &number[18];
	while (a[0]||a[1]) {
		*sp-- = '0'+divide(10L);
		for (q=0; q<2; q++) a[q]=b[q];
	}
	printf("%s",++sp);
	for (q=0; q<2; q++) a[q]=u[q];
}

bomb()
{
	if (pflag) fprintf(stderr,"usage: primes [ number ... ]\n");
	else fprintf(stderr,"usage: factor [ number ... ]\n");
	exit(1);
}

void bailout()
{
	printf("...\n");
	fclose(stdout);
	exit(-1);
}

colonel@sunybcs.UUCP (George Sicherman) (12/08/83)

.TH FACTOR 1 
.UC 4
.SH NAME
factor, primes \- factor a number, generate large primes
.SH SYNOPSIS
.B factor
[
.B \-p
] [ number ... ]
.PP
.B primes
[ number ... ]
.br
.SH DESCRIPTION
When
.I factor
is invoked without an argument,
it waits for a number to be typed in.
If you type a positive number less than
2\(ua56
(about 7.2E16),
it
prints the number's prime factors,
each the proper number of times.
Then it waits for another number.
.PP
If
.I factor
is invoked with arguments,
it factors them as above and then exits.
.PP
Maximum time to factor
is proportional to the square root of
.I n
and occurs when
.I n
is prime or the square of a prime.
.PP
When
.I primes
is invoked,
it waits for a number to be typed in.
If you type a positive number less than
2\(ua56,
it prints all primes greater than or equal
to the number.
.PP
If
.I primes
is invoked with an argument,
it uses the argument as the number.
With several arguments,
.I primes
prints only the least prime
greater than or equal to each,
and then exits.
.PP
With the
.B -p
option,
.I factor
is equivalent to
.I primes.
.PP
.SH BUGS
This program is a home-brewed imitation
of the corresponding program in Unix Version 7.
It is warranted neither for speed nor correctness.
.SH AUTHOR
Col. Sicherman