THE	FOLLOWING UP TO $include IS PFCOMF
	implicit real*8 (a-h,o-z)
	parameter (ntot=20000,nvar=200,nterm=2500,nch=15,nchp=nch+1,
     +   nder=70, ntott=ntot/2, ntermt=nterm/2, mitay=9,
     +   ninch=8000, ninl=300, nchpl=127, nlab=20, ndind=9)
	character*1 a,ha,hc,hn,pv,ap,inpl,pl
	common/ch/a(nchpl),ha(26),hc(26),hn(10),pv(nchp,nvar)
     ,   ,ap(nchpl),inpl(ninch),pl(nchp,nlab)
	common/sc/nb,nd,nv,ns,np,nt,maxe,maxnt,maxite,itetop,ndump
     ,   ,minns,lintot,idflg,line,lines,linret,labels,tnow,tbeg
	common/i1/itb(nterm),ite(nterm)
	common/i2/itbt(ntermt),itet(ntermt)
	common/i3/ivn(ntot)
	common/i4/ivp(ntot)
	common/i5/ifnt(ntott),ifpt(ntott)
	common/i6/ider(3,nder),ittb(nvar),itte(nvar),ipv(nvar)
     ,   ,inpba(ninl),inpea(ninl),linlab(nlab),ipl(nlab)
     ,   ,idt(ndind,nvar),ios(nvar),iosi(nvar)
	common/f1/tm(nterm)
	common/f2/tmt(ntermt)
$include: 'pfcom.for'
	open(unit=7,file='algin',status='old')
	open(unit=8,file='algout',status='new')
	call datim
	tbeg=tnow
	do 4 i=1,nvar
4	ittb(i)=0
	ios(1)=1
	iosi(1)=1
	ipv(1)=1
	pv(1,1)='?'
	ittb(1)=1
	itte(1)=1
	itb(1)=1
	ite(1)=1
	ivn(1)=2
	ivp(1)=1
	nb=7
	nd=0
	nv=1
	ns=1
	np=2
	nt=1
	maxe=9999
	maxnt=0
	maxite=0
	itetop=1
	ndump=0
	minns=9999
	lintot=0
	idflg=1
	ha(1)='a'
	call defc(ha,26)
	hc(1)='A'
	call defc(hc,26)
	hn(1)='0'
	call defc(hn,10)
	ma=ntot
	mb=nvar
	mc=nterm
	md=ninch
	me=ninl
	mg=nder
	mh=nlab
	mi=ndind
	mj=mitay
	write(8,1)ma,mb,mc,md,me,mg,mh,mi,mj
1	format(' version 85/6/21 ',
     +   ' NTOT=',i6,' NVAR=',i4,' NTERM=',i5,' NINCH=',i5,/
     +  ,' NINL=',i6,' NDER=',i3,' NLAB=',i3,' NDIND=',i2
     +  ,' MITAY=',i3)
	call monitr
	call datim
	write(8,90)nd,nv,ns,maxnt,maxite,lintot
90	format(' nu. derivatives=',i3,'  nu. variables=',i3,
     +    '  nu. sums=',i3,/,' max nu. terms=',i5,
     +  '  max total nu. variables=',i6,'  input lines=',i5)
	t=tnow-tbeg
	write(8,91)t
91	format(' time used=',f9.2,' seconds')
	stop
	end
 
 
	subroutine add(is,js,t,ians)
$include: 'pfcom.for'
	nt1=nt+1
	ib=ittb(is)
	ie=itte(is)
	jb=ittb(js)
	je=itte(js)
	i=itb(ie)
	if(ivn(i).ne.2 .or. ivp(i).le.maxe)goto 3
	do 1 it=ib,ie
	i=itb(it)
	if(ivn(i).eq.2 .and. ivp(i).gt.maxe)goto 2
1	continue
2	ie=it-1
3	i=itb(je)
	if(ivn(i).ne.2 .or.ivp(i).le.maxe)goto 7
	do 4 it=jb,je
	i=itb(it)
	if(ivn(i).eq.2 .and. ivp(i).gt.maxe)goto 5
4	continue
5	je=it-1
7	if(ib.gt.ie)goto 70
	if(jb.gt.je)goto 12
	if(icomp(ib,jb))8,60,30
8	jl=jb
	if(jb.eq.je)goto 10
	if(icomp(ib,je))10,16,18
10	call cpynt(jb,je,t)
12	call cpynt(ib,ie,1.d0)
14	call defsum(ians,nt1,nt)
	return
16	save=tm(je)
	tm(je)=tm(ib)/t+tm(je)
	call cpynt(jb,je,t)
	tm(je)=save
	ib=ib+1
	jb=je+1
	goto 7
18	jr=je
20	if(jr-jl.eq.1)goto 28
	jc=(jl+jr)/2
	if(icomp(ib,jc))22,26,24
22	jl=jc
	goto 20
24	jr=jc
	goto 20
26	save=tm(jc)
	tm(jc)=tm(ib)/t+tm(jc)
	call cpynt(jb,jc,t)
	tm(jc)=save
	ib=ib+1
	jb=jc+1
	goto 7
28	call cpynt(jb,jl,t)
	jb=jr
30	il=ib
	if(ib.eq.ie)goto 31
	if(icomp(jb,ie))31,36,38
31	call cpynt(ib,ie,1.d0)
32	call cpynt(jb,je,t)
	goto 14
36	save=tm(ie)
	tm(ie)=tm(ie)+t*tm(jb)
	call cpynt(ib,ie,1.d0)
	tm(ie)=save
	jb=jb+1
	ib=ie+1
	goto 70
38	ir=ie
40	if(ir-il.eq.1)goto 48
	ic=(il+ir)/2
	if(icomp(jb,ic))42,46,44
42	il=ic
	goto 40
44	ir=ic
	goto 40
46	save=tm(ic)
	tm(ic)=tm(ic)+t*tm(jb)
	call cpynt(ib,ic,1.d0)
	tm(ic)=save
	jb=jb+1
	ib=ic+1
	goto 7
48	call cpynt(ib,il,1.d0)
	ib=ir
	goto 8
60	save=tm(ib)
	tm(ib)=tm(ib)+t*tm(jb)
	call cpynt(ib,ib,1.d0)
	tm(ib)=save
	ib=ib+1
	jb=jb+1
	goto 7
70	if(jb.gt.je)goto 14
	goto 32
	end
 
 
	subroutine bomb(i,h)
$include: 'pfcom.for'
	character*8 h
	goto(100,100,30,40,50),i
30	write(8,31)h
31	format(a8)
	goto 100
40	write(8,41)h,h
41	format(a8,' is too small for your problem.',
     +       '  Recompile with bigger ',a8)
	goto 100
50	continue
100	write(8,101)line,ap
101	format('input line number',i5,' was:',/,127a1)
	call datim
	write(8,90)nd,nv,ns,maxnt,maxite,line
90	format(/,'nu. derivatives=',i3,
     +      '   nu. variables=',i3,'nu. sums=',i3,'  max nu. terms=',
     +      i5,/,'  max total nu. variables=',i6,' input lines=',i5)
	if(ndump.ne.0)call dump
	t=tnow-tbeg
	write(8,91)t
91	format('time used=',f9.2,' seconds')
	stop
	return
	end
 
 
	subroutine bkup(j,k)
$include: 'pfcom.for'
	if(j.eq.k)goto 2
	do 1 i=j,k-1
	itb(i)=itb(i+1)
	ite(i)=ite(i+1)
1	tm(i)=tm(i+1)
2	k=k-1
	return
	end
 
 
	subroutine call
$include: 'pfcom.for'
	call nxt(3,nchpl,' ',j1)
	call numl(j1,nchpl,j2,j3,j4)
	if(linlab(j4).eq.-1)call bomb(3,'call    ')
	linret=line
	line=linlab(j4)
	return
	end
 
 
	subroutine coef(is,j3,c)
c  c=coefficient, j3=start of variable (on return)
$include: 'pfcom.for'
	f=1.d0
	sign=1.d0
	call nxtnb(is,nchpl,j2)
	call srch(j2,nchpl,'*',jj)
	if(jj.eq.-1)jj=nchpl
	call srch(j2,nchpl,'/',jk)
	if(jk.eq.-1)jk=nchpl
	call srchch(j2,nchpl,j3)
	j3=min0(jj,jk,j3)
	if(a(j2).eq.'-')goto 4
	if(a(j2).eq.'+')goto 10
	j4=j2
	goto 12
4	sign=-1.d0
10	call nxtnb(j2+1,nchpl,j4)
12	c=sign
	if(a(j4).eq.'.')goto 15
	call intval(j4,i)
	if(i.eq.-1 .and. j4.ne.j3)call bomb(3,'coef1   ')
	if(i.eq.-1)return
15	c=0.d0
	do 30 j1=j4,nchpl
	if(a(j1).eq.'.' .and. f.lt..5d0)call bomb(3,'coef2   ')
	if(a(j1).eq.'.')f=.1d0
	if(a(j1).eq.'.')goto 30
	call intval(j1,i)
	if(i.eq.-1)goto 40
	if(f.gt..5d0)c=10.d0*c+i
	if(f.lt..5d0)c=c+i*f
	if(f.lt..5d0)f=f*.1d0
30	continue
40	c=c*sign
	return
	end
 
 
	subroutine comput
$include: 'pfcom.for'
300	call nxt(1,nchpl,' ',j1)
	call numv(j1,nchpl,j2,j3,j4)
	call nxt(j3,nchpl,'=',j5)
	call srch(j5,nchpl,'/',j6)
	if(j6.eq.-1)goto 320
c  /
	call nxtnb(j5+1,j6-1,j7)
	if(a(1).eq.'d')goto 301
	call numv(j7,j6-1,j8,j9,j10)
	call numv(j6+1,nchpl,j11,j12,j13)
	if(ittb(j10).eq.0 .or. ittb(j13).eq.0)call bomb(3,'not sum ')
	call recip(j13,1)
	nt1=nt+1
	call mult(ittb(j10),itte(j10),ittb(1),itte(1))
	call defsum(j4,nt1,nt)
	call outqu(j4)
	minns=min0(minns,ios(1))
	it=ittb(1)
	itte(1)=it
	ite(it)=itb(it)
	k=itb(it)
	ivn(k)=2
	ivp(k)=0
	call pack
	return
c  compute sj4=dsj10/dfj14
301	call numv(j7+1,j6-1,j8,j9,j10)
	if(a(j7).ne.'d')call bomb(3,'comput1 ')
	if(idflg.ne.0)call dertab
	call nxt(j6,nchpl,'d',j11)
	call numv(j11+1,nchpl,j12,j13,j14)
	call srch(j13,nchpl,'=',j15)
	iord=1
	if(j15.eq.-1)goto 311
	call nxtnb(j15+1,nchpl,j16)
	call intval(j16,j17)
	if(j17.ne.-1)goto 360
	call numv(j16,nchpl,j19,j20,j18)
	if(ittb(j18).eq.0)call bomb(3,'comput14')
	iord=1.000000001d0*tm(ittb(j18))
	if(iord.gt.0)goto 311
	call bomb(3,'comput15')
360	call getexp(j15+1,iord,j16,0)
	if(j16.ne.-1 .or. iord.le.0)call bomb(3,'comput2 ')
311	do 313 ii=1,iord
	if(ii.eq.1)call deriv(j10,j14,1)
	if(ii.ne.1)call deriv(1,j14,1)
	call outqu(1)
	if(tm(ittb(1)).eq.0.)goto 314
