#! /bin/sh
# From factor@src.dec.com Sun Jun  6 11:27:05 1993
# This is a shell archive.  If you save it in file 'foo', unpack it
# by typing 'sh foo'
chmod 755 .
echo Extracting README
cat >README <<'!End!of!file!'
To use this package on Unix, read at least the short instructions.  To use
it on VMS, you'll also need to request VMSprogram#, just the same way you
retrieved this.  You can retrieve some user-contributed programs (a C
version of runfactor, a System V renice, and faster multiplication
routines for Sun-3's) by asking for contrib#.

THE SHORT INSTRUCTIONS:

If you just want to run factoring on a single kind of machine with a single
CPU per machine, and a Berkeley-esque version of Unix (such as SunOS or
Ultrix), you should be able to do the following:

1. Create a directory, and unpack this archive into it by saving it in
a file, say 'factor.shar', and typing 'sh factor.shar'.  If you're
reading this, you've probably done that already.

2. Send mail to factor@src.dec.com containing two lines:

	path# a-mail-path-from-SRC-back-to-you
	task# 10

Replace 10 by 10*the number of machines you will run factoring on.  The path#
and task# lines should be flush against the left margin in order for the
automatic mail reading program to recognize them.

3. Type 'mv mpqssrc/* .; rmdir mpqssrc; make' (you might want to change
the Makefile a little before you do this)

4. When you get mail back from the factor service containing lines to
add to mutable/TASKLIST, come back to this directory, and do it.

5. For each machine on which you want to factor, log into that machine, cd
to this directory, and type 'checkrun &'.

6. Once every few days, check that TASKLIST is non-empty, and that
factoring is running (ps ux may be helpful here).  You can look at the
two log files in this directory for extra entertainment value.  You should
get mail telling you that the TASKLIST is empty, if that happens.


A FEW IMPORTANT WARNINGS:

Some mailers like to play funny games with long lines, or with certain
characters in the text.  In particular, the number to be factored may wind
up on two lines instead of one.  This is bad.  Also, in basic.c, there's
an up-arrow (^) that sometimes gets turned into a tilde (~).  If they
looked the same to you in the line that just went by, they're probably
wrong below, too; there's a comment in the lines just below the problem
area, so you'll know what to do; we want the xor-assign operator in C, not
the complement-assign operator.


THE LONG INSTRUCTIONS:

First, create a directory, and unpack the shell archive in that directory.

Next, decide if you plan to run factoring on a single architecture, or on
multiple architectures, sharing a single TEMPLATE, RESTART (usually about
a megabyte, so sharing it is worthwhile) and TASKLIST.  If a single
architecture is sufficient, type 
'mv mpqssrc/* . ; rmdir mpqssrc; make'.  Otherwise, for each architecture,
create a subdirectory whose name is the name of the architecture (e.g.,
'sun3', 'vax').  In each of these directories, copy the Makefile from
mpqqsrc, and create links (or copy) the files basic.c and mpqs.c to this
directory.  Edit the Makefile for each architecture as desired; we've had
good luck on Sun 3's and on VAXes using the Gnu C compiler 'gcc' instead
of the standard C compiler.  'gcc' probably works on many other machines,
but we haven't tried it.  Make executables for each type of machine by
typing 'make' in the appropriate directory, on the appropriate type of
machine (or by setting things up to use the appropriate cross-compiler).
There's no point in doing this unless the directory you're doing all this
in is remotely visible from all the machines you want to run things on; in
particular, you might put the factoring directory someplace where all the
other machines you want to run factoring on can get to it using NFS.  If
you don't have a remote file system, you can still run factoring on
several machines, but you'll have to set up directories of this sort on
each one, and make individual TASKLISTs, etc. for each one.  This is more
work for you, and more of a headache to keep running, but it will work.

Next, you need to get a range of UIDs to work with.  Send mail to the
factor account here at SRC.  If you're on the internet, you can use
the address 'factor@src.dec.com'.  If you're a UUCP site, you can use
the address 'decwrl!factor'; most major sites talk directly to decwrl,
so you just need to prepend the name of the nearest major site; for
instance, uunet talks to decwrl.  The content of the mail should be
two lines.  The first line should contain the string 'path#' followed
by a mail path that you think might work to get mail from me back to
you.  The second line should contain the string 'task#' followed by
the number of tasks you estimate you will need.  A task is consumed
each time a processor starts helping out; you will need a fresh one
whenever a machine is rebooted.  You can always get more by reexcuting
these instructions.  I recommend using 10 * number of processors as a
starting point.  If all is working, you should get a mail message
telling you what to do to set up your TASKLIST file.  If you don't get
a response within a day or two, try again with a different return address.
If you still don't get a response, send your mail to 'msm@src.dec.com', and
I'll see what I can do.  If nothing else is working, and you really want to
help, you can give me a call (collect) at 1-(415)853-2221.

Once you have a tasklist, and you know how to send mail to the factor account,
you're almost ready to start.  If you can't send mail from the machine you
plan on running on, you should change the 25 on the next-to-last line of
mutable/TEMPLATE to a 0. The results will accumulate in a file on that
machine; see the comments at the start of mpqs.c for details; also read the
logfile.

Now, you're ready to run.  If you were able to use factor@src.dec.com as
the path by which to send mail, all is well.  If not, you need to set the
environment variable MAILTO to a path that worked.  If you need to do this
for all the machines you're using, you might want to edit environ to set
the default value of MAILTO automatically.  Otherwise, set the environment
variable in the appropriate way for your shell.  Then, if you're on a single
machine, you can just type 'checkrun &', and everything ought to work.  You
may have trouble on System V Unix, because of long names for files.  Explore
the shell scripts, and change them, as needed.  You may want to prevent
factoring from running when other things want to use the machine; if so,
read about pauseload and pausetty, below.  Later System V machines changed
the syntax of the ps command; if you need to say ps -p 175 to check on process
175, set the PS environment variable to 'ps -p' before starting.  They also
broke 'ln'; the underlying link(2) call still gets an EEXIST error, but 'touch
foo; touch bar; ln foo bar' doesn't raise an error.  We rely on the atomicity
of ln in the code for locking (see acquire, below); you should make sure
that the environment variable LN is set up correctly to be either /bin/ln
or /etc/link as appropriate.

If you're using multiple machines, but a single architecture, you can
just log into each of them, cd to the appropriate directory, and start
them as in the single machine case.  You should start one of them
first, and wait until it has created a RESTART file before starting
the others.  This will work, unless hostname doesn't distinguish
between the machines.

If you're using multiple architectures, you can do the same thing, but
you should set FACTORPATH to 'vax:.' first, assuming that you're
running it on a vax and you named the subdirectory with vax
exectuables 'vax'.  What setting a variable consists of depends on whether
you're using the C-shell or the Bourne shell.  If you're using the C-shell,
type 'setenv FACTORPATH .:vax'.  If you're using the Bourne shell, type
'FACTORPATH=.:vax; export FACTORPATH'.  They mean the same thing.

One thing that we do here, is that we have a command /bin/machine that prints
the name of the architecture of the machine on which it's run.  If you have
something like that, you might want to change the default value of FACTORPATH
to `/bin/machine`:.

In any case, you can cause the log files to be placed in a different directory
than the factoring programs.  You will need to set FACTORPATH to include
the names of the directories containing the shell scripts, and the
architecture-specific executable.  In my case, that's
'/usr/msm/factordir/vax:/usr/msm/factordir'.  You will also need to set
FACTORDIR; in my case I set it to '/usr/msm/factordir/mutable'.  That refers
to the directory containing RESTART, TEMPLATE, TASKLIST, and a few other
goodies.  Then I create the directory where I want the logs and lockfiles,
cd to it, and type $FACTORDIR/../checkrun&

You might find the log files entertaining.  I don't, but some people
are strange.  You can truncate them periodically; there's nothing much
in there worth saving.  They just help convince you that the program
is doing something other than infinite looping.  You should check on
TASKLIST periodically, to make sure that the sum of the numbers in the
second column isn't really small.  You should probably make sure that
it stays above twice the number of participating machines; send mail
as above to get more tasks.  If you don't want logfiles, set 'control' and
'log' to /dev/null in your environment.


MAKING FACTORING LESS (OR MORE) OF A HOG

By default, factoring tries to stay out of your way: it will run the
factoring program at close to the lowest priority, and will also pause the
program whenever it seems that your machine is in use.

It tries to pause when you're typing:  If you set pausetty to some positive
integer (before starting factoring) then whenever you start using your machine
factoring will notice within (by default) five minutes, and wait until it
discovers that your machine has been idle for pausetty before resuming.  On
some workstations (like Suns) you'll want to set ttydevs to include /dev/mouse
(or /dev/win*, in some versions) so that factoring will pause even if the
user is just using the mouse.  You might want to add /dev/kbd, too.  If you're
using X, and neither /dev/mouse nor /dev/kbd works, you might try using xtime
(included) to fetch the X server's idea of the last event time.  Change the
newsum line in runfactor to:

      newsum="`ls -lu $ttydevs | sum` `xtime`"

and do the same to the oldsum line.  You'll need to make sure that the
DISPLAY environment variable is correctly set.  If you use this, be
careful.  It detects idleness by a probe to the X server which may change
the very thing that it's trying to detect.  This doesn't happen with
current DEC X servers, but the R4 MIT X server does do it.  To see
if you can use this, see if the command "xtime; xtime" produces the same
output on both lines (and don't move the mouse or type anything while it's
running, or else it certainly won't).  If so, you can use this.  If not,
you should hope that /dev/kbd or /dev/mouse works.  You can check those by
typing "ls -lu /dev/kbd", and seeing if the output changes after you type
some more, and similarly for the mouse.

You can make it pause when your machine is busy: you can make factoring pause
whenever the load average exceeds $pauseload (excluding the contribution
of factoring).  It will wait until the load average drops below this
threshhold before resuming.

If you want to pause factoring for a while, you can go to the directory
where your log files are, and run pausefactor from the distribution
directory.  By default it pauses for an hour, but you can tell it how many
seconds to pause, if you want.

If you want to be more aggressive, you can go to the directory where the
logfiles are and touch a file shortkill.machineName, mediumkill.machineName, or
longkill.machineName, where machineName is what running hostname prints out.
This will actually kill the factoring process (but not the shell scripts)
for about 4, 6, or 8 hours, respectively.  If you're bothered by factoring
during your workday, you might try putting an entry in cron, or running
an at job to touch the lonkill file at just about the time you get in to
work.

If you decide to give up on factoring completely, you can create a
file called stop.MyMachineName in the same directory as the log files.
This will permanently terminate all factoring activity on your
machine, if your machine is named MyMachineName.  That string should
actually match the one that follows log in the name of the logfile.

If you want to change the pause parameters, you can kill the runfactor
job on your machine, and start a new one with new values.  If your
machine doesn't crash for the next few hours, this will work fine.  If
it does, the wrong at file may cause the old values to be restored.
You can fix this by using 'at -r' to get rid of the old restart scripts
before running the new copy of checkrun.

If you have a system with user-level cron files, you might want to execute
checkrun from that, and simplify ticktock.  It will make things more
reliable, too.

If you have a multiprocessor, treat it as several different machines.  You'll
need to set MACH to some distinct names for each copy (such as hydra.1,
hydra.2, etc.)  You'll need to increase the pausefactor value enough to
account for the number of copies you're running, or else turn it off.

IF THINGS DON'T WORK

First, check that the following all exist, preferably in the indicated
directories:

awk, cat, chmod, date, echo, expr, hostname, kill, ln, ls, mv, ps, rm,
sed, sh, test (usually found in /bin)

at, sleep, sum, touch (usually found in /usr/bin)

uptime (usually found in /usr/ucb)

If any of these are in directories other than the ones listed, or in
other than your normal search path, you will have to add the
directories in which they can be found to the shell scripts in this
directory or to the variable FACTORPATH, decsribed below.  You can
survive if hostname or uptime is missing; the others will be
a problem.

If hostname is missing, you will have to explicitly set the shell
variable MACH, which is used to distinguish different machines sharing
the same directory.  You can explicitly set it to any reasonably short
string.

If uptime is missing, and cannot be supplied, you should not set
pauseload, described above.

If at doesn't work on your machine, you will need to do a little more work.
If you can change the crontab on the machines you will be factoring on, you
can simply make an entry that periodically (once every hour will be fine)
invokes a shell script that sets up some environment variables and the uid
correctly, and then invokes checkrun, just like the at script would do.  You
can then simplify the ticktock script to just sleep for $TICKSEC.  If you
can't make at work, and you can't change crontab, you will just have to
periodically make sure that factoring is running yourself, and restart it
if it isn't.  You should change ticktock to just sleep for $TICKSEC in that
case as well.  If you use crontab to run factoring, note that the
factoring process will not have a control terminal, and you will have
to change the environ file to use kill -STOP instead of kill -TSTP; this
can be a problem in general under SunOS.

All of the other executables are fundamental to the operation of these
scripts; fortunately, they exist in every version of Unix since
version 7.

If all else fails, send me mail at msm@src.dec.com, or call me collect
at 1-(415)853-2221.

Mark
!End!of!file!
echo Extracting acquire
cat >acquire <<'!End!of!file!'
#! /bin/sh

# A way to acquire a lock named $1 that works across the network on most
# systems; on System V, you should either set LN to /etc/link, or change
# the default below.

OLDLS=""; EQUAL=0;

until touch $1.x; ${LN=ln} $1.x $1 2> /dev/null
do
  LS="`ls -l $1 2> /dev/null`"
  if test "$OLDLS" = "$LS"
  then
    EQUAL=`expr $EQUAL + 1`
    if test $EQUAL -ge 30
    then
      rm -f $1; OLDLS=""; EQUAL=0
    fi
  else
    OLDLS="$LS"; EQUAL=0
  fi
  sleep 30
done
rm -f $1.x
!End!of!file!
chmod 755 acquire
echo Extracting atcheck
cat >atcheck <<'!End!of!file!'
#! /bin/sh
checkrun
!End!of!file!
chmod 755 atcheck
echo Extracting checkrun
cat >checkrun <<'!End!of!file!'
#! /bin/sh

# This file checks to see if the .run file has been touched recently
# If not, it starts up runfactor.  You should make sure that your environment
# contains correct settings for all the variables in environ, and that
# FACTORPATH is set correctly.

if test -z "${FPATH+yes}"
then
  FPATH=/bin:/usr/bin:/usr/ucb:/etc:${FACTORPATH-.}; export FPATH 
fi
PATH=$FPATH; export PATH

. environ

TICKS=`expr 90 '*' $TICKLEN`

run=$OUTDIR/.run.${MACH=`hostname`}
lock=$OUTDIR/.lockfile.$MACH

oldval="`ls -l $run 2> /dev/null`"
sleep $TICKS
acquire $lock
newval="`ls -l $run 2> /dev/null`"
if test "$oldval" = "$newval"
then
  touch $run
  runfactor &
fi
rm -f $lock
!End!of!file!
chmod 755 checkrun
echo Extracting environ
cat >environ <<'!End!of!file!'
#! /bin/sh

# This script contains the code common to checkrun and runfactor listing
# all the useful environment variables you might want to set.

# A mail path back to SRC for orphaned OUT files.  You may want to change
# this to e.g., decsrc::factor or factor%src.dec.com@local.mail.router
MAILTO=${MAILTO-factor@src.dec.com}

# The local command for sending mail, both to you and to SRC
MAILCMD=${MAILCMD-mail}

# The user to whom to send mail when the TASKLIST is empty
MAILUSER=${MAILUSER-$USER}

# Whether to send mail when the TASKLIST is empty; set to "" to not send
DOMAIL=${DOMAIL-yes}

# directory for finding mutable files
FACTORDIR=${FACTORDIR-./mutable}; export FACTORDIR

# directory for output files (and working directory for mpqs)
OUTDIR=${OUTDIR-$FACTORDIR}; export OUTDIR

# location of template file
template=${template-$FACTORDIR/TEMPLATE}

# machine instance name; on a multiprocessor, make sure that you either
# use different directories for each process, or set (and export) MACH
# to distinct values.
MACH=${MACH-`hostname`}; export MACH

# log files -- set to /dev/null to suppress output
control=${control-$OUTDIR/log.$MACH.ctrl}; log=${log-$OUTDIR/log.$MACH}

# number of multiples of TICKLEN to pause when terminal activity detected
# set (and export) this to 0 in your environment to not pause at all
pausetty=${pausetty-1}

# devices examined to see if factoring should pause; other good things to
# inspect might include /dev/mouse and /dev/kbd
ttydevs=${ttydevs-"/dev/tty*"}

# load average above which to pause factoring, exluding factoring.
# on a multiprocessor, increase this enough to let all the copies run
pauseload=${pauseload-0.50}

# A ps command that takes a process id as a parameter; ps -p on System V.
PS=${PS-ps}

# command to suspend and resume a process.  If none (some System V),
# set PAUSECMD to kill and RESUMECMD to echo, or both to "echo" to
# ignore pausing.
PAUSECMD=${PAUSECMD-"kill -TSTP"}
RESUMECMD=${RESUMECMD-"kill -CONT"}

# amount of time to wait between checks of activity, in minutes; must be
# 1, 2, 5, or 10.
TICKLEN=${TICKLEN-5}; export TICKLEN

# A version of the 'link' command that generates an error if the
# file already exists: on System V, you probably want /etc/link
LN={$LN-ln}

!End!of!file!
chmod 755 environ
echo Extracting getnum
cat >getnum <<'!End!of!file!'
#! /bin/sh

# This routine extracts a task range from TASKLIST

lock=$FACTORDIR/.lockfile
file=$FACTORDIR/TASKLIST
ofile="$file".z

acquire $lock
rm -f $ofile
touch $ofile
awk '
	{ if (!go && NF >= 1 && $1 != "#") {
		go = 1;
		print $1;
		if (NF == 2) {
			cnt = $2;
			if (--cnt)
				print ($1 + 1), cnt >>"'"$ofile"'";
		}
	} else
		print $0 >> "'"$ofile"'";}
END	{ if (!go)
		print "";}' $file
mv -f $ofile $file
touch $file
chmod 666 $file
rm -f $lock
!End!of!file!
chmod 755 getnum
mkdir mpqssrc
chmod 755 mpqssrc
echo Extracting mpqssrc/Makefile
cat >mpqssrc/Makefile <<'!End!of!file!'
# We recommend using the optimizer, and gcc if you have it.  On VAXes
# and Suns, the following work well:
#
# CC = gcc
# CFLAGS = -O
#
# On DECstations, SGI boxes, and MIPS computers, at least, we like:
#
# CFLAGS = -O4
#

# If you want to build xtime, so that you can include it in checking for
# activity, you'll want
# PROGS = mpqs xtime
# instead.
PROGS = mpqs

all: $(PROGS)

mpqs: mpqs.c basic.c
	$(CC) $(CFLAGS) -o mpqs mpqs.c basic.c -lm

xtime: xtime.c
	$(CC) $(CFLAGS) -o xtime xtime.c -lX11
!End!of!file!
echo Extracting mpqssrc/basic.c
cat >mpqssrc/basic.c <<'!End!of!file!'
/* Here is basic.c.  Someone who knows gcc can probably make better versions /*
/* of zdiv21 and zsmulmod (like zaddmulp contribution); be careful */
/* because zsmulmod is also referenced in the new mpqs.c.   */

#include <stdio.h>
#include <math.h>
/* NBITS must be even */
#define NBITS 30
#define NBITSH (NBITS>>1)
#define RADIX (1<<NBITS)
#define RADIXM (RADIX-1)
#define RADIXH (1<<(NBITS-1))
#define RADIXROOT (1<<NBITSH)
#define RADIXROOTM (RADIXROOT-1)
#define ZSIZEP 17
#define ZZSIZEP (2*ZSIZEP-1)
#define	ZDECSIZE	291

#define	PRIMBND		5000
#define	INDBND		48

/* a = a+t+b*d % RADIX; t = b*d / RADIX */
/* should be optimized */
 
#if 0
#define		zaddmulp(_a,_b,_d,_t) \
{ register long lb = (_b), ld = (_d), \
  		b1 = (_b) & RADIXROOTM, d1 = (_d) & RADIXROOTM; \
  register long aa = *(_a) + b1 * d1; \
  b1 = b1* ( ld >>= NBITSH ) + d1* ( lb >>= NBITSH ) + (aa >> NBITSH); \
  aa = (aa & RADIXROOTM) + ((b1 & RADIXROOTM) << NBITSH) + *(_t); \
  *(_t) = ld*lb + (b1 >> NBITSH) + (aa >> NBITS); \
  *(_a) = (aa & RADIXM); \
} 
#endif
#if 1
/* 
	This definition of zaddmulp presumes a two's complement 
	machine in which integer overflow is ignored
	and where double precision arithmetic is fast.
	The 0.25 allows round to nearest or round towards zero
	(value being rounded should be integer except for roundoff.)
*/

static double fradix_inv = 1.0/RADIX;	/* Global constant */

#define zaddmulp(_a,_b,_d,_t) \
{ register at = *(_a) + *(_t); \
  register long aa = (at + (_b)*(_d)) & RADIXM; \
  *(_t) = (long)(0.25 + fradix_inv*(  ((double)(at-aa)) \
		                    + ((double)(_b)) * ((double)(_d)) ) ); \
  *(_a) = aa; \
}
#endif

/* a += d*b, optimize this one */
/* a and b not at overlapping addresses */
/* (except for d=1) */
 
#define 	zaddmul(ams,ama,amb) \
{ register long lami, lams = (ams), *lama = (ama), *lamb = (amb); \
  if (lams == 1) { /* Special case in zmod and zdiv */ \
     lams = 0; \
     for (lami=(*lamb++);lami>0;lami--) \
	{  lams += (*lama + *lamb++); \
	  *lama++ = lams & RADIXM; \
	   lams >>= NBITS; } \
     *lama += lams; \
  } else { long lamcarry=0; \
     for (lami=(*lamb++);lami>0;lami--) \
	{ zaddmulp(lama,*lamb,lams,&lamcarry); lama++; lamb++; } \
     *lama += lamcarry; \
  } \
}
/* Be careful, the last lama is unnormalized */

/* a[s:] += (a[s]*zminv)*b so that a[s] = 0 */
/* only in montgomery multiplication zmont */
/* should be optimized. The meaning of s is clear from zmont, */
/* someting like zmtop-i +- 1 */

#define		zmmulp(mmpa,mmpb) \
{ register long *lmmpa = (mmpa), *lmmpb =(mmpb); \
  register long lmmi = (*lmmpa) >> NBITSH, lmmd = (*lmmpa) & RADIXROOTM; \
  long zmmtemp = 0; \
  lmmd=(zminv1*lmmd+(((zminv1*lmmi+zminv2*lmmd)&RADIXROOTM)<<NBITSH))&RADIXM; \
  for (lmmi = *lmmpb++;lmmi>0;lmmi--) \
	{ zaddmulp(lmmpa,*lmmpb,lmmd,&zmmtemp); lmmpa++; *lmmpb++; } \
  if (((*lmmpa += zmmtemp) & RADIX) > 0 ) { (*lmmpa++)&=RADIXM; (*lmmpa)++; } \
}

/* That was somewhat dirty. If you prefer you can replace the */
/* above three macro's (zaddmulp, zaddmul, and zmmulp) by the */
/* following three C-routines */

/*
zaddmulp(a,b,d,t)
long *a,b,d,*t;
{ register long b1 = b & RADIXROOTM, d1 = d & RADIXROOTM;
  register long aa = *a + b1 * d1;
  b1 = b1* ( d >>= NBITSH ) + d1* ( b >>= NBITSH ) + (aa >> NBITSH);
  aa = (aa & RADIXROOTM) + ((b1 & RADIXROOTM) << NBITSH) + *t;
  *t = d*b + (b1 >> NBITSH) + (aa >> NBITS);
  *a = (aa & RADIXM);
}

zaddmul(d,a,b)
long d,*a,*b;
{ register long i;
  long carry=0;
  for (i=(*b++);i>0;i--) zaddmulp(a++,*b++,d,&carry);
  *a += carry;
}

zmmulp(pa,pb)	
long *pa,*pb;
{ register long i = (*pa) >> NBITSH, d = (*pa) & RADIXROOTM;
  long temp = 0;
  d=(zminv1*d+(((zminv1*i+zminv2*d)&RADIXROOTM)<<NBITSH))&RADIXM;
  for (i = *pb++;i>0;i--) zaddmulp(pa++,*pb++,d,&temp);
  if (((*pa += temp) & RADIX) > 0 ) { (*pa++)&=RADIXM; (*pa)++; }
}
*/


/* global variables */

static double epsilon, fradix, fudge; 			/* for long division */
static long zseed[ZSIZEP], zranp[ZSIZEP], zprroot[ZSIZEP];	/* for random */
static long zr[ZSIZEP], zrr[ZSIZEP], zrrr[ZSIZEP], znm[ZSIZEP],
    zmtop, zminv1, zminv2;		/* for montgomery multiplication */
static short int xxprimes[PRIMBND], interval[PRIMBND];	/* to generate primes */
static long pshift = -2;
static long pindex;

/* end global variables */

zstart()		/* set values global variables */
{ double one = (double)1, half = 1/(double)2;
  fradix = (double)RADIX;	epsilon = one;
  fudge = one + epsilon; /* this is sick, but we have to */
  while ( one != fudge ) { /* because the test one != (one+epsilon) */
	epsilon = epsilon * half; /* leads to disasters on some machines */
	fudge = one + epsilon;
	}
  epsilon += epsilon;
  if ((epsilon*RADIX) > one) printf("decrease RADIX and recompile\n");
  else epsilon *= 3;
  fudge = fradix+epsilon*fradix;
  zseed[0] = -1;
}

zsubmul(r,a,b)		/* a[s:] -= d*b, optimize this one */
			/* a and b not at overlapping addresses */
long r,*a,*b;
{ register long rd = RADIX-r, i;
  long carry = RADIX;
for (i=(*b++);i>0;i--) 
	{ zaddmulp(a,*b,rd,&carry); a++;
	  carry += RADIXM-(*b++);
	}
*a += carry-RADIX; /* unnormalized */
}

zhalt(e)	/* error handling routine */
long e;
{ printf("fatal error number %ld occurred\n",e); fflush(stdout);
  /* Call your favorite error routine here */
  exit(0);
}

zcopy(a,b)	/* b = a */
long *a, *b;
{ register long i;
  for (i= *a; i>=0; i--) *b++ = *a++;
}

zintoz(d,a)	/* a = d, d long */
long d,a[];
{ register long i=0;
  a[1]=0;
  while (d > 0)
	{ i++; a[i] = d&RADIXM; d >>= NBITS; }
  if (i>0) a[0] = i;
  else a[0] = 1; 
}

long ztoint(a)	/* returns (long)a, no overflow check */
long *a;
{ register long d, i, sa;
  sa = *a; a += sa; d = *a;
  for (i=sa-1; i>0; i--)
	{ d <<= NBITS; d += *(--a); }
  return(d);
}
  

long zcompare(a,b)	/* a>b:1; a=b:0; a<b:-1 */
long *a, *b;
{ register long sa = *a, sb = *b, i;
if (sa > sb ) return(1);
if (sa < sb ) return(-1);
a += sa; b += sa;
for (i=sa;i>0;i--)
	{ if (*a > *b) return(1);
	  if (*a-- < *b-- ) return(-1);
	}
return(0);
}

zsadd(a,d,b)	/* b = a+d, abs(d)<RADIX, output can be input */
		/* assume result positive */
long a[],d,b[];
{ long x[2];
  x[0] = 1;
  if (d >= 0) { x[1] = d; zadd(a,x,b); }
  else { x[1] = -d; zsub(a,x,b); }
}

zaddls(a,b,c)	/* c = a+b, b shorter than a, output can be input */
long *a,*b,c[];
{ register long sa = *a, sb = *b, carry=0, i, *pc;
pc = &c[0];
for (i=1; i<=sb; i++)
	{ if (( *(++pc) = (*(++a))+(*(++b))+carry) < RADIX ) carry = 0;
	  else	{ *pc -= RADIX; carry = 1; }
	}
for (; i<=sa; i++)
	{ if (( *(++pc) = (*(++a))+carry) < RADIX ) carry = 0;
	  else	{ *pc -= RADIX; carry = 1; }
	}
if ( carry == 1 )
	{ c[0] = sa+1; *(++pc) = 1; }
else c[0] = sa;
}

zadd(a,b,c)	/* c = a+b, output can be input */
long *a,*b,c[];
{ register long sa = *a, carry=0, i, *pc;
pc = &c[0];
if ( sa == *b )
	{ 
	for (i=sa; i>0; i-- )
		{ if (( *(++pc) = (*(++a))+(*(++b))+carry) < RADIX ) carry = 0;
		  else	{ *pc -= RADIX; carry = 1; }
		}
	if ( carry == 1 )
		{ c[0] = sa+1; *(++pc) = 1; }
	else c[0] = sa;
	}
else if ( sa > *b ) zaddls(a,b,c);
else zaddls(b,a,c);
}

zsub(a,b,c)	/* c = a-b for a >= b, output can be input */
long *a,*b,c[];
{ register long sa = *a, sb = *b, carry=0, i, *pc;
pc = &c[0];
for (i=1; i<= sb; i++)
	{ if ((*(++pc) = (*(++a))-(*(++b))-carry) >= 0 ) carry = 0;
	  else { *pc += RADIX; carry = 1; };
	}
for (; i<=sa; i++)
	{ if ((*(++pc) = (*(++a))-carry) >= 0 ) carry = 0;
	  else { *pc += RADIX; carry = 1; };
	}
i = sa;
while ((i>1) && (*pc == 0 )) { i--; pc--; }
c[0] = i;
}

zsmul(a,d,b)	/* b = d*a, d < RADIX, output can be input */
long *a,d,b[];
{ register long sa = *a, i, *pb;
pb = &b[0];
if (d == 0 )
	{ *pb++ = 1; *pb = 0; }
else	{ for (i=sa;i>=0;i--) *pb++ = *a++;
	  sa ++; *pb = 0;
	  zaddmul(d-1,&b[1],&b[0]);
	  while ((sa>1) && (b[sa] == 0)) sa --;
	  b[0] = sa;
	}
}

zmul(a,b,c)	/* c =a*b, output cannot be input */
long *a,*b,c[];
{ register long i, *pc, sc = *a + *b;
  pc = &c[1];
  for (i= sc; i>0; i--) *pc++ = 0;
  pc = &c[1];
  if (*a <= *b ) for (i= *a;i>0;i--) { zaddmul(*(++a),pc++,b); }
  else 	for (i= *b;i>0;i--) { zaddmul(*(++b),pc++,a); }
  while ((sc>1) && (c[sc] == 0)) sc --;
  c[0] = sc;
}

/*
	zdiv21 returns quot, rem so

	quot = (numhigh*RADIX + numlow)/denom;
	rem  = (numhigh*RADIX + numlow)%denom;

Assumes 0 <= numhigh < denom < RADIX and 0 <= numlow < RADIX.
*/

#define zdiv21(numhigh, numlow, denom, quot, rem) {\
    register long lr21, lq21 = (long)(((fradix*(double)(numhigh)) \
		              + (double)(numlow)) / (double)(denom)); \
    if (1) { \
      long lprodhigh = 0, lprodlow = 0; \
      zaddmulp(&lprodlow, lq21, denom, &lprodhigh); \
      lr21 = RADIX*(numhigh - lprodhigh) + (numlow - lprodlow); \
    } else { \
/* Following works in many two's complement architectures. */ \
      lr21 = (numhigh << NBITS) + numlow  - lq21*denom; \
    }\
    if (lr21 < 0) { do {lq21--;} while ((lr21 += denom) < 0); \
    } else      { while (lr21 >= denom) {lr21 -= denom; lq21++;}; \
    }  \
    *quot = lq21; *rem = lr21; \
}

long zsdiv(a,d,b)	/* b = a/d, returns a%ld, output can be input */
long *a,d,*b;
{ 
if (d >= RADIX) {
	long zd[ZSIZEP];
	zintoz(d,zd); zdiv(a,zd,b,zd); return(zd[1]); 
} else {
	register long den = d, sa = *a, carry=0, i;
	if (a[sa] < den && sa > 1) {
	    carry = a[sa--];
        }	
 	if (d < RADIXROOT) {
	    register long rdivd = RADIX/den, rmodd = RADIX - den*rdivd, temp;
  	    for (i=sa; i>0; i--)
		{ temp = a[i] + rmodd*carry;
	  	b[i] = rdivd*carry + temp/d;
	  	carry = temp % d;
		}
	} else {
	    for (i = sa; i > 0; i--) 
	        { long newcarry;
	    	zdiv21(carry, a[i], den, &b[i], &newcarry);
	    	carry = newcarry;
		}
        }
  	while ((sa>1) && (b[sa] == 0)) sa--;
  	b[0] = sa;
  	return(carry);
}
}

zread(a)		/* read a from standard input */
long a[];
{ char in[ZDECSIZE];
  register long d=0;
  scanf("%s",in);
  a[0] = 1; a[1] = 0;
  while (in[d] != '\0')
	{ if (in[d] == '\\') 
		{ scanf("%s",in); d = 0; }
	  else
		{ zsmul(a,(long)10,a);
		  zsadd(a,(long)(in[d++]-'0'),a);
		};
}	}

fzread(fil,a)		/* read a from file fil */
FILE *fil;
long a[];
{ char in[ZDECSIZE];
  register long d=0;
  fscanf(fil,"%s",in);
  a[0] = 1; a[1] = 0;
  while (in[d] != '\0')
	{ if (in[d] == '\\') 
		{ fscanf(fil,"%s",in); d = 0; }
	  else
		{ zsmul(a,(long)10,a);
		  zsadd(a,(long)(in[d++]-'0'),a);
		};
}	}

long zwrite(a)	/* write a on standard output */
long a[];
{ long ca[ZZSIZEP], out[ZDECSIZE], zsdiv();
  register long sa = a[0], i;
  long result;
  for (i=0; i<=sa;i++) ca[i] = a[i];
  i = -1;
  do out[++i] = zsdiv(ca,(long)10,ca) ;
  while (ca[0] > 1 || ca[1] != 0 );
  sa = 0;
  if (i < 0 ) i=0;
  result = i+1;
  while (i >= 0 )
	{ printf("%ld",out[i--]);
	  sa ++;
	  if ((i >= 0) && (sa == 70 ))
		{ sa = 0; printf("\\\n"); }
	}
  return(result);
}

long fzwrite(fil,a,lll,eol,bol)
		/* write a on file fil, at most lll digits per line */
		/* end line with eol, begin next line with bol */
FILE *fil;
long a[],lll;
char *eol, *bol;
{ long ca[ZZSIZEP], out[ZDECSIZE], zsdiv();
  register long sa = a[0], i;
  long result;
  for (i=0; i<=sa;i++) ca[i] = a[i];
  i = -1;
  do out[++i] = zsdiv(ca,(long)10,ca) ;
  while (ca[0] > 1 || ca[1] != 0 );
  sa = 0;
  if (i < 0 ) i=0;
  result = i+1;
  while (i >= 0 )
	{ fprintf(fil,"%ld",out[i--]);
	  sa ++;
	  if ((i >= 0) && (sa == lll ))
		{ sa = 0; fprintf(fil,"%s\n%s",eol,bol); }
	}
  return(result);
}

zdiv(a,b,q,r)	/* q = a/b, r = a%b, output can be input */
long a[],b[],q[],r[];
{ register long i, qq;
  long c[ZZSIZEP], d[ZZSIZEP], sa, sb = b[0], sq, *p, *pc, zsdiv();
  double btopinv, aux;
  p = &b[sb];
  if ((sb == 1) && (*p < RADIX)) zintoz(zsdiv(a,*p,q),r);
  else
	{ sa = a[0]; sq = sa-sb; btopinv = (double)(*p)*fradix;
	  if (sb > 1) btopinv += (*(--p));
	  btopinv *= fradix;
	  if (sb > 2) btopinv += (*(p-1));
	  btopinv = fudge / btopinv;
	  p = &a[1]; pc = &c[0];
	  for (i=sa;i>0;i--) *pc++ = *p++;
	  p = pc; sa = 1-sb;
	  for (i=(-sq);i>0;i--) *p++ = 0;
	  *pc = 0; d[1] = 0; p = &d[sq+1];
	  for (i=sq;i>=0;i--)
		{ aux = fradix*(fradix*(*pc)+(*(pc-1)))+1.0;
		  if (i > sa) aux += (*(pc-2));
		  qq = (long)(btopinv*aux+0.5);
		/* dirty, but safe. On most machines (long)(btopinv*aux) */
		/* is correct, or one too big; on some however it becomes*/
		/* too small. Could change zstart, but +0.5 and a while */
		/* instead of one single if is safer */
		  if (qq > 0)
			{ zsubmul(qq,&c[i],&b[0]); 
			  while (*pc < 0) 
				{ qq --;
				  { zaddmul(1,&c[i],&b[0]); }
				}
			}
		  pc--; *p-- = qq;
		} /* loop on i */
	sb --;
	while ((sb > 0) && (c[sb] == 0)) sb --;
	sb ++;r[0] = sb; p = &r[1]; pc = &c[0];
	for (i=sb;i>0;i--) *p++ = *pc++;
	if (sq < 0)
		{ q[0] = 1; q[1] = d[1]; }
	else
		{ sq ++;
		  while ((sq > 1) && (d[sq] == 0)) sq --;
		  q[0] = sq; p = &q[1]; pc = &d[1];
		  for (i=sq;i>0;i--) *p++ = *pc++;
		}
	} /* non-trivial case */

}

zmod(a,b,r)	/* r = a%b, output can be input */
long a[],b[],r[];
{ register long i, qq;
  long c[ZZSIZEP], sa, sb = b[0], sq, *p, *pc, zsdiv();
  double btopinv, aux;
  p = &b[sb];
  if ((sb == 1) && (*p < RADIX)) zintoz(zsdiv(a,*p,c),r);
  else
	{ sa = a[0]; sq = sa-sb; btopinv = (double)(*p)*fradix;
	  if (sb > 1) btopinv += (*(--p));
	  btopinv *= fradix;
	  if (sb > 2) btopinv += (*(p-1));
	  btopinv = fudge / btopinv;
	  p = &a[1]; pc = &c[0];
	  for (i=sa;i>0;i--) *pc++ = *p++;
	  p = pc; sa = 1-sb;
	  for (i=(-sq);i>0;i--) *p++ = 0;
	  *pc = 0;
	  for (i=sq;i>=0;i--)
		{ aux = fradix*(fradix*(*pc)+(*(pc-1)))+1.0;
		  if (i > sa) aux += (*(pc-2));
		  qq = (long)(btopinv*aux+0.5);
		/* see comment in zdiv */
		  if (qq > 0)
			{ zsubmul(qq,&c[i],&b[0]);
			  while (*pc < 0) 
				{ { zaddmul(1,&c[i],&b[0]); }
				}
			}
		  pc--;
		} /* loop on i */
	  sb --;
	  while ((sb > 0) && (c[sb] == 0)) sb --;
	  sb ++;r[0] = sb; p = &r[1]; pc = &c[0];
	  for (i=sb;i>0;i--) *p++ = *pc++;
	} /* non-trivial case */

}

zaddmod(a,b,n,c)	/* c = (a+b)%n, with 0<=a,b<n */
long a[],b[],n[],c[];
{ long zcompare();
  zadd(a,b,c);
  if (zcompare(c,n) >= 0 ) zsub(c,n,c);
}

zsubmod(a,b,n,c)	/* c = (a-b)%n, with 0<=a,b<n */
long a[],b[],n[],c[];
{ long d[ZSIZEP], zcompare();
  if (zcompare(a,b) >= 0) zsub(a,b,c);
  else { zadd(a,n,d); zsub(d,b,c); }
}

long zsmulmod(a, b, n)	/* return (a*b)%n */
long a,b,n;
{
  register long lqmul = (long)(((double)a)*((double)b) /((double)n));
#if 0
  register long lr;
  long lprodhigh1=0, lprodhigh2=0, lprodlow1=0, lprodlow2=0;
  zaddmulp(&lprodlow1,  a, b, &lprodhigh1);
  zaddmulp(&lprodlow2, lqmul, n, &lprodhigh2);
  lr = RADIX*(lprodhigh1 - lprodhigh2) + (lprodlow1 - lprodlow2); 
#else
/* Many machines compute the following modulo 2^32, which is OK */
  register long lr = a*b - lqmul*n;
#endif
  while (lr >= n) lr -= n;
  while (lr < 0) lr += n;
  return(lr);
}

zmulmod(a,b,n,c)	/* c = (a*b)%n */
long a[],b[],n[],c[];
{ long mem[ZZSIZEP];
  zmul(a,b,mem);	zmod(mem,n,c);
}

zmstart(n)	/* set up global values for Montgomery's multiplication */
long n[];
{ long x[ZZSIZEP], i, zsdiv();
  if ((n[1]&1) == 0) zhalt((long)3);
  zmtop = n[0];
  if (n[zmtop] >= (RADIX >>1)) zmtop ++;
  if (zmtop >= ZSIZEP ) zhalt((long)4);
  zsadd(n,(long)1,zrr); i=zsdiv(zrr,(long)2,x);	/* x = (n+1)/2	*/
  zrr[0]=1; zrr[1]=NBITS; zexp(x,zrr,n,x);	/* x = 2^-NBITS %n	*/
  for (i=x[0]; i>0; i--) x[i+1] = x[i];
  x[1] = 0; x[0] ++;
  zsadd(x,(long)-1,x);
  zdiv(x,n,x,zrr); zminv1 = x[1];	/* zminv1=n[1]^-1 modulo 2^NBITS */
  zminv2 = zminv1 >> NBITSH;	zminv1 &= RADIXROOTM;
  for (i=1; i<=zmtop; i++) x[i] = 0; x[zmtop+1] = 1; x[0] = zmtop+1;
  zdiv(x,n,zrr,zr);		zmulmod(zr,zr,n,zrr);
  zmulmod(zr,zrr,n,zrrr);	zsubmod(n,zr,n,znm);
}

ztom(a,n,b)	/* b is montgomery representation of a */
long a[],n[],b[];
{ long zcompare();
  if (zcompare(n,a) < 0)
	{ zmod(a,n,b); zmont(b,zrr,n,b);}
  else	zmont(a,zrr,n,b);
}

zmtoz(a,n,b)	/* montgomery a is backtransformed to normal b */
long a[],n[],b[];
{ long x[2];
  x[0] = 1; x[1] = 1;
  zmont(a,x,n,b);
}

zmont(a,b,n,c)	/* c is mont product of mont's a and b */
		/* output can be input */
long *a,*b,*n,c[];
{ register long i;
  long x[ZZSIZEP], *px, *pc, zcompare();
  px = &x[0];
  for (i=ZZSIZEP;i>0;i--) *px++ = 0;
  px = &x[0];
  if (*a <= *b ) for (i= *a;i>0;i--) { zaddmul(*(++a),px++,b); }
  else 	for (i= *b;i>0;i--) { zaddmul(*(++b),px++,a); }
  px = &x[0];
  for (i=zmtop; i>0; i--) zmmulp(px++,n);
  pc = &c[1];
  for (i=zmtop; i>0; i--) *pc++ = *px++;
  i = zmtop;
  while ((i > 1) && (*--pc == 0)) i --;
  c[0] = i;
  if (zcompare(n,c) <= 0) zsub(c,n,c);
}

zexp(a,e,n,b)		/* b = (a^e)&n, b and a can be identical, but */
			/* both unequal to n and e */
			/* only if RADIX is power of two */
long a[],e[],n[],b[];
{ register long i,j, cntr, ei;
  long sq[ZZSIZEP], *pae;
  if ((e[0] == 1) && (e[1] == 0 ))
	{ b[0] = 1; b[1] = 1; }
  else
	{ zmod(a,n,sq);
	  b[0] = 1; b[1] = 1; cntr = 0; pae = &e[0];
	  for (i=(*pae); i>0; i--)
		{ for (j=cntr; j>0; j--) zmulmod(sq,sq,n,sq);
		  ei = *(++pae); cntr = NBITS;
		  if (ei > 0)
		  	{ while (ei > 1)
				{ if ( (ei&1) == 1 ) zmulmod(b,sq,n,b);
			  	  ei >>= 1;
			  	  zmulmod(sq,sq,n,sq); cntr--;
				}
			  zmulmod(b,sq,n,b);
			}
		}
	}
}

long zsexp(a,e,n)	/* Single precision version of zexp */
long a,e,n;
{ 
    if (e == 0)
     	return(1); 
    else
	{ register long aa = a%n, ipow2 = 1;
	  long b = aa;
	  while ((ipow2 << 1) <= e) ipow2 <<= 1; 
	  while ((ipow2 >>= 1) != 0) { 
		/* Answer = (b^(2*ipow2))*(aa^(e%(2*ipow2))) mod n */
		b = zsmulmod(b, b, n);
	     	if ((e & ipow2) != 0) b = zsmulmod(b, aa, n);
		}
	  return(b);
	}
}

/* returns the inverse of n modulo the prime p	*/
/* first the old version, undefined if not coprime */
/* then a faster binary version */
/*
long sinv(n,p)
long n,p;
{ register long q=n, r, nq, nr, u=1, w, nw, par=0;
if (n > p) { w = 0; r = p; }
else r = p- (w = p/n) *n;
nr = r;
while (nr != 0) {
	if ( (nr = q - r) < r) nq = 1;
	else if ( ( nr -= r ) < r) nq = 2;
	else nr = q - ( nq = q/r ) *r;
	nw = nq*w+u;	u = w;
	w = nw;		q = r;
	r = nr;		par = 1-par;
	}
if (par == 0) return(u);
else return(p-u);
}
*/

long sodinv(n,p) long n, p; { /* only for odd p>=3 and n>0 */
	register long n1 = n, n2 = p, preg = p, m1 = 1, m2 = 0;
		/* m1*n == n1 mod p,  m2*n == n2 mod p */
                /* n2 is odd, n1 is odd after initialization */
	while (!(n1&1)) {
		n1 >>= 1;
		if (m1&1) m1 += preg;
		m1 >>= 1;
	}
	if (n1 == 1) return m1;
	while (n1 != n2) {
		if (n1 >n2) {
			n1 -= n2; m1 -= m2;
			if (m1 < 0) m1 += preg;
			do {
				n1 >>= 1;
				if (m1&1) m1 += preg;
				m1 >>= 1;
			} while (!(n1 & 1));
			if (n1 == 1) return m1;
		} else {
			n2 -= n1; m2 -= m1;
			if (m2 < 0) m2 += preg;
			do {
				n2 >>= 1;
				if (m2&1) m2 += preg;
				m2 >>= 1;
			} while (!(n2 & 1));
			if (n2 == 1) return m2;
		}
	}
	zhalt(8); /* non-trivial gcd n1 */
}


long sinv (n,p) long n,p; { /* halts if non-trivial gcd or even p */
	if (p < 3 || (!(p & 1)) || n <= 0) zhalt(7);
	return( sodinv(n,p) );
}

long zinv(ain,nin,inv)	/* if 0, then ain*inv%nin=1 */
			/* if 1, then gcd(ain,nin)=inv > 1 */
long ain[],nin[],inv[];
{ long a[ZZSIZEP], n[ZZSIZEP], q[ZZSIZEP], u[ZZSIZEP],
      w[ZZSIZEP], x[ZZSIZEP], y[ZZSIZEP], z[ZZSIZEP], 
      diff, ilo, sa, sn, temp, e,
      try11, try12, try21, try22, got11, got12, got21, got22,
      fast, parity, gotthem, *pin, *p, i;
  double hi, lo, dt, fhi, flo, num, den;
  fhi = 1.0 + epsilon; flo = 1.0 - epsilon;
  pin = &ain[0]; p = &a[0];
  for (i=(*pin); i>=0; i--) *p++ = *pin++;
  pin = &nin[0]; p = &n[0]; 
  for (i=(*pin); i>=0; i--) *p++ = *pin++;
  u[0] = 1; u[1] = 1; w[0] = 1; w[1] = 0;
  while (n[0] > 1 || n[1] > 0)
	{ gotthem = 0; sa = a[0]; sn = n[0]; diff = sa-sn;
 	  if (diff == 0 || diff == 1 )
		{ sa = a[0]; p = &a[sa];
		  num = (double)(*p)*fradix;
		  if (sa > 1) num += (*(--p));
		  num *= fradix;
		  if (sa > 2) num += (*(p-1));
		  sn = n[0]; p = &n[sn];
		  den = (double)(*p)*fradix;
		  if (sn > 1) den += (*(--p));
		  den *= fradix;
		  if (sn > 2) den += (*(p-1));
		  hi = fhi*(num+1.0)/den;
		  lo = flo*num/(den+1.0);
		  if (diff > 0) { hi *= fradix; lo *= fradix; }
		  try11 = 1; try12 = 0; try21 = 0; try22 = 1;
		  parity = 1; fast = 1;
		  while (fast > 0)
			{ parity ^= 1;
/* be careful, some mailers might change the previous line into */
/* open curly bracket parity tilde-equal one semicolon */
/* the correct line is */
/* open curly bracket parity hat-equal one semicolon */
			  if (hi >= fradix )
				{ fast = 0; }
			  else	{ ilo = (long)lo;
				  if (ilo == 0 || ilo < (long)hi)
					{ fast = 0; }
				  else	{ dt = lo; lo = flo/(hi-ilo);
					  if (dt > ilo) hi = fhi/(dt-ilo);
					  else hi = fradix;
					  temp = try11; try11 = try21;
					  if ((RADIX-temp)/ilo < try21)
						fast = 0;
					  else try21 = temp+ilo*try21;
					  temp = try12; try12 = try22;
					  if ((RADIX-temp)/ilo < try22)
						fast = 0;
					  else try22 = temp+ilo*try22;
					  if ((fast > 0) && (parity >0))
						{ gotthem = 1;
						  got11 = try11; got12 = try12;
						  got21 = try21; got22 = try22;
						}
					}
				}
			}
		}
	  if (gotthem > 0)
		{ zsmul(u,got11,x); zsmul(w,got12,y);
		  zsmul(u,got21,z); zsmul(w,got22,w);
		  zadd(x,y,u);     zadd(z,w,w);
		  zsmul(a,got11,x); zsmul(n,got12,y);
		  zsmul(a,got21,z); zsmul(n,got22,n);
		  zsub(x,y,a);     zsub(n,z,n);
		}
	  else	{ zdiv(a,n,q,a); zmul(q,w,x); zadd(u,x,u);
		  if (a[0] > 1 || a[1] > 0)
			{ zdiv(n,a,q,n); zmul(q,u,x); zadd(w,x,w); }
		  else	{ p = &a[0]; pin = &n[0];
			  for (i=(*pin); i>=0; i--) *p++ = *pin++;
			  n[0] = 1; n[1] = 0; zsub(nin,w,u);
			}
		}
	}
  p = &inv[0];
  if (a[0] != 1) e = 1;
  else { if (a[1] == 1) e = 0; else e = 1; }
  if (e == 0) pin = &u[0];
  else 	{ pin = &a[0]; }
  for (i=(*pin); i>=0; i--) *p++ = *pin++;
  return(e);
}

zrstart(s)	/* initalize globals zseed = s, zranp = 2^89-1, */
		/* zprroot = 3^101 mod zranp */
long s;
{ register long i, ones, res;
  long e[2];
  zintoz(s,zseed);
  ones = 89/NBITS; res = 89%NBITS;
  for (i=1; i<=ones; i++) zranp[i] = RADIXM;
  if (res != 0) zranp[++ones] = ( RADIXM >> (NBITS - res) );
  zranp[0] = ones;
  zintoz((long)3,zprroot); e[0] = 1; e[1] = 101;
  zexp(zprroot,e,zranp,zprroot);
}

long zrandom(b) /* returns random in [0,b) */
long b;
{ long zb[ZSIZEP], ztoint();
  if (zseed[0] < 0) zhalt((long)5);
  if (b <= 0) b = 1;
  zmulmod(zseed,zprroot,zranp,zseed);
  zintoz(b,zb);
  zmod(zseed,zb,zb);
  return(ztoint(zb));
}

long zprime(m,t)		/* 0 if m not prime, 1 if m passed t */
			/* pseudo prime tests */
long m[],t;
{ long sm, i, j, s, ok, u[ZSIZEP], m1[ZSIZEP], a[ZSIZEP], zrandom(), zcompare(),
	zsdiv();
 if ((m[0] == 1) && (m[1] <= RADIXROOT))
	{ sm = m[1];
	  if (sm < 2) return(0);
	  if (sm == 2) return(1);
	  if ((sm&1) == 0) return(0);
	  i = 3;
	  while (i*i <= sm )
		{ if ((sm%i) == 0) return(0);
		  i++; i++;
		}
	  return(1);
	}
 else	{ if ((m[1]&1) == 0) return(0);
	  zsadd(m,(long)-1,m1); i=zsdiv(m1,(long)2,u); s = 0;
	  while ((u[1]&1) == 0) { i=zsdiv(u,(long)2,u); s ++; }
	  if (m[1] <= RADIXROOT ) sm = RADIXM-2;
	  else sm = m[1];
	  for (i=t; i>0; i--)
		{ zintoz(2+zrandom(sm),a);
		  zexp(a,u,m,a);
		  if ( (a[0] != 1 || a[1] != 1) && (zcompare(m1,a) != 0))
			{ j = s; ok = 0;
			  while ((j > 0) && (ok == 0))
				{ zmulmod(a,a,m,a);
				  if (zcompare(a,m1) == 0) ok = 1;
				  j--;
				}
			  if (ok == 0) return(0);
			}
		}
	  return(1);
	}
}

zpstart()
{ long i, j, jstep, jstart;
  short int *p;
  if (pshift == -2)
  	{ p = &xxprimes[0]; for (i=PRIMBND; i>0; i--) *p++ = 1;
  	  jstep = 1; jstart = -1;
  	  for (i=0; i<=INDBND; i++)
		  { jstep += 2; jstart += 2*(jstep-1);
  	  	    if (xxprimes[i] > 0)
			  { for (j=jstart; j<PRIMBND; j+=jstep) xxprimes[j] = 0;
			  }
		  }
	}
  pshift = 0; pindex = -1;
}

zpshift()
{ long i, j, jstep, jstart, ibound;
  short int *p, *pi;
  pi = &interval[0];
  if (pshift == 0)
	{ p = &xxprimes[0]; pi = &interval[0];
	  for (i=PRIMBND; i>0; i--) *pi++ = *p++;
	}
 else	{ for (i=PRIMBND; i>0; i--) *pi++ = 1;
	  jstep = 3; ibound = pshift+2*PRIMBND+1;
	  for (i=0; jstep*jstep<=ibound; i++)
		{ if (xxprimes[i] > 0)
			{ jstart = (pshift+2)/jstep+1;
			  if ((jstart&1) == 0) jstart++;
			  if (jstart <= jstep) jstart = jstep;
			  jstart = (jstart*jstep-pshift-3)/2;
			  for (j=jstart;j<PRIMBND; j+=jstep)
				interval[j] = 0;
			}
		  jstep += 2;
		}
	}
}

long zpnext()	/* returns next prime number, starting at 2 */
		/* restart with pshift = -1 (get 2 first) */
		/* or with zpstart() (get 3 first) */
{ if (pshift < 0) { zpstart(); return(2); }
  if (pindex < 0) { pindex = 0; zpshift(); return(3); }
  while ((++pindex) < PRIMBND )
	{ if (interval[pindex] > 0) return(pshift+2*pindex+3); }
  pshift += 2*PRIMBND;
  if (pshift > 2*PRIMBND*(2*PRIMBND+1)) { zpstart(); return(2); }
  zpshift();
  pindex = -1;
  while ((++pindex) < PRIMBND )		/* assuming interval non-empty */
	{ if (interval[pindex] > 0) return(pshift+2*pindex+3); }
/*NOTREACHED*/	/* to shut up lint */
}

long isqrt(n)
long n;
{ long a, ndiva, newa;
  if (n <= 0) return (0);
  if (n <= 3) return (1);
  if (n <= 8) return (2);
  a = n / 2;
  while (1) {
	ndiva = n/a; newa = (ndiva+a)/2;
	if (newa-ndiva <= 1) 
		if (newa*newa <= n) return(newa);
		else return(ndiva);
	a = newa;
	}
}

zsqrt(n,r)
long n[],r[];
{ long a[ZZSIZEP], ndiva[ZZSIZEP], diff[ZZSIZEP], i, isqrt(), zcompare(),
	zsdiv();
  if (n[0] == 1)
	zintoz(isqrt(n[1]),r);
  else	{
	a[0]=(n[0]+1)/2;
	for (i=1;i<a[0];i++) a[i]=0;
	a[a[0]] = isqrt(n[n[0]]) + 1;
	if ((n[0]&1) ==0) a[a[0]] <<= NBITSH;
	if (a[a[0]]&RADIX) {
		a[a[0]] = 0; a[0] ++; a[a[0]] = 1;
	}
  	while (1) { 
		zdiv(n,a,ndiva,r);
		zadd(a,ndiva,r);
		i=zsdiv(r,(long)2,r);
		if (zcompare(r,ndiva) <= 0) break;
		zsub(r,ndiva,diff);
		if ((diff[0] == 1) && (diff[1] <= 1)) {
			zmul(r,r,diff);
			if (zcompare(diff,n) > 0) zcopy(ndiva,r);
			break;
			}
		zcopy(r,a);
		}
	}
}

!End!of!file!
echo Extracting mpqssrc/mpqs.c
cat >mpqssrc/mpqs.c <<'!End!of!file!'
/* improved ppmpqs */
/* mpqs restartfile [infile logfile [outputdirectory/ [mail-cmd]]] */
#include <math.h>
#include <stdio.h>
#include <signal.h>
#include <time.h>

/* NBITS must be even */
#define NBITS 30
#define NBITSH (NBITS>>1)
#define RADIX (1<<NBITS)
#define RADIXM (RADIX-1)
#define RADIXH (1<<(NBITS-1))
#define RADIXROOT (1<<NBITSH)
#define RADIXROOTM (RADIXROOT-1)
#define ZSIZEP 17
#define ZZSIZEP (2*ZSIZEP-1)

#define LINELENGTH 50
#define STARTER 1
#define SLAVE 2
#define FIREFLY 3
#ifndef CHECK 
#define CHECK 1000        /* check every CHECK polynomials */
#endif
#ifndef REPORT 
#define REPORT 25         /* report every REPORT polynomials */
#endif
#ifndef MAXPOL
#define MAXPOL 100000     /* who will ever reach that in one run */
#endif
#define MAXSMALLPRIME 97
#define MAXSMALLPRIMEP 98
#ifndef MAXSIEVE
#define MAXSIEVE 500000
#endif
#define MINSIEVE 100000
#define FIRST_GE_SENTINEL 4
#define SMALLCORRECT 2
#ifdef RISC
#define BYTSOF1 1
#else
#define BYTSOF1 0x01010101
#endif
#define MAXFACNB 200

#ifndef BUILD_RESTART
#define BUILD_RESTART 1 /* 1 to write restart file, 0 to omit */
	/* When creating RESTART, USEHEX reduces file size  */
	/* by 40% but the resultant file is harder to read. */
#endif
#ifdef USEHEX
#define SPNUM " %lx"
#define NOSPNUM "%lx"
#else
#define SPNUM " %ld"
#define NOSPNUM "%ld"
#endif

#ifdef CONTRIBUTION_3
/*  
   Contribution 3 has assembly language versions of selected
   mpqs.c and basic.c routines for 68020-based machines.
*/ 
#define ASM_bigp_siever 1
#define ASM_lowp_siever 1
#define ASM_medp_siever 1
#define ASM_sieveinit   1
#endif
#ifdef CONTRIBUTION_8
/*
   Contribution 8 has assembly language versions of selected
   mpqs.c and basic.c routines for SPARC (i.e., SUN 4) systems.
*/
#define ASM_bigp_siever 1
#define ASM_first_ge    1
#define ASM_lowp_siever 1
#define ASM_medp_siever 1
#define ASM_sieveinit   1
#endif

/* Globals */

#ifdef RISC
typedef long uchar;
#else
typedef unsigned char uchar;
#endif

typedef long *row;
typedef uchar *urow;
extern long zsdiv(), zpnext(), zsmulmod();

FILE *mailit;
FILE *restart;

long n[ZSIZEP], kn[ZSIZEP], sqkn[ZSIZEP], mail_batchcount,
        multiplier, sizefacbase, lastprime, sievelength, sievesize, maxbig,
        startd[ZSIZEP], altstartd[ZSIZEP],
        unique, how, outputcounter = 0, notallocated = 1, identification,
	nbsmallps, nbsmallpsp, nextlog[30], logcnt, closeit,
	ful = 0, maxful = 0, minful = 100000, sumful = 0,
        par = 0, maxpar = 0, minpar = 100000, sumpar = 0,
        ppr = 0, maxppr = 0, minppr = 100000, sumppr = 0,
        bad = 0, maxbad = 0, minbad = 100000, sumbad = 0,
        pol = 0, numfac = 0, interrupted = 0, started = 0;

char mail_execstring[200], *mail_cmd = NULL;
char nameoutput[200], hostname[100], *restartfn;
char parid[20], fulid[20], pprid[20], ppridc[20], paridc[20], fulidc[20];
char default_restart[] = "RESTART";
char nothing[] = "", backsl[] = "\\";
char *file_prefix = "";

long luchar = sizeof(uchar);
 
uchar smallsievebnd, sievebound, firstlog;
double smlogp[MAXSMALLPRIMEP];
 
long smallprimes[] = {0,0,128,243,0,125,0,343,0,0,0,121,0,169,
                        0,0,0,289,0,361,0,0,0,529,0,0,0,0,0,841,0,961,
0,0,0,0,0,1369,0,0,0,1681,0,1849,0,0,0,2209,0,0,0,0,0,2809,0,
0,0,0,0,3481,0,3721,0,0,0,0,0,4489,0,0,0,5041,0,5329,0,0,0,0,0,6241,
0,0,0,6889,0,0,0,0,0,7921,0,0,0,0,0,0,0,9409 };
typedef struct {long r1; long r2;} roots_t;
roots_t *roots;
row primbase, sqrtkn;
urow sieve;
 
/* End Globals */
 
long smalldiv()
{ long p, q[ZSIZEP];
zpstart();
p = zpnext();
while (p <= lastprime) {
        if  (zsdiv(n,p,q) == 0) return(p);
        p = zpnext();
        }
return(0);
}
 
long findmultiplier(t)        /* returns selected multiplier for n, where   */
                        /* t primes will be considered per candidate        */
long t;
{ long longp[2], p, longexp[2], knep[ZSIZEP], try, besttry = 1,
      nmod4, i, j, nbcand = 42;
  double dp, tryval, besttryval = 1.0;
  static long cand[] = { 1,3,5,7,11,13,15,17,19,21,23,29,31,33,35,37,
                        39,41,43,47,51,53,55,57,59,61,65,67,69,71,73,79,
                        83,85,89,91,93,95,97,101,103,105 };
/* that is probably enough! */
nmod4 = n[1] % 4; longp[0] = 1; longexp[0] = 1;
for (j=0;j<nbcand;j++) {
        if (((try = cand[j])*nmod4 % 4) != 1) continue;
        zpstart();                        /* so that next prime will be 3 */
        i = 0;
        tryval = -0.5*log((double)try);
        zsmul(n,try,kn);
        if ((kn[1]%8) == 1) { i++; tryval += 1.38629; } /* log(4) */
        p = zpnext();
        while (i < t) {
                longexp[1] = p/2; longp[1] = p;
                zmod(kn,longp,knep);
                zexp(knep,longexp,longp,knep);
                if (knep[1] == 1) {
                        i++; dp = (double)p; dp = (log(dp))/dp;
                        if (try%p == 0) tryval += dp;
                        else tryval += 2.0*dp;
                        }
                p = zpnext();
                }
        if (tryval > besttryval) { besttry = try; besttryval = tryval; }
        }
return(besttry);
}
 
int jacobi(m, n)	/* Compute Jacobi symbol m over n */
long m, n;
{
    register long m2 = abs(n);
    register long m1;
    register int jcbi = 1;
 
    if ((m2 & 1) == 0) {
	fprintf(stderr, "jacobi - denominator even %ld %ld\n", m, n);
	exit(0);
	}
 
    m1 = (m >= 0) ? (m % m2) : (m2 - ((-m) % m2));
 
    while (1) {
	switch (m1 & 15) {
	case 0:
	    if (m1 == 0)
		return ((m2 == 1) ? jcbi : 0);
 
	    m1 >>= 4;		/* Divide numerator by 16 */
	    continue;           /* continue while loop */
 
	case 8:
	    m1 >>= 2;		/* Divide numerator by 4 */
	    continue;           /* continue while loop */
 
	case 4:
	    m1 >>= 2;		/* Divide numerator by 4 */
	    break;
 
	case 1:
	case 5:
	case 9:
	case 13:
	    break;
 
	case 2:
	case 10:
	    m1 >>= 1;
	    if ((m2 & 7) == 3 || (m2 & 7) == 5)
		jcbi = -jcbi;
	    break;
 
	case 12:
	    m1 >>= 2;		/* Divide numerator by 4, */
				/* fall through to 3 mod 4 cases */
	case 3:
	case 7:
	case 11:
	case 15:
	    if ((m2 & 2) != 0)	/* negate if  m1 == m2 == 3 mod 4 */
		jcbi = -jcbi;
	    break;
 
	case 6:
	case 14:
	    m1 >>= 1;
	    if ((m2 & 4) != 0)
		jcbi = -jcbi;
	    break;
	}                       /* end switch */
 
	{			/* m1 and m2 are odd; invert them */
				/* First divide m2 by 4, mod m1 */
	    register long m3;
 
	    m3 = ((((m2 + m1) & 3) == 0) ? m2 + m1 : m2 - m1) >> 2;
	    if (m3 >= m1) m3 %= m1;
	    m2 = m1;
	    m1 = m3;
	}
    }           /* end while(1) */
}
 
long root(u,p,lit)        /* returns square root of u modulo prime p */
/* dirty, could make this more efficient, needed only in initialization */
long u,p,lit;
{ long e, q, n, v, w, zsmulmod(), zsexp();
  if (u <= 1) return(u);
  if (lit) {  /* separate code for prime powers (smallprimes) */
          for (n=0; n<p; n++) if ((((n*n)-u) % p) ==0 ) return(n) ;
        return(0);
        }
  e = 0; q = (p-1)/2;
  while ((q&1) == 0) { e++; q >>= 1; }  /* p = 1 + q*(2^(e+1)), q odd */
  w = zsexp(u, (q-1)/2, p);
  v = zsmulmod(u, w, p);  /* u^((q+1)/2) */
  if (e != 0) {
     long r, yroots[NBITS];	/* yroots[k] = (2^k) root of -1 */
     w = zsmulmod(v, w, p);  /* u^q */
     for (n=5; ; n+=6) {
	long k;
  	if (n >= 1000) {
	    printf("root - can't find non-residue for %ld \n", p);
	    exit(0);
        }
        if (jacobi(n, p) != -1) continue;
	yroots[e] = zsexp(n, q, p);
        for (k=e; k>0; k--) {
	    yroots[k-1] = zsmulmod(yroots[k], yroots[k], p);
	}
	if (yroots[0] == p-1) break;
     }
     r = e+1;
     while (w != 1) {	/* v^2 == u*w, w^(2^(r-1)) == 1 */
          long k = 1, z = w;
	  while (z != p-1) {	
	       z = zsmulmod(z, z, p);
               if (k++ >= r) {
		  printf("Error in root, perhaps a non-residue:\n");
		  printf("p=%ld, u=%ld, v=%ld, w=%ld, r=%ld, z=%ld, n=%ld\n",
		 	p, u, v, w, r, z, n);
	 	  exit(0);
	       }
	  }
	  v = zsmulmod(v, yroots[k], p);
          w = zsmulmod(w, yroots[k-1], p);
          r=k;
/*
	  printf("root - u = %ld, p=%ld, v=%ld, w=%ld, r=%ld\n",
		u, p, v, w, r);
*/
     } /* while w */
  } /* if e */
  if (v > p/2) v = p-v;
  if (zsmulmod(v,v,p) != u) {
	printf("Error in root: %ld %ld %ld \n", u, p, v);
	exit(0);
  }
  return(v);
}
 
createfactorbase()
{ long lit, p, knep[ZSIZEP], knmodp, i = 1;
  long zsexp(), zsdiv(), zpnext();
uchar addlog;
long exp_addlog;
zpstart();        /* so that next prime will be 3 */
if ((kn[1]%8) == 1) {
        primbase[1] = 2;
        sqrtkn[1] = root(kn[1],smallprimes[2],(long)1);
        i = 2; nbsmallps = 1;
	}
else nbsmallps = 0;
while (i <= sizefacbase) {
	lit=0;
	p = zpnext();
        knmodp = zsdiv(kn, p, knep);
        if (jacobi(knmodp, p) != -1 || p == multiplier) {
                primbase[i] = p;
                if (p <= MAXSMALLPRIME) {
			lit =1;
			nbsmallps ++;
                        p = smallprimes[p];
                        knmodp = zsdiv(kn, p, knep);
		}	
                sqrtkn[i] = root(knmodp,p,lit);
                i++;
                }
        }
nbsmallpsp = nbsmallps + 1;
lastprime = primbase[sizefacbase];
for (i=1; i<=nbsmallps; i++)
	smlogp[i]=(log((double)primbase[i])/1.0);
logcnt = 0;
if (i<=sizefacbase) {
	firstlog = (uchar)(log((double)primbase[i])/1.0);
	addlog = firstlog;
	exp_addlog = (long)(exp((double)(addlog+1))/1.0);
}
i++;
for (; i<=sizefacbase; i++) {
	if (primbase[i] > exp_addlog ) {
		nextlog[logcnt++] = primbase[i];
		addlog ++; /* step > 1 possible error, but unlikely */
				/* to affect many logs, at most a few */
		exp_addlog = (long)(exp((double)(addlog+1))/1.0);
	}
}
nextlog[logcnt] = lastprime+1;
}
 
void findsievebound()
/* better to determine experimentally */
/* this one gives a rather low result (which is better, */
/* because it gives more information. In the course of the */
/* computation the bound should be adapted.        */
{ double lkn, loglast;
lkn = log((double)kn[kn[0]]) + (double)(kn[0]-1)*log((double)RADIX);
loglast = 2 * log((double)maxbig);
sievebound=(uchar)(((lkn+log((double)sievelength))/2.0-loglast)/1.0-4.0);
smallsievebnd= sievebound - (uchar)SMALLCORRECT;
}
 
long dumpdata() {
/* dump in decimal notation, and not using fwrite, so */
/* that different machines can read each others restart file */
/* slow, but saves more time that it wastes */
/* put on separate lines, to make it easier to change RESTART */
/* using some editor (long line can cause problems) */
register long i;
if (fprintf(restart,"%ld",sizefacbase) < 0) return(1);
if (fprintf(restart," %ld\n",identification) < 0) return(1);
fzwrite(restart,n,1000,backsl,nothing); fprintf(restart,"\n");
if (fprintf(restart,"%ld\n",multiplier) < 0) return(1);
fzwrite(restart,sqkn,1000,backsl,nothing); fprintf(restart,"\n");
if (fprintf(restart,"%ld",sievelength) < 0) return(1);
if (fprintf(restart," %ld",(long)sievebound) < 0) return(1);
if (fprintf(restart," %ld\n",maxbig) < 0) return(1);
if (fprintf(restart," %ld\n",nbsmallps) < 0) return(1);
for (i=1; i<=nbsmallps; i++) {
	if (fprintf(restart," %ld",(long)smlogp[i]) < 0) return(1);
        if (i%5 == 0) fprintf(restart,"\n");
        }
if (fprintf(restart," %ld\n",(long)firstlog) < 0) return(1);
if (fprintf(restart," %ld\n",logcnt) < 0) return(1);
for (i=0; i<=logcnt; i++) {
	if (fprintf(restart," %ld",nextlog[i]) < 0) return(1);
        if (i%5 == 0) fprintf(restart,"\n");
        }
if (i%5 != 1) fprintf(restart,"\n");
fzwrite(restart,startd,1000,backsl,nothing); fprintf(restart,"\n");
for (i=1; i<=sizefacbase; i++) {
#ifdef USEHEX
        long gap = (i <= 2) ? primbase[i]
                            : ((primbase[i]-primbase[i-1]) >> 1);
#else
        long gap = primbase[i];
#endif
        if (fprintf(restart," %ld",gap) < 0) return(1);
        if (i%5 == 0) fprintf(restart,"\n");
        }
for (i=1; i<=sizefacbase; i++) {
        if (fprintf(restart,SPNUM,sqrtkn[i]) < 0) return(1);
        if (i%5 == 0) fprintf(restart,"\n");
        }
fclose(restart);
return(0);
}
 
makestartd() {
	long zcompare();
        if (unique < 0) {
       		altstartd[1] = 0; altstartd[2] = -unique;
       		altstartd[0] = 2;
		if (zcompare(altstartd,startd) >= 0) {
                    printf("negative unique process ID not allowed, exit\n");
               		exit(0);
		}
		zsub(startd,altstartd,startd);
		startd[1] |= 3;
        }
	else {
		altstartd[1] = 0; altstartd[2] = unique;
		if (unique == 0) altstartd[0] = 1;
		else altstartd[0] = 2;
		zadd(altstartd,startd,startd);
        	startd[1] |= 3;
	}
/* no overlap is almost guaranteed in this way, unless machines are */
/* VERY fast, and run for a LONG time. resulting values are close */
/* enough to optimal for our range of interest (10^90 - 10^110) */
}
 
getmemory() {
long bytefacbase = (sizefacbase+1) * sizeof(long);
if (((primbase = (long*)malloc(bytefacbase))==NULL)||
        ((sqrtkn=(long*)malloc(bytefacbase))==NULL) ||
        ((roots=(roots_t*)malloc((sizefacbase+1)*sizeof(roots_t)))==NULL)) {
                printf("unable to get memory, exit\n");
                exit(0);
        }
        notallocated = 0;
#if 0
	printf("primbase = %lx, sqrtkn = %lx, roots = %lx\n",
		    (long)primbase, (long)sqrtkn, (long)roots);
#endif
}
 
long getandalloc() {
register long i;
long dummy, nn[ZSIZEP];
#ifndef VMS
if (how == FIREFLY) {
        if ((restart=fopen("./RESTART","r"))==NULL)
                return(1);
        if (fscanf(restart,"%ld",&sizefacbase) == EOF) goto wrong;
        if (fscanf(restart,"%ld",&dummy) == EOF) goto wrong;
        if (dummy != identification) goto wrong;
        fzread(restart,n);
 
}
else
#endif /* VMS */
if (how == SLAVE) {
        printf("Attempt to read restart file %s\n", restartfn);
	fflush(stdout);
#ifdef VMS
        restart=fopen(restartfn,"r","mbf=4","mbc=32","rop=RAH");
#else
        restart=fopen(restartfn,"r");
#endif /* VMS */
        if (restart==NULL) return(1);
        if (fscanf(restart,"%ld",&sizefacbase) == EOF) goto wrong;
        if (fscanf(restart,"%ld",&dummy) == EOF) goto wrong;
        if (dummy != identification) goto wrong;
        fzread(restart,nn);
        if (zcompare(n,nn) != 0) goto wrong;
}
else {
        printf("impossible call to getandalloc, exit\n");
        exit(0);
}
 
/*successfully opened file, and contains the right n */
 
getmemory();
if (fscanf(restart,"%ld",&multiplier)==EOF) goto wrong;
fzread(restart,sqkn);
if (fscanf(restart,"%ld",&sievelength)==EOF) goto wrong;
if (fscanf(restart,"%ld",&dummy)==EOF) goto wrong;
sievebound = (uchar)dummy;
smallsievebnd= sievebound - (uchar)SMALLCORRECT;
if (fscanf(restart,"%ld",&maxbig)==EOF) goto wrong;
if (fscanf(restart,"%ld\n",&nbsmallps)==EOF) goto wrong;
nbsmallpsp= nbsmallps+1;
for (i=1; i<=nbsmallps; i++) {
	if (fscanf(restart,"%ld",&dummy)==EOF) goto wrong;
        }
if (fscanf(restart,"%ld\n",&dummy)==EOF) goto wrong;
firstlog = (uchar)dummy;
if (fscanf(restart,"%ld\n",&logcnt)==EOF) goto wrong;
for (i=0; i<=logcnt; i++) {
	if (fscanf(restart,"%ld",&nextlog[i])==EOF) goto wrong;
        }
fzread(restart,startd);
for (i=1;i<=sizefacbase;i++) {
        if (fscanf(restart,"%ld",&primbase[i]) == EOF) goto wrong;
#ifdef USEHEX
 	if (i > 2) primbase[i] = primbase[i-1] + 2*primbase[i];
#endif
 	}
for (i=1;i<=sizefacbase;i++)
        if (fscanf(restart,NOSPNUM,&sqrtkn[i]) == EOF) goto wrong;
 
fclose(restart);
 
for (i=1; i<=nbsmallps; i++)
	smlogp[i]=(log((double)primbase[i])/1.0);
 
lastprime = primbase[sizefacbase];
zsmul(n,multiplier,kn);
if ((multiplier>1) && (multiplier <= MAXSMALLPRIME))
        smallprimes[multiplier] = multiplier;
makestartd();
return(0);
 
wrong:
fclose(restart);
if (how == SLAVE) {
        unlink(restartfn);
        printf("Attempt failed, restart file removed, creating %s\n",
            restartfn);
        fflush(stdout);
}
else printf("Exit: unsuccessful restart\n");
return(1);
}
 

/* set up the "mail_execstring" variable */
setupmailer()
{
        char template_mail_cmd[200];

        printf("the results will be written to the file %s\n",nameoutput);
        scanf("%ld",&mail_batchcount);

	/* the default mail_cmd can be overridden by argv[5] */
	if (mail_cmd == (char*)NULL) {
	    /* no command given as param, let's try TEMPLATE file */
	    mail_cmd = template_mail_cmd;
	    gets(mail_cmd);		/* eat NL */
	    if (gets(mail_cmd) == NULL) {
		/* as last resort, use compiled-in command */
		strcpy(mail_cmd,
#ifdef VMS
			"MAIL/NOSELF %s.; decsrc::factor" /* DEC internal */
#else
		        "mail factor@src.dec.com < %s"    /* Internet */
#endif /*VMS*/
			);
	    }
	}
        if (mail_batchcount)
	{
	    printf("As soon as it contains %ld relations,\n",
		   mail_batchcount);
	    printf("this file will be mailed using the command:\n");
	    sprintf(mail_execstring, mail_cmd, nameoutput);
        } else {
	    printf("Please process %s every now and then:\n",nameoutput);
#ifdef VMS
	    printf("   rename %s MAILOUT\n", nameoutput);
#else
	    printf("   mv %s MAILOUT\n", nameoutput);
#endif /*VMS*/
	    printf(" (wait a moment, maybe process is writing to it)\n");
	    sprintf(mail_execstring, mail_cmd, "MAILOUT");
	}
	printf("   %s\n",mail_execstring);
        fflush(stdout);
}
 
setup () {
    long bytefacbase, dummy, mem1[ZSIZEP], mem2[ZSIZEP];
#if BUILD_RESTART
    char restartname[150];
#endif
 
    zstart(); /* set up long integer arithmetic */
 
    if (how != STARTER) {
        scanf("%ld %ld",&identification,&unique);
        sprintf(parid,"par%ld#",identification);
        sprintf(pprid,"ppr%ld#",identification);
        sprintf(fulid,"ful%ld#",identification);
        sprintf(paridc,"par%ldc#",identification);
        sprintf(ppridc,"ppr%ldc#",identification);
        sprintf(fulidc,"ful%ldc#",identification);
        sprintf(nameoutput,"%sOUT%07ld",file_prefix,unique);
        zrstart(abs(unique)+1);
    }
 
#ifndef VMS
    if (how == FIREFLY) {
        gethostname(hostname,100); /* for debugging*/
        if (getandalloc()) exit(0);
        mail_batchcount = 0;
        started = 1;
        return;
    }
#endif /*not VMS*/
 
    if (how == SLAVE) {
        zread(n);
        if (getandalloc() == 0) {
                printf("successful restart from file %s\n", restartfn);
                scanf("%ld",&dummy); /* multiplier does not change */
                scanf("%ld",&dummy); /* sizefacbase does not change */
                scanf("%ld",&sievelength); /* sievelength may be changed */
                scanf("%ld",&dummy);
                sievebound = (uchar)dummy; /* sievebound may be changed */
                scanf("%ld",&maxbig); /* maxbig may be changed */
                zread(startd);
                makestartd();
        }
        else {
                scanf("%ld %ld",&multiplier,&sizefacbase);
                zsmul(n,multiplier,kn); zsqrt(kn,sqkn);
                if ((multiplier>1) && (multiplier <= MAXSMALLPRIME))
                        smallprimes[multiplier] = multiplier;
                if (notallocated) getmemory();
                scanf("%ld %ld %ld",&sievelength,&dummy,&maxbig);
                sievebound = (uchar)dummy;
 
                printf("start creation new restart data\n");
                fflush(stdout);
                createfactorbase();
                zread(startd);
#if BUILD_RESTART
#ifdef VMS
		/* VMS has versioning, don't need to make unique */
                strcpy(restartname, restartfn);
#else
                sprintf(restartname,"%s%07ld",restartfn, unique);
#endif /*VMS*/
                if ((restart = fopen(restartname,"w")) == NULL) {
                        printf("Cannot create restart file %s, non fatal\n",
			       restartname);
                }
                else {
                        if (dumpdata()) {
                                printf(
				  "Cannot dump restart data, non fatal\n");
                                unlink(restartname);
                        }
                        else {
#ifndef VMS
                                link(restartname, restartfn);
                                unlink(restartname);
#endif /*VMS*/
                                printf(
				  "Created restart file %s\n", restartfn);
                        }
                }
                fflush(stdout);
#endif /* BUILD_RESTART */
                makestartd();
        }
        smallsievebnd= sievebound - (uchar)SMALLCORRECT;
        printf("multiplier: %ld\n",multiplier);
        printf("sizefacbase: %ld\n",sizefacbase);
        printf("sievelength: %ld\n",sievelength);
        printf("sievebound: %ld\n",(long)sievebound);
        printf("maxbig: %ld\n",maxbig);
        printf("startd:"); zwrite(startd);
        printf("\n");
        setupmailer();
        started = 1;
        return;
    }
 
#ifndef VMS
    if (how == STARTER) { /*includes FIREFLYSTART */
        /* reads n, n's ID, sizefacbase, and maxbig */
        /* builds RESTART, prints restartvalues, and stops */
 
        printf("start up, create restart file in this directory\n");
        printf("give following data:\n");
        printf("        number to be factored,\n");
        printf("        its unique identification,\n");
        printf("        size of factor base,\n");
        printf("        and bound for largest prime in partials.\n");
        fflush(stdout);
 
        zread(n);
        scanf("%ld",&identification);
        if ((n[1] & 1) == 0) {
                printf("input is even, exit\n"); exit(0);
        }
 
        scanf("%ld %ld",&sizefacbase,&maxbig);
        printf("begin computing restart data for factorization of\n");
        zwrite(n); printf("\n");
        printf("identification: %ld\n",identification);
        printf("size factor base: %ld\n",sizefacbase);
        printf("primes in partial factorizations at most: %ld\n",maxbig);
        fflush(stdout);
 
        multiplier = findmultiplier((long)100);
        printf("multiplier: %ld\n",multiplier); fflush(stdout);
        zsmul(n,multiplier,kn); zsqrt(kn,sqkn);
        if ((kn[1] % 4) != 1) {
                printf("wrong multiplier selected, exit\n");
                exit(0);
        }
        if ((multiplier>1) && (multiplier <= MAXSMALLPRIME))
                smallprimes[multiplier] = multiplier;
        bytefacbase = (sizefacbase+1) * sizeof(long);
        if (((primbase = (long*)malloc(bytefacbase)) == NULL) ||
                ((sqrtkn = (long*)malloc(bytefacbase)) == NULL)) {
                printf("unable to get memory, exit\n");
                exit(0);
        }
 
        createfactorbase(); /* puts data in primbase, r1, and log info */
                        /* also sets lastprime */
        printf("factor base created\n"); fflush(stdout);
        if ((dummy = smalldiv()) > 0 ) {
                printf("small factor: %ld, exit\n",dummy);
                exit(0);
        }
        printf("largest prime in factor base: %ld\n",lastprime);
        if ((sievelength = 5*lastprime) > 400)
                sievelength = 400*((sievelength+200)/400);
        if (sievelength > MINSIEVE)
                sievelength = MAXSIEVE*((sievelength+MAXSIEVE/2)/MAXSIEVE);
        printf("sieve length: %ld\n",sievelength);
 
        findsievebound();
        printf("sievebound: %ld\n",(long)sievebound);
 
        zsdiv(kn,(long)2,mem1);                zsqrt(mem1,mem2);
        zintoz(sievelength,mem1);        zdiv(mem2,mem1,mem1,mem2);
        zsqrt(mem1,startd);                startd[1] |= 3;
        if ((startd[0] > 1 ) || (startd[1] > 4))
		zsadd(startd,(long)-4,startd);
        printf("polynomial generator starts at: ");
                zwrite(startd);
        printf("\n");
 
        if ((restart = fopen(restartfn,"w")) == NULL) {
                printf("Could not open restart file, exit\n");
                exit(0);
        }
        if (dumpdata()) {
                printf("Could not write restart file, exit\n");
                unlink(restartfn);
                exit(0);
        }
        printf("restart file created, done\n");
        printf("\n");
        printf("For a firefly mpqs restart do the following:\n");
        printf("- copy the file RESTART to your home directory ");
        printf("on jumbo, bigtop, srcf0d and srcf34.\n");
        printf("- give the i-th processor the following input:\n");
        printf("--------- ><8 --------- \n");
        printf("%ld %ld i\n", FIREFLY, identification);
        printf("--------- ><8 --------- \n");
        printf("The i should be unique per processor to guarantee that\n");
        printf("different processors process different curves/intervals,\n");
        printf("the %ld means it is a firefly restart, and\n",FIREFLY);
        printf("the %ld is the identification for n.\n",identification);
        printf("Fireflies write results to standard output;\n");
        printf("relations are on one line, begin with ful, par, or ppr,\n");
        printf("immediately followed by an identification number\n");
        printf("(unique per number to be factored) and a colon.\n");
        printf("\n");
        printf("Give other machines the following input:\n");
        printf("--------- ><8 --------- \n");
        printf("%ld %ld i\n",SLAVE,identification);
        zwrite(n); printf("\n");
        printf("%ld %ld %ld %ld %ld\n",multiplier,sizefacbase,
                        sievelength,(long)sievebound,maxbig);
        zwrite(startd);
        printf("\n0/1\n");
        printf("address\n");
        printf("\n");
        printf("--------- ><8 --------- \n");
        printf(
	  "The second to last line contains a zero or a one; a one if\n");
        printf("the results should be mailed to address (last line), a\n");
        printf("zero if the owner of that machine will mail the\n");
        printf(
	  "results. The mpqs-process will create a file RESTART which\n");
        printf("can be used by other processes running in the same\n");
        printf("directory; this file will not be created if it already\n");
        printf("exists for this n (with the same identification).\n");
        printf(" The i on the first line should be an integer and at\n");
        printf("least zero. It should be different from any i given to\n");
        printf("the Fireflies, and different from any i given to other\n");
        printf("processes. Never use any i more than once, even on the\n");
        printf("same machine, to avoid duplicating curves/relations.\n");
        exit(0);
        }
#endif /* not VMS */

/* how is not STARTER, SLAVE or FIREFLY: */
        printf("wrong restart value, exit\n");
        exit(0);
}
 
findnextpolynomial(a,b,d2inv)
long a[], b[], d2inv[];
/* after trying all kinds of nice methods, like the method described */
/* in Carl Pomerance's paper, and Carl's incremental QS, I decided to go*/
/* back to the old straightforward way of generating polynomials. */
/* Reason: fancier methods are only slightly faster, produce fewer */
/* relations (?), and cause all kinds of unpredictable troubles on */
/* other machines (file manipulations, fseek, etc.) */
/* Another advantage of this way of generating polynomials is that */
/* it is much easier to provide different machines with different */
/* inputs that will not 'overlap', and that are not too far from */
/* the optimal values. */
/* Still I would like to have a method that is faster than this one, */
/* but with similar low memory requirements, and, like this method, */
/* without any file manipulations. */
{ long t1[ZSIZEP], t2[ZSIZEP], t3[ZZSIZEP], t2kn[ZZSIZEP], t4[ZZSIZEP],
      ro1[2], ro2[2], pp[2];
register long i, p, ap, bp, bpc;
long sodinv();
long rem1, mod1 = 1033715657;		/* 13*23*37*41*43*53 */
long rem2, mod2 = 1035563755;		/*  5*17*29*61*67*97 */
long rem3, mod3 = 1036405977;		/*  3*19*47*59*79*83 */
long rem4, mod4 = 1039058713;		/*  7*11*31*67*73*89 */
 
i = 1;
zcopy(startd,altstartd);
interrupted = 1;
startd[1] |= 3; /* to be sure */
rem1 = zsdiv(startd, mod1, t1); rem2 = zsdiv(startd, mod2, t1);
rem3 = zsdiv(startd, mod3, t1); rem4 = zsdiv(startd, mod4, t1);
while (i > 0) {
        zsadd(startd,(long)4,startd);
	rem1 += 4; rem2 += 4; rem3 += 4; rem4 += 4;
        if (    jacobi(rem1,mod1) == 0 || jacobi(rem2,mod2) == 0
             || jacobi(rem3,mod3) == 0 || jacobi(rem4,mod4) == 0) continue;
			/* Reject d if divisible by odd prime < 100 */
	zsdiv(startd,(long)4,t1);	/* t1 = (d-3)/4 */
	zexp(kn,t1,startd,t2); 		/* t2 = 1/sqrt(kn) mod d if d good */
	zmul(t2,kn,t2kn); 		/* t2kn = t2*kn = kn^((d+1)/4) mod d */
	zmul(t2,t2kn,t3); zdiv(t3,startd,t4,t1); /* kn^((d-1)/2) mod d */
        if ((t1[0] != 1) || (t1[1] != 1)) continue; /* nonresidue or nonprime*/
        if (zprime(startd, (long)1) > 0) i = 0;	 /* insurance */
        }
interrupted = 0;
zmod(t4,kn,d2inv); if (d2inv[1]&1) zadd(d2inv, kn, d2inv);
zsdiv(d2inv,(long)2,d2inv); zsub(kn,d2inv,d2inv); /* d2inv = 1/2d mod kn */

zmul(startd,startd,a); zmod(t2kn,a,t2kn); zmod(t3,a,t3);
/* Compute b = t2*kn*(t2^2*kn - 3)/2 = sqrt(kn) mod a */
zsadd(t3,(long)-3,t3);        zmulmod(t2kn,t3,a,b);
if (b[1] & 1) zadd(a,b,b);      zsdiv(b,(long)2,b); 
if ((b[1] & 1) == 0) zsub(a,b,b);		/* Odd square root */

zadd(sqkn,b,t1);        zdiv(t1,a,t1,t2);
zsdiv(t1,(long)2,t1);                roots[0].r1 = -t1[1];
zsub(sqkn,b,t1);        zdiv(t1,a,t1,t2);
zsdiv(t1,(long)2,t1);                roots[0].r2 = t1[1];
 
ro1[0] = 1; ro2[0] = 1; pp[0] = 1;
for (i=nbsmallpsp; i<=sizefacbase; i++) {
	pp[1] = (p = primbase[i]);
	bp = zsdiv(b,p,t1);
	ap = zsdiv(startd,p,t1);
	ap = zsmulmod(ap,ap,p);	/* a mod p (a = d^2) */
	if (ap != 0) {
                bpc = bp;
                if ((bp += sqrtkn[i]) >= p) bp -= p;
                ap = sodinv(2*ap,p);
                roots[i].r1 = zsmulmod(p-ap, bp, p);
                if ((bp = sqrtkn[i]-bpc) < 0) bp += p;
                roots[i].r2 = zsmulmod(ap, bp, p);
        } else {
                zmul(b,b,t1);
                if (zcompare(kn,t1) >= 0) zsub(kn,t1,t1);
                else { zsub(t1,kn,t1); bp = p-bp; }
                zsdiv(t1,(long)4,t1);                zdiv(t1,a,t2,t1);
                ro1[1] = zsdiv(t2,p,t1);
                ro2[1] = sodinv(bp,p);        zmulmod(ro1,ro2,pp,ro1);
                roots[i].r2 = (roots[i].r1 = ro1[1]);
	}
}
}
 
checkroot(a,b)
/* This one could be done much nicer, but it is not important */
/* if CHECK has a reasonably large value */
long a[], b[];
{ long i, misser, r[ZSIZEP], c[ZSIZEP], p, minc;
  long zsmulmod(), zsdiv(), sodinv();
misser = 0;
for (i=1; i<=sizefacbase; i++) {
	p = primbase[i];     if (p < MAXSMALLPRIME) p = smallprimes[p];
        if (zsdiv(kn,p,r) != zsmulmod(sqrtkn[i], sqrtkn[i], p))
                { printf("Wrong at primbase[%ld]=%ld\n",i,p);
                fflush(stdout);
                misser ++; }
 
        }
if (misser == 0) printf("squareroots correct\n");
else printf("%ld mistakes in the squareroots\n",misser);
fflush(stdout);
 
zmul(b,b,r);
if (zcompare(kn,r) >= 0) { zsub(kn,r,r); minc = -1; }
else { zsub(r,kn,r); minc = 1; printf("b is very large!\n"); }
if (zsdiv(r,(long)4,r) != 0)  {
        printf("kn-b^2 is not 0 modulo 4, error, exit\n");
        exit(0);
        }
zdiv(r,a,c,r);
if ((r[0] != 1 ) || (r[1] != 0)) {
        printf("kn-^2 not divisible by a, error, exit\n");
        exit(0);
        }
 
misser = 0;
for (i=nbsmallpsp; i<=sizefacbase; i++) {
	long ap, bp, cp, res1, res2;
	long root1 = roots[i].r1, root2 = roots[i].r2;
        p = primbase[i];        if (p <= MAXSMALLPRIME) p = smallprimes[p];
        ap = zsdiv(a,p,r);
        bp = zsdiv(b,p,r);
	cp = minc*zsdiv(c,p,r);
	res1 = zsmulmod(ap, root1, p) + bp;
	res2 = zsmulmod(ap, root2, p) + bp;
	if (res1 >= p) res1 -= p;
	if (res2 >= p) res2 -= p;
	res1 = zsmulmod(res1, root1, p) + cp;
	res2 = zsmulmod(res2, root2, p) + cp;
        if (res1 != 0 && res1 != p) {
                misser ++;
                printf("wrong first root at i=%ld p=%ld\n",i,p);
                printf("ap=%ld, bp=%ld, cp=%ld, root=%ld\n", ap,bp,cp,root1);
                printf("sqrtkn mod p = %ld, (2a)^-1 mod p = %ld\n",sqrtkn[i],
                        sodinv(2*ap,p));
                fflush(stdout);
                }
        if (root1 == root2) printf("double root at p=%ld\n",p);
        if (res2 != 0 && res2 != p) {
                misser ++;
                printf("wrong second root at i=%ld p=%ld\n",i,p);
                printf("ap=%ld, bp=%ld, cp=%ld, root=%ld\n",ap,bp,cp,root2);
                printf("sqrtkn mod p = %ld, (2a)^-1 mod p = %ld\n",
                        sqrtkn[i],sodinv(2*ap,p));
                fflush(stdout);
                }
        }
if (misser == 0) printf("roots correct\n");
else {
        printf("%ld mistakes in the roots, exit\n",misser);
        exit(0);
        }
fflush(stdout);
}
 
#ifndef RISC
#ifndef ASM_sieveinit
sieveinit(lng, inval) /* use that sizeof(long) = 4 */
long lng, inval;
{
register long val = inval, i = (lng + 3) >> 2, *psiev = (long*)sieve;
while (i >= 100) {	/* Do 100 longwords at at time */
	i -= 100;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val;
        }
while (i > 0) {*psiev++ = val; i--;}
}
#endif /*not RISC*/
#endif /* ASM_sieveinit */

double todoub(a) long a[]; { /* of length at most 2 */
	if (a[0]==2) return((double)RADIX*(double)a[2]+(double)a[1]);
	if (a[0]==1) return((double)a[1]);
	return((double)0);
}
 
long gcdbnd(n,m,lw) long n[],m[],lw; {
	/* return 0 if the gcd <= lw, where lw < RADIX */
	/* return 0 if gcd >= RADIX, return gcd otherwise */
	/* only if n[0] and m[0] <=2, n>m, and n odd */
	long a[3], b[3], c[3];
	a[0]=n[0];a[1]=n[1];a[2]=n[2];
	b[0]=m[0];b[1]=m[1];b[2]=m[2];
	while ((b[1]&1)==0) {
		b[1] >>= 1;
		if (b[0]>1) {
			if (b[2]&1) b[1] += RADIXH;
			if ((b[2] >>= 1) == 0) b[0] = 1;
		}
	} /* so that b is odd */
	if (b[0] == 1 && b[1] <= lw) return(0);
	zsub(a,b,c);
	while (c[1]!=0 || c[0] != 1) {
		do {
			c[1] >>= 1;
			if (c[0]>1) {
				if (c[2]&1) c[1] += RADIXH;
				if ((c[2] >>= 1) == 0) c[0] = 1;
			}
		} while ((c[1]&1) == 0);
		if (c[0] == 1 && c[1] <= lw) return(0);
		if (zcompare(b,c)>=0) {
			a[0] = b[0]; a[1] = b[1]; a[2] = b[2];
			b[0] = c[0]; b[1] = c[1]; b[2] = c[2];
		} else {
			a[0] = c[0]; a[1] = c[1]; a[2] = c[2];
		}
		zsub(a,b,c);
	}
	if (a[0] > 1) return(0);
	return(a[1]);
}
 
long pollardrho(n,res,lw) long n[], res[],lw; {
	register long i=0,j=100,adder=1;
	long xi[ZSIZEP], x2i[ZSIZEP], dif[ZSIZEP];
	zmstart(n);
	zintoz(1,dif);
	zintoz(1,xi);
	ztom(xi,n,xi); zcopy(xi,x2i);
	while (i<100000) {
		for (;i<j;i++) {
			zmont(xi,xi,n,xi); zsadd(xi,adder,xi);
			zmont(x2i,x2i,n,x2i); zsadd(x2i,adder,x2i);
			zmont(x2i,x2i,n,x2i); zsadd(x2i,adder,x2i);
			if (zcompare(xi,x2i)>=0) zsub(xi,x2i,res);
			else zsub(x2i,xi,res);
			if ((res[1]!=0)||(res[0]!=1)) zmont(dif,res,n,dif);
			else {
				zintoz(1,dif);
				zintoz(1,xi);
				ztom(xi,n,xi); zcopy(xi,x2i);
				adder += 2;
			}
		}
	if ((dif[1]!=0)||(dif[0]!=1)) {
                        if (res[1]=gcdbnd(n,dif,lw)) {
                                res[0] = 1;
				zdiv(n,res,n,xi);
				if (xi[1]!=0 || xi[0]!=1) return 0;
				if (zcompare(res,n)>0) {
					zcopy(res,dif);
					zcopy(n,res);
					zcopy(dif,n);
				}
				return i;
			}
		} else zintoz(1,dif);
		j=i+60;
	}
	return(0);
}
 
dumpfac(w,wid,widc,ind,fac,bp,bp2,qi)
long *w, ind, fac[], bp, bp2,qi[];
char wid[], widc[]; {
	register long i;
        int status;
 
	if (how == FIREFLY) {
		mailit = stdout; closeit = 0;
	} else if ((mailit = fopen(nameoutput,"a+")) == NULL) {
		printf("Cannot open file %s.\n",nameoutput);
		printf("results to standard output\n");
		mailit = stdout; closeit = 0;
	}
	
    	(*w) ++;
	if (how != FIREFLY )
	fprintf(mailit,"%c%c%c %ld from %d\n",wid[0],wid[1],wid[2],*w,unique);
	fprintf(mailit,"%s%ld %ld %ld",wid,bp,bp2,ind);
    	for (i=1;i<=ind;i++) {
            	if (how != FIREFLY)
            	if (i%10 ==0) fprintf(mailit,"\n%s",widc);
            	fprintf(mailit," %ld",fac[i]);
	}
    	if (how != FIREFLY) {
		fprintf(mailit,"\n%s 0 ",widc);
    		fzwrite(mailit,qi,LINELENGTH,nothing,widc);
	} else {
		fprintf(mailit," 0 ");
		fzwrite(mailit,qi,1000,nothing,nothing);
	}
	
	fprintf(mailit,"\n");
 	fflush(mailit);
	if (closeit) fclose(mailit);
	if (++outputcounter == mail_batchcount) {
		outputcounter = 0;
		send_mpqs_mail();
	}
	numfac += ind/2;
}
 
sievget(a,b,d2inv,sign,i,rshift,value)
long a[],b[],d2inv[],sign,i,rshift;	/* sign no longer used */
uchar value;
{
    long factors[MAXFACNB], qisq[ZSIZEP], t1[ZSIZEP], t2[ZSIZEP],
        index = 1;
/*
	Set P(X) = ((2*a*x + b)^2 - kn)/4*d^2, where a = d^2
	(see findnextpolynomial).  sievget is called when
	siever detects that P(i) is highly composite.
	This routine computes P(i) explicitly and possibly
	outputs a relation containing qisq = (2*a*i + b)/2d
	and the factors of P(i) == qisq mod kn.
*/

    closeit = 1;
    zsmul(a,2*(i > 0 ? i : -i),t1);
    if (i > 0) zadd(t1,b,t1);
    else if (zcompare(t1,b) >= 0) zsub(t1,b,t1);
    else zsub(b,t1,t1);		/* t1 = +- (2*a*i + b) */
    zmulmod(t1,d2inv,kn,qisq);	
    zmulmod(qisq,qisq,kn,t1);
    if (i>=roots[0].r1 && i<=roots[0].r2)
        { factors[1]= -1; zsub(kn,t1,t1); }
    else
	factors[1]=1;
    {
	register long j;
	double contrib = 0.0;

	for (j=1;j<=nbsmallps;j++) {
		long pj = primbase[j];
		if (zsdiv(t1,pj,t2) == 0) {
			long cntr = 1;
			zcopy(t2,t1); cntr=1;
			while (zsdiv(t1,pj,t2) == 0) {
				cntr ++; zcopy(t2,t1);
			}
			contrib += (double)cntr*smlogp[j];
			factors[++index] = cntr;
			factors[++index] = j;
		}
	} /* for j */
	if ((double)value+contrib < (double)sievebound) return;
    }
    {
	register long j;
	register long offset = i - rshift;	/* negative */
	register roots_t *loc_roots = roots;
	register row      loc_primbase = primbase;
	register long     loc_sizefacbase = sizefacbase;

	for (j = nbsmallpsp; j <= loc_sizefacbase; j++) {
		register long imodp = offset;
		register long pj = loc_primbase[j];
/*
		Test whether roots[j].r1 - offset or
		roots[j].r1 + roots[j].r2 - offset
		is divisible by pj.  Recall offset < 0 and
		0 <= roots[j].r1 <= roots[j].r1 + roots[j].r2 < p .
		Ensure that first argument to % is nonnegative.	
*/
		if (imodp <= -(pj << 3)) {
			imodp = pj - 1 - (-1 - imodp)%pj;
		} else {	/* quotient at most 8, often much less */
			do {} while ((imodp += pj) < 0);
		}
		imodp -= loc_roots[j].r1;
        	if (imodp == 0 || imodp == loc_roots[j].r2) {
                	if (zsdiv(t1,pj,t2) == 0) {
				long cntr = 1;
				zcopy(t2,t1);
				while (zsdiv(t1,pj,t2) == 0) {
					cntr ++; zcopy(t2,t1);
				}
				if (index >= MAXFACNB-2) {
				    printf("Increase MAXFACNB\n");
				    return; 
				}
				factors[++index] = cntr;
				factors[++index] = j;
			} else {
				printf("Sievget - bad root, p = %ld\n", pj);
			}
		} /* if imodp */
	} /* for j */
    }
    {
	long ival= (long)value;
	double t1value;
	int goodluck;

	if (t1[0] > 2) {
		goodluck = 0;		/* Cofactor too large */
	} else if (t1[0] == 1 && t1[1] <= maxbig) {
		long pj = t1[1];
		if (pj == 1) {
			dumpfac(&ful,fulid,fulidc,index,factors,1,1,qisq);
                	if (ival>maxful) maxful = ival;
               		if (ival<minful) minful = ival;
               		sumful += ival;
		} else { /* pj */
			/* believe pj is prime, even */
			/* if it isn't. who cares, it */
			/* doesn't matter */
 
			dumpfac(&par,parid,paridc,index,factors,1,pj,qisq);
               		if (ival>maxpar) maxpar = ival;
               		if (ival<minpar) minpar = ival;
               		sumpar += ival;
		}
		goodluck = 1;
	} else if (   (t1value = todoub(t1))
		    < (double)lastprime*(double)lastprime) {
		goodluck = 0;		/* Cofactor must be prime */
	} else if (t1value > (double)maxbig*(double)maxbig) {
		goodluck = 0;		/* Cannot have two factors < maxbig */
	} else if (zprime(t1,(long)1)) { /* then we believe t1 is prime */
		goodluck = 0;
	} else { /* now t1 has precisely two factors */
			/* is it pp or worthless? */
			/* let's first try to factor it */
		if (pollardrho(t1,t2,lastprime)) {
			/* biggest in t1, smallest in t2 */
			if (   t1[0]>1 || t2[0]>1
			    || t1[1]>=maxbig || t2[1]>=maxbig) {
				goodluck = 0;
			} else {	/* now they are both small enough */
		    		if (t2[1] <= lastprime) {
					printf("Sievget - missed f = %ld\n",
						t2[1]);
				}
				/* the pp case */
				dumpfac(&ppr,pprid,ppridc,index,
					factors,t2[1],t1[1],qisq);
                		if (ival>maxppr) maxppr = ival;
				if (ival<minppr) minppr = ival;
				sumppr += ival;
				goodluck = 1;
			}
		} else {	/* pollardrho */
			printf("Missed pollardrho factorization: ");
			zwrite(t1); printf("\n"); fflush(stdout);
			goodluck = 0;
		}
	} /* end of ifs over t1 */
 
	if (!goodluck) {
#if MAXPOL <= 100
		printf("                 unfactored quotient %3ld: ", ival);
		zwrite(t1); printf("\n"); fflush(stdout);
#endif
		bad ++;
		if (ival>maxbad) maxbad = ival;
		if (ival<minbad) minbad = ival;
		sumbad += ival;
	}
    }
} /* sievget */
 
#ifdef ASM_first_ge
long first_ge();
#else
long first_ge(array, start, lower_bound)
/* Return index of first element in array which is >= lower_bound. */
/* Should be optimized for case where very few elements are >= lower_bound. */
long start;
uchar array[], lower_bound;
{
    register uchar *parray = &array[start];
    register uchar  bound  = lower_bound;
    if (sizeof(long) == 4*sizeof(uchar) & FIRST_GE_SENTINEL >= 4) {
/* Advance to longword boundary (N.B. array was allocated by malloc) */
	do {
	    if (*parray >= bound) return (parray-array);
        } while ((++parray - array)&3);
/* Check four bytes at a time. */
/* If something interesting found, complete search one byte at a time. */
/* This attempts to cater to both big-endian and little-endian machines. */
	{ register long long_bound = (0x01010101)*(bound - 1);
	  register long *plarray = (long*)(parray);
 
	  while (
	    (((long_bound - plarray[0])
	    | (long_bound - plarray[1])
	    | (long_bound - plarray[2])
	    | (long_bound - plarray[3]))
	    & 0x80808080) == 0)
		plarray += 4;
 
	  parray = (uchar*)(plarray);
/*
	printf("first %ld\n",(long)(parray-1-array));
	  printf("plarray = %lx, parray = %x %x %x %x, bound=%x\n",
	    plarray[-1], parray[0], parray[1], parray[2], parray[3],bound);
	  fflush(stdout);
*/
	}
    }
    while (*parray++ < bound) {};	/* null body */
/*
	printf("final %ld\n",(long)(parray-1-array));
	  fflush(stdout);
*/
 
    return (parray-1-array);
}
#endif 	/* ASM_first_ge */
 
/* sievget_check checks for elements of sieve which are >= smallsievebnd */
#define sievget_check(direction,shifter) {i = -1; \
	sieve[rsl] = smallsievebnd; \
        while ((i = first_ge(sieve, i+1, smallsievebnd)) < rsl) {\
            sievget(a,b,d2inv,(long)(direction),i+ishift,\
			(long)(shifter),sieve[i]);}}
 
