#!/bin/sh
# This is a shell archive (produced by shar 3.50)
# To extract the files from this archive, save it to a file, remove
# everything above the "!/bin/sh" line above, and type "sh file_name".
#

# Source directory /usr/cp2/emclean/src
#
# existing files will NOT be overwritten unless -c is specified
#
# This shar contains:
# length  mode       name
# ------ ---------- ------------------------------------------
#   2521 -rw------- readme
#    673 -rw------- makefile
#  11171 -rw------- sd.c
#  47619 -rw------- czero.c
#
# ============= readme ==============
if test -f 'readme' -a X"$1" != X"-c"; then
        echo 'x - skipping readme (File already exists)'
else
echo 'x - extracting readme (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'readme' &&
X
X Here is an alpa version of my lapack-J stuff.
X Spellings and usages are very likely to change.
X I hope to get around to including programs for
X least squares, the generalized eigenvalue
X problem, and QL.
X
X But for now the following are implemented
X (Actually balance needs a touch up )
X
X  Real and Complex version of :
X
/* 30!:0 SVD factorization                           */
/* 30!:1 Matrix norms                                */
/* 30!:2 LU factorization                            */
/* 30!:3 eigvect - eigval                            */
/* 30!:4 balance a matrix                            */
/* 30!:5 H and P in A=PHP' (hessenberg)              */
/* 30!:6 condition numbers                           */
/* 30!:7 QR factorization                            */
/* 30!:8 Cholesky factorization                      */
/* 30!:9 H and Z in A =ZHZ' (schur)                  */
X
X To compile this with J you must
X compile or have libF77.a and libI77.a from
X att.research.com in netlib.
X
X You also must get create the libraries in rblas.tar.Z and
X rlap-0.30.tar.Z from RLaB. RLaB can be found at :
X
X       evans.ee.adfa.oz.au     /pub/RLaB
X
X       mizar.docs.uu.se        /pub/gnu/RLaB
X
X       csi.jpl.nasa.gov        pub/matlab/RLaB
X
X Speaking of RLaB, you might say I hope to provide all of the
X functionality of RLaB - and more. That's my goal. RLaB
X has very good online documentation and I am focusing on an Emacs
X j-interaction-mode.el which I will make J the easy choice
X should anyone compare. With Smillie's introductory
X paper, greate doc, more functionality, the rank operator,
X more applications (CDFLIB is next) J will be the envy of
X RLaB in few months :-)
X
X Oh yes. Lets get back to installation.
X
X Then you include the following line in x.c in the J source :
X
X if(30==p){AF*f1,*f2;
X  ASSERT(jla(q,&f1,&f2),EVDOMAIN);
X  R CDERIV(CIBEAM, f1,f2, RMAXL,RMAXL,RMAXL);}
X
X czero is not actually a lapack hack, so for czero just include
X it in the directory in which you are compiling J. (I want to
X put a TOMS 386 GCD program on as a front end to czero, but
X I hope to get to that later.)
X
X For now I'd perfer not to ask Dr.{*filter*}ey to put this at watserv
X as, I repeat, this is an *ALPHA* version.  I'm simply eager to have
X someone else aside from myself experiment with the code. (I hope
X someone else is interested!
X
X I've included my makefile. If you do not have jk.c the keyed file
X stuff I posted then you will have to delete jk.o from the OBJ line.
X Similarly you may not have ft.o, the file for the fourier transform
X stuff.
X
X Emmett
SHAR_EOF
chmod 0600 readme ||
echo 'restore of readme failed'
Wc_c="`wc -c < 'readme'`"
test 2521 -eq "$Wc_c" ||
        echo 'readme: original size 2521, current size' "$Wc_c"
fi
# ============= makefile ==============
if test -f 'makefile' -a X"$1" != X"-c"; then
        echo 'x - skipping makefile (File already exists)'
else
echo 'x - extracting makefile (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'makefile' &&
# J Makefile for Unix
X
CC     = cc
CFLAGS = -g
HDR    = a.h io.h j.h jc.h je.h js.h jt.h p.h v.h x.h f2c.h
OBJ    = sd.o ft.o jk.o czero.o cb.o j.o a.o ai.o ap.o au.o \
X       c.o cc.o cd.o cf.o cg.o cp.o cr.o ct.o cx.o \
X       f.o i.o io.o  k.o  m.o  p.o pc.o pv.o \
X       r.o rl.o rt.o s.o t.o ta.o u.o ut.o \
X       v.o ve.o vg.o vi.o vm.o vp.o vs.o vx.o vz.o \
X       w.o wn.o x.o xa.o xf.o xl.o  
X
OB2    = vh.o vb.o xs.o
X
j      : $(OBJ) $(OB2)
X       $(CC) -g $(OBJ) $(OB2) \
X    -L./lib/lapack-C -lClapack\
X    -L./lib/blas-C -lCblas\
X    -L./lib/fftpack-C -lCfftpack\
X    -L./lib/libF77  -lF77\
X    -L./lib/libI77  -lI77\
X     -lm -o j
X
$(OBJ) : $(HDR)
X

X       $(CC) -g -c $<
X
SHAR_EOF
chmod 0600 makefile ||
echo 'restore of makefile failed'
Wc_c="`wc -c < 'makefile'`"
test 673 -eq "$Wc_c" ||
        echo 'makefile: original size 673, current size' "$Wc_c"
fi
# ============= sd.c ==============
if test -f 'sd.c' -a X"$1" != X"-c"; then
        echo 'x - skipping sd.c (File already exists)'
else
echo 'x - extracting sd.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'sd.c' &&
#include "j.h"
#include <stdio.h>
#define max(a,b) ((a) >= (b) ? (a) : (b))
#define min(a,b) ((a) <= (b) ? (a) : (b))
#define NUMERIC  (BOOL+INT+FL+CMPX)
X
static A jsvd();          /* (0) SVD factorization                           */
static A jnorm();         /* (1) Matrix norms                                */
static A jlu();           /* (2) LU factorization                            */
static A jeev();          /* (3) eigvect - eigval                            */
static A jbal();          /* (4) balance a matrix                            */
static A jhess();         /* (5) H and P in A=PHP' (hessenberg)              */
static A jrcond();        /* (6) condition numbers                           */
static A jqr();           /* (7) QR factorization                            */
static A jchol();         /* (8) Cholesky factorization                      */
static A jschur();        /* (9) H and Z in A =ZHZ' (schur)                  */
static A jzero();         /* (10) napack complex zero finder                 */
X
void (*select)();
static F1(jschur){PROLOG;
X  A wr,wi,ww,vs,rwork,bwork,work;I n,lwork,sdim=0,info=0;
X  F1RANK(2,jschur,0);RZ(w);
X  if(AN(w)==0)EPILOG(link(mtv,link(mtv,mtv)));
X  if(AN(w)==1)w=reshape(v2(1,1),raze(w));
X  ASSERT((AR(w)==2)&&(AS(w)[0]==AS(w)[1])&&NUMERIC&AT(w),EVDOMAIN);
X  if((BOOL+INT)&AT(w))w=cvt(FL,w);
X  n=AS(w)[0];n=AS(w)[1];w=cant1(w);
X  if(FL&AT(w)){
X   lwork=max(1,3*n);
X   GA(wr,FL,n,1,0);GA(wi,FL,n,1,0);
X   GA(work,FL,lwork,1,0);GA(bwork,BOOL,2,1,0);
X   GA(vs,FL,n*n,2,0);AS(vs)[0]=AS(vs)[1]=n;
X   dgees_("V","N",select,&n,AV(w),&n,&sdim,AV(wr),AV(wi),
X          AV(vs),&n,AV(work),&lwork,AV(bwork),&info,1,1);}
X  else {
X   lwork=max(1,3*n);
X   GA(ww,CMPX,n,1,0);
X   GA(work,CMPX,lwork,1,0);
X   GA(rwork,FL,n,1,0);
X   GA(bwork,BOOL,2,1,0);
X   GA(vs,CMPX,n*n,2,0);AS(vs)[0]=AS(vs)[1]=n;
X   zgees_("V","N",select,&n,AV(w),&n,&sdim,AV(ww),
X          AV(vs),&n,AV(work),&lwork,AV(rwork),AV(bwork),&info,1,1);
X   }
X  EPILOG(link(cant1(w),link(cant1(vs),sc(info))));

X
static F1(jchol){PROLOG;I n,info;
X  F1RANK(2,jchol,0);RZ(w);
X  if(AR(w)<2)w=reshape(v2(1,1),raze(w));
X  ASSERT((AR(w)==2)&&(AS(w)[0]==AS(w)[1])&&NUMERIC&AT(w),EVDOMAIN);
X  if((BOOL+INT)&AT(w))w=cvt(FL,w);n=AS(w)[0];w=cant1(w);
X  if(FL&AT(w))
X   dpotrf_("U",&n,AV(w),&n,&info,1);
X  else
X   zpotrf_("U",&n,AV(w),&n,&info,1);
X EPILOG(link(cant1(tymes(df2(iota(sc(n)),iota(sc(n)),slash(ds(CGE))),w)),sc(info)));
X
static F1(jzero){PROLOG; A work,b,z; I n,lw,nz=0;
X   F1RANK(1,jzero,0);RZ(w);
X   ASSERT(NUMERIC&AT(w),EVDOMAIN);
X   if(AN(w)<=1) EPILOG(mtv);
X   b=df1(ne(zero,w),bsdot(slash(ds(CGT))));
X   b=df1(b,bsdot(slash(ds(COR))));
X   w=repeat(b,w);w=divide(w,tail(w));
X   n=i0(tally(w))-1;
X   lw=max(1000,6*n);
X   GA(work,CMPX,lw,1,0);
X   z=reshape(sc(1000),zero); z=cvt(CMPX,z);
X   w=cvt(CMPX,w);
X   czero1_(AV(z),AV(w),&n,AV(work),&lw,&nz);
X   EPILOG(take(sc(nz),z));
X
static F1(jbal){PROLOG;
X  A a,scale;I m,n,ilo,ihi,info=0;
X  F1RANK(2,jbal,0);RZ(w);
X  if(AN(w)==0)EPILOG(link(mtv,link(mtv,mtv)));
X  if(AN(w)==1)w=reshape(v2(1,1),raze(w));
X  ASSERT((AS(w)[0]==AS(w)[1])&&(AR(w)==2)&&NUMERIC&AT(w),EVDOMAIN);
X  if((BOOL+INT)&AT(w))w=cvt(FL,w);
X  n=AS(w)[0];w=cant1(w);a=scc('S');
X   GA(scale,FL,n,1,0);
X if(FL&AT(w))
X   dgebal_(AV(a), &n, AV(w),&n,&ilo,&ihi, AV(scale),&info,1);
X  else
X   zgebal_(AV(a), &n, AV(w),&n,&ilo,&ihi, AV(scale),&info,1);
X  EPILOG(link(cant1(w),link(table(scale),sc(info))));
X
static F1(jqr){PROLOG;
X  A r,tau,work;I k,m,n,lwork,rinfo=0,qinfo=0;
X  F1RANK(2,jqr,0);RZ(w);
X  if(AN(w)==0)EPILOG(over(link(mtv,mtv),link(mtv,mtv)));
X  if(AN(w)==1)w=reshape(v2(1,1),raze(w));
X  ASSERT(2==AR(w)&&NUMERIC&AT(w),EVDOMAIN);
X  if((BOOL+INT)&AT(w))w=cvt(FL,w);
X  m=AS(w)[0];n=AS(w)[1];k=min(m,n);
X  ASSERT(m>=n,EVLENGTH);w=cant1(w);
X  if(FL&AT(w)){
X   GA(tau,FL,min(n,m),1,0);GA(work,FL,n,1,0);
X   dgeqrf_(&m, &n, AV(w), &m, AV(tau), AV(work), &n, &rinfo);
X   r=reshape(shape(w),w);
X   dorgqr_(&m, &n, &n, AV(w), &m, AV(tau), AV(work), &n, &qinfo);}
X  else {
X   GA(tau,CMPX,min(n,m),1,0);GA(work,CMPX,n,1,0);
X   zgeqrf_(&m, &n, AV(w), &m, AV(tau), AV(work), &n, &rinfo);
X   r=reshape(shape(w),w);
X   zungqr_(&m, &n, &n, AV(w), &m, AV(tau),AV(work), &n, &qinfo);
X   }
X   if(m>n)r=cant1(take(sc(n),cant1(r)));
X   r=tymes(df2(iota(sc(n)),iota(sc(n)),slash(ds(CGE))),r);
X   EPILOG(link(cant1(w),link(cant1(r),link(sc(qinfo),sc(rinfo)))));
X  /* EPILOG(link(cant1(w),cant1(r))); */
X
static F2(jrcond){PROLOG;
X  A wnorm,work,iwork,rwork,ipiv;I info,m,n; D rcond;C norm;
X  F2RANK(0,2,jrcond,0);RZ(a&&w);
X  ASSERT((1==AN(a))&&(CHAR&AT(a)),EVDOMAIN);norm=*(C*)AV(a);
X  ASSERT('1'==norm||'0'==norm||'I'==norm,EVDOMAIN);
X  ASSERT(NUMERIC&AT(w)&&AS(w)[0]==AS(w)[1],EVDOMAIN);
X  if(AN(w)==0)EPILOG(link(mtv,link(mtv,mtv)));
X  if(AN(w)==1)w=reshape(v2(1,1),raze(w));
X  m=AS(w)[0];n=AS(w)[1];
X  if((BOOL+INT)&AT(w))w=cvt(FL,w);
X  if('0'==norm)a=scc('1');wnorm=jnorm(a,w);w=cant1(w);
X  if(FL&AT(w)) {
X   GA(work,FL,4*n,1,0);GA(iwork,INT,INT*n,1,0);GA(ipiv,INT,INT*min(m,n),1,0);
X   dgetrf_(&m,&n,AV(w),&m,AV(ipiv),&info);
X   dgecon_(AV(a),&n,AV(w),&m,AV(wnorm),&rcond,AV(work),
X           AV(iwork),&info,1);}
X  else {
X   GA(work,CMPX,2*n,1,0);GA(rwork,FL,2*n*FL,1,0);GA(ipiv,INT,INT*min(m,n),1,0);
X   zgetrf_(&m, &n, AV(w), &m, ipiv, &info);
X   zgecon_(AV(a),&n,AV(w),&n,AV(wnorm),&rcond,AV(work),
X            AV(rwork),&info,1);
X  }
X  EPILOG(link(scf(rcond),link(wnorm,sc(info))));
X  /* EPILOG(scf(rcond)); */
X
static F1(jhess){PROLOG;
X  A w1,w2,ww,h,tau,work;
X  I n,ilo,ihi,info=0;
X  F1RANK(2,jhess,0);RZ(w);
X  if(AN(w)==0)EPILOG(link(mtv,link(mtv,mtv)));
X  if(AN(w)==1)w=reshape(v2(1,1),raze(w));
X  ASSERT((AR(w)==2)&&(AS(w)[0]==AS(w)[1])&&NUMERIC&AT(w),EVDOMAIN);
X  if((BOOL+INT)&AT(w))w=cvt(FL,w);n=AS(w)[0];
X  ilo=1;ihi=n;w=cant1(w);
X  if(FL&AT(w)){
X   GA(tau,FL,n,1,0);GA(work,FL,FL*n,1,0);
X   dgehrd_(&n,&ilo,&ihi,AV(w),&n,AV(tau),AV(work),&n,&info);
X   h=tymes(df2(increm(iota(sc(n))),iota(sc(n)),slash(ds(CGE))),w);
X   dorghr_(&n,&ilo,&ihi,AV(w),&n,AV(tau),AV(work),&n,&info);}
X  else {
X   GA(tau,CMPX,n,1,0);GA(work,FL,CMPX*n,1,0);
X   zgehrd_(&n,&ilo,&ihi,AV(w),&n,AV(tau),AV(work),&n,&info);
X   h=tymes(df2(increm(iota(sc(n))),iota(sc(n)),slash(ds(CGE))),w);
X   zunghr_(&n,&ilo,&n,AV(w),&n,AV(tau),AV(work),&n,&info);}
X  EPILOG(link(cant1(h),link(cant1(w),sc(info))));
X
static F1(jeev){PROLOG;
X  A o,vl,vr,wr,wi,work,rwork,b,t,z; I i,info,lwork,m,n,*ps;
X  F1RANK(2,jeev,0);RZ(w);
X  if(AN(w)==0)EPILOG(link(mtv,link(mtv,mtv)));
X  if(AN(w)==1)w=reshape(v2(1,1),raze(w));
X  ASSERT(NUMERIC&AT(w)&&(AS(w)[0]==AS(w)[1]),EVDOMAIN);
X  if((BOOL+INT)&AT(w))w=cvt(FL,w);w=cant1(w);n=AS(w)[0];
X  if(FL&AT(w)) {
X   lwork = 4*n;
X   GA(wr,FL,n,1,0);ps=(I*)AS(wr);*ps=n;
X   GA(wi,FL,n,1,0);ps=(I*)AS(wi);*ps=n;
X   GA(vl,FL,1,0,0);
X   GA(vr,FL,n*n,2,0);ps=(I*)AS(vr);*ps++=n;*ps=n;
X   GA(work,FL,lwork,2,0); }
X  else {
X   lwork = 2*n;
X   GA(z,CMPX,n,1,0);ps=(I*)AS(z);*ps=n;
X   GA(vl,FL,1,0,0);
X   GA(vr,CMPX,n*n,2,0);ps=(I*)AS(vr);*ps++=n;*ps=n;
X   GA(work,CMPX,lwork,1,0);
X   GA(rwork,CMPX,2*n,1,0); }
X  if(FL&AT(w)) {
X   dgeev_("N","V",&n,AV(w),&n,AV(wr),AV(wi),
X          AV(vl),&n,AV(vr),&n,AV(work),&lwork,&info,1,1); }
X  else {
X   zgeev_("N","V",&n,AV(w),&n,AV(z),AV(vl),
X          &n,AV(vr),&n,AV(work),&lwork,AV(rwork),&info,1,1); }
X  if(CMPX&AT(w))
X   EPILOG(link(table(z),link(cant1(vr),sc(info))));
X  if(0==i0(eps(one,signum(mag(wi)))))  
X   EPILOG(link(table(wr),link(cant1(vr),sc(info))));
X  b=raze(ne(zero,mag(wi)));z=reshape(v2(n,0),mtv);
X  for(i=0;i<n;i++){
X   if(0==i0(from(sc(i),b)))
X    z=overr(z,from(sc(i),vr));
X   else{
X    t=repeat(repeat(sc(n),sc(2)),cant1(from(v2(i,i+1),vr)));
X    t=tymes(t,reshape(v2(2*n,2),over(v2(1,1),v2(1,-1))));
X    z=overr(z,reshape(v2(n,2),df1(cant1(t),slash(ds(CJDOT)))));
X    i++; }
X  }
X  EPILOG(link(table(plus(wr,tymes(a0j1,wi))),link(z,sc(info))));
X
static F1(jlu){PROLOG; A ind,p,v,b,u;
X  I i,m,n,nn,info=0;F1RANK(2,jlu,0);RZ(w);
X  ASSERT(NUMERIC&AT(w),EVDOMAIN);
X  if((BOOL+INT)&AT(w))w=cvt(FL,w);
X  if(AN(w)==0)EPILOG(over(link(mtv,mtv),link(mtv,mtv)));
X  if(AN(w)==1)w=reshape(v2(1,1),raze(w));
X  ASSERT(AR(w)==2,EVDOMAIN);
X  m=AS(w)[0];n=AS(w)[1];w=cant1(w);
X  p=reshape(sc(min(m,n)),zero);p=cvt(INT,p);
X  if(FL&AT(w))dgetrf_(&m,&n,AV(w),&m,AV(p),&info);
X  else        zgetrf_(&m,&n,AV(w),&m,AV(p),&info);
X  w=cant1(w);
X  u=tymes(w,le(tymes(iota(sc(m)),increm(sc(n))),iota(v2(m,n))));
X  w=tymes(w,gt(tymes(iota(sc(m)),increm(sc(n))),iota(v2(m,n))));
X  ind=take(sc(min(n,m)),tymes(iota(sc(m)),increm(sc(n))));
X  w=df2(reshape(sc(min(n,m)),one),w,rbrace(ind));
X  if(m<n) w=df2(sc(m),w,qq(ds(CTAKE)),sc(-1));else u=take(sc(n),u);
X  b=v=iota(sc(m));nn=i0(tally(p));p=decrem(p);
X  for(i=0;i<nn;i++){
X  ind=over(from(sc(i),p),sc(i));
X  b=from(df2(reverse(ind),v,rbrace(ind)),b);}
X  p=minv(from(b,df2(v,v,slash(ds(CEQ)))));
X  EPILOG(link(p,link(w,link(u),sc(info))));
X
static F2(jnorm){PROLOG; A work,msg=mtv; I l,n,m;
X  D norm, dlange_(),zlange_(); F2RANK(0,2,jnorm,0);RZ(a&&w);
X  ASSERT(CHAR&AT(a)&&NUMERIC&AT(w),EVDOMAIN);
X  ASSERT(i0(eps(a,str(11,"Mm1OoIiFfEe"))),EVDOMAIN);
X  if(AR(w)<2)w=reshape(v2(1,1),raze(w));
X  ASSERT(2==AR(w),EVDOMAIN);
X  if((BOOL+INT)&AT(w))w=cvt(FL,w);
X  if('I'==*(C*)AV(a))
X   GA(work,FL,m,1,0)
X  else
X   GA(work,FL,1,1,0);
X  m=AS(w)[0];n=AS(w)[1];w=cant1(w);
X  if(FL&AT(w)){
X   norm=dlange_(AV(a),&m,&n,AV(w),&m,AV(work),1);}
X  else{
X   norm=zlange_(AV(a),&m,&n,AV(w),&m,AV(work),1);}
X  EPILOG(scf(norm));
static F1(jsvd){ PROLOG; A u,v,s,rwork,work;
X  I m,n,lwork,info=0; F1RANK(2,jsvd,0);
X  ASSERT(NUMERIC&AT(w),EVDOMAIN);
X  if((BOOL+INT)&AT(w))w=cvt(FL,w);
X  if(AN(w)==0)EPILOG(over(link(mtv,mtv),link(mtv,mtv)));
X  if(AN(w)==1)w=reshape(v2(1,1),raze(w));
X  ASSERT(2==AR(w),EVDOMAIN);
X  m=AS(w)[0];n=AS(w)[1];w=cant1(w);
X  if(FL&AT(w)){
X   lwork=max(3*min(m,n)+max(m,n),5*min(m,n)-4);
X   GA(u,FL,(m*m),2,0);AS(u)[0]=AS(u)[1]=m;
X   GA(s,FL,min(m,n),2,0);AS(s)[0]=min(m,n);AS(s)[1]=1;
X   GA(v,FL,(n*n),2,0);AS(v)[0]=AS(v)[1]=n;
X   GA(work,FL,lwork,2,0);
X   dgesvd_("A","A",&m,&n,AV(w),&m,AV(s),AV(u),
X           &m,AV(v),&n,AV(work),&lwork,&info,1,1);}
X  else {
X   lwork=2*min(m,n)+max(m,n);
X   GA(u, CMPX,(m*m),2,0);AS(u)[0]=AS(u)[1]=m;
X   GA(s,FL,min(m,n),2,0);AS(s)[0]=min(m,n);AS(s)[1]=1;
X   GA(v,CMPX,(n*n),2,0);AS(v)[0]=AS(v)[1]=n;
X   GA(rwork,CMPX,5*max(m,n),2,0);
X   GA( work,CMPX,lwork,2,0);
X   zgesvd_("A","A",&m,&n,AV(w),&m,AV(s),AV(u),
X           &m,AV(v),&n,AV(work),&lwork,AV(rwork),&info,1,1);}
X   EPILOG(link(cant1(u),link(s,link(v,sc(info)))));
X
C jla (k, f1, f2) I k;
AF * f1, *f2;
{
X  switch (k) {
X    case 0: *f1 = jsvd    ; *f2 = NULL     ; R 1;                      
X    case 1: *f1 = NULL    ; *f2 = jnorm    ; R 1;                      
X    case 2: *f1 = jlu     ; *f2 = NULL     ; R 1;      
X    case 3: *f1 = jeev    ; *f2 = NULL     ; R 1;      
X    case 4: *f1 = jbal    ; *f2 = NULL     ; R 1;      
X    case 5: *f1 = jhess   ; *f2 = NULL     ; R 1;      
X    case 6: *f1 = NULL    ; *f2 = jrcond   ; R 1;      
X    case 7: *f1 = jqr     ; *f2 = NULL     ; R 1;      
X    case 8: *f1 = jchol   ; *f2 = NULL     ; R 1;      
X    case 9: *f1 = jschur  ; *f2 = NULL     ; R 1;
X    case 10: *f1 = jzero  ; *f2 = NULL     ; R 1;
X    default:
X      ASSERT (0, EVNONCE);    
X  }
SHAR_EOF
chmod 0600 sd.c ||
echo 'restore of sd.c failed'
Wc_c="`wc -c < 'sd.c'`"
test 11171 -eq "$Wc_c" ||
        echo 'sd.c: original size 11171, current size' "$Wc_c"
fi
# ============= czero.c ==============
if test -f 'czero.c' -a X"$1" != X"-c"; then
        echo 'x - skipping czero.c (File already exists)'
else
echo 'x - extracting czero.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'czero.c' &&
/* file.f -- translated by f2c (version of 23 April 1993  18:34:30).
X   You must link the resulting object file with the libraries:
X       -lF77 -lI77 -lm -lc   (in that order)
*/
X
#include "f2c.h"
X
/* Table of constant values */
X
static doublecomplex c_b21 = {1.,0.};
static integer c__2 = 2;
static integer c__9 = 9;
static integer c__1 = 1;
static integer c__1000 = 1000;
static integer c__7 = 7;
static real c_b142 = (float)8.;
static real c_b143 = (float)7.;
static real c_b144 = (float)59.;
static real c_b145 = (float)56.;
static real c_b146 = (float)834.;
static real c_b147 = (float)-19.;
static real c_b148 = (float)-26161.;
static real c_b149 = (float)2382.;
static real c_b150 = (float)-159372.;
static real c_b151 = (float)-161959.;
static real c_b152 = (float)-949059.;
static real c_b153 = (float)-1102644.;
static real c_b154 = (float)-12620270.;
static real c_b155 = (float)1210779.;
static real c_b156 = (float)152570664.;
static real c_b157 = (float)-30182594.;
static real c_b158 = (float)590670352.;
static real c_b159 = (float)787660776.;
static real c_b160 = (float)3031287952.;
static real c_b161 = (float)4148487952.;
static real c_b162 = (float)32217708704.;
static real c_b163 = (float)-8128629584.;
static real c_b164 = (float)-188996316672.;
static real c_b165 = (float)43166023136.;
static real c_b166 = (float)-321964024320.;
static real c_b167 = (float)-560530933248.;
static real c_b168 = (float)-1179344517120.;
static real c_b169 = (float)-1674091261440.;
static real c_b170 = (float)-5319936921600.;
static real c_b171 = (float)2194678886400.;
static integer c__3 = 3;
X
/* czero1.m4 - last revised K M Briggs 1993 April 9 */
/* To regenerate fortran versions: */
/* m4 -DX file.m4 >file.f */
/* where X=S, or Q for single, double, or quadruple. */
/* Jenkins-Traub algorithm to compute zeros of a complex polynomial. */
/* A version of czero.f from NAPACK by Bill Hager from netlib, */
/* Adapted by K M Briggs. */
/* but note changed calling sequence. */
/* For use with Bailey's transmp. */
X
/* Changes by KMB: */
/* 1. Check w is large enough */
/* 2. Generic names for intrinsics */
/* 3. CMP directives */
/* 4. cz1test subroutine for self-check */
/* 5. Fix bug in linear case */
/* Set desired MP decimal places here... */
/* MP+ PRECISION LEVEL 100 */
/* Need TYPE ERROR OFF as work array w used as real and complex... */
/* MP+ TYPE ERROR OFF */
/* Test program... */
/*--- DO NOT CHANGE BELOW HERE ---------------------------------------------*/
/* Subroutine */ int czero1_(z, p, nd, w, lw, nz)
doublecomplex *z, *p;
integer *nd;
doublecomplex *w;
integer *lw, *nz;
{
X    /* System generated locals */
X    integer i__1, i__2, i__3;
X    doublereal d__1;
X    doublecomplex z__1, z__2, z__3;
X
X    /* Builtin functions */
X    double z_abs();
X    void z_div();
X    double log();
X    integer s_wsfe(), e_wsfe();
X
X    /* Local variables */
X    extern /* Subroutine */ int dmag_();
X    static doublereal a, b, c, d, e;
X    static integer f, g, h, i, j, k, l, m, n, o, q;
X    static doublecomplex r, s;
X    static doublereal t;
X    static doublecomplex u, v, x;
X    static doublereal y, theta;
X    static integer n1, o1, o2, o3, o4;
X    static doublereal s1, t1, t2, v1, y0;
X    extern /* Subroutine */ int cc_(), ff_();
X    static integer nc;
X    static doublereal dm, ct;
X    extern /* Subroutine */ int jt_();
X    static doublereal st;
X    extern doublereal dpreal_();
X    extern /* Double Complex */ int dpcmpl_();
X    extern /* Subroutine */ int dpcssn_(), cef_(), qad_(), rad_(), rfn_(),
X           err_();
X
X    /* Fortran I/O blocks */
X    static cilist io___39 = { 0, 6, 0, "(\"WARNING: czero1 fails to separate\
X roots.\")", 0 };
X
X
/*      ________________________________________________________ */
/*     |                                                        | */
/*     | use Jenkins-Traub algorithm to compute zeros of a monic| */
/*     |  polynomial p(1) + p(2)x +...+ p(nd)x**(nd-1) + x**nd  | */
/*     |                                                        | */
/*     |    input:                                              | */
/*     |                                                        | */
/*     |         p     --the first nd coefficients of polynomial| */
/*     |                 (complex real*8 array)       | */
/*     |                                                        | */
/*     |         nd    --degree of polynomial                   | */
/*     |                                                        | */
/*     |         w     --work array (required length depends on | */
/*     |                 the difficulty in computing zeros. in  | */
/*     |                 many cases, between 4nd+32 and 6nd     | */
/*     |                 complex real*8 entries are   | */
/*     |                 sufficient)                            | */
/*     |         lw   --the CDP size of w                       | */
/*     |                                                        | */
/*     |                                                        | */
/*     |    output:                                             | */
/*     |                                                        | */
/*     |         z     --zeros (complex real*8 array) | */
/*     |                                                        | */
/*     |         nz    --number of computed zeros               | */
/*     |                (-1 indicates lw too small)             | */
/*     |    note:                                               | */
/*     |                                                        | */
/*     |      this algorithm can generate underflows since the  | */
/*     |      function values tend to zero as the iterations    | */
/*     |      approach a root. thus the appropriate compiler    | */
/*     |      options must be invoked so that an underflow is   | */
/*     |      replaced by zero and the execution continues. for | */
/*     |      ibm fortran, try the following statement:         | */
/*     |               call errset(208,256,-1,-1,0)             | */
/*     |                                                        | */
/*     |    subroutines: cc,cef,dmag,err,ff,jt,qad,rad,rfn      | */
/*     |________________________________________________________| */
X
/* Original by: */
/* Bill Hager */
/* Mathematics Department */
/* University of Florida */
/* Gainesville, FL 32611 */
X

X
/* MP+ MULTIP COMPLEX p,w,z,r,s,u,v,x */
/* MP+ MULTIP REAL a,b,c,d,e,s1,t,t1,t2,v1,y,y0,theta,ct,st,dm */
/* DP pi... */
/* DP machine epsilon... */
X
X    /* Parameter adjustments */
X    --w;
X    --p;
X    --z;
X
X    /* Function Body */
X    n = *nd;
X    *nz = n;
/*     is poly constant ? */
X    if (n == 0) {
X       return 0;
X    }
X    i__1 = n;
X    for (i = 1; i <= i__1; ++i) {
X       i__2 = i;
X       i__3 = i;
X       z[i__2].r = p[i__3].r, z[i__2].i = p[i__3].i;
/* L10: */
X    }
L20:
X    if (n == 1) {
/*       poly is linear */
X       z__1.r = -p[1].r, z__1.i = -p[1].i;
X       z[1].r = z__1.r, z[1].i = z__1.i;
X       return 0;
X    }
X    if (n > 2) {
X       goto L30;
X    }
X    qad_(&z[2], &z[1]);
X    return 0;
/*     -------------------------------- */
/*     |    remove zeros at origin    | */
/*     -------------------------------- */
L30:
X    dmag_(&z[1], &dm);
X    if (dm != 0.) {
X       goto L50;
X    }
X    i__1 = n;
X    for (i = 2; i <= i__1; ++i) {
X       i__2 = i - 1;
X       i__3 = i;
X       z[i__2].r = z[i__3].r, z[i__2].i = z[i__3].i;
/* L40: */
X    }
X    i__1 = n;
X    z[i__1].r = 0., z[i__1].i = 0.;
X    --n;
X    goto L20;
L50:
X    *nz = n;
/*     --------------------------------- */
/*     |    compute machine epsilon    | */
/*     --------------------------------- */
/* Changed to fixed 2.22d-16 by KMB */
/* Next line to allow MP conversion... */
X    e = 2.22e-16;
X    m = 5;
X    x.r = 1., x.i = 0.;
X    y = 65536.;
X    y0 = 1. / y;
L70:
X    n1 = n + 1;
X    o1 = n1 + n;
X    o2 = o1 + n;
X    o3 = o2 + n1 / 2;
X    o4 = o3 + n1 / 2;
X    o = o4 - 1;
X    if (o1 + n > *lw || o2 + n > *lw) {
/*       w too small for cc */
X       *nz = -1;
X       return 0;
X    }
X    cc_(&z[1], &w[o1], &w[o2], &n, &y, &y0);
X    g = o4 + n1;
X    nc = n << 5;
X    k = 10;
X    l = 8;
X    q = 0;
X    j = 0;
X    f = 0;
L80:
X    ++f;
X    q += j;
X    if (f < 6) {
X       goto L90;
X    }
X    if (l > nc) {
X       goto L370;
X    }
L90:
X    l += l;
/*     -------------------------------------------------------------- */
/*     |    estimate radius of smallest circle containing a zero    | */
/*     -------------------------------------------------------------- */
X    if (n1 + n > *lw) {
/* ! w too small for rad */
X       *nz = -1;
X       return 0;
X    }
X    rad_(&z[1], &w[1], &w[n1], &k, &n, &a);
X    if (o1 + n > *lw || o2 + n > *lw || o3 + n > *lw || o4 + n > *lw) {
/*     ! w too small for cef */
X       *nz = -1;
X       return 0;
X    }
X    cef_(&w[o4], &z[1], &w[o1], &w[o2], &w[o3], &a, &x, &n, &n1, &t, &t1, &y,
X           &y0);
X    t = n * t * e;
X    j = 1;
X    i = n1;
X    h = i;
L100:
X    i /= 2;
X    j += j;
X    if (i > 0) {
X       goto L100;
X    }
X    if (h + h == j) {
X       goto L120;
X    }
X    h = j;
X    j = h + o;
X    if (j > *lw) {
/*      ! w too small */
X       *nz = -1;
X       return 0;
X    }
X    i__1 = j;
X    for (i = g; i <= i__1; ++i) {
X       i__2 = i;
X       w[i__2].r = 0., w[i__2].i = 0.;
/* L110: */
X    }
L120:
X    j = l + o4;
X    if (o4 + (l << 1) > *lw || j + l > *lw) {
/*      ! w too small for ff */
X       *nz = -1;
X       return 0;
X    }
X    ff_(&w[o4], &w[j], &h, &l);
X    dmag_(&w[o4], &c);
X    j = 0;
X    h = l - 1;
/*     ---------------------------------- */
/*     |    determine starting guess    | */
/*     ---------------------------------- */
X    i__1 = h;
X    for (i = 1; i <= i__1; ++i) {
X       dmag_(&w[o4 + i], &d);
X       if (d > c) {
X           goto L130;
X       }
X       j = i;
X       c = d;
L130:
X       ;
X    }
X    theta = j * 2. * 3.14159265358979323846 / l;
X    dpcssn_(&theta, &ct, &st);
X    z__2.r = a * x.r, z__2.i = a * x.i;
X    dpcmpl_(&z__3, &ct, &st);
X    z__1.r = z__2.r * z__3.r - z__2.i * z__3.i, z__1.i = z__2.r * z__3.i +
X           z__2.i * z__3.r;
X    s.r = z__1.r, s.i = z__1.i;
X    s1 = z_abs(&s);
X    if (s1 < 1.) {
X       s1 = 1.;
X    }
X    theta = 3.14159265358979323846 / l / 2.;
X    dpcssn_(&theta, &ct, &st);
X    dpcmpl_(&z__2, &ct, &st);
X    z__1.r = x.r * z__2.r - x.i * z__2.i, z__1.i = x.r * z__2.i + x.i *
X           z__2.r;
X    x.r = z__1.r, x.i = z__1.i;
X    d__1 = z_abs(&x);
X    z__1.r = x.r / d__1, z__1.i = x.i / d__1;
X    x.r = z__1.r, x.i = z__1.i;
/*     ---------------------------------- */
/*     |    jenkins-traub iterations    | */
/*     ---------------------------------- */
X    i__1 = m;
X    for (i = 1; i <= i__1; ++i) {
X       jt_(&z[1], &w[1], &n, &v, &s, &s1);
/* L140: */
X    }
X    c = 0.;
X    j = 0;
L150:
X    jt_(&z[1], &w[1], &n, &v, &s, &s1);
X    dmag_(&v, &v1);
X    if (v1 == 0.) {
X       goto L240;
X    }
X    if (s1 > 1.) {
X       goto L170;
X    }
X    u.r = 1., u.i = 0.;
X    i = n;
L160:
X    --i;
X    i__1 = i;
X    z__2.r = u.r * s.r - u.i * s.i, z__2.i = u.r * s.i + u.i * s.r;
X    z__1.r = w[i__1].r + z__2.r, z__1.i = w[i__1].i + z__2.i;
X    u.r = z__1.r, u.i = z__1.i;
X    if (i > 1) {
X       goto L160;
X    }
X    dmag_(&u, &dm);
X    if (dm == 0.) {
X       goto L80;
X    }
X    goto L190;
L170:
X    u.r = w[1].r, u.i = w[1].i;
X    z_div(&z__1, &c_b21, &s);
X    r.r = z__1.r, r.i = z__1.i;
X    i__1 = n;
X    for (i = 2; i <= i__1; ++i) {
X       i__2 = i;
X       z__2.r = u.r * r.r - u.i * r.i, z__2.i = u.r * r.i + u.i * r.r;
X       z__1.r = w[i__2].r + z__2.r, z__1.i = w[i__2].i + z__2.i;
X       u.r = z__1.r, u.i = z__1.i;
/* L180: */
X    }
X    dmag_(&u, &dm);
X    if (dm == 0.) {
X       goto L80;
X    }
X    z__1.r = u.r * r.r - u.i * r.i, z__1.i = u.r * r.i + u.i * r.r;
X    u.r = z__1.r, u.i = z__1.i;
L190:
X    z_div(&z__2, &v, &u);
X    z__1.r = s.r - z__2.r, z__1.i = s.i - z__2.i;
X    r.r = z__1.r, r.i = z__1.i;
X    ++j;
X    z__1.r = r.r - s.r, z__1.i = r.i - s.i;
X    dmag_(&z__1, &d);
X    s.r = r.r, s.i = r.i;
X    if (s1 <= a) {
X       goto L200;
X    }
X    err_(&z[1], &n, &s1, &t2);
X    t2 = n * t2 * e;
X    if (v1 <= t2) {
X       goto L240;
X    }
X    goto L210;
L200:
X    if (n * log(t1 / s1) + log(t / v1) >= 0.) {
X       goto L230;
X    }
L210:
X    s1 = z_abs(&s);
X    if (s1 < 1.) {
X       s1 = 1.;
X    }
X    if (j < 3) {
X       goto L220;
X    }
X    if (d * .5 > b + c) {
X       goto L80;
X    }
X    if (j < m) {
X       goto L220;
X    }
X    if (d * 4 > b + c) {
X       goto L80;
X    }
L220:
X    b = c;
X    c = d;
X    goto L150;
L230:
X    err_(&z[1], &n, &s1, &t2);
X    t2 = n * e * t2;
X    if (v1 > t2) {
X       goto L210;
X    }
L240:
X    i__1 = n;
X    z__1.r = z[i__1].r + s.r, z__1.i = z[i__1].i + s.i;
X    r.r = z__1.r, r.i = z__1.i;
X    i = n;
L250:
X    --i;
X    i__1 = i;
X    z__2.r = r.r * s.r - r.i * s.i, z__2.i = r.r * s.i + r.i * s.r;
X    z__1.r = z[i__1].r + z__2.r, z__1.i = z[i__1].i + z__2.i;
X    v.r = z__1.r, v.i = z__1.i;
X    i__1 = i;
X    z[i__1].r = r.r, z[i__1].i = r.i;
X    r.r = v.r, r.i = v.i;
X    if (i > 1) {
X       goto L250;
X    }
X    i__1 = n;
X    z[i__1].r = s.r, z[i__1].i = s.i;
X    q += j;
X    --n;
X    if (n > 2) {
X       goto L70;
X    }
X    qad_(&z[2], &z[1]);
X    nc = *nd;
X    n = 1;
L260:
X    l = n;
X    n = *nd;
X    m = n + n;
X    j = 1;
X    k = n - 1;
X    if (max(n,m) + k > *lw) {
/*      ! w too small */
X       *nz = -1;
X       return 0;
X    }
X    i__1 = k;
X    for (i = 2; i <= i__1; ++i) {
X       i__2 = j;
X       i__3 = i;
X       d__1 = (doublereal) j;
X       z__1.r = d__1 * p[i__3].r, z__1.i = d__1 * p[i__3].i;
X       w[i__2].r = z__1.r, w[i__2].i = z__1.i;
X       i__2 = j + n;
X       d__1 = (doublereal) (j * i);
X       i__3 = i + 1;
X       z__1.r = d__1 * p[i__3].r, z__1.i = d__1 * p[i__3].i;
X       w[i__2].r = z__1.r, w[i__2].i = z__1.i;
X       dmag_(&p[j], &w[j + m]);
X       j = i;
/* L270: */
X    }
X    if (m + n > *lw) {
/*      ! w too small */
X       *nz = -1;
X       return 0;
X    }
X    i__1 = k;
X    i__2 = n;
X    d__1 = (doublereal) k;
X    z__1.r = d__1 * p[i__2].r, z__1.i = d__1 * p[i__2].i;
X    w[i__1].r = z__1.r, w[i__1].i = z__1.i;
X    i__1 = n;
X    w[i__1].r = (doublereal) n, w[i__1].i = 0.;
X    i__1 = m - 1;
X    d__1 = (doublereal) (n * k);
X    w[i__1].r = d__1, w[i__1].i = 0.;
X    i__1 = m;
X    w[i__1].r = 0., w[i__1].i = 0.;
X    dmag_(&p[k], &w[m + k]);
X    dmag_(&p[n], &w[m + n]);
X    o1 = m + 1;
X    n1 = n + 1;
/*     -------------------------- */
/*     |    refine the zeros    | */
/*     -------------------------- */
X    i__1 = *nz;
X    for (k = l; k <= i__1; ++k) {
X       i__2 = k;
X       s.r = z[i__2].r, s.i = z[i__2].i;
X       for (j = 1; j <= 3; ++j) {
X           rfn_(&p[1], &w[1], &w[n1], &w[o1], &n, &s, &v, &c, &t);
X           t = n * t * e;
X           if (c <= t) {
X               goto L330;
X           }
/* L280: */
X       }
X       a = z_abs(&s);
X       if (a > 1.) {
X           goto L300;
X       }
X       v.r = 1., v.i = 0.;
X       t = 1.;
X       i = n;
L290:
X       z__2.r = v.r * s.r - v.i * s.i, z__2.i = v.r * s.i + v.i * s.r;
X       i__2 = i;
X       z__1.r = z__2.r + p[i__2].r, z__1.i = z__2.i + p[i__2].i;
X       v.r = z__1.r, v.i = z__1.i;
X       t = t * a + dpreal_(&w[i + m]);
X       --i;
X       if (i > 0) {
X           goto L290;
X       }
X       goto L320;
L300:
X       z_div(&z__1, &c_b21, &s);
X       r.r = z__1.r, r.i = z__1.i;
X       a = 1. / a;
X       v.r = 0., v.i = 0.;
X       t = 0.;
X       i__2 = n;
X       for (i = 1; i <= i__2; ++i) {
X           z__2.r = v.r * r.r - v.i * r.i, z__2.i = v.r * r.i + v.i * r.r;
X           i__3 = i;
X           z__1.r = z__2.r + p[i__3].r, z__1.i = z__2.i + p[i__3].i;
X           v.r = z__1.r, v.i = z__1.i;
X           t = t * a + dpreal_(&w[i + m]);
/* L310: */
X       }
X       z__2.r = v.r * r.r - v.i * r.i, z__2.i = v.r * r.i + v.i * r.r;
X       z__1.r = z__2.r + 1., z__1.i = z__2.i;
X       v.r = z__1.r, v.i = z__1.i;
X       t = t * a + 1.;
L320:
X       dmag_(&v, &a);
X       if (a > n * t * e) {
X           goto L340;
X       }
L330:
X       i__2 = k;
X       z[i__2].r = s.r, z[i__2].i = s.i;
L340:
X       ;
X    }
X    i__1 = n;
X    for (i = l; i <= i__1; ++i) {
X       i__2 = n1 - i;
X       i__3 = i;
X       w[i__2].r = z[i__3].r, w[i__2].i = z[i__3].i;
/* L350: */
X    }
X    l = n1 - l;
X    i__1 = l;
X    for (i = 1; i <= i__1; ++i) {
X       i__2 = i;
X       i__3 = i;
X       z[i__2].r = w[i__3].r, z[i__2].i = w[i__3].i;
/* L360: */
X    }
X    *nz = nc;
X    return 0;
L370:
X    nc = *nd - n + 1;
X    s_wsfe(&io___39);
X    e_wsfe();
/*     write(6,*) 'after',q,'iterations, we stop while computing' */
/*     write(6,*) 'root number',nc */
X    i__1 = n;
X    z[i__1].r = s.r, z[i__1].i = s.i;
X    goto L260;
X
/* % */
/*     ------------------------------ */
/*     |  roots of a quadratic      | */
/*     ------------------------------ */
/* Subroutine */ int qad_(b, c)
doublecomplex *b, *c;
{
X    /* System generated locals */
X    doublereal d__1;
X    doublecomplex z__1, z__2, z__3, z__4, z__5, z__6;
X
X    /* Builtin functions */
X    double sqrt();
X    void z_sqrt(), pow_zi(), z_div();
X
X    /* Local variables */
X    extern /* Subroutine */ int dmag_();
X    static doublecomplex r, s;
X    static doublereal t, dm;
X
/* MP+ MULTIP COMPLEX b,c,r,s */
/* MP+ MULTIP REAL t,dm */
X
X    dmag_(b, &t);
X    d__1 = 4.;
X    z__1.r = d__1 * c->r, z__1.i = d__1 * c->i;
X    r.r = z__1.r, r.i = z__1.i;
X    dmag_(&r, &dm);
X    if (t > sqrt(dm)) {
X       goto L10;
X    }
X    z__2.r = b->r * b->r - b->i * b->i, z__2.i = b->r * b->i + b->i * b->r;
X    z__1.r = z__2.r - r.r, z__1.i = z__2.i - r.i;
X    s.r = z__1.r, s.i = z__1.i;
X    z_sqrt(&z__1, &s);
X    r.r = z__1.r, r.i = z__1.i;
X    goto L20;
L10:
X    d__1 = 1. / t;
X    s.r = d__1, s.i = 0.;
X    z__4.r = s.r * b->r - s.i * b->i, z__4.i = s.r * b->i + s.i * b->r;
X    pow_zi(&z__3, &z__4, &c__2);
X    z__6.r = r.r * s.r - r.i * s.i, z__6.i = r.r * s.i + r.i * s.r;
X    z__5.r = z__6.r * s.r - z__6.i * s.i, z__5.i = z__6.r * s.i + z__6.i *
X           s.r;
X    z__2.r = z__3.r - z__5.r, z__2.i = z__3.i - z__5.i;
X    z_sqrt(&z__1, &z__2);
X    r.r = z__1.r, r.i = z__1.i;
X    z__1.r = t * r.r, z__1.i = t * r.i;
X    r.r = z__1.r, r.i = z__1.i;
L20:
X    z__1.r = r.r + b->r, z__1.i = r.i + b->i;
X    s.r = z__1.r, s.i = z__1.i;
X    z__1.r = r.r - b->r, z__1.i = r.i - b->i;
X    r.r = z__1.r, r.i = z__1.i;
X    dmag_(&r, &t);
X    dmag_(&s, &dm);
X    if (t >= dm) {
X       goto L30;
X    }
X    z__3.r = c->r + c->r, z__3.i = c->i + c->i;
X    z__2.r = -z__3.r, z__2.i = -z__3.i;
X    z_div(&z__1, &z__2, &s);
X    c->r = z__1.r, c->i = z__1.i;
X    z__1.r = s.r * -.5, z__1.i = s.i * -.5;
X    b->r = z__1.r, b->i = z__1.i;
X    return 0;
L30:
X    if (t == 0.) {
X       return 0;
X    }
X    z__2.r = c->r + c->r, z__2.i = c->i + c->i;
X    z_div(&z__1, &z__2, &r);
X    c->r = z__1.r, c->i = z__1.i;
X    z__1.r = r.r * .5, z__1.i = r.i * .5;
X    b->r = z__1.r, b->i = z__1.i;
X    return 0;
X
/* % */
/*     --------------------------------- */
/*     |    jenkins-traub iteration    | */
/*     --------------------------------- */
/* Subroutine */ int jt_(p, q, n, v, s, s1)
doublecomplex *p, *q;
integer *n;
doublecomplex *v, *s;
doublereal *s1;
{
X    /* System generated locals */
X    integer i__1, i__2, i__3;
X    doublecomplex z__1, z__2, z__3, z__4;
X
X    /* Builtin functions */
X    void z_div();
X
X    /* Local variables */
X    extern /* Subroutine */ int dmag_();
X    static integer i, m;
X    static doublecomplex r, t, u, w;
X    static doublereal dm;
X
/* MP+ MULTIP COMPLEX p,q,r,s,t,u,v,w */
/* MP+ MULTIP REAL dm,s1 */
X    /* Parameter adjustments */
X    --q;
X    --p;
X
X    /* Function Body */
X    if (*s1 > 1.) {
X       goto L70;
X    }
X    u.r = 0., u.i = 0.;
X    v->r = 1., v->i = 0.;
X    i = *n;
L10:
X    i__1 = i;
X    z__2.r = v->r * s->r - v->i * s->i, z__2.i = v->r * s->i + v->i * s->r;
X    z__1.r = p[i__1].r + z__2.r, z__1.i = p[i__1].i + z__2.i;
X    v->r = z__1.r, v->i = z__1.i;
X    i__1 = i;
X    z__2.r = u.r * s->r - u.i * s->i, z__2.i = u.r * s->i + u.i * s->r;
X    z__1.r = q[i__1].r + z__2.r, z__1.i = q[i__1].i + z__2.i;
X    u.r = z__1.r, u.i = z__1.i;
X    --i;
X    if (i > 0) {
X       goto L10;
X    }
X    dmag_(&u, &dm);
X    if (dm != 0.) {
X       goto L50;
X    }
X    m = *n;
L20:
X    i__1 = m;
X    r.r = q[i__1].r, r.i = q[i__1].i;
X    i__1 = m;
X    q[i__1].r = 0., q[i__1].i = 0.;
X    i = m;
L30:
X    --i;
X    i__1 = i;
X    z__2.r = r.r * s->r - r.i * s->i, z__2.i = r.r * s->i + r.i * s->r;
X    z__1.r = q[i__1].r + z__2.r, z__1.i = q[i__1].i + z__2.i;
X    u.r = z__1.r, u.i = z__1.i;
X    i__1 = i;
X    q[i__1].r = r.r, q[i__1].i = r.i;
X    r.r = u.r, r.i = u.i;
X    if (i > 1) {
X       goto L30;
X    }
X    u.r = 0., u.i = 0.;
X    i = m;
L40:
X    i__1 = i;
X    z__2.r = u.r * s->r - u.i * s->i, z__2.i = u.r * s->i + u.i * s->r;
X    z__1.r = q[i__1].r + z__2.r, z__1.i = q[i__1].i + z__2.i;
X    u.r = z__1.r, u.i = z__1.i;
X    --i;
X    if (i > 0) {
X       goto L40;
X    }
X    --m;
X    dmag_(&u, &dm);
X    if (dm == 0.) {
X       goto L20;
X    }
L50:
X    z_div(&z__1, v, &u);
X    r.r = z__1.r, r.i = z__1.i;
X    i__1 = *n;
X    i__2 = *n;
X    z__3.r = r.r * q[i__2].r - r.i * q[i__2].i, z__3.i = r.r * q[i__2].i +
X           r.i * q[i__2].r;
X    z__2.r = p[i__1].r - z__3.r, z__2.i = p[i__1].i - z__3.i;
X    z__1.r = z__2.r + s->r, z__1.i = z__2.i + s->i;
X    t.r = z__1.r, t.i = z__1.i;
X    i__1 = *n;
X    q[i__1].r = 1., q[i__1].i = 0.;
X    i = *n;
L60:
X    --i;
X    i__1 = i;
X    i__2 = i;
X    z__3.r = r.r * q[i__2].r - r.i * q[i__2].i, z__3.i = r.r * q[i__2].i +
X           r.i * q[i__2].r;
X    z__2.r = p[i__1].r - z__3.r, z__2.i = p[i__1].i - z__3.i;
X    z__4.r = t.r * s->r - t.i * s->i, z__4.i = t.r * s->i + t.i * s->r;
X    z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
X    u.r = z__1.r, u.i = z__1.i;
X    i__1 = i;
X    q[i__1].r = t.r, q[i__1].i = t.i;
X    t.r = u.r, t.i = u.i;
X    if (i > 1) {
X       goto L60;
X    }
X    return 0;
L70:
X    z_div(&z__1, &c_b21, s);
X    w.r = z__1.r, w.i = z__1.i;
X    u.r = q[1].r, u.i = q[1].i;
X    v->r = p[1].r, v->i = p[1].i;
X    i__1 = *n;
X    for (i = 2; i <= i__1; ++i) {
X       z__2.r = u.r * w.r - u.i * w.i, z__2.i = u.r * w.i + u.i * w.r;
X       i__2 = i;
X       z__1.r = z__2.r + q[i__2].r, z__1.i = z__2.i + q[i__2].i;
X       u.r = z__1.r, u.i = z__1.i;
X       z__2.r = v->r * w.r - v->i * w.i, z__2.i = v->r * w.i + v->i * w.r;
X       i__2 = i;
X       z__1.r = z__2.r + p[i__2].r, z__1.i = z__2.i + p[i__2].i;
X       v->r = z__1.r, v->i = z__1.i;
/* L80: */
X    }
X    z__1.r = u.r * w.r - u.i * w.i, z__1.i = u.r * w.i + u.i * w.r;
X    u.r = z__1.r, u.i = z__1.i;
X    z__2.r = v->r * w.r - v->i * w.i, z__2.i = v->r * w.i + v->i * w.r;
X    z__1.r = z__2.r + 1., z__1.i = z__2.i;
X    v->r = z__1.r, v->i = z__1.i;
X    dmag_(&u, &dm);
X    if (dm != 0.) {
X       goto L120;
X    }
X    m = *n;
L90:
X    i__1 = m;
X    r.r = q[i__1].r, r.i = q[i__1].i;
X    i__1 = m;
X    q[i__1].r = 0., q[i__1].i = 0.;
X    i = m;
L100:
X    --i;
X    i__1 = i;
X    z__2.r = r.r * s->r - r.i * s->i, z__2.i = r.r * s->i + r.i * s->r;
X    z__1.r = q[i__1].r + z__2.r, z__1.i = q[i__1].i + z__2.i;
X    u.r = z__1.r, u.i = z__1.i;
X    i__1 = i;
X    q[i__1].r = r.r, q[i__1].i = r.i;
X    r.r = u.r, r.i = u.i;
X    if (i > 1) {
X       goto L100;
X    }
X    u.r = 0., u.i = 0.;
X    --m;
X    i__1 = *n;
X    for (i = 1; i <= i__1; ++i) {
X       z__2.r = u.r * w.r - u.i * w.i, z__2.i = u.r * w.i + u.i * w.r;
X       i__2 = i;
X       z__1.r = z__2.r + q[i__2].r, z__1.i = z__2.i + q[i__2].i;
X       u.r = z__1.r, u.i = z__1.i;
/* L110: */
X    }
X    z__1.r = u.r * w.r - u.i * w.i, z__1.i = u.r * w.i + u.i * w.r;
X    u.r = z__1.r, u.i = z__1.i;
X    dmag_(&u, &dm);
X    if (dm == 0.) {
X       goto L90;
X    }
L120:
X    z_div(&z__1, v, &u);
X    r.r = z__1.r, r.i = z__1.i;
X    z__2.r = r.r * q[1].r - r.i * q[1].i, z__2.i = r.r * q[1].i + r.i * q[1]
X           .r;
X    z__1.r = p[1].r - z__2.r, z__1.i = p[1].i - z__2.i;
X    u.r = z__1.r, u.i = z__1.i;
X    z__2.r = -w.r, z__2.i = -w.i;
X    z__1.r = z__2.r * u.r - z__2.i * u.i, z__1.i = z__2.r * u.i + z__2.i *
X           u.r;
X    q[1].r = z__1.r, q[1].i = z__1.i;
X    i__1 = *n;
X    for (i = 2; i <= i__1; ++i) {
X       i__2 = i;
X       i__3 = i;
X       z__3.r = r.r * q[i__3].r - r.i * q[i__3].i, z__3.i = r.r * q[i__3].i
X               + r.i * q[i__3].r;
X       z__2.r = p[i__2].r - z__3.r, z__2.i = p[i__2].i - z__3.i;
X       z__4.r = u.r * w.r - u.i * w.i, z__4.i = u.r * w.i + u.i * w.r;
X       z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
X       u.r = z__1.r, u.i = z__1.i;
X       i__2 = i;
X       z__2.r = -w.r, z__2.i = -w.i;
X       z__1.r = z__2.r * u.r - z__2.i * u.i, z__1.i = z__2.r * u.i + z__2.i *
X                u.r;
X       q[i__2].r = z__1.r, q[i__2].i = z__1.i;
/* L130: */
X    }
X    i__1 = *n;
X    q[i__1].r = 1., q[i__1].i = 0.;
X    return 0;
X
/* % */
/*     ----------------------- */
/*     |    a newton-step    | */
/*     ----------------------- */
/* Subroutine */ int rfn_(p, q, r, z, n, s, v, b, c)
doublecomplex *p, *q, *r, *z;
integer *n;
doublecomplex *s, *v;
doublereal *b, *c;
{
X    /* System generated locals */
X    integer i__1, i__2;
X    doublecomplex z__1, z__2;
X
X    /* Builtin functions */
X    double z_abs();
X    void z_div();
X
X    /* Local variables */
X    extern /* Subroutine */ int dmag_();
X    static doublereal a;
X    static integer i;
X    static doublecomplex m, t, u, w;
X    extern doublereal dpreal_();
X
/* MP+ MULTIP COMPLEX p,q,r,z,m,s,t,u,v,w */
/* MP+ MULTIP REAL a,b,c */
X    /* Parameter adjustments */
X    --z;
X    --r;
X    --q;
X    --p;
X
X    /* Function Body */
X    a = z_abs(s);
X    if (a > 1.) {
X       goto L20;
X    }
X    v->r = 1., v->i = 0.;
X    u.r = 0., u.i = 0.;
X    w.r = u.r, w.i = u.i;
X    *c = 1.;
X    i = *n;
L10:
X    z__2.r = v->r * s->r - v->i * s->i, z__2.i = v->r * s->i + v->i * s->r;
X    i__1 = i;
X    z__1.r = z__2.r + p[i__1].r, z__1.i = z__2.i + p[i__1].i;
X    v->r = z__1.r, v->i = z__1.i;
X    z__2.r = u.r * s->r - u.i * s->i, z__2.i = u.r * s->i + u.i * s->r;
X    i__1 = i;
X    z__1.r = z__2.r + q[i__1].r, z__1.i = z__2.i + q[i__1].i;
X    u.r = z__1.r, u.i = z__1.i;
X    z__2.r = w.r * s->r - w.i * s->i, z__2.i = w.r * s->i + w.i * s->r;
X    i__1 = i;
X    z__1.r = z__2.r + r[i__1].r, z__1.i = z__2.i + r[i__1].i;
X    w.r = z__1.r, w.i = z__1.i;
X    *c = *c * a + dpreal_(&z[i]);
X    --i;
X    if (i > 0) {
X       goto L10;
X    }
X    goto L40;
L20:
X    z_div(&z__1, &c_b21, s);
X    t.r = z__1.r, t.i = z__1.i;
X    a = 1. / a;
X    v->r = p[1].r, v->i = p[1].i;
X    u.r = q[1].r, u.i = q[1].i;
X    w.r = r[1].r, w.i = r[1].i;
X    *c = dpreal_(&z[1]);
X    i__1 = *n;
X    for (i = 2; i <= i__1; ++i) {
X       z__2.r = v->r * t.r - v->i * t.i, z__2.i = v->r * t.i + v->i * t.r;
X       i__2 = i;
X       z__1.r = z__2.r + p[i__2].r, z__1.i = z__2.i + p[i__2].i;
X       v->r = z__1.r, v->i = z__1.i;
X       z__2.r = u.r * t.r - u.i * t.i, z__2.i = u.r * t.i + u.i * t.r;
X       i__2 = i;
X       z__1.r = z__2.r + q[i__2].r, z__1.i = z__2.i + q[i__2].i;
X       u.r = z__1.r, u.i = z__1.i;
X       z__2.r = w.r * t.r - w.i * t.i, z__2.i = w.r * t.i + w.i * t.r;
X       i__2 = i;
X       z__1.r = z__2.r + r[i__2].r, z__1.i = z__2.i + r[i__2].i;
X       w.r = z__1.r, w.i = z__1.i;
X       *c = *c * a + dpreal_(&z[i]);
/* L30: */
X    }
X    z__2.r = v->r * t.r - v->i * t.i, z__2.i = v->r * t.i + v->i * t.r;
X    z__1.r = z__2.r + 1., z__1.i = z__2.i;
X    v->r = z__1.r, v->i = z__1.i;
X    z__1.r = u.r * t.r - u.i * t.i, z__1.i = u.r * t.i + u.i * t.r;
X    u.r = z__1.r, u.i = z__1.i;
X    z__1.r = w.r * t.r - w.i * t.i, z__1.i = w.r * t.i + w.i * t.r;
X    w.r = z__1.r, w.i = z__1.i;
X    *c = *c * a + 1.;
L40:
X    dmag_(v, b);
X    if (*b == 0.) {
X       return 0;
X    }
X    dmag_(&u, &a);
X    if (a == 0.) {
X       return 0;
X    }
X    z_div(&z__1, v, &u);
X    t.r = z__1.r, t.i = z__1.i;
X    dmag_(&w, &a);
X    if (a == 0.) {
X       goto L50;
X    }
X    z_div(&z__1, &u, &w);
X    u.r = z__1.r, u.i = z__1.i;
X    z__1.r = u.r - t.r, z__1.i = u.i - t.i;
X    m.r = z__1.r, m.i = z__1.i;
X    dmag_(&m, &a);
X    if (a == 0.) {
X       return 0;
X    }
X    z_div(&z__1, &u, &m);
X    m.r = z__1.r, m.i = z__1.i;
X    z__2.r = m.r * t.r - m.i * t.i, z__2.i = m.r * t.i + m.i * t.r;
X    z__1.r = s->r - z__2.r, z__1.i = s->i - z__2.i;
X    s->r = z__1.r, s->i = z__1.i;
X    return 0;
L50:
X    z__1.r = s->r - t.r, z__1.i = s->i - t.i;
X    s->r = z__1.r, s->i = z__1.i;
X    return 0;
X
/* % */
/*     -------------------------------------------------------------- */
/*     |    estimate radius of smallest circle containing a zero    | */
/*     -------------------------------------------------------------- */
/* Subroutine */ int rad_(p, q, r, k, n, a)
doublecomplex *p, *q, *r;
integer *k, *n;
doublereal *a;
{
X    /* Initialized data */
X
X    static integer j = 1;
X
X    /* System generated locals */
X    integer i__1, i__2, i__3, i__4;
X    doublereal d__1, d__2;
X    doublecomplex z__1, z__2;
X
X    /* Builtin functions */
X    double z_abs();
X    void z_div();
X    double pow_dd();
X
X    /* Local variables */
X    extern /* Subroutine */ int dmag_();
X    static doublereal c, d, f;
X    static integer i, l, m;
X    static doublereal s;
X    static doublecomplex t;
X    extern doublereal dpnrt_();
X    static doublereal dm;
X
/* MP+ MULTIP COMPLEX p,q,r,t */
/* MP+ MULTIP REAL a,c,d,dm,f,s */
X    /* Parameter adjustments */
X    --r;
X    --q;
X    --p;
X
X    /* Function Body */
/*     print*,'rad: j=',j,' k=',k,' n=',n */
X    c = *n * z_abs(&p[1]);
X    f = 1. / *n;
X    if (j == *k) {
X       i__1 = *n;
X       for (i = 1; i <= i__1; ++i) {
X           i__2 = i;
X           i__3 = i;
X           q[i__2].r = r[i__3].r, q[i__2].i = r[i__3].i;
/* L30: */
X       }
X    } else {
X       j = 1;
X       l = 1;
X       i__1 = *n;
X       for (i = 2; i <= i__1; ++i) {
X           i__2 = l;
X           i__3 = i;
X           d__1 = (doublereal) l;
X           z__2.r = d__1 * p[i__3].r, z__2.i = d__1 * p[i__3].i;
X           d__2 = (doublereal) (*n);
X           z__1.r = z__2.r / d__2, z__1.i = z__2.i / d__2;
X           q[i__2].r = z__1.r, q[i__2].i = z__1.i;
X           l = i;
/* L10: */
X       }
X       i__1 = *n;
X       q[i__1].r = 1., q[i__1].i = 0.;
/*       c = n*abs(p(1)) */
/*       f = 1.0d0/n */
X       d__1 = z_abs(&p[1]);
X       *a = dpnrt_(&d__1, n);
X    }
/* L40: */
X    l = j + *k;
L50:
X    dmag_(&q[1], &dm);
X    if (dm == 0.) {
X       goto L90;
X    }
X    ++j;
X    z_div(&z__1, &p[1], &q[1]);
X    t.r = z__1.r, t.i = z__1.i;
X    i__1 = *n;
X    for (i = 2; i <= i__1; ++i) {
X       i__2 = i - 1;
X       i__3 = i;
X       i__4 = i;
X       z__2.r = t.r * q[i__4].r - t.i * q[i__4].i, z__2.i = t.r * q[i__4].i
X               + t.i * q[i__4].r;
X       z__1.r = p[i__3].r - z__2.r, z__1.i = p[i__3].i - z__2.i;
X       q[i__2].r = z__1.r, q[i__2].i = z__1.i;
/* L60: */
X    }
X    s = 1. / j;
X    d__1 = 1. - s;
X    d__2 = z_abs(&t);
X    f = pow_dd(&f, &d__1) * pow_dd(&d__2, &s);
X    d__1 = c / z_abs(&q[1]);
X    d = f * pow_dd(&d__1, &s);
X    if (d < *a) {
X       *a = d;
X    }
L70:
X    if (j < l) {
X       goto L50;
X    }
X    *k = j;
X    i__1 = *n;
X    for (i = 1; i <= i__1; ++i) {
X       i__2 = i;
X       i__3 = i;
X       r[i__2].r = q[i__3].r, r[i__2].i = q[i__3].i;
/* L80: */
X    }
X    return 0;
L90:
X    m = j;
L100:
X    ++j;
X    i__1 = *n;
X    for (i = 2; i <= i__1; ++i) {
X       i__2 = i - 1;
X       i__3 = i;
X       q[i__2].r = q[i__3].r, q[i__2].i = q[i__3].i;
/* L110: */
X    }
X    i__1 = *n;
X    q[i__1].r = 0., q[i__1].i = 0.;
X    dmag_(&q[1], &dm);
X    if (dm == 0.) {
X       goto L100;
X    }
X    s = 1. / j;
X    d__1 = m * s;
X    f = pow_dd(&f, &d__1);
X    d__1 = c / z_abs(&q[1]);
X    d = f * pow_dd(&d__1, &s);
X    if (d < *a) {
X       *a = d;
X    }
X    z_div(&z__1, &p[1], &q[1]);
X    t.r = z__1.r, t.i = z__1.i;
X    i__1 = *n;
X    for (i = 2; i <= i__1; ++i) {
X       i__2 = i - 1;
X       i__3 = i;
X       i__4 = i;
X       z__2.r = t.r * q[i__4].r - t.i * q[i__4].i, z__2.i = t.r * q[i__4].i
X               + t.i * q[i__4].r;
X       z__1.r = p[i__3].r - z__2.r, z__1.i = p[i__3].i - z__2.i;
X       q[i__2].r = z__1.r, q[i__2].i = z__1.i;
/* L120: */
X    }
X    i__1 = *n;
X    q[i__1].r = 1., q[i__1].i = 0.;
X    d__1 = z_abs(&t);
X    f *= pow_dd(&d__1, &s);
X    goto L70;
X
/* % */
/*     -------------------------------- */
/*     |    fast fourier transform    | */
/*     -------------------------------- */
/* Subroutine */ int ff_(a, b, n, l)
doublecomplex *a, *b;
integer *n, *l;
{
X    /* System generated locals */
X    integer i__1, i__2, i__3, i__4, i__5, i__6;
X    doublecomplex z__1, z__2;
X
X    /* Builtin functions */
X    void z_sqrt();
X
X    /* Local variables */
X    static integer maxi, i, j, k, m;
X    static doublecomplex s, t, w;
X    static doublereal dppic, theta;
X    static integer i1, i2, l1, l2, m0, m1, il;
X    static doublereal ct, st;
X    extern /* Double Complex */ int dpcmpl_();
X    static integer maxipj;
X    extern /* Subroutine */ int dpcssn_();
X    static integer maxipi1, maxipi2;
X
/* MP+ MULTIP COMPLEX a,b,s,t,w */
/* MP+ MULTIP REAL theta,ct,st */
X    /* Parameter adjustments */
X    --b;
X    --a;
X
X    /* Function Body */
X    dppic = 3.14159265358979323846;
X    maxi = 0;
X    maxipj = 0;
X    maxipi1 = 0;
X    maxipi2 = 0;
X    m = *n;
X    if (*l >= m) {
X       goto L20;
X    }
X    --m;
X    i__1 = m;
X    i__2 = *l;
X    for (j = *l; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
X       i__3 = *l;
X       for (i = 1; i <= i__3; ++i) {
X           i__4 = i;
X           i__5 = i;
X           i__6 = i + j;
X           z__1.r = a[i__5].r + a[i__6].r, z__1.i = a[i__5].i + a[i__6].i;
X           a[i__4].r = z__1.r, a[i__4].i = z__1.i;
/*  maxi = max(maxi,i) */
/* L10: */
X       }
X    }
X    m = *l;
L20:
X    theta = m * dppic / *l;
X    dpcssn_(&theta, &ct, &st);
X    dpcmpl_(&z__1, &ct, &st);
X    w.r = z__1.r, w.i = z__1.i;
X    l1 = *l - 1;
X    l2 = *l / 2;
X    if (m == *l) {
X       goto L40;
X    }
X    i__3 = l1;
X    i__2 = m;
X    for (j = m; i__2 < 0 ? j >= i__3 : j <= i__3; j += i__2) {
X       i__1 = m;
X       for (i = 1; i <= i__1; ++i) {
X           i__4 = i + j;
X           i__5 = i;
X           a[i__4].r = a[i__5].r, a[i__4].i = a[i__5].i;
/*  maxipj = max(maxipj,i+j) */
/* L30: */
X       }
X    }
L40:
X    m0 = m;
X    m /= 2;
X    m1 = m - 1;
X    if (m == 0) {
X       return 0;
X    }
X    s.r = 1., s.i = 0.;
X    t.r = -1., t.i = 0.;
X    i1 = 0;
X    i__1 = *l;
X    i__2 = m0;
X    for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
X       il = j + m1;
X       i2 = i1 + l2;
X       i__3 = il;
X       for (i = j; i <= i__3; ++i) {
X           k = i + m;
X           i__4 = i + i1;
X           i__5 = i;
X           i__6 = k;
X           z__2.r = s.r * a[i__6].r - s.i * a[i__6].i, z__2.i = s.r * a[i__6]
X                   .i + s.i * a[i__6].r;
X           z__1.r = a[i__5].r + z__2.r, z__1.i = a[i__5].i + z__2.i;
X           b[i__4].r = z__1.r, b[i__4].i = z__1.i;
X           i__4 = i + i2;
X           i__5 = i;
X           i__6 = k;
X           z__2.r = t.r * a[i__6].r - t.i * a[i__6].i, z__2.i = t.r * a[i__6]
X                   .i + t.i * a[i__6].r;
X           z__1.r = a[i__5].r + z__2.r, z__1.i = a[i__5].i + z__2.i;
X           b[i__4].r = z__1.r, b[i__4].i = z__1.i;
/*  maxipi1 = max(maxipi1,i+i1) */
/*  maxipi2 = max(maxipi2,i+i2) */
/* L50: */
X       }
X       z__1.r = w.r * s.r - w.i * s.i, z__1.i = w.r * s.i + w.i * s.r;
X       s.r = z__1.r, s.i = z__1.i;
X       z__1.r = w.r * t.r - w.i * t.i, z__1.i = w.r * t.i + w.i * t.r;
X       t.r = z__1.r, t.i = z__1.i;
X       i1 -= m;
/* L60: */
X    }
X    i__2 = *l;
X    for (i = 1; i <= i__2; ++i) {
X       i__1 = i;
X       i__3 = i;
X       a[i__1].r = b[i__3].r, a[i__1].i = b[i__3].i;
/* L70: */
X    }
X    z_sqrt(&z__1, &w);
X    s.r = z__1.r, s.i = z__1.i;
X    goto L40;
X
/* % */
/* Subroutine */ int dmag_(z, d)
doublecomplex *z;
doublereal *d;
{
X    /* Builtin functions */
X    double z_abs();
X
/* MP+ MULTIP COMPLEX z */
/* MP+ MULTIP REAL d */
X    *d = z_abs(z);
X    return 0;
X
/* % */
/* Subroutine */ int cc_(z, f, e, n, y, y0)
doublecomplex *z, *f;
doublereal *e;
integer *n;
doublereal *y, *y0;
{
X    /* System generated locals */
X    integer i__1, i__2;
X    doublecomplex z__1;
X
X    /* Local variables */
X    extern /* Subroutine */ int dmag_();
X    static integer i;
X    static doublereal s, t;
X    static doublecomplex x;
X
/* MP+ MULTIP COMPLEX f,z,x */
/* MP+ MULTIP REAL e,s,t,y,y0 */
X    /* Parameter adjustments */
X    --e;
X    --f;
X    --z;
X
X    /* Function Body */
X    i__1 = *n;
X    for (i = 1; i <= i__1; ++i) {
X       s = 0.;
X       i__2 = i;
X       x.r = z[i__2].r, x.i = z[i__2].i;
X       dmag_(&x, &t);
X       if (t > 1.) {
X           goto L20;
X       }
X       if (t > *y0) {
X           goto L30;
X       }
X       if (t == 0.) {
X           goto L30;
X       }
L10:
X       t *= *y;
X       z__1.r = *y * x.r, z__1.i = *y * x.i;
X       x.r = z__1.r, x.i = z__1.i;
X       s += -1.;
X       if (t <= *y0) {
X           goto L10;
X       }
X       goto L30;
L20:
X       t *= *y0;
X       z__1.r = *y0 * x.r, z__1.i = *y0 * x.i;
X       x.r = z__1.r, x.i = z__1.i;
X       s += 1.;
X       if (t > 1.) {
X           goto L20;
X       }
L30:
X       e[i] = s;
X       i__2 = i;
X       f[i__2].r = x.r, f[i__2].i = x.i;
/* L40: */
X    }
X    return 0;
X
/* % */
/* Subroutine */ int cef_(w, z, f, e, d, r, x, n, n1, t, t1, y, y0)
doublecomplex *w, *z, *f;
doublereal *e, *d, *r;
doublecomplex *x;
integer *n, *n1;
doublereal *t, *t1, *y, *y0;
{
X    /* System generated locals */
X    integer i__1, i__2, i__3;
X    doublereal d__1;
X    doublecomplex z__1;
X
X    /* Builtin functions */
X    double pow_dd();
X
X    /* Local variables */
X    extern /* Subroutine */ int dmag_();
X    static doublereal a, b;
X    static doublecomplex c;
X    static integer i;
X    static doublecomplex p;
X    static doublereal q, s, dm;
X
/* MP+ MULTIP COMPLEX f,w,z,c,x,p */
/* MP+ MULTIP REAL d,dm,e,a,b,q,r,s,t,t1,y,y0 */
X
X    /* Parameter adjustments */
X    --d;
X    --e;
X    --f;
X    --z;
X    --w;
X
X    /* Function Body */
X    p.r = 1., p.i = 0.;
X    q = 0.;
X    s = e[1];
X    z__1.r = *r * x->r, z__1.i = *r * x->i;
X    c.r = z__1.r, c.i = z__1.i;
X    if (*r > 1.) {
X       goto L50;
X    }
X    *t = 1.;
X    *t1 = 1.;
X    i__1 = *n;
X    for (i = 1; i <= i__1; ++i) {
X       a = e[i] + q;
X       if (a > s) {
X           s = a;
X       }
X       i__2 = i;
X       i__3 = i;
X       z__1.r = p.r * f[i__3].r - p.i * f[i__3].i, z__1.i = p.r * f[i__3].i
X               + p.i * f[i__3].r;
X       w[i__2].r = z__1.r, w[i__2].i = z__1.i;
X       d[i] = a;
X       dmag_(&z[*n1 - i], &dm);
X       *t = *r * *t + dm;
X       z__1.r = p.r * c.r - p.i * c.i, z__1.i = p.r * c.i + p.i * c.r;
X       p.r = z__1.r, p.i = z__1.i;
X       dmag_(&p, &dm);
X       if (dm > *y0) {
X           goto L20;
X       }
L10:
X       z__1.r = *y * p.r, z__1.i = *y * p.i;
X       p.r = z__1.r, p.i = z__1.i;
X       q += -1.;
X       dmag_(&p, &dm);
X       if (dm <= *y0) {
X           goto L10;
X       }
L20:
X       ;
X    }
L30:
X    if (q > s) {
X       s = q;
X    }
X    i__1 = *n;
X    for (i = 1; i <= i__1; ++i) {
X       a = d[i] - s;
X       i__2 = i;
X       i__3 = i;
X       d__1 = pow_dd(y, &a);
X       z__1.r = d__1 * w[i__3].r, z__1.i = d__1 * w[i__3].i;
X       w[i__2].r = z__1.r, w[i__2].i = z__1.i;
/* L40: */
X    }
X    a = q - s;
X    i__1 = *n1;
X    d__1 = pow_dd(y, &a);
X    z__1.r = d__1 * p.r, z__1.i = d__1 * p.i;
X    w[i__1].r = z__1.r, w[i__1].i = z__1.i;
X    return 0;
L50:
X    *t1 = *r;
X    b = 1. / *r;
X    *t = 0.;
X    i__1 = *n;
X    for (i = 1; i <= i__1; ++i) {
X       a = e[i] + q;
X       if (a > s) {
X           s = a;
X       }
X       i__2 = i;
X       i__3 = i;
X       z__1.r = p.r * f[i__3].r - p.i * f[i__3].i, z__1.i = p.r * f[i__3].i
X               + p.i * f[i__3].r;
X       w[i__2].r = z__1.r, w[i__2].i = z__1.i;
X       d[i] = a;
X       dmag_(&z[i], &dm);
X       *t = b * *t + dm;
X       z__1.r = p.r * c.r - p.i * c.i, z__1.i = p.r * c.i + p.i * c.r;
X       p.r = z__1.r, p.i = z__1.i;
X       dmag_(&p, &dm);
X       if (dm <= 1.) {
X           goto L70;
X       }
L60:
X       z__1.r = *y0 * p.r, z__1.i = *y0 * p.i;
X       p.r = z__1.r, p.i = z__1.i;
X       q += 1.;
X       dmag_(&p, &dm);
X       if (dm > 1.) {
X           goto L60;
X       }
L70:
X       ;
X    }
X    *t = b * *t + 1.;
X    goto L30;
X
/* % */
/* Subroutine */ int err_(z, n, s, t)
doublecomplex *z;
integer *n;
doublereal *s, *t;
{
X    /* System generated locals */
X    integer i__1;
X
X    /* Local variables */
X    extern /* Subroutine */ int dmag_();
X    static integer i;
X    static doublereal r, dm;
X
/* MP+ MULTIP COMPLEX z */
/* MP+ MULTIP REAL dm,r,s,t */
X
X    /* Parameter adjustments */
X    --z;
X
X    /* Function Body */
X    if (*s > 1.) {
X       goto L20;
X    }
X    *t = 1.;
X    i = *n;
L10:
X    dmag_(&z[i], &dm);
X    *t = *t * *s + dm;
X    --i;
X    if (i > 0) {
X       goto L10;
X    }
X    return 0;
L20:
X    r = 1. / *s;
X    dmag_(&z[1], t);
X    i__1 = *n;
X    for (i = 2; i <= i__1; ++i) {
X       dmag_(&z[i], &dm);
X       *t = *t * r + dm;
/* L30: */
X    }
X    *t = *t * r + 1.;
X    return 0;
X
X
doublereal dpreal_(a)
doublecomplex *a;
{
X    /* System generated locals */
X    doublereal ret_val;
X
X    ret_val = a->r;
X    return ret_val;
X
X
doublereal dpimag_(a)
doublecomplex *a;
{
X    /* System generated locals */
X    doublereal ret_val;
X
X    /* Builtin functions */
X    double d_imag();
X
X    ret_val = d_imag(a);
X    return ret_val;
X
X
doublereal dpnrt_(a, n)
doublereal *a;
integer *n;
{
X    /* System generated locals */
X    doublereal ret_val, d__1;
X
X    /* Builtin functions */
X    double pow_dd();
X
X    d__1 = 1. / *n;
X    ret_val = pow_dd(a, &d__1);
X    return ret_val;
X
X
/* Subroutine */ int dpcssn_(t, ct, st)
doublereal *t, *ct, *st;
{
X    /* Builtin functions */
X    double cos(), sin();
X
X    *ct = cos(*t);
X    *st = sin(*t);
X    return 0;
X
X
/* Double Complex */ int dpcmpl_( ret_val, a, b)
doublecomplex * ret_val;
doublereal *a, *b;
{
X    /* System generated locals */
X    doublecomplex z__1, z__2;
X
/* Not all compilers like the next line... */
/*     dpcmpl = (a,b) */
/* So do it a safe way... */
X    z__2.r = *b * (float)0., z__2.i = *b * (float)1.;
X    z__1.r = *a + z__2.r, z__1.i = z__2.i;
X     ret_val->r = z__1.r,  ret_val->i = z__1.i;
X    return ;
X
X
/* Subroutine */ int cz1test_()
{
X    /* System generated locals */
X    integer i__1;
X    doublecomplex z__1;
X
X    /* Builtin functions */
X    integer s_wsle(), do_lio(), e_wsle();
X    /* Subroutine */ int s_stop();
X    integer s_wsfe(), e_wsfe();
X
X    /* Local variables */
X    static integer i;
X    static doublecomplex p[20], w[1000], z[20];
X    extern /* Subroutine */ int czero1_();
X    static integer nd, nz;
X    extern /* Double Complex */ int dpcmpl_();
X
X    /* Fortran I/O blocks */
X    static cilist io___106 = { 0, 6, 0, 0, 0 };
X    static cilist io___111 = { 0, 6, 0, 0, 0 };
X    static cilist io___112 = { 0, 6, 0, "(/)", 0 };
X    static cilist io___113 = { 0, 6, 0, "(/)", 0 };
X    static cilist io___114 = { 0, 6, 0, 0, 0 };
X    static cilist io___115 = { 0, 6, 0, 0, 0 };
X    static cilist io___116 = { 0, 6, 0, "(/)", 0 };
X    static cilist io___117 = { 0, 6, 0, 0, 0 };
X    static cilist io___118 = { 0, 6, 0, 0, 0 };
X    static cilist io___119 = { 0, 6, 0, 0, 0 };
X    static cilist io___120 = { 0, 6, 0, 0, 0 };
X    static cilist io___121 = { 0, 6, 0, "(/)", 0 };
X    static cilist io___122 = { 0, 6, 0, 0, 0 };
X    static cilist io___123 = { 0, 6, 0, 0, 0 };
X    static cilist io___124 = { 0, 6, 0, 0, 0 };
X    static cilist io___125 = { 0, 6, 0, "(/)", 0 };
X
X
/* Workspace size... */
/* MP+ MULTIP COMPLEX z,p,w */
X
/* TEST 1 */
X
X    nd = 4;
X    p[0].r = (float)24., p[0].i = (float)0.;
X    p[1].r = (float)-20., p[1].i = (float)18.;
X    p[2].r = (float)10., p[2].i = (float)-15.;
X    p[3].r = (float)-5., p[3].i = (float)3.;
X    s_wsle(&io___106);
X    do_lio(&c__9, &c__1, "TEST 1: zeroes should be i, 2, 3, -4i...", 40L);
X    e_wsle();
X    czero1_(z, p, &nd, w, &c__1000, &nz);
X    if (nz <= 0) {
X       s_stop("czero1 failed!", 14L);
X    }
X    i__1 = nz;
X    for (i = 1; i <= i__1; ++i) {
/*        write(6,'(1x,1pd20.12,'' + '',1pd20.12,''i'')') z(i) */
/*         ! for MP version */
X       s_wsle(&io___111);
X       do_lio(&c__7, &c__1, (char *)&z[i - 1], (ftnlen)sizeof(doublecomplex))
X               ;
X       e_wsle();
X       s_wsfe(&io___112);
X       e_wsfe();
/* L1: */
X    }
X
/* TEST 2 */
X
/*     10       9         8          7           6           5            
4*/
/*    x   - 55 x  + 1320 x  - 18150 x  + 157773 x  - 902055 x  + 3416930 x
*/
X
/*             3             2 */
/*  - 8409500 x  + 12753576 x  - */
/*        10628640 x + 3628800 */
X    nd = 10;
X    p[9].r = -55., p[9].i = 0.;
X    p[8].r = 1320., p[8].i = 0.;
X    p[7].r = -18150., p[7].i = 0.;
X    p[6].r = 157773., p[6].i = 0.;
X    p[5].r = -902055., p[5].i = 0.;
X    p[4].r = 3416930., p[4].i = 0.;
X    p[3].r = -8409500., p[3].i = 0.;
X    p[2].r = 12753576., p[2].i = 0.;
X    p[1].r = -10628640., p[1].i = 0.;
X    p[0].r = 3628800., p[0].i = 0.;
X    s_wsfe(&io___113);
X    e_wsfe();
X    s_wsle(&io___114);
X    do_lio(&c__9, &c__1, "TEST 2: zeroes should be 1, 2, 3, 4, ..., 10.", 45L)
X           ;
X    e_wsle();
X    czero1_(z, p, &nd, w, &c__1000, &nz);
X    if (nz <= 0) {
X       s_stop("czero1 failed!", 14L);
X    }
X    i__1 = nz;
X    for (i = 1; i <= i__1; ++i) {
/*        write(6,'(1x,1pd20.12,'' + '',1pd20.12,''i'')') z(i) */
/*         ! for MP version */
X       s_wsle(&io___115);
X       do_lio(&c__7, &c__1, (char *)&z[i - 1], (ftnlen)sizeof(doublecomplex))
X               ;
X       e_wsle();
X       s_wsfe(&io___116);
X       e_wsfe();
/* L2: */
X    }
X
/* TEST 3 */
X
X    s_wsle(&io___117);
X    do_lio(&c__9, &c__1, "TEST 3: zeroes should be 1,2i,-3,-4i,5,6i,-7,-8i,",
X           49L);
X    e_wsle();
X    s_wsle(&io___118);
X    do_lio(&c__9, &c__1, "        8,9i,-10,-10i,11,12i,-13,-14i", 37L);
X    e_wsle();
X    nd = 16;
X    dpcmpl_(&z__1, &c_b142, &c_b143);
X    p[15].r = z__1.r, p[15].i = z__1.i;
X    dpcmpl_(&z__1, &c_b144, &c_b145);
X    p[14].r = z__1.r, p[14].i = z__1.i;
X    dpcmpl_(&z__1, &c_b146, &c_b147);
X    p[13].r = z__1.r, p[13].i = z__1.i;
X    dpcmpl_(&z__1, &c_b148, &c_b149);
X    p[12].r = z__1.r, p[12].i = z__1.i;
X    dpcmpl_(&z__1, &c_b150, &c_b151);
X    p[11].r = z__1.r, p[11].i = z__1.i;
X    dpcmpl_(&z__1, &c_b152, &c_b153);
X    p[10].r = z__1.r, p[10].i = z__1.i;
X    dpcmpl_(&z__1, &c_b154, &c_b155);
X    p[9].r = z__1.r, p[9].i = z__1.i;
X    dpcmpl_(&z__1, &c_b156, &c_b157);
X    p[8].r = z__1.r, p[8].i = z__1.i;
X    dpcmpl_(&z__1, &c_b158, &c_b159);
X    p[7].r = z__1.r, p[7].i = z__1.i;
X    dpcmpl_(&z__1, &c_b160, &c_b161);
X    p[6].r = z__1.r, p[6].i = z__1.i;
X    dpcmpl_(&z__1, &c_b162, &c_b163);
X    p[5].r = z__1.r, p[5].i = z__1.i;
X    dpcmpl_(&z__1, &c_b164, &c_b165);
X    p[4].r = z__1.r, p[4].i = z__1.i;
X    dpcmpl_(&z__1, &c_b166, &c_b167);
X    p[3].r = z__1.r, p[3].i = z__1.i;
X    dpcmpl_(&z__1, &c_b168, &c_b169);
X    p[2].r = z__1.r, p[2].i = z__1.i;
X    dpcmpl_(&z__1, &c_b170, &c_b171);
X    p[1].r = z__1.r, p[1].i = z__1.i;
X    p[0].r = (float)6.974263296e12, p[0].i = (float)0.;
X    czero1_(z, p, &nd, w, &c__1000, &nz);
X    if (nz <= 0) {
X       s_wsle(&io___119);
X       do_lio(&c__9, &c__1, "nz=", 3L);
X       do_lio(&c__3, &c__1, (char *)&nz, (ftnlen)sizeof(integer));
X       e_wsle();
X       s_stop("czero1 failed!", 14L);
X    }
X    i__1 = nz;
X    for (i = 1; i <= i__1; ++i) {
/*        write(6,'(1x,1pd20.12,'' + '',1pd20.12,''i'')') z(i) */
/*         ! for MP version */
X       s_wsle(&io___120);
X       do_lio(&c__7, &c__1, (char *)&z[i - 1], (ftnlen)sizeof(doublecomplex))
X               ;
X       e_wsle();
X       s_wsfe(&io___121);
X       e_wsfe();
/* L3: */
X    }
X
/* TEST 4 */
X
/*  PRODUCT(x-i,i,1,20)= */
X
/*x**20-210*x**19+20615*x**18-1256850*x**17+53327946*x**16-1672280820*x**1
5+4017*/
/*1771630*x**14-756111184500*x**13+11310276995381*x**12-135585182899530*x
**11+13*/
/*07535010540395*x**10-10142299865511450*x**9+63030812099294896*x**8-31133
364316*/
/*1390640*x**7+1206647803780373360*x**6-3599979517947607200*x**5+803781182
264505*/
/*1776*x**4-12870931245150988800*x**3+13803759753640704000*x**2-8752948036
761600*/
/* 000*x+2432902008176640000 */
X
X    s_wsle(&io___122);
X    do_lio(&c__9, &c__1, "TEST 4: zeroes should be i=1, 2, 3, ... , 20", 44L);
X    e_wsle();
X    nd = 20;
X    p[19].r = (float)-210., p[19].i = (float)0.;
X    p[18].r = (float)20615., p[18].i = (float)0.;
X    p[17].r = (float)-1256850., p[17].i = (float)0.;
X    p[16].r = (float)53327946., p[16].i = (float)0.;
X    p[15].r = (float)-1672280820., p[15].i = (float)0.;
X    p[14].r = (float)40171771630., p[14].i = (float)0.;
X    p[13].r = (float)-756111184500., p[13].i = (float)0.;
X    p[12].r = (float)11310276995381., p[12].i = (float)0.;
X    p[11].r = (float)-135585182899530., p[11].i = (float)0.;
X    p[10].r = (float)1307535010540395., p[10].i = (float)0.;
X    p[9].r = (float)-10142299865511450., p[9].i = (float)0.;
X    p[8].r = (float)63030812099294896., p[8].i = (float)0.;
X    p[7].r = (float)-311333643161390640., p[7].i = (float)0.;
X    p[6].r = (float)1206647803780373400., p[6].i = (float)0.;
X    p[5].r = (float)-3599979517947607200., p[5].i = (float)0.;
X    p[4].r = (float)8037811822645051800., p[4].i = (float)0.;
X    p[3].r = (float)-1.2870931245150989e19, p[3].i = (float)0.;
X    p[2].r = (float)1.3803759753640704e19, p[2].i = (float)0.;
X    p[1].r = (float)-8.7529480367616e18, p[1].i = (float)0.;
X    p[0].r = (float)2.43290200817664e18, p[0].i = (float)0.;
X    czero1_(z, p, &nd, w, &c__1000, &nz);
X    if (nz <= 0) {
X       s_wsle(&io___123);
X       do_lio(&c__9, &c__1, "nz=", 3L);
X       do_lio(&c__3, &c__1, (char *)&nz, (ftnlen)sizeof(integer));
X       e_wsle();
X       s_stop("czero1 failed!", 14L);
X    }
X    i__1 = nz;
X    for (i = 1; i <= i__1; ++i) {
/*        write(6,'(1x,1pd20.12,'' + '',1pd20.12,''i'')') z(i) */
/*         ! for MP version */
X       s_wsle(&io___124);
X       do_lio(&c__7, &c__1, (char *)&z[i - 1], (ftnlen)sizeof(doublecomplex))
X               ;
X       e_wsle();
X       s_wsfe(&io___125);
X       e_wsfe();
/* L4: */
X    }
X    return 0;
X
SHAR_EOF
chmod 0600 czero.c ||
echo 'restore of czero.c failed'
Wc_c="`wc -c < 'czero.c'`"
test 47619 -eq "$Wc_c" ||
        echo 'czero.c: original size 47619, current size' "$Wc_c"
fi
exit 0