313	continue
314	call swap(j4)
	return
c  not /
320	call srch(j5,nchpl,'+',j6)
	if(a(1).eq.'d')call bomb(3,'comput3 ')
	if(j6.ne.-1)goto 321
	call srch(j5,nchpl,'-',j6)
	tmm=-1.d0
	if(j6.ne.-1)goto 322
	call srch(j5,nchpl,'*',j6)
	if(j6.ne.-1)goto 339
	call srch(j5,nchpl,'^',j6)
	if(j6.ne.-1)goto 339
c  compute sj4=sj12
	call numv(j5,nchpl,j6,j11,j12)
	if(ittb(j12).eq.0)call bomb(3,'comput4 ')
	call copys(j12,j4)
	call outqu(j4)
	return
c  compute sj4=sj9  + or -  sj12
321	tmm=1.d0
322	call numv(j5+1,j6-1,j7,j8,j9)
	call numv(j6+1,nchpl,j10,j11,j12)
	call add(j9,j12,tmm,j4)
	call outqu(j4)
	return
c  compute sj4=sj9*sj12
339	call numv(j5+1,j6-1,j7,j8,j9)
	if(a(j6).eq.'^')goto 350
	if(a(j6+1).eq.'*')goto 349
	call numv(j6+1,nchpl,j10,j11,j12)
	nt1=nt+1
	call mult(ittb(j9),itte(j9),ittb(j12),itte(j12))
	call defsum(j4,nt1,nt)
	call outqu(j4)
	return
349	j6=j6+1
350	call numv(j6+1,nchpl,j10,j11,j12)
	call expnd(j9,j12,j4)
	return
	end
 
 
	subroutine copy(a,b,n)
	character*1 a,b
	dimension a(2),b(2)
	do 1 i=1,n
1	b(i)=a(i)
	return
	end
 
	subroutine copys(j1,j2)
$include: 'pfcom.for'
	nt1=nt+1
	call cpynt(ittb(j1),itte(j1),1.d0)
	call defsum(j2,nt1,nt)
	return
	end
 
 
	subroutine count
$include: 'pfcom.for'
	call nxt(3,nchpl,' ',j1)
	call numv(j1,nchpl,j3,j4,j5)
	do 12 i=nchpl,j4+1,-1
	if(a(i).ne.' ')goto 14
12	continue
	call bomb(3,'count1  ')
14	do 16 j=i,j4+1,-1
	if(a(j).eq.' ')goto 18
16	continue
18	call numv(j,i,j6,j7,j8)
	if(ittb(j8).eq.0)call bomb(3,'count2  ')
	if(ittb(j5).ne.0)goto 20
	k=itetop+1
	call incnt
	itb(nt)=k
	call nite(k)
	tm(nt)=itte(j8)+1-ittb(j8)
	ivn(k)=2
	ivp(k)=0
	call defsum(j5,nt,nt)
	return
20	ia=ittb(j5)
	ib=itte(j5)
	k=itb(ia)
	l=ite(ia)
	ivn(k)=2
	ivp(k)=0
	itte(j5)=ia
	ite(ia)=k
	tm(ia)=itte(j8)+1-ittb(j8)
	if(ia.ne.ib .or. k.ne.l)minns=min0(minns,ios(j5))
	call pack
	return
	end
 
 
	subroutine cpynt(ia,ib,tmm)
c  copy terms ia-ib to nt+1,..., coefficients are multiplied by tmm
$include: 'pfcom.for'
	if(ib.lt.ia)return
	if(ia.lt.1 .or. ib.gt.nt)call bomb(3,'cpynt2  ')
	do 15 i=ia,ib
	ja=itb(i)
	if(ivn(ja).eq.2 .and. ivp(ja).gt.maxe)goto 15
	k=itetop
	call incnt
	tm(nt)=tm(i)*tmm
	itb(nt)=k+1
	do 10 j=ja,ite(i)
	k=k+1
	ivn(k)=ivn(j)
10	ivp(k)=ivp(j)
	call nite(k)
15	continue
16	return
	end
 
 
	subroutine datim
$include: 'pfcom.for'
	call date(iy,mon,id)
	call time(ih,min,is,if)
	tnow=is+60.*(min+60.*ih)+if/100.
	write(8,1)iy,mon,id,ih,min,is
1	format(' job run  ',i4,'/',i2,'/',i2,i9,':',i2,':',i2)
	return
	end
 
 
	subroutine datime
$include: 'pfcom.for'
	call time(ih,min,is,if)
	tnow=is+60.*(min+60.*ih)+if/100.
	return
	end
 
 
	subroutine defc(a,n)
	character*1 a(2)
	do 1 i=2,n
	j=a(i-1)
1	a(i)=j+1
	return
	end
 
 
	subroutine defsum(is,ittbv,ittev)
$include: 'pfcom.for'
	if(ittev.ne.nt .or. ittev+1.ne.ittbv)goto 1
	call niltrm
	ittev=nt
1	if(ittbv.gt.ittev)call bomb(3,'defsum1 ')
	if(ittb(is).eq.0)goto 20
	j=ios(is)
	minns=min0(minns,j)
	if(ittbv.ge.ittb(is) .and. ittev.le.itte(is))goto 8
	if(j.eq.ns .and. ittbv.ge.ittb(is))goto 8
	if(ittbv.gt.itte(iosi(ns)))goto 4
	call bomb(3,'defsum2 ')
4	ios(is)=ns
	do 6 i=j,ns-1
	k=iosi(i+1)
	iosi(i)=k
6	ios(k)=i
	iosi(ns)=is
8	ittb(is)=ittbv
	itte(is)=ittev
	call pack
	return
20	if(ittbv.le.itte(iosi(ns)))call bomb(3,'defsum3 ')
	ns=ns+1
	ios(is)=ns
	iosi(ns)=is
	ittb(is)=ittbv
	itte(is)=ittev
	minns=min0(minns,ns)
	call pack
	return
	end
 
 
	subroutine delet
$include: 'pfcom.for'
c  delete sj4
	call nxt(4,nchpl,' ',j1)
1	call numv(j1,nchpl,j2,j3,j4)
	call delsum(j4)
	if(j3.eq.nchpl)goto 10
	call srchnb(j3+1,nchpl,j1)
	if(j1.ne.-1)goto 1
10	return
	end
 
 
	subroutine delsum(is)
$include: 'pfcom.for'
	if(ittb(is).eq.0)return
	j=ios(is)
	minns=min0(minns,j)
	if(j.eq.0 .or. ns.eq.0)call bomb(3,'delsum1 ')
	ns=ns-1
	if(j.gt.ns)goto 10
	do 2 i=j,ns
	k=iosi(i+1)
	iosi(i)=k
2	ios(k)=i
10	iosi(ns+1)=0
	ios(is)=0
	ittb(is)=0
	itte(is)=0
	call pack
	return
	end
 
 
	subroutine deriv(j10,ind,jans)
c  d(sum j10)/d(var ind) is put in sum j4
c  ider(1,nd) is dependent variable in derivative
c  ider(2,nd) is independent variable in derivative
c  ider(3,nd) is variable equal to derivative
$include: 'pfcom.for'
	character*1 bl
	common/i/if
	nt1=nt+1
	if(ittb(j10).lt.1 .or. itte(j10).gt.nt)call bomb(3,'deriv2  ')
	do 90 it=ittb(j10),itte(j10)
	ia=itb(it)
	ib=ite(it)
	do 85 ic=ia,ib
	ie=ivn(ic)
	if(ie.ne.ind)goto 40
	call incnt
	itb(nt)=itetop+1
	tm(nt)=tm(it)*ivp(ic)
	if=itetop
	do 36 i=ia,ib
	if(i.eq.ic)goto 32
	if=if+1
	ivn(if)=ivn(i)
	ivp(if)=ivp(i)
	goto 36
32	if(ivp(i).eq.1 .and. ia.ne.ib)goto 36
	if=if+1
	ivn(if)=ivn(i)
	ivp(if)=ivp(i)-1
36	continue
	call nite(if)
	if(if.lt.itb(nt))call bomb(3,'deriv3  ')
	goto 85
40	if(nd.eq.0)goto 85
	i1=1
41	j1=idt(i1,ie)
	if(j1.eq.0)goto 85
	k1=ider(2,j1)
	level=1
	if(k1.eq.ind)goto 50
	i2=1
42	j2=idt(i2,k1)
	if(j2.eq.0)goto 77
	k2=ider(2,j2)
	level=2
	if(k2.eq.ind)goto 50
	i3=1
43	j3=idt(i3,k2)
	if(j3.eq.0)goto 76
	k3=ider(2,j3)
	level=3
	if(k3.eq.ind)goto 50
	i4=1
44	j4=idt(i4,k3)
	if(j4.eq.0)goto 75
	k4=ider(2,j4)
	level=4
	if(k4.eq.ind)goto 50
	i5=1
45	j5=idt(i5,k4)
	if(j5.eq.0)goto 74
	k5=ider(2,j5)
	level=5
	if(k5.eq.ind)goto 50
	l1=ipv(ie)
	l2=ipv(k1)
	l3=ipv(k2)
	l4=ipv(k3)
	l5=ipv(k4)
	l6=ipv(k5)
	bl=' '
	write(8,49)(pv(i,ie),i=1,l1),bl,(pv(i,k1),i=1,l2),bl
     +  ,(pv(i,k2),i=1,l3),bl,(pv(i,k3),i=1,l4),bl,(pv(i,k4),i=1,l5)
     +  ,bl,(pv(i,k5),i=1,l6)
49	format('chain too long: ',90a1)
	call bomb(3,'deriv4  ')
50	call incnt
	if=itetop
	itb(nt)=if+1
	tm(nt)=tm(it)*ivp(ic)
	do 55 i=ia,ib
	if(i.eq.ic)goto 52
	if=if+1
	ivn(if)=ivn(i)
	ivp(if)=ivp(i)
	goto 55
52	if(ivp(i).eq.1)goto 55
	if=if+1
	ivn(if)=ivn(i)
	ivp(if)=ivp(i)-1
55	continue
	call deriv2(j1)
	if(level.eq.1)goto 77
	call deriv2(j2)
	if(level.eq.2)goto 76
	call deriv2(j3)
	if(level.eq.3)goto 75
	call deriv2(j4)
	if(level.eq.4)goto 74
	call deriv2(j5)
	i5=i5+1
	if(i5.le.ndind)goto 45
74	i4=i4+1
	if(i4.le.ndind)goto 44
75	i3=i3+1
	if(i3.le.ndind)goto 43
76	i2=i2+1
	if(i2.le.ndind)goto 42
77	i1=i1+1
	if(i1.le.ndind)goto 41
85	continue
90	continue
	call defsum(jans,nt1,nt)
	call order(jans)
	return
	end
 
 
	subroutine deriv2(j)
$include: 'pfcom.for'
	common/i/if
	if=if+1
	call nite(if)
	ivn(if)=ider(3,j)
	ivp(if)=1
	return
	end
 
 
	subroutine dertab
$include: 'pfcom.for'
	idflg=0
	do 302 i=1,ndind
	do 302 j=1,nvar
302	idt(i,j)=0
	do 304 i=1,nd
	idepv=ider(1,i)
	do 303 j=1,ndind
	if(idt(j,idepv).eq.0)goto 304
303	continue
	call bomb(4,'ndind   ')
304	idt(j,idepv)=i
	return
	end
 
 
	subroutine dfine