/*
	Procedures bigp_siever, lowp_siever, medp_siever
	all have four arguments (ibeg, iend, rsl, logp).
	Each sieves for all primes p = primbase[i],
	where ibeg <= i <= iend, by adding logp
	to all appropriate entries of the sieve.
	All assume the new roots are roots[i].r1 and
	roots[i].r1 + roots[i].r2, where

		0 <= roots[i].r1 <= roots[i].r1 + roots[i].r2 <= p-1

	The difference is the range of values of p
	assumed by each procedure:

		bigp_siever -- p >= rsl          (at most one hit per root)
		lowp_siever -- p <= rsl/2        (at least two hits per root)
		medp_siever -- rsl/2 <= p <= rsl (one or two hits per root)
*/


#ifndef ASM_bigp_siever
bigp_siever(ibeg, iend, rsl, logp)
long ibeg;	/* lower bound on i for sieving */
long iend;	/* upper bound on i for sieving */
long rsl;	/* width of sieve interval */
uchar logp;	/* logarithm of these primes */
{
    register long loc_rsl = rsl, i;
    register uchar loc_logp = logp;
    register urow sieve_rsl = sieve + loc_rsl;
    register roots_t *rootsi = roots + ibeg;
    register row pi = primbase + ibeg;

    for (i = ibeg - iend - 1; i < 0; i++) {
	register long root1 = rootsi->r1 - loc_rsl;
				/* root relative to end of sieve interval */
#if 0
	if (   rootsi->r1< 0 || rootsi->r2 < 0
	    || rootsi->r1 + rootsi->r2 >= *pi || *pi < loc_rsl) {
	    printf("Bad r1 = %ld, r2 = %ld, bigp = %ld\n",
		rootsi->r1, rootsi->r2, *pi);
	}
#endif
	if (root1 < 0) {	
	    register long root2 = root1 + rootsi->r2;

	    sieve_rsl[root1] += loc_logp;
	    root1 += *pi;	     /* Advance smaller root by p */
	    if (root2 < 0) {
		sieve_rsl[root2] += loc_logp;
					     /* Difference unchanged */
	    } else {
		rootsi->r2 = root1 - root2;  /* New difference */
		root1 = root2;		     /* New smaller root */
	    }
	}
	rootsi->r1 = root1;
	rootsi++;
	pi++;
    }	/* for i */
}
#endif	/* ASM_bigp_siever */

