alternative factor(1) - SOURCE

[George Sicherman]

/*
 *	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);
}