$include: 'pfcom.for'
	call nxt(4,nchpl,' ',j1)
	call nxtch(j1,nchpl,j2)
	call numv(j2,nchpl,j4,j5,j6)
	ittb6=nt+1
	call srch(j5,nchpl,'=',j7)
	if(j7.eq.-1)goto 410
	call term(j7+1)
	goto 412
410	call read(ib)
	if(np.ge.2)call prntap
	if(ib.eq.1)goto 412
	call term(1)
	goto 410
412	call defsum(1,ittb6,nt)
	call order(1)
	call outqu(1)
	call swap(j6)
	return
	end
 
	subroutine dump
$include: 'pfcom.for'
	character*1 bl
	data bl/' '/
	write(8,1)nd,nv,ns,np,nt,line,maxe,maxnt,maxite
1	format(/,'nd=',i3,' nv=',i3,' ns=',i3,' np=',i2,' nt=',i5,
     +     ' line=',i3,' maxe=',i2,' maxnt=',i5,' maxite=',i6)
	do 10 i=1,nv
	j=ipv(i)
	if(j.eq.nch)goto 3
	jp=j+1
	do 4 jj=jp,nch
4	pv(jj,i)=bl
3	write(8,5)i,(pv(k,i),k=1,nch),ittb(i),itte(i),ios(i),iosi(i)
5	format('i=',i5,'   pv=',15a1,'   ittb=',i5,'   itte=',i5
     +      ,'  ios=',i3,'  iosi=',i3)
10	continue
	write(8,11)
11	format(' ')
	write(8,13)
13	format('       i     tm     itb   ite')
	do 20 i=1,nt,4
	j=i+1
	k=i+2
	l=i+3
	write(8,21)i,tm(i),itb(i),ite(i),j,tm(j),itb(j),ite(j)
     +          ,k,tm(k),itb(k),ite(k),l,tm(l),itb(l),ite(l)
21	format(4(i9,f8.2,2i6))
20	continue
	write(8,11)
	do 22 i=1,nd
	j=ider(3,i)
	l=ipv(j)
22	write(8,24)ider(1,i),ider(2,i),(pv(m,j),m=1,l)
24	format('the derivative of',i4,' wrt',i4,' is ',15a1)
	write(8,11)
	write(8,40)(i,ivn(i),ivp(i),i=1,itetop)
40	format(10(1x,i5,2i3))
	return
	end
 
 
	subroutine expnd(jb,jp,jans)
$include: 'pfcom.for'
	ia=ittb(jb)
	if(ia.eq.0)call bomb(3,'expnd1  ')
	if(tm(ia).ne.1.d0)call bomb(3,'expnd2  ')
	ic=itb(ia)
	if(ivn(ic).ne.2 .or. ivp(ic).ne.0)call bomb(3,'expnd3  ')
	ic=ia+1
	id=itte(jb)
5	call cpynt(ia,ia,1.d0)
	if(ic.le.id)goto 10
6	call defsum(jans,nt,nt)
	return
10	ig=nt
	ie=ittb(jp)
	tp=tm(ie)
	if(tp.eq.0.d0)goto 6
	jprvb=ig
	jprve=ig
	tprv=1.d0
	do 20 kp=1,nb
	if(dabs(tp).lt.1.d-11)goto 25
	tprv=tprv*tp/kp
	tmt(kp)=tprv
	tp=tp-1.d0
	itbt(kp)=nt+1
	call mult(jprvb,jprve,ic,id)
	itet(kp)=nt
	if(itbt(kp).gt.nt)goto 25
	jprvb=itbt(kp)
20	jprve=nt
	kp=nb+1
25	kp=kp-1
	do 40 i=1,kp
	ih=itbt(i)
	ii=itet(i)
	tp=tmt(i)
	do 30 il=ih,ii
30	tm(il)=tm(il)*tp
40	continue
	call defsum(jans,ig,ii)
	call order(jans)
	return
	end
 
 
	subroutine extrct
$include: 'pfcom.for'
	call nxt(2,nchpl,' ',j1)
	call srch(j1,nchpl,'-',j6)
	if(j6.ne.-1)goto 300
	num=1
	do 1 i=2,10
	if(a(j1-1).eq.hn(i))num=i-1
1	continue
	call numv(j1,nchpl,j3,j4,j5)
	call nxt(j4,nchpl,'=',j6)
	call nxt(j6,nchpl,'m',j7)
	call numv(j7+1,nchpl,j8,j9,j10)
	k=ittb(j10)
	if(k.eq.0)call bomb(3,'extrct1 ')
	kk=tm(k)+.1
	if(kk.gt.0)goto 15
12	write(8,14)tm(k),(a(i),i=j8,j9)
14	format('bad term number  ',d12.3,5x,15a1)
	call bomb(3,'extrct2 ')
15	call nxt(j9+1,nchpl,'o',j11)
	call numv(j11+2,nchpl,j12,j13,j14)
	kl=itte(j14)+1-ittb(j14)
	if(kk.gt.kl)goto 12
	km=ittb(j14)+kk-1
	nt1=nt+1
	ktop=min0(km+num-1,itte(j14))
	call cpynt(km,ktop,1.d0)
	call defsum(j5,nt1,nt)
	call outqu(j5)
	return
300	call numv(j1,nchpl,j2,j3,j4)
	call nxt(j3,j6,'=',j5)
	call nxt(j5,j6,'t',j9)
	call nxt(j9,j6,' ',j10)
	call nxtnb(j10,j6,j11)
	call srch(j11,j6,' ',j12)
	if(j12.eq.-1)j12=j6
	call getint(j11,j12-1,j13)
	call srchnb(j6+1,nchpl,j8)
	call srch(j8,nchpl,' ',j14)
	call getint(j8,j14-1,j15)
	call nxt(j14,nchpl,'o',j19)
	if(a(j19+1).ne.'f' .and. a(j19+2).ne.' ')
     +    call bomb(3,'extrct5 ')
330	call numv(j19+2,nchpl,j16,j17,j18)
	if(ittb(j18).eq.0)call bomb(3,'extrct6')
	if(j15.gt.itte(j18)-ittb(j18)+1)j15=itte(j18)-ittb(j18)+1
	if(j15.lt.j13)call bomb(3,'extrct7 ')
	nt1=nt+1
	call cpynt(ittb(j18)+j13-1,ittb(j18)+j15-1,1.d0)
	call defsum(j4,nt1,nt)
	call outqu(j4)
	return
	end
 
 
	subroutine fmtc(it,if,il,posch)
$include: 'pfcom.for'
	character*1 aa,pa,posch
	common/pr1/pa(232),aa(22)
	character*127 cdum
	common/fmt/cdum
	il=22
	tt=tm(it)
	if(tt.ne.0.)goto 3
	aa(9)='0'
	if=9
	il=9
	return
3	ab=tm(it)*64.d0
	if(ab.gt.0.d0)j=ab+.5d0
	if(ab.lt.0.d0)j=ab-.5d0
	t=j/64.d0
	ac=dmax1(1.d0,dabs(tm(it)))
	if(dabs((t-tt)/ac).lt.1.d-11)tm(it)=t
	tt=tm(it)
	if(dabs(tt).gt.9999.9d0)goto 16
6	format(22a1)
14	if(dabs(tt).gt.1.d0)goto 20
	if(dabs(tt).gt..001d0)goto 24
16	write(cdum,18)tt
18	format(d22.13)
	read(cdum,6)aa
	ig=22
181	ig=ig-1
	if(aa(ig).ne.'E' .and. aa(ig).ne.'e')goto 181
	ih=ig
182	ih=ih-1
	if(aa(ih).eq.'0')goto 182
	do 183 i=ig,22
183	aa(ih+1+i-ig)=aa(i)
	il=ih+1+22-ig
	goto 32
20	write(cdum,22)tt
22	format(f22.12)
	goto 28
24	write(cdum,26)tt
26	format(f22.16)
28	read(cdum,6)aa
30	if(aa(il).ne.'0')goto 32
	il=il-1
	goto 30
32	do 34 if=1,il
	if(aa(if).ne.' ')goto 36
34	continue
36	if(tt.gt.0.d0)if=if-1
	if(tt.gt.0.d0)aa(if)=posch
	return
	end
 
 
	subroutine fort(it,ib)
$include: 'pfcom.for'
	character*1 aa,pa
	common/pr1/pa(232),aa(22)
	character*127 cdum
	common/fmt/cdum
	if(np.eq.0)return
	j=ib
	call fmtc(it,if,il,'+')
	do 42 i=if,il
	ib=ib+1
42	a(ib)=aa(i)
	ic=itb(it)
	id=ite(it)
	do 60 iv=ic,id
	ip=iabs(ivp(iv))
	in=ivn(iv)
	ie=ipv(in)
	iq=0
	if(ip.eq.1)goto 47
	write(cdum,43)ip
43	format(i20)
	read(cdum,6)aa
6	format(22a1)
	il=20
	do 44 if=1,20
	if(aa(if).ne.' ')goto 45
44	continue
45	iq=il+1-if+2
47	if(iq+ie+1.le.72)goto 52
	j=0
	write(8,50)(a(ich),ich=1,ib)
50	format(80a1)
	do 51 i=1,80
51	a(i)=' '
	a(6)='.'
	ib=11
52	ib=ib+1
	a(ib)='*'
	if(ivp(iv).lt.0)a(ib)='/'
	do 53 i=1,ie
	ib=ib+1
53	a(ib)=pv(i,in)
	if(iq.eq.0)goto 60
	ib=ib+1
	a(ib)='*'
	ib=ib+1
	a(ib)='*'
	do 54 i=if,il
	ib=ib+1
54	a(ib)=aa(i)
60	continue
	k=0
	if(a(j+2).eq.'1' .and. a(j+3).eq.'.' .and. a(j+4).eq.'*')k=1
	if(k.eq.0)goto 61
	a(j+2)=' '
	a(j+3)=a(j+1)
	a(j+4)=' '
	a(j+1)=' '
61	write(8,50)(a(ich),ich=1,ib)
	return
	end
 
 
	subroutine fortrn
$include: 'pfcom.for'
0	call nxt(1,nchpl,' ',j1)
	call numv(j1+1,nchpl,j2,j3,j4)
	if(np.eq.0 .or. ittb(j4).eq.0)goto 10
	do 1 i=1,nchpl
1	a(i)=' '
	ia=ipv(j4)
	ib=6
	do 2 i=1,ia
	ib=ib+1
2	a(ib)=pv(i,j4)
	ib=ib+1
	a(ib)='='
	ic=ittb(j4)
	id=itte(j4)
	j=0
	do 4 it=ic,id
	call fort(it,ib)
	do 3 i=1,nchpl
3	a(i)=' '
	a(6)='.'
	ib=7
	j=j+1
	if(j.lt.10)goto 4
	a(6)=' '
	call copy(pv(1,j4),a(7),ia)
	a(7+ia)='='
	call copy(pv(1,j4),a(8+ia),ia)
	ib=7+ia+ia
	j=0
4	continue
10	return
	end
 
 
	subroutine getint(j1,j2,j3)
c  on return j3=integer value of a(j1)-a(j2)
	j3=0
	if(j2.lt.j1)call bomb(3,'getint1 ')
	do 1 i=j1,j2
	call intval(i,j)
1	j3=10*j3+j
	return
	end
 
 
	subroutine getexp(j1,ie,k,if)