#ifndef ASM_lowp_siever
lowp_siever(ibeg, iend, rsl, logp)
long ibeg;	/* lower bound on i for sieving */
long iend;	/* upper bound on i for sieving */
long rsl;	/* width of sieve interval */
uchar logp;	/* logarithm of these primes */
{
             long loc_rsl = rsl;
    register long i;
    register urow sieve_rsl = sieve + loc_rsl, sieve1;
    register uchar loc_logp = logp;
    register roots_t *rootsi = roots + ibeg;
    register row pi = primbase + ibeg;

    for (i = ibeg - iend - 1; i < 0; i++) {
	register long rootdif = rootsi->r2;
	register long root2 = rootsi->r1 - loc_rsl + rootdif;
			/* larger root relative to end of sieve interval. */
			/* Will always be negative initially */
	register long p = *pi++;

#if 0
	if (   rootsi->r1 < 0 || rootdif < 0
	    || rootsi->r1 + rootdif >= p || p+p > loc_rsl) {
	    printf("Bad r1 = %ld, r2 = %ld, lowp = %ld\n",
			rootsi->r1, rootdif, p);
	}
#endif
 
	sieve1 = sieve_rsl - rootdif;
			/* Note  sieve1  = sieve + (rsl-p) + (p-rootdif)
					>= sieve,
			   as required by ANSI C. */
	do {
	    sieve1[root2] += loc_logp; sieve_rsl[root2] += loc_logp;
	    root2 += p;	
	    sieve1[root2] += loc_logp; sieve_rsl[root2] += loc_logp;
	    root2 += p;	
	} while (root2 < -p);
	if (root2 < 0) {
	    sieve1[root2] += loc_logp; sieve_rsl[root2] += loc_logp;
	    root2 += p;	
	}	
	root2 -= rootdif;		/* Replace root2 by smaller root */
	if (root2 < 0) {		/* Extra hit for smaller root */
	    sieve_rsl[root2] += loc_logp;
	    root2 += rootdif;		/* Restore previous value */
	    rootsi->r2 = p - rootdif;  	/* New difference */
	}
	rootsi->r1 = root2;		/* New smaller root */
	rootsi++;
    }	/* for i */
}
#endif	/* ASM_lowp_siever */

