[net.sources] Prime factorization, sigma function, phi function

spaf@gatech.UUCP (07/26/84)

The following program has been sitting around in my directory
for a while.  I've cleaned it up in case anybody out there
has any use for it.  It can be used to get the unique prime
factorization of any 32 bit unsigned integer.  It can also be
used to obtain the sigma function (sum of all divisors) or phi
function (number of integers less than and relatively prime to).

Enjoy.


: to unbundle, "sh" this file -- DO NOT use csh
:
echo x - Makefile
sed 's/^X//' >Makefile <<'+FUNKY+STUFF+'
X#  Makefile for "genprimes" and "factor"
X# 
X#  	Eugene Spafford	
X# 	School of Information and Computer Science
X# 	Georgia Institute of Technology
X# 	Atlanta, Georgia  30332
X# 
X# 	...gatech!spaf		spaf%gatech@CSnet-Relay
X# 
X# 
X#   Copyright (c) 1984 by Eugene Spafford and the Georgia Tech
X# 		Research Institute
X# 	The author hereby gives permission for anyone to
X# 	copy or use this software provided that it is not 
X# 	sold, traded or otherwise exchanged for profit,
X# 	and that this notice and the author's address both
X# 	appear in any copies.
X# 
X# 
X# 
X#    Modification History:
X# 	07/25/84	EHS	Release coding
X# 
X# 
X
X#  Change the following definition to something reasonable
X#   on your system (like in your local library)
X
XTABLE= \"/usr/lib/primetable\"
X
X
X
Xall:	genprimes factor  clean
X
Xclean:	
X	rm -f *.o a.out
X
X
Xgenprimes: factor.h genprimes.c
X	cc -O genprimes.c -DPTAB=$(TABLE) 
X	a.out
X
Xlint:	
X	lint -abchp genprimes.c
X	lint -abchp factor.c
X
Xfactor:	factor.h factor.c
X	cc -O -DPTAB=$(TABLE) factor.c -lm -o factor
+FUNKY+STUFF+
lf -l Makefile
echo x - factor.1
sed 's/^X//' >factor.1 <<'+FUNKY+STUFF+'
X.TH FACTOR 1
X.SH NAME
Xfactor \- find prime factorization, phi, and sigma values
X.SH SYNOPSIS
X.B factor
X[
X.B \-v
X] [
X.B \-s
X] [
X.B \-p
X]
Xpositive integers....
X.SH DESCRIPTION
X.I Factor
Xuses a table of pre-calculated prime numbers to reduce the input
Xto a unique prime factorization.  If the user requests only that
Xfactorization, then that factorization is printed.  Otherwise, it
Xis used to calculate the sigma function (sum of all divisors), or
Xthe phi function (all numbers less than the input which are relatively
Xprime to the input).
X
XOf interest to number theorists, the sigma function is a definition
Xof perfect numbers.  Any input whose sigma function value is equal
Xto twice the input is perfect.  For example, sigma(6) = 12 and
X6 is the first perfect number.
X
XThe phi function is also of interest for various reasons.  One
Xinteresting feature is that phi(p) where p is prime is precisely
Xequal to p-1 iff p is prime.
X
XThe program will accept any valid 32 bit, unsigned integer. Output is
Xdefined for all valid inputs (i.e., the output does not have to
Xfit into a 32 bit unsigned integer).
X.PP
XOptions:
X.TP 6
X.B \-v
XLabel each output line with the input value for which the result is
Xbeing computed.
X.TP 6
X.B \-s
XCalculate the sigma function value for all of the input values.
XMutually exclusive with the "-p" option.
X.TP 6
X.B  \-p
XCalculate the phi function for the input values.  Mutually exclusive
Xwith the "-s" option.
X.SH "SEE ALSO"
Xgenprimes(8)
X.SH DIAGNOSTICS
XUnknown options will cause the program to end after printing a usage
Xmessage to stderr.  Improper format of input (non-numeric characters
Xor improper use of signs) will cause an error message to stderr with
Xthe offending input, but processing of other values will continue.
XValues which are too large for 32 bits will also be flagged as
Xnon-fatal errors.
X.SH FILES
X/usr/lib/primetable
X.SH AUTHOR
XEugene Spafford (...gatech!spaf)
X.SH BUGS
XIt requires a significant amount of time to read in the table of
Xprime numbers.  Currently, the program assumes the whole read
Xcan be done in one block -- if not, the code needs modification.
+FUNKY+STUFF+
lf -l factor.1
echo x - factor.c
sed 's/^X//' >factor.c <<'+FUNKY+STUFF+'
X/* factor --- factor integers, and calculate phi and sigma functions
X *
X * 	Eugene Spafford	
X *	School of Information and Computer Science
X *	Georgia Institute of Technology
X *	Atlanta, Georgia  30332
X *
X *	...gatech!spaf		spaf%gatech@CSnet-Relay
X *
X *
X *  Copyright (c) 1984 by Eugene Spafford and the Georgia Tech
X *		Research Institute
X *	The author hereby gives permission for anyone to
X *	copy or use this software provided that it is not 
X *	sold, traded or otherwise exchanged for profit,
X *	and that this notice and the author's address and name
X *	appear in any copies.
X *
X *
X *   Description:
X *	This program allows the user to find the unique prime
X *	factorization of a 32 bit integer.  The user may also
X *	find the appropriate value of the Euler phi function
X *	(the number of relatively prime values less than the
X *	number), and the sigma function (the sum of all factors
X *	of the number).
X *
X *	factor <number> <number> <number>
X *	will give the prime factorizations, one per line.
X *
X *	factor -s <number> <number> <number>
X *	gives the sigma function value
X *
X *	factor -p <number> <number> <number>
X *	gives the phi function value.
X *
X *	The -v flag is provided to label output lines with the
X *	value being operated upon.  This may be used alone or
X *	combined with other flags, as in:
X *
X *	factor -s -v
X *	factor -pv
X *	factor -v
X *
X *	The "genprimes" program must first have been run, and the table
X *	of prime numbers built for use by this program.
X *
X *
X *   Modification History:
X *	07/25/84	EHS	Release coding
X *
X */
X
X#include <stdio.h>
X#include <ctype.h>
X#include "factor.h"
X
X#define ONE (double) 1.0
X
X
Xmain(argc, argv)
Xint argc;
Xchar **argv;
X{
X	unsigned value;
X	int fd, flag = 0, verify = 0, typen, i;
X	int numw = TABSIZ*sizeof(int);
X	extern int read(), open(), check(), atou(), strcmp();
X
X
X	if (argc < 2) 		/* No arguments results in a usage message */
X		error("", "", 1);
X
X	/* Scan for valid options.
X	 *
X	 *   -p  causes calculation of the phi function for arguments
X	 *   -s  causes calculation of the sigma function for arguments
X	 *   -v  causes verbose mode; each line is tagged with the argument
X	 *
X	 *   if neither -s nor -p are given, the unique prime factorization
X	 *   is printed.
X	 */
X	for (i = 1; --argc && argv[i][0] == '-'; i++) {
X		int j;
X		for (j = 1; argv[i][j]; j++)
X			switch (argv[i][j]) {
X			case 's':
X				flag = 1;
X				break;
X			case 'p':
X				flag = 2;
X				break;
X			case 'v':
X				verify = 1;
X				break;
X			default:
X				error("Unknown option: %s\n", argv[i], 1);
X			}
X	}
X
X	/*  Read in our table of prime numbers */
X
X	if ((fd = open(PTAB, 0)) < 0) 
X		error("Trying to open %s\n", PTAB, 2);
X	else if (read(fd, divisors, numw) != numw) 
X		error("Trying to read %s\n", PTAB, 2);
X	close(fd);
X
X	/* Now, for each number provided, calculate our values.
X	 *  Valid inputs can contain leading blanks or space, or one
X	 *  leading + sign.  The values must fit into an unsigned 32
X	 *  bit integer variable.  Other embedded characters are not 
X	 *  allowed, and neither are trailing spaces.
X	 */
X	for (; argc--; i++) {
X		typen = atou(argv[i], &value);
X		if (typen == -1)
X			error("Not a valid positive integer %s\n", argv[i], 0);
X		else if (typen == -2)
X			error("Out-of-range: %s\n", argv[i], 0);
X		else if (typen == -3)
X			error("Cannot evaluate negative value: %s\n", argv[i], 0);
X		else if (value < 2)
X			error("not defined for %s\n", argv[i], 0);
X		else {
X			if (verify)
X				printf("%u: ", value);
X			check(value, flag);
X			printf ("\n");
X		}
X	}
X}
X
X
Xint check(try, flag)		/* routine which does the dirty work */
Xregister unsigned try;
Xint flag;
X{
X	register int iy, ix;
X	register unsigned partial;
X	double factors[FACTAB], powers[FACTAB];
X	double sigma = ONE;
X	extern double pow();
X
X
X	if (flag) 
X		for (ix = 0; ix < FACTAB; ix++)
X			factors[ix] = powers[ix] = 0.0;
X
X	iy = -1, ix = 0;
X	do {
X		partial = try/divisors[ix];
X		if (partial < divisors[ix]) 
X			break;
X		else if (try == partial*divisors[ix]) {
X			try = partial;
X			if (NOT flag)
X				printf ("%u   ", divisors[ix]);
X			else {
X				if (iy < 0 || (int) factors[iy] != divisors[ix])
X					factors[++iy] = (double) divisors[ix];
X				powers[iy]++;
X			}
X		}
X		else
X			ix++;
X	} while (ix < TABSIZ);
X
X	if (NOT flag) {
X		if (try > 1)
X			printf("%u", try);
X		return;
X	}
X
X	if (try > 1) {
X		if (iy < 0 || (int) factors[iy] != try)
X			factors[++iy] = (double) try;
X		powers[iy]++;
X	}
X	for (sigma = ONE, ix = 0; ix <= iy; ix++)  
X		sigma *= (flag == 1 ?
X		(pow(factors[ix], ++powers[ix])-ONE)/(factors[ix]-ONE) :
X		pow(factors[ix], --powers[ix])*(factors[ix]-ONE));
X
X	printf("%.0f", sigma);
X}
X
X
X
Xint atou(p, value)	/* decode arguments into unsigned integers */
Xregister char *p;
Xregister unsigned *value;
X{
X	register unsigned new;
X	int sign = 0;
X
X
X	*value = 0;
X	while (!isdigit(*p)) {
X		switch(*p++) {
X		case ' ':
X		case '\t':
X			if (sign)
X				return -1;
X			continue;
X		case '-':
X			return -3;	/* we don't allow negative values */
X		case '+':
X			if (sign++) 	/* we don't allow multiple signs */
X				return -1;
X			continue;
X		default:
X			return -1;
X		}
X	}
X
X	while(isdigit(*p)) {
X		new = *value*10 + *p++ - '0';
X		if (new/10 != *value) 		/* Check for overflow */
X			return -2;
X		else
X			*value = new;
X	}
X	return (*p) ? -1 : 1;	/* Anything other than EOS is an error */
X}
X
X
Xerror(format, arg, die)		/* Standard error message formatting */
Xchar *format, *arg;
Xint die;
X{
X	fprintf(stderr, format, arg);
X	if (die == 1)
X		fprintf(stderr, "usage: factor [-p] [-s] [-v] <numbers>\n");
X	if (die)
X		exit(1);
X}
+FUNKY+STUFF+
lf -l factor.c
echo x - factor.h
sed 's/^X//' >factor.h <<'+FUNKY+STUFF+'
X/* factor.h --- Common definitions for "genprimes" and "factor"
X *
X * 	Eugene Spafford	
X *	School of Information and Computer Science
X *	Georgia Institute of Technology
X *	Atlanta, Georgia  30332
X *
X *	...gatech!spaf		spaf%gatech@CSnet-Relay
X *
X *
X *  Copyright (c) 1984 by Eugene Spafford and the Georgia Tech
X *		Research Institute
X *	The author hereby gives permission for anyone to
X *	copy or use this software provided that it is not 
X *	sold, traded or otherwise exchanged for profit,
X *	and that this notice and the author's address and name
X *	appear in any copies.
X *
X *
X *   Modification History:
X *	07/25/84	EHS	Release coding
X *
X */
X
X/*  The following number is carefully chosen to minimze table size.
X *  Any integer that will fit in 32 unsigned bits can be factored
X *  using a sieve of the first 6543 prime numbers.  Trust me.
X */ 
X
X#define TABSIZ 6543
X
X/*  The next value delimits the number of unique prime factors for
X *  any unsigned 32-bit value.  You can prove it by constructing
X *  the example case of 2*3*5*7*11*13*17*19*23.  If you multiply
X *  that by 29, you exceed 32 bits of precision.
X */
X
X#define FACTAB 9
X
X/*  Next, we have a wired in pathname for where to place your
X *  table or prime numbers.  You can obviously change this with
X *  the command line D= option to cc
X */ 
X
X#ifndef PTAB
X#define PTAB "/u/clouds/spaf/stuff/primetable"
X#endif
X
X/*  Where we store the prime number themselves...
X */
X
Xunsigned divisors[TABSIZ];
X
X#define NOT !
X
+FUNKY+STUFF+
lf -l factor.h
echo x - genprimes.1
sed 's/^X//' >genprimes.1 <<'+FUNKY+STUFF+'
X.TH GENPRIMES 8
X.SH NAME
Xgenprimes \- generate a table of prime numbers
X.SH SYNOPSIS
X.B genprimes
X.SH DESCRIPTION
X.I Genprimes
Xbuilds an array of the first 6543 prime numbers and then writes them
Xto a library file.  The library file may then be used by other
Xprogram which might need those prime numbers.  The first 6543
Xprime numbers are sufficient to factor any unsigned 32 bit quantity.
X
XThe algorithm involved is a variation of the classic sieve method.
X.SH "SEE ALSO"
Xfactor(1)
X.SH DIAGNOSTICS
XThe program complains if it cannot open the destination file or if
Xit cannot write the whole array with a single write.
X.SH FILES
X/usr/lib/primetable
X.SH AUTHOR
XEugene Spafford (...gatech!spaf)
X.SH BUGS
XThe algorithm could be more efficient, but it is not expected
Xthat the user will run this program often.
+FUNKY+STUFF+
lf -l genprimes.1
echo x - genprimes.c
sed 's/^X//' >genprimes.c <<'+FUNKY+STUFF+'
X/* genprimes --- generate the table of primes needed by "factor"
X *
X * 	Eugene Spafford	
X *	School of Information and Computer Science
X *	Georgia Institute of Technology
X *	Atlanta, Georgia  30332
X *
X *	...gatech!spaf		spaf%gatech@CSnet-Relay
X *
X *
X *  Copyright (c) 1984 by Eugene Spafford and the Georgia Tech
X *		Research Institute
X *	The author hereby gives permission for anyone to
X *	copy or use this software provided that it is not 
X *	sold, traded or otherwise exchanged for profit,
X *	and that this notice and the author's address and name
X *	appear in any copies.
X *
X *
X *   Description:
X *	This program simply builds a large array of the first
X *	TABSIZ prime numbers and then writes them to a library
X *	file.  The algorithm is a relatively unsoophisitcated
X *	sieve.
X *
X *   Modification History:
X *	07/25/84	EHS	Release coding
X *
X */
X
X#include <stdio.h>
X#include "factor.h"
X
Xint index = 1;
X
X
Xmain()
X{
X	int fd;
X	unsigned try = 1;
X	int numw = TABSIZ*sizeof(int);
X
X
X	if ((fd = creat(PTAB, 0644)) < 0) 
X		perror("Creating PTAB");
X
X	divisors[0] = 2;
X	while (index < TABSIZ) {
X		try += 2;		/* No sense trying even numbers! */
X		if (check(try)) 
X			divisors[index++] = try;
X	}
X
X	if (write(fd, divisors, numw) != numw)
X		perror ("Writing PTAB");
X	close(fd);
X	exit(0);
X}
X
X
Xint check(try)
Xregister unsigned try;
X{
X	register int ix, partial;
X
X
X	for (ix = 1; ix < index; ix++)  {
X		partial = try/divisors[ix];
X		if (partial < divisors[ix])
X			return 1;		/* Faster than a sqrt test */
X		else if (try == divisors[ix]*partial)
X			return 0 ;
X	}
X	return 1;
X}
+FUNKY+STUFF+
lf -l genprimes.c
echo x - mkpath.1
sed 's/^X//' >mkpath.1 <<'+FUNKY+STUFF+'
X.TH MKPATH 1 
X.SH NAME
Xmkpath \- make a Usenet path data base
X.SH SYNOPSIS
X.B mkpath
X[
X.B \-2
X] [
X.B \-m
X] [
X.B \-lsite
X] [
X.B \-dsite
X] [
X.B \-dsite1!site2
X] [
X.B \-asite1!site2
X]
Xmapfiles...
X.SH DESCRIPTION
X.I Mkpath
Xcreates a list of USENET sites and the UUCP routing path to each site.
XThe input is a list of USENET sites and what other sites each site is
Xconnected to.  The input data should look like:
X.nf
X
X	Name: (name of site)
X	(any other fields)
X	News: (site names separated by blanks, commas, or newlines)
X	(more site names)
X	Mail: (still more site names)
X	(even more site names)
X	(any other fields)
X
X.fi
XMark and Karen Horton just happen to send Usenet connection information
Xin this format once a month to the news group net.news.map.
XContinuation lines must start with either a blank or tab.
XThe output is sent to stdout, and looks like:
X.nf
X
X	(name of site)	(printf string of how to get there)
X
X.fi
XThe routes to all possible sites are listed.
X
X.PP
XOptions:
X.TP 6
X.B \-2
XInterpret all connections as two-way links.
X.TP 6
X.B \-m
XUse the Mail: field in the connection information.  Without this option,
Xonly the News: field is used.
X.TP 6
X.B  \-lsite
XSet the originating system name to
X.I site.
XIf this option is not used,
X.I mkpath
Xgets the system name from one of several places, depending upon
Xcompile-time options.
X.TP 6
X.B  \-dsite
XDelete all references to
X.I site
Xbefore building the UUCP paths.
X.TP 6
X.B \-dsite1!site2
XDelete the connection from
X.I site1
Xto
X.I site2
Xbefore building the UUCP paths.
X.TP 6
X.B \-asite1!site2
XAdd the connection from
X.I site1
Xto
X.I site2
Xbefore building the UUCP paths.
X.dt
X.SH "SEE ALSO"
Xnmail(1)
X.SH DIAGNOSTICS
XExit code 1 is returned for an invalid command line, 2 for
Xan inaccessible map input, 3 for starting system not in
Xthe input data, and 4 for not enough memory.
X.SH AUTHOR
XMike Mitchell  (duke!mcnc!ncsu!mcm)
X.SH ACKNOWLEDGEMENTS
XMark and Karen Horton	(cbosgd!mark)
X.br
XJohn Carl Zeigler	(duke!mcnc!ncsu!jcz)
XBruce Parker		(psuvax!parker)
X.SH BUGS
X.I Mkpath assumes every site is connected via UUCP and exchanges mail.
+FUNKY+STUFF+
lf -l mkpath.1
-- 
Off the Wall of Gene Spafford
The Clouds Project, School of ICS, Georgia Tech, Atlanta GA 30332
Phone:	(404) 894-6169, (404) 894-6170 [messages]
CSNet:	Spaf @ GATech		ARPA:	Spaf%GATech.CSNet @ CSNet-Relay.ARPA
uucp:	...!{akgua,allegra,hplabs,ihnp4,masscomp,ut-ngp}!gatech!spaf
	...!{rlgvax,sb1,uf-cgrl,unmvax,ut-sally}!gatech!spaf