c  begin looking at a(j1), exponent put in ie, a(k) is next nonblank
c     (or k=-1 if only blanks are found)
$include: 'pfcom.for'
	ie=1
	call srchnb(j1,nchpl,k)
	if(k.eq.-1)return
	if(a(k).eq.'*'.and. a(k+1).ne.'*')return
	ise=1
	if(a(k).eq.'*'.and. a(k+1).eq.'*')goto 35
	if(a(k).eq.'/')return
	call srchch(j1,nchpl,j2)
	if(j2.eq.k .and. if.ne.0)return
	if(j2.ne.k)goto 9
32	call numv(j2,nchpl,j3,j4,j5)
	if(ittb(j5).eq.0)call bomb(3,'getexp10')
	tmittb=tm(ittb(j5))
	ie=tmittb+.5d0
	if(tmittb.lt.0.d0)ie=tmittb-.5d0
	re=ie
	if(dabs(re-tmittb).gt.1.d-11)call bomb(3,'getexp11')
	k=-1
	if(j4.ge.nchpl-1)return
	call srchnb(j4+1,nchpl,k)
	return
9	if(a(k).eq.'+')goto 20
	if(a(k).eq.'-')goto 10
	if(a(k).eq.'^')goto 30
	j3=k
1	me=0
5	call intval(j3,j4)
	me=10*me+j4
	j3=j3+1
	if(j3.gt.nchpl)goto 70
	if(a(j3).eq.' ')goto 50
	if(a(j3).eq.'*')goto 60
	if(a(j3).eq.'/')goto 60
	goto 5
10	ise=-1
20	if(k+1.gt.nchpl)call bomb(3,'getexp2 ')
	call nxtnb(k+1,nchpl,j3)
	if(j3.eq.-1)call bomb(3,'getexp4 ')
	goto 1
30	call nxtnb(k+1,nchpl,j3)
31	k=j3
	j2=j3
	call srchch(j3,nchpl,j4)
	if(j3.eq.j4)goto 32
	if(a(k).eq.'-')goto 10
	if(a(k).eq.'+')goto 20
	goto 1
35	if(k.ge.nchpl-1)call bomb(3,'getexp6 ')
	call nxtnb(k+2,nchpl,j3)
	if(j3.eq.-1)call bomb(3,'getexp8 ')
	goto 31
40	k=k+1
	goto 30
50	ie=me*ise
	call srchnb(j3,nchpl,k)
	return
60	k=j3
	if(j3.eq.nchpl)call bomb(3,'getexp9 ')
	ie=me*ise
	return
70	ie=me*ise
	k=-1
	return
	end
 
 
	function icomp(ia,ib)
$include: 'pfcom.for'
c  if ia should precede ib icomp=1
c  if ib should precede ia icomp=-1
c  if terms are the same except possibly coefficients, icomp=0
	ja=itb(ia)
	jb=itb(ib)
	ka=ite(ia)
	kb=ite(ib)
1	if(ivn(ja)-ivn(jb))2,8,9
2	if(ivp(ja))3,10,4
4	icomp=-1
	return
5	if(jb.gt.kb)goto 7
3	icomp=1
	return
6	ja=ja+1
11	jb=jb+1
	if(ja.gt.ka)goto 5
	if(jb.gt.kb)goto 4
	goto 1
7	icomp=0
	return
8	if(ivp(ja)-ivp(jb))3,6,4
9	if(ivp(jb))4,11,3
10	ja=ja+1
	if(ja.gt.ka)goto 3
	goto 1
	end
 
 
	subroutine iftest(go)
$include: 'pfcom.for'
	character*1 aj1,aj2
	go=0.
	call nxt(2,nchpl,' ',j1)
	call nxt(j1,nchpl,'.',j)
	call numv(j1,j-1,j3,j4,j5)
	aj1=a(j+1)
	aj2=a(j+2)
	if(aj1.eq.'n')goto 20
	if(aj1.eq.'l' .or. aj1.eq.'g')goto 22
	if(aj1.ne.'e')call bomb(3,'iftest1 ')
	if(aj2.ne.'q')call bomb(3,'iftest2 ')
4	if(a(j+3).ne.'.')call bomb(3,'iftest3 ')
	call numv(j+4,nchpl,j6,j7,j8)
	if(j5.ne.j8)goto 10
	if(aj1.eq.'n')return
5	call nxtch(j7+1,nchpl,j9)
	do 6 j=j9,nchpl
	a(j+1-j9)=a(j)
6	a(j)=' '
	go=1.d0
	call copy(a,ap,nchpl)
	return
10	k1=itb(ittb(j5))
	k2=ite(itte(j5))
	k3=itb(ittb(j8))
	k4=ite(itte(j8))
	t5=tm(ittb(j5))
	t8=tm(ittb(j8))
	if(aj1.eq.'l')goto 23
	if(aj1.eq.'g')goto 25
	if(k2-k1.ne.k4-k3)goto 30
	k=k3
	do 11 i=k1,k2
	if(ivn(i).ne.ivn(k))goto 30
	if(ivp(i).ne.ivp(k))goto 30
11	k=k+1
	k1=ittb(j5)
	k2=itte(j5)
	k3=ittb(j8)
	k4=itte(j8)
	if(k2-k1.ne.k4-k3)goto 30
	k=k3
	do 12 i=k1,k2
	den=dabs(tm(i))+dabs(tm(k))
	if(den.lt.1.d-11)goto 12
	if(dabs(tm(i)-tm(k))/den.gt.1.d-11)goto 30
12	k=k+1
	if(aj1.eq.'n')return
	goto 5
20	if(a(j+2).ne.'e')call bomb(3,'iftest4 ')
	goto 4
22	if(aj2.eq.'e' .or. aj2.eq.'t')goto 4
	call bomb(5,'iftest5 ')
23	if(aj2.eq.'e')goto 24
	if(t5.lt.t8)goto 5
	return
24	if(t5.le.t8)goto 5
	return
25	if(aj2.eq.'e')goto 26
	if(t5.gt.t8)goto 5
	return
26	if(t5.ge.t8)goto 5
	return
30	if(aj1.eq.'n')goto 5
	return
	end
 
 
	subroutine incnt
$include: 'pfcom.for'
	nt=nt+1
	maxnt=max0(nt,maxnt)
	if(maxnt.gt.nterm)call bomb(4,'nterm   ')
	return
	end
 
 
	subroutine intval(j1,i)
c  on return i=integer value of a(j1)
$include: 'pfcom.for'
	do 1 j=1,10
	if(a(j1).eq.hn(j))goto 2
1	continue
	i=-1
	return
2	i=j-1
	return
	end
 
 
	subroutine monitr
$include: 'pfcom.for'
	call readin
	ntp=0
	tprev=tnow
10	if(ha(1).ne.'a')call bomb(3,'monitr1 ')
	if(ittb(2).ne.0)call bomb(3,'useorder')
	if(ha(1).eq.'a')goto 11
	if(ns.lt.2)goto 14
	do 20 ms=2,ns
	is=iosi(ms-1)
	js=iosi(ms)
	if(is.eq.0 .or. js.eq.0)call bomb(3,'mon1    ')
	if(ios(is).ne.ms-1 .or. ios(js).ne.ms)call bomb(3,'mon2    ')
	if(ittb(is).eq.0 .or. ittb(js).eq.0)call bomb(3,'mon3    ')
	ia=itte(is)
	ib=ittb(js)
	if(ite(ia)+1.ne.itb(ib))call bomb(3,'mon4    ')
	do 16 it=ib,itte(js)
	if(itb(it).gt.ite(it))call bomb(3,'mon5    ')
	if(dabs(tm(it)).gt.1.d-11)goto 15
	if(itb(it).ne.ite(it))call bomb(3,'mon6    ')
	if(ivn(itb(it)).ne.2)call bomb(3,'mon7    ')
	if(itb(it).gt.ite(it))call bomb(3,'mon8    ')
15	do 17 i=itb(it),ite(it)
	if(ivn(i).eq.1)call bomb(3,'mon9    ')
	if(ittb(ivn(i)).ne.0)call bomb(3,'mon10   ')
	if(i.eq.itb(it))goto 17
	if(ivp(i).eq.0)call bomb(3,'mon11   ')
	if(ivn(i-1).ge.ivn(i))call bomb(3,'mon12   ')
17	continue
16	continue
20	continue
14	if(ns.eq.0)call bomb(3,'mon13   ')
	ic=iosi(ns)
	if(itte(ic).ne.nt)call bomb(3,'mon14   ')
	if(ite(itte(ic)).ne.itetop)call bomb(3,'mon15   ')
11	call read(ib)
	if(ib.eq.1 .and. np.ge.1)call prntap
	if(ib.eq.1)goto 10
	call datime
	t=tnow-tprev
	ttot=tnow-tbeg
	tprev=tnow
	ntmntp=nt-ntp
	if(np.ge.3)write(8,3)t,ttot,ntmntp,nt,itetop,maxnt,maxite
3	format(10x,'t=',f6.2,' total t=',f8.1,' ntmntp=',i5,' nt=',i5,
     +       ' itetop=',i6,' maxnt=',i5,' maxite=',i6)
	if(np.ge.2)call prntap
5	ntp=nt
	call copy(a,ap,nchpl)
	if(a(1).eq.'c')goto 300
	if(a(1).eq.'d')goto 400
	if(a(1).eq.'e')goto 500
	if(a(1).eq.'f')goto 600
	if(a(1).eq.'g')goto 700
	if(a(1).eq.'i')goto 900
	if(a(1).eq.'l')goto 10
	if(a(1).eq.'m')goto 1300
	if(a(1).eq.'n')goto 1400
	if(a(1).eq.'o')goto 1500
	if(a(1).eq.'p')goto 1600
	if(a(1).eq.'r')goto 1800
	if(a(1).eq.'s')goto 1900
	if(a(1).eq.'t')goto 2000
	if(a(1).eq.'z')goto 2600
	call bomb(3,'bad cmnd')
300	if(a(2).eq.'o')goto 310
	if(a(2).ne.'a')call bomb(3,'monitr2 ')
	call call
	goto 10
310	if(a(3).eq.'m')goto 320
	if(a(3).ne.'u')call bomb(3,'monitr3 ')
	call count
	goto 10
320	call comput
	call pack
	goto 10
400	if(a(2).eq.'u')goto 470
	if(a(2).eq.'i')goto 320
	if(a(2).ne.'e')call bomb(3,'monitr4 ')
	if(a(3).eq.'l')goto 450
	if(a(3).eq.'f')goto 440
	call bomb(3,'monitr5 ')
440	call dfine
	goto 10
450	call delet
	call pack
	goto 10
470	call dump
	goto 10
500	if(a(2).eq.'n')return
	if(a(2).ne.'x')call bomb(3,'monitr6 ')
	call extrct
	goto 10
600	if(a(2).eq.'l')goto 610
	call fortrn
	goto 10
610	call readin
	goto 10
700	call nxt(1,nchpl,' ',j9)
	call numl(j9,nchpl,j10,j11,j12)
	if(linlab(j12).eq.-1)write(8,710)
	if(linlab(j12).eq.-1)call bomb(3,'monitr7 ')
710	format('label not found')
	line=linlab(j12)-1
	goto 10
900	call iftest(go)
	if(go.eq.1.d0)goto 5
	goto 10
1300	call srch(1,nchpl,'=',j)
	call srchnb(j+1,nchpl,k)
	call getexp(k,maxe,j,0)
	goto 10
1400	call nxt(3,nchpl,'=',j1)
	call getexp(j1+1,j3,j2,0)
	if(a(2).eq.'d')goto 1405
	if(a(2).eq.'p')goto 1410
	if(a(2).eq.'b')goto 1415
	call bomb(3,'monitr8 ')