#ifndef ASM_medp_siever
medp_siever(ibeg, iend, rsl, logp)
long ibeg;	/* lower bound on i for sieving */
long iend;	/* upper bound on i for sieving */
long rsl;	/* width of sieve interval */
uchar logp;	/* logarithm of these primes */
{
    register long loc_rsl = rsl, i;
    register uchar loc_logp = logp;
    register urow sieve_rsl = sieve + loc_rsl;
    register roots_t *rootsi = roots + ibeg;
    register row pi = primbase + ibeg;

    for (i = ibeg - iend - 1; i < 0; i++) {
	register long root1 = rootsi->r1 - loc_rsl;
	register long root2 = root1 + rootsi->r2;
	register long p = *pi++;
				/* roots relative to end of sieve interval */
/*
		This code resembles that in bigp_siever, except we
		sieve at least once per root.
*/

#if 0
	if (   rootsi->r1 < 0 || rootsi->r2 < 0
	    || rootsi->r1 + rootsi->r2 >= p
	    || p+p < loc_rsl || p > loc_rsl) {
	    printf("Bad r1 = %ld, r2 = %ld, medp = %ld\n",
		rootsi->r1, rootsi->r2, p);
	}
#endif
	sieve_rsl[root1] += loc_logp;
	sieve_rsl[root2] += loc_logp;
	root1 += p;
	if (root1 < 0) {	
	    sieve_rsl[root1] += loc_logp;
	    root1 += p;	     		/* Advance smaller root by p */
	    root2 += p;	
	    if (root2 < 0) {
		sieve_rsl[root2] += loc_logp;
					     /* Difference unchanged */
	    } else {
		rootsi->r2 = root1 - root2;  /* New difference */
		root1 = root2;		     /* New smaller root */
	    }
	}
	rootsi->r1 = root1;
	rootsi++;
    }	/* for i */
}
#endif /* ASM_medp_siever */