1405	if(a(3).eq.'u')goto 1406
	call bomb(3,'monitr9 ')
1406	ndump=j3
	goto 10
1410	np=j3
	goto 10
1415	nb=j3
	goto 10
1500	if(nv.ne.1)call bomb(3,'putfirst')
	call ordset
	goto 10
1600	if(a(3).eq.'i')goto 1620
	if(a(3).ne.'o')call bomb(3,'monitr10')
1605	call read(ib)
	if(np.ge.2)call prntap
	if(a(1).ne.'r' .or. a(2).ne.'e' .or. a(3).ne.'t')goto 1605
	goto 10
1620	call nxt(4,nchpl,' ',j1)
	call numv(j1,nchpl,j2,j3,j4)
	if(np.eq.1)write(8,1)
	if(np.eq.1)call prntap
1	format(130a1)
	call print(j4)
	goto 10
1800	if(a(2).ne.'e' .or. a(3).ne.'t')call bomb(3,'return? ')
	line=linret
	goto 10
1900	call sub
	call pack
	goto 10
2000	if(a(2).eq.'a')goto 2010
	call thderv
	goto 10
2010	call taylor
	goto 10
2600	call numv(2,nchpl,j1,j2,j3)
	call zsub(j3)
	call pack
	goto 10
	end
 
 
	subroutine mult(ia,ib,ic,id)
c  multiply terms ia-ib by terms ic-id
c  result goes to nt+1,.....
$include: 'pfcom.for'
	common/nwt/iitb,iite,jitb,jite,ntb,if
	if(ia.eq.0 .or. ib.eq.0 .or. ic.eq.0 .or. id.eq.0)
     +   call bomb(3,'mult1   ')
	ntb=nt+1
	if=0
	do 20 i=ia,ib
	iitb=itb(i)
	iite=ite(i)
	do 10 j=ic,id
	jitb=itb(j)
	jite=ite(j)
	if(ivn(iitb)+ivn(jitb).eq.4 .and. ivp(iitb)+ivp(jitb).gt.maxe)
     +                                                       goto 20
	call sub2(tm(i)*tm(j))
10	continue
20	continue
	return
	end
 
 
	subroutine niltrm
$include: 'pfcom.for'
	k=itetop+1
	call incnt
	tm(nt)=0.d0
	itb(nt)=k
	call nite(k)
	ivn(k)=2
	ivp(k)=0
	return
	end
 
 
	subroutine nite(j)
$include: 'pfcom.for'
	itetop=j
	ite(nt)=j
	maxite=max0(maxite,j)
	if(maxite.gt.ntot)call bomb(4,'ntot    ')
	return
	end
 
 
	subroutine numl(j1,j2,j3,j4,j5)
c  search j1-j2  find name in j3-j4, j5 is index of label
$include: 'pfcom.for'
	character*1 c
	if(j2.lt.j1)call bomb(3,'numl    ')
	call nxtch(j1,j2,j3)
	do 10 j=j3,j2
	c=a(j)
	do 5 k=1,26
	if(c.eq.ha(k) .or. c.eq.hc(k))goto 10
5	continue
	do 6 k=1,10
	if(c.eq.hn(k))goto 10
6	continue
	if(c.eq.' ' .or. c.eq.'+' .or. c.eq.'-')goto 12
	if(c.eq.'*' .or. c.eq.'/' .or. c.eq.'=' .or. c.eq.'^')goto 12
10	continue
	j4=j2
	goto 14
12	j4=j-1
14	l=j4-j3+1
	do 18 i=1,labels
	if(ipl(i).ne.l)goto 18
	do 16 j=1,l
	if(pl(j,i).ne.a(j+j3-1))goto 18
16	continue
	j5=i
	return
18	continue
	labels=labels+1
	if(labels.gt.nlab)call bomb(4,'nlab    ')
	j5=labels
	ipl(j5)=l
	call copy(a(j3),pl(1,j5),l)
	return
	end
 
 
	subroutine numv(j1,j2,j3,j4,j5)
c  search j1-j2  find name in j3-j4, j5 is index of variable
$include: 'pfcom.for'
	character*1 c
	if(j2.lt.j1)call bomb(3,'numv    ')
	call nxtch(j1,j2,j3)
	do 10 j=j3,j2
	c=a(j)
	do 5 k=1,26
	if(c.eq.ha(k) .or. c.eq.hc(k))goto 10
5	continue
	do 6 k=1,10
	if(c.eq.hn(k))goto 10
6	continue
	if(c.eq.' ' .or. c.eq.'+' .or. c.eq.'-')goto 12
	if(c.eq.'*' .or. c.eq.'/' .or. c.eq.'=' .or. c.eq.'^')goto 12
10	continue
	j4=j2
	goto 14
12	j4=j-1
14	l=j4-j3+1
	do 28 ii=ns,1,-1
	i=iosi(ii)
	if(ipv(i).ne.l)goto 28
	do 26 j=1,l
	if(pv(j,i).ne.a(j+j3-1))goto 28
26	continue
	j5=i
	return
28	continue
	do 18 i=1,nv
	if(ittb(i).ne.0)goto 18
	if(ipv(i).ne.l)goto 18
	do 16 j=1,l
	if(pv(j,i).ne.a(j+j3-1))goto 18
16	continue
	j5=i
	return
18	continue
	nv=nv+1
	j5=nv
	if(nv.gt.nvar)call bomb(4,'nvar    ')
	ipv(nv)=l
	call copy(a(j3),pv(1,nv),l)
	return
	end
 
 
	subroutine nxt(j1,j2,h,j3)
	character*1 h
	call srch(j1,j2,h,j3)
	if(j3.ne.-1)return
	write(8,1)h,j1,j2
1	format(' nxt did not find ',a1,' in ',2i4)
	call bomb(3,'nxt     ')
	return
	end
 
 
	subroutine nxtch(j1,j2,j3)
	call srchch(j1,j2,j3)
	if(j3.eq.-1)call bomb(3,'nxtch   ')
	return
	end
 
 
	subroutine nxtnb(j1,j2,j3)
	call srchnb(j1,j2,j3)
	if(j3.eq.-1)call bomb(3,'nxtnb   ')
	end
 
 
	subroutine order(is)
$include: 'pfcom.for'
	ina=ittb(is)
	inb=itte(is)
	if(inb.lt.ina .or. ina.eq.0)call bomb(3,'order1  ')
	it=ina
2	ia=itb(it)
	if(dabs(tm(it)).le.1.d-11)goto 41
	ib=ite(it)
3	if(ib.eq.ia)goto 25
	i=ia
4	if(ivp(i).eq.0)goto 20
	i=i+1
	if(i.le.ib)goto 4
	if(ia.eq.ib)goto 25
23	i=ia+1
5	if(ivn(i)-ivn(i-1))10,7,6
6	i=i+1
	if(i.le.ib)goto 5
	goto 40
7	ivp(i-1)=ivp(i-1)+ivp(i)
	if(i.eq.ib)goto 9
	ib=ib-1
	do 8 j=i,ib
	ivn(j)=ivn(j+1)
8	ivp(j)=ivp(j+1)
	if(ivp(i-1).ne.0)goto 5
	goto 3
9	ib=ib-1
	if(ivp(ib).ne.0)goto 40
	goto 3
10	if(i.eq.ia+1)goto 19
	ic=i-2
11	if(ivn(i)-ivn(ic))12,13,16
12	if(ic.eq.ia)goto 18
	ic=ic-1
	goto 11
13	ivp(ic)=ivp(ic)+ivp(i)
	if(i.eq.ib)goto 15
	ib=ib-1
	do 14 j=i,ib
	ivn(j)=ivn(j+1)
14	ivp(j)=ivp(j+1)
	if(ivp(ic).eq.0)goto 3
	goto 5
15	ib=ib-1
	if(ivp(ic).ne.0)goto 40
	goto 3
16	ip=ivp(i)
	in=ivn(i)
	do 17 j=i,ic+2,-1
	ivn(j)=ivn(j-1)
17	ivp(j)=ivp(j-1)
	ivn(ic+1)=in
	ivp(ic+1)=ip
	goto 6
18	ic=ia-1
	goto 16
19	ip=ivp(i)
	in=ivn(i)
	ivp(i)=ivp(ia)
	ivn(i)=ivn(ia)
	ivp(ia)=ip
	ivn(ia)=in
	goto 6
20	if(i.eq.ib)goto 22
	ib=ib-1
	do 21 j=i,ib
	ivn(j)=ivn(j+1)
21	ivp(j)=ivp(j+1)
	goto 4
22	ib=ib-1
	if(ib.gt.ia)goto 23
	ib=ia
25	if(ivp(ia).eq.0)ivn(ia)=2
40	ite(it)=ib
	if(ivn(ia).ne.2 .or. ivp(ia).le.maxe)goto 52
41	call bkup(it,inb)
	goto 53
52	it=it+1
53	if(it.le.inb)goto 2
c  sort terms ina-inb
	if(ina.lt.inb)goto 70
	if(ina.eq.inb)goto 82
51	inb=ina
	tm(ina)=0.
	k=itb(ina)
	ite(ina)=k
	ivn(k)=2
	ivp(k)=0
	goto 82
70	it=ina
68	it=it+1
69	if(it.gt.inb)goto 80
	if(dabs(tm(it)).le.1.d-11)goto 75
	if(icomp(ina,it))71,73,74
71	call tzap(it,ina)
	goto 68
73	tm(ina)=tm(ina)+tm(it)
75	call bkup(it,inb)
	goto 69
74	if(icomp(it-1,it))62,61,68
61	tm(it-1)=tm(it-1)+tm(it)
	call bkup(it,inb)
	goto 69
62	la=ina
	lb=it-1
63	if(lb-la.eq.1)goto 64
	lc=(la+lb)/2
	if(icomp(lc,it))65,66,67
64	call tzap(it,lb)
	goto 68
65	lb=lc
	goto 63
66	tm(lc)=tm(lc)+tm(it)
	call bkup(it,inb)
	goto 69
67	la=lc
	goto 63
80	if(ina.gt.inb)goto 51
82	call defsum(is,ina,inb)
	return
	end
 
 
	subroutine ordset
$include: 'pfcom.for'
	call nxt(2,nchpl,' ',j1)
	call srchnb(j1,nchpl,j2)
	if(j2.eq.-1)goto 20
	if(a(j2).eq.';')goto 20
	ntsav=nt
	itesav=itetop
	call term(j1)
	nt=ntsav
	itetop=itesav
	return
20	if=1
	do 10 il=1,300
	call read(ib)
	if(np.ge.2)call prntap
	if(ib.eq.1)return
	is=1
2	call srchch(is,nchpl,j1)
	if(j1.eq.-1)goto 10
	call numv(j1,nchpl,j2,j3,j4)
	if=if+1
	if(j4.lt.if)call bomb(3,'ordset1 ')
	if(j4.eq.if)goto 8
	do 5 i=nv,if,-1
	ipv(i)=ipv(i-1)
	do 4 ic=1,nch
4	pv(ic,i)=pv(ic,i-1)
5	continue
	call copy(a(j2),pv(1,if),j3-j2+1)
	ipv(if)=j3-j2+1
8	is=j3+1
	if(is.le.nchpl)goto 2
10	continue
	return
	end
 
 
	subroutine outqu(j5)
$include: 'pfcom.for'
1	do 6 it=ittb(j5),itte(j5)
	do 5 iv=itb(it),ite(it)
	if(ittb(ivn(iv)).ne.0)goto 9