long prim_atmost(p)
long p;
{
/*
	Return the largest index i such that primbase[i] <= p.
	Use a binary search.  We treat primbase[0] as -infinity
	and primbase[sizefacbase+1] as +infinity.
*/
    register long mn = 0, mx = sizefacbase;
    while (mx != mn) {
		/* primbase[mn] <= p < primbase[mx+1] */
	register long mid = (mn + mx + 1) >> 1;
	if (primbase[mid] > p) {
	    mx = mid - 1;
	} else {
	    mn = mid;	
	}
    } /* while */
    return mn;
}

#ifndef ASM_sieveroot_adjust
sieveroot_adjust(amt)
long amt;	/* assumed nonnegative */
{
/*
	findnextpolynomial stored the roots of a*x^2 + b*x + c mod p
	into r1 and r2, in arbitrary order.  That is also the
	convention used in checkroot.  Add (amt % p) to these
	roots, keeping them in [0, p-1).  Then store the smaller root
	in r1 and the difference in r2 (a convention used in
	bigp_siever, lowp_siever, medp_siever, and sievget).

*/
    register row          pi       = &primbase[nbsmallpsp];
    register roots_t *rootsi       =    &roots[nbsmallpsp];
    register roots_t *rootsi_bound =    &roots[sizefacbase];

    for (; rootsi <= rootsi_bound; rootsi++) {
        register long p     = *pi++;
        register long srem  = amt % p;
        register long root1 = rootsi->r1 + srem;
        register long root2 = rootsi->r2 + srem;

        if (root1 >= p) root1 -= p;	/* Ensure in [0, p-1) */
        if (root2 >= p) root2 -= p;

        if (root1 < root2) {	
	    rootsi->r1 = root1;		/* smaller root */
	    rootsi->r2 = root2-root1;	/* difference */
        } else {
	    rootsi->r1 = root2;
	    rootsi->r2 = root1-root2;
        }
    } /* for rootsi */
}
#endif	/* ASM_sieveroot_adjust */