5	continue
6	continue
	return
9	jv=ivn(iv)
	call subs(jv,jv,j5)
	goto 1
	end
 
 
	subroutine pack
$include: 'pfcom.for'
	if(minns.le.ns)goto 9
	is=iosi(ns)
	nt=itte(is)
	itetop=ite(nt)
	minns=ns+1
	return
9	if(minns.le.0)call bomb(3,'pack4   ')
	js=minns
	kt=0
	kv=0
	if(js.eq.1)goto 13
	kt=itte(iosi(js-1))
	kv=ite(kt)
13	do 30 ms=js,ns
	is=iosi(ms)
	if(is.le.0)call bomb(3,'pack6   ')
	ia=ittb(is)
	ib=itte(is)
	lv=0
	lt=0
	do 14 it=ia,ib
	if(dabs(tm(it)).lt.1.d-11)goto 14
	lt=lt+1
	if(lt.gt.ntermt)call bomb(4,'ntermt  ')
	tmt(lt)=tm(it)
	itbt(lt)=lv+1
	do 12 i=itb(it),ite(it)
	if(ivp(i).ne.0)goto 31
	do 32 j=itb(it),ite(it)
	if(ivp(j).ne.0)goto 12
32	continue
	lv=lv+1
	ifnt(lv)=2
	ifpt(lv)=0
	goto 33
31	lv=lv+1
	ifnt(lv)=ivn(i)
	ifpt(lv)=ivp(i)
12	continue
33	itet(lt)=lv
	if(lv.gt.ntott)call bomb(4,'ntott   ')
14	continue
	if(lv.ne.0)goto 15
	kt=kt+1
	ittb(is)=kt
	itte(is)=kt
	kv=kv+1
	tm(kt)=0.d0
	itb(kt)=kv
	ite(kt)=kv
	ivn(kv)=2
	ivp(kv)=0
	goto 30
15	ittb(is)=kt+1
	do 18 it=1,lt
	kt=kt+1
	tm(kt)=tmt(it)
	itb(kt)=kv+1
	do 16 i=itbt(it),itet(it)
	kv=kv+1
	ivn(kv)=ifnt(i)
16	ivp(kv)=ifpt(i)
18	ite(kt)=kv
	itte(is)=kt
30	continue
	nt=kt
	itetop=kv
	minns=ns+1
	return
	end
 
 
	subroutine print(ii)
$include: 'pfcom.for'
	common/pr1/pa(232),aa(22)
	character*1 aa,pa,bl,sp,sm,cf
	data bl,sp,sm,cf/' ','+','-','^'/
	if(np.eq.0)return
	ia=ittb(ii)
	if(ia.eq.0)write(8,2)
	if(ia.eq.0)return
2	format('there is nothing to print')
	ib=itte(ii)
	nst=21
	do 40 it=ia,ib
	n=nst
	do 4 i=1,232
4	pa(i)=' '
	ka=itb(it)
	kb=ite(it)
	do 30 k=ka,kb
	n=n+1
	pa(n)=bl
	if=ivn(k)
	jp=ivp(k)
	ifl=ipv(if)
	do 20 j=1,ifl
	n=n+1
20	pa(n)=pv(j,if)
	if(jp.eq.1)goto 30
	n=n+1
	if(jp.ge.0)pa(n)=cf
	if(jp.lt.0)pa(n)=sm
	iap=iabs(jp)
	n=n+1
	j=mod(iap,10)
	pa(n)=hn(j+1)
	if(iap.le.9)goto 30
	j=iap/10
	if(j.ge.10)call bomb(3,'print2  ')
	n=n+1
	pa(n)=pa(n-1)
	pa(n-1)=hn(j+1)
30	continue
	if(n.gt.232)call bomb(3,'print3  ')
	call fmtc(it,if,il,' ')
	if(il-if.eq.2 .and. aa(if+1).eq.'1' .and. aa(il).eq.'.')il=if
	do 32 i=if,il
32	pa(nst-il+i)=aa(i)
	j=it+1-ia
33	if(j.ge.1000)j=j-1000
	if(j.ge.1000)goto 33
	write(8,34)j,(pa(k),k=1,n)
34	format(i3,2x,127a1,/,27x,105a1)
40	continue
	return
	end
 
 
	subroutine prntap
$include: 'pfcom.for'
	do 1 n=nchpl,1,-1
	if(ap(n).ne.' ')goto 2
1	continue
	n=1
2	write(8,3)(ap(i),i=1,n)
3	format(131a1)
	return
	end
 
 
	subroutine read(ib)
$include: 'pfcom.for'
c  ib=0 on return if line has nonblanks
c  ib=1 on return if line is all blanks
	line=line+1
	ib=inpba(line)
	ie=inpea(line)
	il=ie-ib+1
	do 1 i=1,nchpl
1	a(i)=' '
	call copy(inpl(ib),a,il)
	call copy(a,ap,nchpl)
	i=1
	if(a(1).eq.'*')goto 8
	do 6 i=1,nchpl
	if(a(i).eq.';')goto 8
6	continue
	goto 14
8	do 10 j=i,nchpl
10	a(j)=' '
14	ib=0
	do 16 i=1,nchpl
	if(a(i).ne.' ')goto 18
16	continue
17	ib=1
	return
18	if(i.eq.1)return
	do 20 j=i,nchpl
	a(j+1-i)=a(j)
20	a(j)=' '
	return
	end
 
 
	subroutine readin
$include: 'pfcom.for'
	do 1 i=1,nlab
1	linlab(i)=-1
	line=0
	lines=0
	labels=0
	iinpl=1
5	read(7,6)a
6	format(127a1)
	do 7 k=nchpl,1,-1
	if(a(k).ne.' ')goto 10
7	continue
	k=1
10	lines=lines+1
	lintot=lintot+1
	if(lines.gt.ninl)call bomb(4,'ninl    ')
	if(k.eq.1)goto 15
	call srchnb(1,k,m)
	if(a(m).eq.'p' .and. a(m+1).eq.'r' .and. a(m+2).eq.'o')goto 11
	if(a(m).ne.'l' .or. a(m+1).ne.'a' .or. a(m+2).ne.'b')goto 15
11	call nxt(m+2,nchpl,' ',j1)
	call nxtch(j1,nchpl,j3)
	call nxt(j3,nchpl,' ',j4)
	call numl(j1,nchpl,j3,j4,j5)
	if(linlab(j5).ne.-1)write(8,12)(a(i),i=j3,j4)
12	format('duplicate label: ',15a1)
	if(linlab(j5).ne.-1)call bomb(3,'readin1 ')
	linlab(labels)=lines
15	inpba(lines)=iinpl
	inpea(lines)=iinpl+k-1
	if(inpea(lines).gt.ninch)call bomb(4,'ninch   ')
	call copy(a,inpl(iinpl),k)
	iinpl=iinpl+k
	if(k.eq.1)goto 5
	if(a(m).eq.'f' .and. a(m+1).eq.'l' .and. a(m+2).eq.'u')goto 16
	if(a(m).ne.'e' .or. a(m+1).ne.'n' .or. a(m+2).ne.'d')goto 5
	if(a(m+3).ne.' ' .or. a(m+4).ne.'i' .or. a(m+5).ne.'n')goto 5
16	write(8,20)lines,labels,iinpl
20	format('number of lines in input file=',i5,' labels='
     +      ,i2,'  characters=',i6)
	return
	end
 
 
	subroutine recip(j13,j4)
$include: 'pfcom.for'
	ia=ittb(j13)
	iz=itte(j13)
	k=itb(ia)
	if(ivn(k).eq.2 .and. ivp(k).ne.0 .and. ia.ne.iz)then
		write(8,2)
		call bomb(3,'recip1  ')
2	format('this simple program cannot process the denominator.',
     +     '  denominators should have 1 term without variable 1.')
	endif
	ja=nt+1
	call cpynt(ia,ia,1.d0)
	tm(ja)=1.d0/tm(ja)
	do 4 i=itb(ja),ite(ja)
4	ivp(i)=-ivp(i)
	if(iz.eq.ia)call defsum(j4,ja,nt)
	if(iz.eq.ia)return
	if(maxe.gt.2)write(8,6)
	if(maxe.gt.2)call bomb(3,'recip2  ')
6	format('this simple program cannot divide if maxe is',
     +      ' greater than 2')
	ib=ia+1
	k=itb(ib)
	if(ivn(k).ne.2 .or. ivp(k).lt.1)write(8,2)
	if(ivn(k).ne.2 .or. ivp(k).lt.1)call bomb(3,'recip3  ')
	if(ivp(k).eq.2)goto 14
	ic=ib
	if(ib.eq.iz)goto 10
8	k=itb(ic+1)
	if(ic.eq.iz .or. ivp(k).eq.2)goto 10
	ic=ic+1
	goto 8
10	jb=nt+1
	call cpynt(ib,ic,1.d0)
	jc=nt
	if(ic.eq.iz)goto 12
11	jd=nt+1
	call cpynt(ic+1,iz,1.d0)
	je=nt
	goto 16
12	jd=nt+1
	je=jd
	itb(jd)=itetop+1
	call incnt
	call nite(itb(jd))
	k=itetop
	ivn(k)=2
	ivp(k)=2
	tm(jd)=0.d0
	goto 16
14	jb=nt+1
	jc=jb
	itb(jb)=itetop+1
	call incnt
	call nite(itb(jb))
	k=itetop
	ivn(k)=2
	ivp(k)=1
	tm(jb)=0.d0
	ic=ib-1
	goto 11
16	jf=nt+1
	call mult(ja,ja,jb,jc)
	jg=nt
	if(maxe.eq.1)goto 18
	jh=nt+1
	call mult(jf,jg,jf,jg)
	ji=nt
	jj=nt+1
	call mult(ja,ja,jd,je)
	jk=nt
	jl=nt+1
	call cpynt(jf,jg,-1.d0)
	call cpynt(jh,ji,1.d0)
	call cpynt(jj,jk,-1.d0)
	jm=nt
	call cpynt(ja,ja,1.d0)
	call mult(ja,ja,jl,jm)
	nt1=jm+1
	goto 24
18	do 19 i=jf,jg
19	tm(i)=-tm(i)
	call cpynt(ja,ja,1.d0)
	call mult(ja,ja,jf,jg)
	nt1=jg+1
24	call defsum(j4,nt1,nt)
	call order(j4)
	return
	end
 
 
	subroutine swap(j4)
$include: 'pfcom.for'
	if=ittb(j4)
	ig=ios(j4)
	if(if.ne.0)minns=min0(minns,ios(j4))
	ios(j4)=ios(1)
	iosi(ios(j4))=j4
	ittb(j4)=ittb(1)
	itte(j4)=itte(1)
	if(if.eq.0)goto 5
	ios(1)=ig
	iosi(ig)=1
	ittb(1)=if
	itte(1)=if
	ite(if)=itb(if)
	call pack
	return
5	call niltrm
	ittb(1)=0
	call defsum(1,nt,nt)
	return
	end
 
 
	subroutine srch(j1,j2,h,j3)
$include: 'pfcom.for'
	character*1 h
	do 1 j3=j1,j2
	if(a(j3).eq.h)return
1	continue
	j3=-1
	return
	end
 
 
	subroutine srchch(j1,j2,j3)
$include: 'pfcom.for'
	do 2 j3=j1,j2
	do 1 i=1,26
	if(a(j3).eq.ha(i))return
	if(a(j3).eq.hc(i))return
1	continue
2	continue
	j3=-1
	end
 
 
	subroutine srchnb(j1,j2,j3)