siever(a,b,d2inv)
long a[], b[], d2inv[];
{
register long i;
long value = (primbase[1] == 2 ? BYTSOF1 : 0);
long ishift;
long logpstart[30];	/* Should be global */
uchar lastlog;		/* Should be global */
long rsl_prev = 0, rsl_crossing, rsl2_crossing;
	/* primbase[rsl_crossing ] <= rsl   <= primbase[rsl_crossing+1 ] */
	/* primbase[rsl2_crossing] <= rsl/2 <= primbase[rsl2_crossing+1] */

/* Find where the logarithm changes.  This processing should be moved
   elsewhere (with lastlog, logpstart being made global)
   rather than be repeated for each polynomial.
*/
{
    register uchar logp = firstlog;

    logpstart[logp] = nbsmallpsp;
    while (1) {
	long logbr = nextlog[logp-firstlog];
	if (logbr > primbase[sizefacbase]) break;
	logpstart[++logp] = prim_atmost(logbr-1) + 1;
    }
    lastlog = logp;
    logpstart[lastlog+1] = sizefacbase + 1;
}

sieveroot_adjust(sievelength);	/* Add sievelength %p to roots */
				/* (negative of initial ishift) */

for (ishift = -sievelength; ishift < sievelength; ishift += sievesize) {
	long  rshift =  (   ishift + sievesize < sievelength
			  ? ishift + sievesize : sievelength);
				/* right end of interval, plus 1 */
	long  rsl = rshift - ishift;
				/* width of current sieve interval */
	register uchar logp;
#ifdef RISC
        for (i=0;i<rsl;i++) sieve[i] = value;
#else
        sieveinit(rsl,value);
#endif

	if (rsl != rsl_prev) {	/* Update rsl_crossing, rsl2_crossing  */
	    rsl_prev = rsl;
	    rsl_crossing  = prim_atmost(rsl     );
	    rsl2_crossing = prim_atmost(rsl >> 1);
	}

        for (logp = firstlog; logp <= lastlog; logp++) {
	    register long ibeg = logpstart[logp];
	    long iendbig = logpstart[logp+1] - 1;
	    long iendlow = (iendbig <= rsl2_crossing ? iendbig : rsl2_crossing);
	    long iendmed = (iendbig <= rsl_crossing  ? iendbig : rsl_crossing );

	    /* Sieve over primes  p <= rsl/2 with this logarithm */

	    if (ibeg <= iendlow) {
		lowp_siever(ibeg, iendlow, rsl, logp);
		ibeg = iendlow + 1;
	    }

	    /* Sieve over primes p with rsl/2 < p <= rsl and this logarithm */

	    if (ibeg <= iendmed) {
		medp_siever(ibeg, iendmed, rsl, logp);
		ibeg = iendmed + 1;
	    }

	    /* Sieve over primes p > rsl with this logarithm */

	    if (ibeg <= iendbig) {
		bigp_siever(ibeg, iendbig, rsl, logp);
	    }
        }   /* for logp */
        sievget_check(1, rshift);
        } /* for ishift */
}