$include: 'pfcom.for'
	do 1 j3=j1,j2
	if(a(j3).ne.' ')return
1	continue
	j3=-1
	return
	end
 
 
	subroutine sub
$include: 'pfcom.for'
	character*3 aj4,aj8
	call nxt(1,nchpl,' ',j1)
	call nxtnb(j1,nchpl,j2)
	call nxtch(j1,nchpl,j3)
	if(j2.ne.j3)goto 10
	call numv(j1,nchpl,j2,j3,j4)
	aj4='sum'
	if(ittb(j4).eq.0)aj4='var'
1	call nxt(j3+1,nchpl,'f',j5)
	if(a(j5+1).ne.'o')call bomb(3,'readin2 ')
	if(a(j5+2).ne.'r')call bomb(3,'readin4 ')
	if(a(j5+3).ne.' ')call bomb(3,'readin6 ')
	call numv(j5+3,nchpl,j6,j7,j8)
	aj8='var'
	if(ittb(j8).gt.0)aj8='sum'
	call nxt(j7+1,nchpl,'i',j9)
	if(a(j9+1).ne.'n')call bomb(3,'readin8 ')
	if(a(j9+2).ne.' ')call bomb(3,'readin10')
	call numv(j9+2,nchpl,j10,j11,j12)
	minns=min0(minns,ios(j12))
	if(aj4.eq.'num')goto 35
	if(aj4.eq.'var')goto 30
	if(aj8.eq.'sum')goto 50
	call subs(j4,j8,j12)
	return
10	call coef(j1,j4,c)
	aj4='num'
	j3=j4-1
	goto 1
30	if(aj8.eq.'sum')goto 40
	call subv(j4,j8,j12)
	return
35	if(aj8.eq.'sum')call bomb(3,'readin12')
	call subn(c,j8,j12)
	return
40	call subr(j8,j4,j12)
	return
50	a(71)='A'
	a(72)='?'
	a(73)='0'
	call numv(71,73,j15,j16,j17)
	call subr(j8,j17,j12)
	call subs(j4,j17,j12)
	return
	end
 
 
	subroutine subn(c,j8,j12)
$include: 'pfcom.for'
	if(ittb(j12).eq.0)call bomb(3,'subn1   ')
	if(ittb(j8).ne.0)call bomb(3,'subn2   ')
	do 8 it=ittb(j12),itte(j12)
	do 5 i=itb(it),ite(it)
	if(ivn(i).ne.j8)goto 5
	tm(it)=tm(it)*c**ivp(i)
	ivp(i)=0
	if(ite(it).eq.itb(it))goto 5
	ite(it)=ite(it)-1
	if(i.gt.ite(it))goto 5
	do 1 j=i,ite(it)
	ivn(j)=ivn(j+1)
1	ivp(j)=ivp(j+1)
5	continue
8	continue
	call order(j12)
	return
	end
 
 
	subroutine subr(ia,ib,ic)
c  replace pattern=first term in sum ia by variable ib in sum ic
$include: 'pfcom.for'
	if(ittb(ia).eq.0)call bomb(3,'subr1   ')
	if(ittb(ib).ne.0)call bomb(3,'subr2   ')
	nt1=nt+1
	j1=ittb(ia)
	j2=itb(j1)
	j3=ite(j1)
	j4=ittb(ic)
	j5=itte(ic)
	do 30 it=j4,j5
	j6=itb(it)
	j7=ite(it)
	minrat=9999
	do 5 j8=j2,j3
	jfp=ivn(j8)
1	jfm=ivn(j6)
	if(jfm.gt.jfp)goto 22
	irat=ivp(j6)/ivp(j8)
	if(jfm.eq.jfp .and. irat.ge.1)goto 4
	jt=it+1-j4
	if(j2.eq.j3 .and. np.ge.2 .and. jfm.eq.jfp .and. irat.le.-1)
     +  write(8,2)jt,(pv(i,ic),i=1,nch),(pv(i,jfp),i=1,nch),
     +            (pv(i,ia),i=1,nch)
2	format('term',i4,' of ',15a1,' has exponent of ',15a1,
     +       ' with different sign from ',15a1)
	j6=j6+1
	if(j6.le.j7)goto 1
	goto 22
4	minrat=min0(minrat,ivp(j6)/ivp(j8))
5	continue
	call incnt
	tm(nt)=tm(it)/tm(j1)
	itb(nt)=itetop+1
	j6=itb(it)
	k=itb(nt)
	do 20 j=j6,j7
	do 12 j8=j2,j3
	if(ivn(j8).eq.ivn(j))goto 16
12	continue
	ivn(k)=ivn(j)
	ivp(k)=ivp(j)
	goto 20
16	ivn(k)=ivn(j)
	ivp(k)=ivp(j)-minrat*ivp(j8)
20	k=k+1
	ivn(k)=ib
	ivp(k)=minrat
	call nite(k)
	goto 30
22	call cpynt(it,it,1.d0)
30	continue
	call defsum(ic,nt1,nt)
	call order(ic)
	return
	end
 
 
	subroutine subv(j4,j8,j12)
$include: 'pfcom.for'
	if(ittb(j12).eq.0)call bomb(3,'subv1   ')
	if(ittb(j4).ne.0 .or. ittb(j8).ne.0)call bomb(3,'subv2   ')
	do 8 it=ittb(j12),itte(j12)
	do 5 i=itb(it),ite(it)
	if(ivn(i).eq.j8)ivn(i)=j4
5	continue
8	continue
	call order(j12)
	return
	end
 
 
	subroutine subs(in1,in2,in3)
c  substitute sum=in1 for variable=in2 in sum=in3
c  result goes in nt+1,...
$include: 'pfcom.for'
	common/nwt/iitb,iite,jitb,jite,ntb,if
c  dimension of ittbs ittes is max jp
	dimension ittbs(30),ittes(30)
	jp=1
	ia=in2
	ib=ittb(in3)
	ic=itte(in3)
	id=ittb(in1)
	ie=itte(in1)
	if(ib.eq.0 .or. id.eq.0)call bomb(3,'subs1   ')
	ntb=nt+1
	lev1=0
	do 30 it=ib,ic
	iitb=itb(it)
	iite=ite(it)
	if(ivn(iitb).eq.2)lev1=ivp(iitb)
	if(lev1.gt.maxe)goto 31
	do 5 if=iitb,iite
	if(ivn(if).ne.ia)goto 5
	ifp=ivp(if)
	lev2=0
	if(ifp.eq.1)goto 10
	if(id.eq.ie)goto 70
	if(ifp.gt.0)goto 17
	itpr=it-ib+1
	if(np.ge.2)write(8,1)(pv(i,ia),i=1,nch),
     +       (pv(i,in3),i=1,nch),itpr
1	format(2x,15a1,' occurs to a negative power in ',15a1,
     +      '  term ',i4)
5	continue
c  here copy term it to nt+1
	if=0
	jitb=0
	call sub2(tm(it))
	goto 30
c  here copy term it with substitution
10	do 68 jt=id,ie
	jitb=itb(jt)
	jite=ite(jt)
	if(ivn(jitb).eq.2)lev2=ivp(jitb)
	if(lev1+lev2.gt.maxe)goto 30
	call sub2(tm(it)*tm(jt))
68	continue
	goto 30
70	if(ivn(itb(id)).eq.2)lev2=ivp(itb(id))*ifp
	if(lev1+lev2.gt.maxe)goto 30
	k=0
	do 71 i=itb(id),ite(id)
	k=k+1
	ifnt(k)=ivn(i)
71	ifpt(k)=ivp(i)*ifp
	jitb=1
	jite=k
	call sub3(tm(it)*tm(id)**ifp)
	goto 30
c  here get in1**ifp
17	jd=id
	je=ie
	ifpsav=ifp
	ntbsav=ntb
	iitbsv=iitb
	iitesv=iite
	ifsav=if
	ntbmu=nt
	itesav=itetop
c  test whether this power is in memory
	if(jp.eq.1)goto 80
	if(ifp.lt.jp)goto 27
	if(ittbs(jp).eq.0)goto 30
	if(ifp.eq.jp)goto 27
	jd=nt+1
	k=itetop
	do 86 i=ittbs(jp),ittes(jp)
	call incnt
	tm(nt)=tmt(i)
	itb(nt)=k+1
	do 84 j=itbt(i),itet(i)
	k=k+1
	ivn(k)=ifnt(j)
84	ivp(k)=ifpt(j)
	call nite(k)
86	continue
	je=nt
	goto 18
80	ks=0
	nts=0
18	jpb=nt+1
	ntbtmu=nt
	call mult(id,ie,jd,je)
	if(nt.eq.ntbtmu)then
		jp=jp+1
		ittbs(jp)=0
		nt=ntbmu
		ntb=ntbsav
		itetop=itesav
		goto 30
		endif
	jpe=nt
	jp=jp+1
	if(jp.gt.30)call bomb(4,'max jp  ')
	jd=jpb
	je=jpe
	ittbs(jp)=nts+1
	do 21 i=jpb,jpe
	nts=nts+1
	tmt(nts)=tm(i)
	itbt(nts)=ks+1
	do 20 j=itb(i),ite(i)
	ks=ks+1
	ifnt(ks)=ivn(j)
20	ifpt(ks)=ivp(j)
	itet(nts)=ks
	if(ks.gt.ntott)call bomb(4,'ntott   ')
	if(nts.gt.ntermt)call bomb(4,'ntermt  ')
21	continue
	ittes(jp)=nts
	if(jp.ne.ifp)goto 18
	nt=ntbmu
	itetop=itesav
	if=ifsav
	ifp=ifpsav
	iitb=iitbsv
	iite=iitesv
	ntb=ntbsav
27	do 54 jt=ittbs(ifp),ittes(ifp)
	jitb=itbt(jt)
	jite=itet(jt)
	if(ifnt(jitb).eq.2)lev2=ifpt(jitb)
	if(lev1+lev2.gt.maxe)goto 30
	call sub3(tm(it)*tmt(jt))
54	continue
30	continue
31	call defsum(in3,ntb,nt)
	return
	end
 
 
	subroutine sub2(coef)
$include: 'pfcom.for'
	common/nwt/iitb,iite,jitb,jite,ntb,if
	if(dabs(coef).lt.1.d-11)return
	call incnt
	k=itetop
	itb(nt)=k+1
	tm(nt)=coef
	ii=iitb
	if(jitb.eq.0)goto 15
	jj=jitb
11	if(ii.gt.iite)goto 13
	if(ii.eq.if)goto 12
	if(jj.gt.jite)goto 15
	if(ivn(ii).lt.ivn(jj))goto 16
	if(ivn(ii).eq.ivn(jj))goto 57
	if(ivp(jj).eq.0)goto 18
	k=k+1
	ivn(k)=ivn(jj)
	ivp(k)=ivp(jj)
18	jj=jj+1
	goto 11
12	ii=ii+1
	goto 11
13	if(jj.gt.jite)goto 60
	if(ivp(jj).eq.0)goto 17
	k=k+1
	ivn(k)=ivn(jj)
	ivp(k)=ivp(jj)
17	jj=jj+1
	goto 13
14	if(ii.gt.iite)goto 60
15	if(ivp(ii).eq.0)goto 9
	k=k+1
	ivn(k)=ivn(ii)
	ivp(k)=ivp(ii)
9	ii=ii+1
	if(ii.eq.if)goto 9
	goto 14
16	if(ivp(ii).eq.0)goto 12
	k=k+1
	ivn(k)=ivn(ii)
	ivp(k)=ivp(ii)
	ii=ii+1
	goto 11
57	if(ivp(ii)+ivp(jj).eq.0)goto 58
	k=k+1
	ivn(k)=ivn(ii)
	ivp(k)=ivp(ii)+ivp(jj)
58	ii=ii+1
	jj=jj+1
	goto 11
60	if(k.gt.itetop)goto 56
	k=itetop+1
	ivn(k)=2
	ivp(k)=0
56	itesav=itetop
	call nite(k)
	if(nt.eq.ntb)goto 68
	if(icomp(ntb,nt))31,33,34
31	call tzap(nt,ntb)
	goto 68
33	tm(ntb)=tm(ntb)+tm(nt)
	nt=nt-1
	itetop=itesav
	goto 68
34	if(icomp(nt-1,nt))62,61,68
61	tm(nt-1)=tm(nt-1)+tm(nt)
	nt=nt-1
	itetop=itesav
	goto 68
62	la=ntb
	lb=nt-1
63	if(lb-la.eq.1)goto 64
	lc=(la+lb)/2
	if(icomp(lc,nt))65,66,67
64	call tzap(nt,lb)
	goto 68
65	lb=lc
	goto 63
66	tm(lc)=tm(lc)+tm(nt)
	nt=nt-1
	itetop=itesav
	goto 68
67	la=lc
	goto 63
68	return
	end
 
 
	subroutine sub3(coef)
$include: 'pfcom.for'
	common/nwt/iitb,iite,jitb,jite,ntb,if
	if(dabs(coef).lt.1.d-11)return
	call incnt
	k=itetop
	itb(nt)=k+1
	tm(nt)=coef
	ii=iitb
	if(jitb.eq.0)goto 15
	jj=jitb
11	if(ii.gt.iite)goto 13
	if(ii.eq.if)goto 12
	if(jj.gt.jite)goto 15
	if(ivn(ii).lt.ifnt(jj))goto 16
	if(ivn(ii).eq.ifnt(jj))goto 57
	if(ifpt(jj).eq.0)goto 18
	k=k+1
	ivn(k)=ifnt(jj)
	ivp(k)=ifpt(jj)
18	jj=jj+1
	goto 11
12	ii=ii+1
	goto 11
13	if(jj.gt.jite)goto 60
	if(ifpt(jj).eq.0)goto 17
	k=k+1
	ivn(k)=ifnt(jj)
	ivp(k)=ifpt(jj)
17	jj=jj+1
	goto 13
14	if(ii.gt.iite)goto 60
15	if(ivp(ii).eq.0)goto 9
	k=k+1
	ivn(k)=ivn(ii)
	ivp(k)=ivp(ii)
9	ii=ii+1
	if(ii.eq.if)goto 9
	goto 14
16	if(ivp(ii).eq.0)goto 12
	k=k+1
	ivn(k)=ivn(ii)
	ivp(k)=ivp(ii)
	ii=ii+1
	goto 11
57	if(ivp(ii)+ifpt(jj).eq.0)goto 58
	k=k+1
	ivn(k)=ivn(ii)
	ivp(k)=ivp(ii)+ifpt(jj)
58	ii=ii+1
	jj=jj+1
	goto 11
60	if(k.gt.itetop)goto 56
	k=itetop+1
	ivn(k)=2
	ivp(k)=0
56	itesav=itetop
	call nite(k)
	if(nt.eq.ntb)goto 68
	if(icomp(ntb,nt))31,33,34
31	call tzap(nt,ntb)
	goto 68
33	tm(ntb)=tm(ntb)+tm(nt)
	nt=nt-1
	itetop=itesav
	goto 68
34	if(icomp(nt-1,nt))62,61,68
61	tm(nt-1)=tm(nt-1)+tm(nt)
	nt=nt-1
	itetop=itesav
	goto 68
62	la=ntb
	lb=nt-1
63	if(lb-la.eq.1)goto 64
	lc=(la+lb)/2
	if(icomp(lc,nt))65,66,67
64	call tzap(nt,lb)
	goto 68
65	lb=lc
	goto 63
66	tm(lc)=tm(lc)+tm(nt)
	nt=nt-1
	itetop=itesav
	goto 68
67	la=lc
	goto 63
68	return
	end
 
 
	subroutine switch(j,k,is)
$include: 'pfcom.for'
	do 2 it=ittb(is),itte(is)
	do 1 i=itb(it),ite(it)
	if(ivn(i).eq.j)ivn(i)=k
1	continue
2	continue
	return
	end
 
 
	subroutine taylor
$include: 'pfcom.for'
	character*4 star
	dimension ltay(mitay)
	data iorder,ifirst/1,0/
	if(idflg.ne.0)call dertab
	call nxt(2,nchpl,' ',j1)
	call srch(j1,nchpl,'=',j2)
	j4=j1
	if(j2.eq.-1)goto 5
	call getexp(j2+1,iorder,j3,0)
	if(iorder.lt.1)call bomb(3,' taylor1')
	return
5	call nxt(2,nchpl,' ',j13)
	star='no'
	if(a(j13-1).eq.'*')star='yes'
10	call nxt(j4,nchpl,'i',j3)
	if(j3.eq.-1)call bomb(3,' taylor2')
	if(a(j3-1).eq.' ' .and. a(j3+1).eq.'n' .and. a(j3+2).eq.' ')
     +       goto 12
	j4=j3+1
	if(j4.lt.nchpl)goto 10
	call bomb(3,' taylor3')
12	do 14 j5=j3-1,j1,-1
	if(a(j5).ne.' ')goto 16
14	continue
	call bomb(3,' taylor4')
16	do 18 j6=j5,j1,-1
	if(a(j6).eq.' ')goto 20
18	continue
	call bomb(3,' taylor5')
20	call numv(j6,j5,j7,j8,j9)
	if(ittb(j9).eq.0)call bomb(3,' taylor6')
	itesav=itetop
	call term(j3+2)
	itay=itetop-itesav
	if(itay.gt.mitay)call bomb(4,'mitay   ')
	do 21 i=1,itay
21	ltay(i)=ivn(itesav+i)
	nt=nt-1
	itetop=itesav
	call niltrm
	call defsum(1,nt,nt)
	if(ifirst.ne.0)goto 22
	a(1)='A'
	a(2)='?'
	a(3)='0'
	call numv(1,3,j14,j15,j10)
	a(3)='1'
	call numv(1,3,j14,j15,j11)
	a(3)='2'
	call numv(1,3,j14,j15,j12)
22	ifirst=1
	nv1=nv+1
	do 34 it=ittb(j9),itte(j9)
	call cpynt(it,it,1.d0)
	call defsum(j10,nt,nt)
	do 32 ivv=1,itay
	iv=ltay(ivv)
	if(star.eq.'yes')call switch(iv,nv1,j10)
	call copys(j10,j11)
	call subn(0.d0,iv,j10)
	f=1.d0
	do 30 io=1,iorder
	call deriv(j11,iv,j12)
	call outqu(j12)
	if(tm(ittb(j12)).eq.0.d0)goto 31
	call copys(j12,j11)
	call subn(0.d0,iv,j12)
	f=f*io
	call incnt
	k=itetop+1
	tm(nt)=1.d0/f
	itb(nt)=k
	ite(nt)=k
	ivn(k)=iv
	ivp(k)=io
	call nite(k)
	nt1=nt+1
	call mult(nt,nt,ittb(j12),itte(j12))
	call defsum(j12,nt1,nt)
	call pack
30	call add(j12,j10,1.d0,j10)
31	if(star.eq.'yes')call switch(nv1,iv,j10)
32	continue
34	call add(j10,1,1.d0,1)
	call swap(j9)
	call order(j9)
	call delsum(j10)
	return
	end
 
 
	subroutine term(ist)
c  each term is on a separate line, free format with coefficient first
$include: 'pfcom.for'
	call incnt
	call coef(ist,j1,c)
	tm(nt)=c
	itb(nt)=itetop+1
	if(j1.eq.-1)then
		k=itb(nt)
		ivn(k)=2
		ivp(k)=0
		goto 10
		endif
	k=itetop
5	ise=1
	if(a(j1).eq.'*')j1=j1+1
	if(a(j1).eq.'/')ise=-1
	if(a(j1).eq.'/')j1=j1+1
	call numv(j1,nchpl,j2,j3,j4)
	ie=1
	j5=-1
	if(j3.eq.nchpl-1 .and. a(nchpl).ne.' ')call bomb(3,'term1   ')
	if(j3.le.nchpl-2)call getexp(j3+1,ie,j5,1)
	k=k+1
	ivn(k)=j4
	ivp(k)=ie*ise
	j1=j5
	if(j5.ne.-1)goto 5
10	call nite(k)
	return
	end
 
 
	subroutine thderv
$include: 'pfcom.for'
	idflg=1
	call nxt(2,nchpl,'o',j1)
	call nxt(j1,nchpl,' ',j2)
	call numv(j2,nchpl,j3,j4,j5)
	call nxt(j4+1,nchpl,'w',j6)
	call nxt(j6,nchpl,' ',j7)
	call numv(j7,nchpl,j8,j9,j10)
	call srch(j9+1,nchpl,'=',j11)
	if(j11.ne.-1)goto 1
	call srch(j9+1,nchpl,'i',j12)
	if(j12.eq.-1)call bomb(3,'thderv1 ')
	call nxt(j12,nchpl,' ',j11)
1	call srchnb(j11+1,nchpl,j12)
	if(j12.eq.-1)call bomb(3,'thderv2 ')
	ig=0
	if(a(j12).eq.'0')ig=1
	if(ig.eq.1)goto 8
	call numv(j11+1,nchpl,j12,j13,j14)
8	if=0
	do 10 i=1,nd
	if(ider(1,i).eq.j5 .and. ider(2,i).eq.j10)if=1
	if(ider(1,i).eq.j5 .and. ider(2,i).eq.j10)in=i
10	continue
	if(if.eq.0 .and. ig.eq.1)goto 50
	if(if.eq.0)goto 14
	ja=ider(3,in)
	jb=ipv(ja)
	if(np.ge.2)write(8,11)(pv(i,ja),i=1,jb)
11	format('REDEFINITION: previously the  derivative was = ',15a1)
	if(ig.eq.0)goto 18
	nd=nd-1
	if(in.gt.nd)goto 50
	do 28 ih=in,nd
	ider(1,ih)=ider(1,ih+1)
	ider(2,ih)=ider(2,ih+1)
28	ider(3,ih)=ider(3,ih+1)
	goto 50
14	nd=nd+1
	in=nd
	if(nd.gt.nder)call bomb(4,'nder    ')
	ider(1,in)=j5
	ider(2,in)=j10
18	ider(3,in)=j14
50	return
	end
 
 
	subroutine tzap(ka,kb)
c  zap term at ka to position kb
$include: 'pfcom.for'
	if(ka.eq.kb)return
	tmka=tm(ka)
	itbka=itb(ka)
	iteka=ite(ka)
	if(ka.gt.kb)goto 5
	do 1 i=ka,kb-1
	tm(i)=tm(i+1)
	itb(i)=itb(i+1)
1	ite(i)=ite(i+1)
	goto 8
5	do 6 i=ka,kb+1,-1
	tm(i)=tm(i-1)
	itb(i)=itb(i-1)
6	ite(i)=ite(i-1)
8	tm(kb)=tmka
	itb(kb)=itbka
	ite(kb)=iteka
	return
	end
 
 
 
	subroutine zsub(j3)
	return
	end