report()
{
/*	
	This report is produced at intervals specified by REPORT.
	The output has one long line, resembling:

	S40; F1=137; P2(126,113); p4(105,102.25,100);
	B148(118,102.82,99); L13.86;12765946021730812111390339

	S40  		    - 40 polynomials have been sieved.
	F=137 		    - 1 full found by sievget, with ival = 137
			      (ival is total logarithm from sieving).
	P2(126,113) 	    - Two partials found, with ival = 126 and 113.
	p4(105,102.25,100)  - Four double partials have been found,
		  	      with largest ival = 105, smallest = 100,
			      average = 102.25.
	B148(118,102.82,99) - 148 false reports to sievget occurred,
		  	      with largest ival = 118, smallest = 99,
			      average = 102.82.
	L13.86 		    - Average value of index on dumpfac call
			      (divisors from factor base per relation) = 13.86.
	12765946021730812111390339 - Latest startd for findnextpolynomial.
*/

printf("S%ld",pol);
if (ful > 0) {
        if (ful == 1) printf("; F1=%ld",maxful);
        else if (ful == 2) printf("; F2(%ld,%ld)",maxful,minful);
        else printf("; F%ld(%ld,%.2f,%ld)",ful,maxful,
                        (double)sumful/ful,minful);
        }
if (par > 0) {
        if (par == 1) printf("; P1=%ld",maxpar);
        else if (par == 2) printf("; P2(%ld,%ld)",maxpar,minpar);
        else printf("; P%ld(%ld,%.2f,%ld)",par,maxpar,
                        (double)sumpar/par,minpar);
        }
if (ppr > 0) {
        if (ppr == 1) printf("; p1=%ld",maxppr);
        else if (ppr == 2) printf("; p2(%ld,%ld)",maxppr,minppr);
        else printf("; p%ld(%ld,%.2f,%ld)",ppr,maxppr,
                        (double)sumppr/ppr,minppr);
        }
if (bad > 0) {
        if (bad == 1) printf("; B1=%ld",maxbad);
        else if (bad == 2) printf("; B2(%ld,%ld)",maxbad,minbad);
        else printf("; B%ld(%ld,%.2f,%ld)",bad,maxbad,
                        (double)sumbad/bad,minbad);
        }
if ((ful+par+ppr) > 0)
	printf("; L%.2f",(double)numfac/(ful+par+ppr));
if (started > 0) {
        printf(";");
        if (interrupted > 0)
                zwrite(altstartd);
        else
                zwrite(startd);
        }
printf("\n");
fflush(stdout);
}
 
void csaver(sig) {
    if (outputcounter &&( how == SLAVE )) {
	if (mail_batchcount) {
	    printf("mailing...\n");
	    send_mpqs_mail();
	}
	else printf("please mail %s to the factor account\n",nameoutput);
    }
    printf("Exit: ");
    report();
    exit(0);
}
 
send_mpqs_mail()
{
    int status;

    printf("executing mail command '%s'\n", mail_execstring);
    status = system(mail_execstring);
#ifdef VMS
    if ((status & 1) == 0)
#else
    if (status)
#endif
	printf("mail failed '%s' (status = %d)\n", mail_execstring, status);
    else	    /* delete only if mail worked */
	unlink(nameoutput);
    fflush(stdout);
}

#ifdef VMS
#include <descrip.h>
#include <lib$routines.h>

unlink(filename)
    char* filename;
{
    int status;
    struct  dsc$descriptor_s delete_desc;
    char file_spec[100];

    strcpy(file_spec, filename);
    strcat(file_spec, ".;");
    delete_desc.dsc$b_dtype  = DSC$K_DTYPE_T;
    delete_desc.dsc$b_class  = DSC$K_CLASS_S;
    delete_desc.dsc$a_pointer = file_spec;
    delete_desc.dsc$w_length = strlen(file_spec);

    status = lib$delete_file(&delete_desc);
    if ((status & 01) != 1)
    {
	printf("VMS delete failed 'lib$delete_file(%s)' (status = %d)\n",
	       file_spec, status);
	fflush(stdout);
    }
}
#endif /*VMS*/


void setupsignals () {
        signal(SIGTERM, csaver);
        signal(SIGINT, csaver);
}
 
iterate()
{
long a[ZSIZEP], b[ZSIZEP], d2inv[ZSIZEP];
 
sievesize = MAXSIEVE;
if (sievesize > sievelength) sievesize = sievelength;
if (sievesize < 400)  {
        if ((sieve = (uchar*)malloc(luchar*(400+FIRST_GE_SENTINEL)))==NULL) {
                printf("unable to get memory for sieve, exit\n");
                exit(0);
                }
        }
else {
        while ((sieve=(uchar*)malloc(
			luchar*(sievesize+FIRST_GE_SENTINEL)))==NULL) {
                if ((sievesize -= MINSIEVE) < 0) {
                      printf("unable to get memory for sieve, exit\n");
                        exit(0);
                        }
                }
        }
 
for (pol=1; pol<=MAXPOL; pol++)  {
        findnextpolynomial(a,b,d2inv);
        if (how == SLAVE)
                if (((pol%CHECK) == 0) /* || (pol == 1 ) */ ) checkroot(a,b);
        siever(a,b,d2inv);
        if ((pol%REPORT) == 0) report();
        }
csaver(0);
}
 

main (argc, argv) char *argv[]; int argc; {
       time_t nowtime = time(NULL); 
#ifdef SYSV
	setpgrp2(0,0);
#endif
        if (argc >= 4) {
            restartfn = argv[1];
            if (strcmp(argv[2], "-") && !freopen(argv[2], "r", stdin)) {
		fprintf(stderr, "Can't open input file '%s'\n",  argv[2]);
		exit(1);
	    }
            if (!freopen(argv[3], "a", stdout)) {
		fprintf(stderr, "Can't open output file '%s'\n",  argv[3]);
		exit(1);
	    }
        }
	printf("%s: started %s\n", ctime(&nowtime), argv[0]);
        if (argc == 1)
            restartfn = default_restart;
        else {
            restartfn = argv[1];
            printf("Restart file is %s\n", restartfn);
        } 
        if (argc >= 5)
            file_prefix = argv[4];
        if (argc >= 6)
            mail_cmd = argv[5];
        fflush(stdout);
        setupsignals();
        scanf("%ld",&how);
        setup();
        iterate();
}

!End!of!file!
echo Extracting mpqssrc/xtime.c
cat >mpqssrc/xtime.c <<'!End!of!file!'
#include <X11/Xlib.h>
#include <X11/Xatom.h>
#include <stdio.h>

main() {
    XSetWindowAttributes xwa;
    Display *dpy = XOpenDisplay(NULL);
    xwa.event_mask = PropertyChangeMask;
    if (dpy) {
        XEvent ev;
        Window w = XCreateWindow(dpy,DefaultRootWindow(dpy),0,0,1,1,0,
            CopyFromParent,CopyFromParent,CopyFromParent,CWEventMask,&xwa);
        XChangeProperty(dpy,w,XA_PRIMARY,XA_ATOM,8,PropModeAppend,"",0);
        for (; XNextEvent(dpy,&ev), ev.xany.type != PropertyNotify; ) ;
        printf("%d\n",ev.xproperty.time);
    }
}
!End!of!file!
mkdir mutable
chmod 777 mutable
echo Extracting mutable/TASKLIST
cat >mutable/TASKLIST <<'!End!of!file!'
# This file should contain pairs of starting task and count
# Send mail to factor@src.dec.com asking for task# nnn to get some tasks

!End!of!file!
chmod 666 mutable/TASKLIST
echo Extracting mutable/TEMPLATE
cat >mutable/TEMPLATE <<'!End!of!file!'
2 2443 UID
1558841406043389331552847078537194140231730515145858751188381329215625\
20451700092630317221590417115548161
17 50000 6500000 88 100000000
74833927072460629735319
25
mail factor@src.dec.com < %s
!End!of!file!
echo Extracting pausefactor
cat >pausefactor <<'!End!of!file!'
#!/bin/sh

# pausefactor [n] causes factoring to pause on machine $MACH for n
# (default 3600) seconds.  To stun factoring more seriously, touch
# shortkill.$MACH, mediumkill.$MACH, or longkill.$MACH.

if test -z "${FPATH+yes}"
then
  FPATH=/bin:/usr/bin:/usr/ucb:/etc:${FACTORPATH-.}; export FPATH 
fi
PATH=$FPATH; export PATH

. environ

pause=$OUTDIR/.pause.$MACH; run=$OUTDIR/.run.$MACH

trap "rm -f $pause; $RESUMECMD $pid >/dev/null 2>&1; exit 0" 1
if test -s $run
then
  pid=`sed 1d $run`
  touch $pause
  $PAUSECMD $pid >/dev/null 2>&1
  sleep ${1-3600}
  $RESUMECMD $pid >/dev/null 2>&1
  rm -f $pause
fi
!End!of!file!
chmod 755 pausefactor
echo Extracting runfactor
cat >runfactor <<'!End!of!file!'
#! /bin/sh

# This shell script loops picking up tasks, touching the .run file so
# that checkrun won't start another copy, and firing off the mpqs program
# whenever a new copy is needed.  It also decides whether to pause
# factoring based on excessive load, terminal activity, and on killing
# factoring based on the existence of a {short,medium,long}kill file or
# a stop file.  Make sure that your FACTORPATH, and all variables named
# in environ have the right values before you start.

trap "" 1

if test -z "${FPATH+yes}"
then
  FPATH=/bin:/usr/bin:/usr/ucb:/etc:${FACTORPATH-.}; export FPATH 
fi
PATH=$FPATH; export PATH

. environ

# The goal of the next line is to change 2 to 200, 2.1 to 210, and 
# 2.103 to 210.  If you can't make this work, just get rid of these
# and use 210 if you want a pauseload of 2.10. 
pauseload=`echo 0 k 100 $pauseload '*' 1 / p | dc`

docheck=`witch atcheck`; tick=`witch ticktock`
run=$OUTDIR/.run.$MACH; dont=$OUTDIR/stop.$MACH
skill=$OUTDIR/shortkill.$MACH; mkill=$OUTDIR/mediumkill.$MACH;
lkill=$OUTDIR/longkill.$MACH; restart=$FACTORDIR/RESTART
pause=$OUTDIR/.pause.$MACH

rm -f $pause $skill $mkill $lkill

TICKSEC=`expr 60 '*' $TICKLEN`; export TICKSEC

if test -s $run
then
  OLDPID=`sed 1d $run`
  $RESUMECMD $OLDPID > /dev/null 2>&1
  kill $OLDPID > /dev/null 2>&1
fi

while test -r $template -a '!' -r $dont
do
  if test -s $run
  then
    oldout=`sed 1q $run`
    test -r $oldout && $MAILCMD $MAILTO < $oldout && rm -f $oldout
  fi

  teller=1

  . $tick

  domail="$DOMAIL"
  while num=`getnum`; test -z "$num" -a -r $template -a '!' -r $dont
  do
    test -n "$domail" && echo "Out of factoring tasks" | $MAILCMD $MAILUSER
    domail=""
    . $tick
  done

  rm -f $run
  echo $num | awk '{printf("'$OUTDIR'/OUT%07d\n", $1);}' > $run

  echo `date`: Starting mpqs with input $num >> $control
  if test -r $template
  then
    oldtemp="`sum $template 2>&1`"
    MAILER="$MAILCMD $MAILTO < `cat $run`"
    sed -e "s/UID/$num/" $template | \
        nice -19 mpqs $restart - $log $OUTDIR/ "$MAILER" &
    pid=$!; echo $pid >> $run; counter=${pausetty-1}
    shouldrun=1; isrun=$shouldrun; oldsum="`ls -lu $ttydevs | sum`"
  fi
  while 
    ps=`${PS} $pid | wc -l`
    if test  -r $template
    then
      newtemp="`sum $template 2>&1`"
    fi
    test $ps -eq 2 -a "$newtemp" = "$oldtemp" -a '!' -r $dont
  do
    . $tick
    why="No reason"
    shouldrun=1
    if test -f $pause
    then
      shouldrun=0
      why="$pause file exists"
    fi
    if test $pausetty -ne 0 -a $shouldrun -ne 0
    then
      newsum="`ls -lu $ttydevs | sum`"
      if test "$newsum" = "$oldsum"
      then
        shouldrun=`expr $counter '>=' $pausetty`
        if test shouldrun -eq 0
        then counter=`expr $counter + 1`
        fi
        if test $isrun -ne $shouldrun
        then why="Terminals idle for $pausetty tick-periods"
        fi
      else
        counter=0; oldsum="$newsum"; shouldrun=0
        if test $isrun -ne $shouldrun
        then why="Terminal activity"
        fi
      fi
    fi
    if test $pauseload -ne 0 -a $shouldrun -ne 0
    then
      load=`uptime | \
awk ' {i = 1; while (i < NF && $i != "average:") i++; print $(i+1); } ' | \
sed -e 's/,//' | sed -e 's/\.//'`
      shouldrun=`expr $shouldrun '&' $load '<=' \
'(' $pauseload '+' '(' '100' '*' $isrun ')' ')'`
      if test $shouldrun -ne $isrun
      then why="Load average is $load"
      fi
    fi
    if test $isrun -ne $shouldrun
    then
      if test $shouldrun -ne 0
      then 
	$RESUMECMD $pid
        echo `date`: $RESUMECMD $pid '(mpqs)--'$why >> $control 
      else 
        $PAUSECMD $pid
	echo `date`: $PAUSECMD $pid '(mpqs)--'$why >> $control
      fi
      isrun=$shouldrun
    fi
    sltel=0
    if test -r $skill
    then
	sltel=45
    fi
    if test -r $mkill
    then
	sltel=90
    fi
    if test -r $lkill
    then
	sltel=135
    fi

    if test $sltel -ne 0
    then
	teller=0
	$RESUMECMD $pid >/dev/null 2>&1
	kill $pid >/dev/null 2>&1
	echo `date`: user kill >> $control
	while test $sltel -ne 0
	do
	  . $tick
	  sltel=`expr $sltel - 1 `
          if test '!' -r $mkill -a '!' -r $lkill -a '!' -r $skill
	  then
		sltel=0
	  fi
	done
	rm -f $mkill $lkill $skill
    fi

  done
  $RESUMECMD $pid >/dev/null 2>&1
  kill $pid >/dev/null 2>&1
  counter=15
  while test $counter -ne 0
  do
    . $tick
    counter=`expr $counter - 1`
  done
done
!End!of!file!
chmod 755 runfactor
echo Extracting ticktock
cat >ticktock <<'!End!of!file!'
#! /bin/sh

# included as a subroutine in runfactor; touches the .run file, keeps
# checkrun in the at queue; sleeps for a few minutes
# If you have user-level cron files, you can run checkrun out of your
# cron file once an hour or so, and then you can cut this file down to:
#    touch $run
#    sleep $TICKSEC

touch $run
if test ${now-x} -eq ${now400-x}
then
# if you have trouble with formatted date, use the following instead:
# now=`date | awk '{ print $4; }' | awk -F: '{print $1 $2;}'`
  now=`date +%H%M`
  now110=`expr '(' $now + 110 ')' '%' 2400`
  if test `expr $now110 '%' 100` -ge 60
  then
    now110=`expr '(' $now110 '+' 40 ')' '%' 2400`
  fi
  at $now110 < $docheck >/dev/null 2>&1
  now400=`expr '(' $now + 400 ')' '%' 2400`
  at $now400 < $docheck >/dev/null 2>&1
elif test ${now-x} -eq ${now110-x}
then
  now110=`expr '(' $now + 110 ')' '%' 2400`
  if test `expr $now110 '%' 100` -ge 60
  then
    now110=`expr '(' $now110 '+' 40 ')' '%' 2400`
  fi
  at $now110 < $docheck >/dev/null 2>&1
fi
sleep $TICKSEC
now=`expr $now '+' $TICKLEN`
if test `expr $now '%' 100` -ge 60
then
  now=`expr '(' $now '+' 40 ')' '%' 2400`
fi
!End!of!file!
chmod 755 ticktock
echo Extracting witch
cat >witch <<'!End!of!file!'
#! /bin/sh

# A version of 'which' that doesn't rerun your .cshrc and trash PATH.

opath=$PATH
PATH=/bin:/usr/bin

case $# in
0)	echo 'Usage: witch command' 1>&2; exit 2
esac

for i in `echo $opath | sed 's/^:/.:/
			     s/::/:.:/g
			     s/:$/:./
			     s/:/ /g'`
do
	if test -f $i/$1
	then
		echo $i/$1
		exit 0
	fi
done
echo "No $1 in $opath" 1>&2
exit 1
!End!of!file!
chmod 755 witch
echo Archive complete
exit 0

