$storage:2
$include: 'picom.for'
	open(unit=7,file='algin',status='old')
	open(unit=8,file='algout',status='new')
	call datim
	itbeg=itime
	do 4 i=1,nvar
4	isb(i)=0
	ios(1)=1
	iosi(1)=1
	ipv(1)=1
	pv(1,1)='?'
	isb(1)=1
	ise(1)=1
	itb(1)=1
	ite(1)=1
	ivn(1)=2
	ivp(1)=1
	nd=0
	nv=1
	ns=1
	np=2
	nt=1
	nb=7
	maxe=9999
	maxnt=0
	maxite=0
	itetop=1
	ndump=0
	minns=9999
	lintot=0
	idflg=1
	call defh
	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/20 ',
     +   ' NTOT=',i6,' NVAR=',i4,' NTERM=',i5,' NINCH=',i5
     +   ,/,' NINL=',i6,' NDER=',i3,' NLAB=',i3,' NDIND=',i2
     +   ,' MITAY=',i2)
	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=(itime-itbeg)/100.
	write(8,91)t
91	format(' time used=',f9.2,' seconds')
	stop
	end
 
 
	subroutine addrat(n1,i1,n2,i2,n3,i3)
c  n3/i3=n1/i1 + n2/i2
	implicit integer*4 (i-n)
	i=i1*i2
	n=n1*i2+n2*i1
	call gcdout(n,i)
	n3=n
	i3=i
	return
	end
 
 
	subroutine adds(is,js,iin,iid,ians)
$include: 'picom.for'
	integer*4 in1,id1,in,id,inn,idd,insav,idsav
	in=iin
	id=iid
	nt1=nt+1
	ib=isb(is)
	ie=ise(is)
	jb=isb(js)
	je=ise(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 copyt(jb,je,iin,iid)
12	call copyt(ib,ie,1,1)
14	call defsum(ians,nt1,nt)
	return
16	insav=itn(je)
	idsav=itd(je)
	idd=id
	if(in.lt.0)idd=-idd
	inn=iabs(in)
	call mulrat(itn(ib),itd(ib),idd,inn,in1,id1)
	call addrat(in1,id1,itn(je),itd(je),itn(je),itd(je))
	call copyt(jb,je,iin,iid)
	itn(je)=insav
	itd(je)=idsav
	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	insav=itn(jc)
	idsav=itd(jc)
	idd=id
	if(in.lt.0)idd=-idd
	inn=iabs(in)
	call mulrat(itn(ib),itd(ib),idd,inn,in1,id1)
	call addrat(in1,id1,itn(jc),itd(jc),itn(jc),itd(jc))
	call copyt(jb,jc,iin,iid)
	itn(jc)=insav
	itd(jc)=idsav
	ib=ib+1
	jb=jc+1
	goto 7
28	call copyt(jb,jl,iin,iid)
	jb=jr
30	il=ib
	if(ib.eq.ie)goto 31
	if(icomp(jb,ie))31,36,38
31	call copyt(ib,ie,1,1)
32	call copyt(jb,je,iin,iid)
	goto 14
36	insav=itn(ie)
	idsav=itd(ie)
	call mulrat(in,id,itn(jb),itd(jb),in1,id1)
	call addrat(in1,id1,itn(ie),itd(ie),itn(ie),itd(ie))
	call copyt(ib,ie,1,1)
	itn(ie)=insav
	itd(ie)=idsav
	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	insav=itn(ic)
	idsav=itd(ic)
	call mulrat(in,id,itn(jb),itd(jb),in1,id1)
	call addrat(in1,id1,itn(ic),itd(ic),itn(ic),itd(ic))
	call copyt(ib,ic,1,1)
	itn(ic)=insav
	itd(ic)=idsav
	jb=jb+1
	ib=ic+1
	goto 7
48	call copyt(ib,il,1,1)
	ib=ir
	goto 8
60	insav=itn(ib)
	idsav=itd(ib)
	call mulrat(in,id,itn(jb),itd(jb),in1,id1)
	call addrat(in1,id1,itn(ib),itd(ib),itn(ib),itd(ib))
	call copyt(ib,ib,1,1)
	itn(ib)=insav
	itd(ib)=idsav
	ib=ib+1
	jb=jb+1
	goto 7
70	if(jb.gt.je)goto 14
	goto 32
	end
 
 
	subroutine bomb(i,h)
$include: 'picom.for'
	character*8 h
	goto(100,20,30,40,50),i
20	write(8,21)
21	format(' parameter error')
	goto 100
30	write(8,31)h
31	format(1x,a8)
	goto 100
40	write(8,41)h,h
41	format(1x,a8,' is too small for your problem.',
     +       '  Recompile with bigger ',a8)
50	continue
100	write(8,101)line,ap
101	format(' BOMB: input line nu.',i5,' was:',/,1x,127a1)
	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)
	call datim
	if(ndump.ne.0)call dump
	t=(itime-itbeg)/100.
	write(8,91)t
91	format(' time used=',f9.2,' seconds')
	stop
	end
 
 
	subroutine bkup(j,k)
$include: 'picom.for'
	if(j.eq.k)goto 2
	do 1 i=j,k-1
	itb(i)=itb(i+1)
	ite(i)=ite(i+1)
	itn(i)=itn(i+1)
1	itd(i)=itd(i+1)
2	k=k-1
	return
	end
 
 
	subroutine call
$include: 'picom.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,jv,in,id)
c is=start, on return jv points beyond coef
$include: 'picom.for'
	integer*4 in,id,ii,ij
	id=1
	isgn=1
	call srchnb(is,nchpl,j1)
	if(j1.ne.-1)goto 5
1	jv=-1
2	in=isgn
	return
5	isgn=-1
	if(a(j1).eq.'-')goto 30
	isgn=1
	if(a(j1).eq.'+')goto 30
	j2=j1
7	jv=j2
	call intval(j2,i)
	if(i.eq.-1)goto 2
	call nxtnnu(j2,nchpl,j3)
	call getint(j2,j3-1,in2,in)
	j5=j3
	if(a(j3).eq.'.')goto 10
	if(a(j3).eq.'/')goto 25
	if(a(j3).eq.' ')goto 20
	goto 15
10	call intval(j3+1,i)
	if(i.ne.-1)goto 16
	j5=j3+1
	if(a(j5).eq.'/')goto 25
	if(a(j5).eq.' ')goto 20
15	jv=j5
	in=in*isgn
	return
16	call nxtnnu(j3+1,nchpl,j7)
	call getint(j3+1,j7-1,in2,ii)
	ij=10**(j7-j3-1)
	in=ij*in+ii
	id=ij
	call gcdout(in,id)
	goto 26
20	call srchnb(j5,nchpl,j1)
	j5=j1
	if(j5.eq.-1)goto 15
	if(a(j5).eq.'/')goto 25
	goto 15
25	call nxtnb(j5+1,nchpl,j6)
	call intval(j6,i)
	if(i.eq.-1)goto 15
	call nxtnnu(j6,nchpl,j7)
	call getint(j6,j7-1,in2,id)
26	if(a(j7).eq.'.')j7=j7+1
	call srchnb(j7,nchpl,jv)
	in=in*isgn
	return
30	call srchnb(j1+1,nchpl,j2)
	if(j2.eq.-1)goto 1
	goto 7
	end
 
 
	subroutine comput
$include: 'picom.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(isb(j10).eq.0 .or. isb(j13).eq.0)call bomb(3,'not sum ')
	call recip(j13,1)
	nt1=nt+1
	call mult(isb(j10),ise(j10),isb(1),ise(1))
	call defsum(j4,nt1,nt)
	call outqu(j4)
	minns=min0(minns,ios(1))
	it=isb(1)
	ise(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(isb(j18).eq.0)call bomb(3,'comput14')
	iord=itn(isb(j18))
	if(itd(isb(j18)).ne.1)call bomb(3,'comput15')
	if(iord.gt.0)goto 311
	call bomb(3,'comput16')
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(itn(isb(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 ')
	id12=1
	if(j6.ne.-1)goto 321
	call srch(j5,nchpl,'-',j6)
	in12=-1
	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(isb(j12).eq.0)call bomb(3,'comput4 ')
	call copys(j12,j4)
	call outqu(j4)
	return
c  compute sj4=sj9  + or -  sj12
321	in12=1
322	call numv(j5+1,j6-1,j7,j8,j9)
	call numv(j6+1,nchpl,j10,j11,j12)
	call adds(j9,j12,in12,id12,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(isb(j9),ise(j9),isb(j12),ise(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 copyc(a,b,n)
	integer*2 i,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: 'picom.for'
	nt1=nt+1
	call copyt(isb(j1),ise(j1),1,1)
	call defsum(j2,nt1,nt)
	return
	end
 
 
	subroutine count
$include: 'picom.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(isb(j8).eq.0)call bomb(3,'count2  ')
	if(isb(j5).ne.0)goto 20
	k=itetop+1
	call incnt
	itb(nt)=k
	call nite(k)
	itn(nt)=ise(j8)+1-isb(j8)
	itd(nt)=1
	ivn(k)=2
	ivp(k)=0
	call defsum(j5,nt,nt)
	return
20	ia=isb(j5)
	ib=ise(j5)
	k=itb(ia)
	l=ite(ia)
	ivn(k)=2
	ivp(k)=0
	ise(j5)=ia
	ite(ia)=k
	itn(ia)=ise(j8)+1-isb(j8)
	itd(ia)=1
	if(ia.ne.ib .or. k.ne.l)minns=min0(minns,ios(j5))
	call pack
	return
	end
 
 
	subroutine copyt(ia,ib,in,id)
c  copy terms ia-ib to nt+1,..., coefficients are multiplied by in/id
$include: 'picom.for'
	integer*4 in4,id4
	in4=in
	id4=id
	if(ib.lt.ia)return
	if(ia.lt.1 .or. ib.gt.nt)call bomb(3,'copyt2  ')
	icop=0
	if(id.eq.1 .and. in.eq.-1)icop=-1
	if(id.eq.1 .and. in.eq. 1)icop= 1
	do 15 i=ia,ib
	ja=itb(i)
	if(ivn(ja).eq.2 .and. ivp(ja).gt.maxe)goto 15
	k=itetop
	call incnt
	itn(nt)=itn(i)*icop
	itd(nt)=itd(i)
	if(icop.eq.0)call
     +     mulrat(itn(i),itd(i),in4,id4,itn(nt),itd(nt))
	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: 'picom.for'
	integer*4 mult
	call date(iy,mon,id)
	call time(ih,min,is,if)
	mult=60
	itime=100*(is+mult*(min+mult*ih))+if
	write(8,1)iy,mon,id,ih,min,is
1	format(' job run  ',i4,'/',i2,'/',i2,i9,':',i2,':',i2)
	return
	end
 
 
	subroutine datime
$include: 'picom.for'
	integer*4 mult
	call time(ih,min,is,if)
	mult=60
	itime=100*(is+mult*(min+mult*ih))+if
	return
	end
 
 
	subroutine defh
$include: 'picom.for'
	character*1 ga(26),gc(26),gn(10)
	data ga/'a','b','c','d','e','f','g','h','i','j','k','l','m'
     ,       ,'n','o','p','q','r','s','t','u','v','w','x','y','z'/
	data gc/'A','B','C','D','E','F','G','H','I','J','K','L','M'
     ,       ,'N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/
	data gn/'0','1','2','3','4','5','6','7','8','9'/
	do 10 i=1,10
10	hn(i)=gn(i)
	do 20 i=1,26
	ha(i)=ga(i)
20	hc(i)=gc(i)
	return
	end
 
 
	subroutine defsum(is,ittbv,ittev)
$include: 'picom.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(isb(is).eq.0)goto 20
	j=ios(is)
	minns=min0(minns,j)
	if(ittbv.ge.isb(is) .and. ittev.le.ise(is))goto 8
	if(j.eq.ns .and. ittbv.ge.isb(is))goto 8
	if(ittbv.gt.ise(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	isb(is)=ittbv
	ise(is)=ittev
	call pack
	return
20	if(ittbv.le.ise(iosi(ns)))call bomb(3,'defsum3 ')
	ns=ns+1
	ios(is)=ns
	iosi(ns)=is
	isb(is)=ittbv
	ise(is)=ittev
	minns=min0(minns,ns)
	call pack
	return
	end
 
 
	subroutine delet
$include: 'picom.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: 'picom.for'
	if(isb(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
	isb(is)=0
	ise(is)=0
	call pack
	return
	end
 
 
	subroutine deriv(j10,ind,jans)
c  d(sum j10)/d(var ind) is put in sum jans
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: 'picom.for'
	integer*4 id1,ivpic
	character*1 bl
	common/i/if
	data id1/1/
	nt1=nt+1
	if(isb(j10).lt.1 .or. ise(j10).gt.nt)call bomb(3,'deriv2  ')
	do 90 it=isb(j10),ise(j10)
	ia=itb(it)
	ib=ite(it)
	do 85 ic=ia,ib
	ie=ivn(ic)
	ivpic=ivp(ic)
	if(ie.ne.ind)goto 40
	call incnt
	itb(nt)=itetop+1
	call mulrat(itn(it),itd(it),ivpic,id1,itn(nt),itd(nt))
	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=idtab(i1,ie)
	if(j1.eq.0)goto 85
	k1=ider(2,j1)
	level=1
	if(k1.eq.ind)goto 50
	i2=1
42	j2=idtab(i2,k1)
	if(j2.eq.0)goto 77
	k2=ider(2,j2)
	level=2
	if(k2.eq.ind)goto 50
	i3=1
43	j3=idtab(i3,k2)
	if(j3.eq.0)goto 76
	k3=ider(2,j3)
	level=3
	if(k3.eq.ind)goto 50
	i4=1
44	j4=idtab(i4,k3)
	if(j4.eq.0)goto 75
	k4=ider(2,j4)
	level=4
	if(k4.eq.ind)goto 50
	i5=1
45	j5=idtab(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
	call mulrat(itn(it),itd(it),ivpic,id1,itn(nt),itd(nt))
	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: 'picom.for'
	common/i/if
	if=if+1
	call nite(if)
	ivn(if)=ider(3,j)
	ivp(if)=1
	return
	end
 
 
	subroutine dertab
$include: 'picom.for'
	idflg=0
	do 302 i=1,ndind
	do 302 j=1,nvar
302	idtab(i,j)=0
	do 304 i=1,nd
	idepv=ider(1,i)
	do 303 j=1,ndind
	if(idtab(j,idepv).eq.0)goto 304
303	continue
	call bomb(4,'NDIND   ')
304	idtab(j,idepv)=i
	return
	end
 
 
	subroutine dfine
$include: 'picom.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
1	format(1x,127a1)
	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: 'picom.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=',i4,' 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),isb(i),ise(i),ios(i),iosi(i)
5	format(' i=',i5,'   pv=',15a1,'   isb=',i5,'   ise=',i5
     +      ,'  ios=',i3,'  iosi=',i3)
10	continue
	write(8,11)
11	format(' ')
	write(8,13)
13	format('        i  itn  itd  itb  ite')
	do 20 i=1,nt,4
	j=i+1
	k=i+2
	l=i+3
	write(8,21)i,itn(i),itd(i),itb(i),ite(i),j,itn(j),itd(j),itb(j),
     +ite(j),k,itn(k),itd(k),itb(k),ite(k),l,itn(l),itd(l),itb(l),ite(l)
21	format(1x,4(i8,4i5))
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: 'picom.for'
	integer*4 kn,kd,ij,ik,jprvn,jprvd
	ia=isb(jb)
	if(ia.eq.0)call bomb(3,'expnd1  ')
	if(itn(ia).ne.1 .or. itd(ia).ne.1)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=ise(jb)
5	call copyt(ia,ia,1,1)
	if(ic.le.id)goto 10
6	call defsum(jans,nt,nt)
	return
10	ig=nt
	ie=isb(jp)
	kn=itn(ie)
	kd=itd(ie)
	if(kn.eq.0)goto 6
	jprvb=ig
	jprve=ig
	jprvn=1
	jprvd=1
	do 20 kp=1,nb
	if(kn.eq.0)goto 25
	call mulrat(jprvn,jprvd,kn,kd*kp,jprvn,jprvd)
	jtn(kp)=jprvn
	jtd(kp)=jprvd
	kn=kn-kd
	jtb(kp)=nt+1
	call mult(jprvb,jprve,ic,id)
	jte(kp)=nt
	if(jtb(kp).gt.nt)goto 25
	jprvb=jtb(kp)
20	jprve=nt
	kp=nb+1
25	kp=kp-1
	do 40 i=1,kp
	ih=jtb(i)
	ii=jte(i)
	ij=jtn(i)
	ik=jtd(i)
	do 30 il=ih,ii
30	call mulrat(itn(il),itd(il),ij,ik,itn(il),itd(il))
40	continue
	call defsum(jans,ig,ii)
	call order(jans)
	return
	end
 
 
	subroutine extrct
$include: 'picom.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=isb(j10)
	if(k.eq.0)call bomb(3,'extrct1 ')
	kk=itn(k)
	if(kk.gt.0 .and. itd(k).eq.1)goto 15
12	write(8,14)itn(k),itd(k),(a(i),i=j8,j9)
14	format(' bad term number  ',2i12,5x, 15a1)
	call bomb(3,'extrct2 ')
15	call nxt(j9+1,nchpl,'o',j11)
	call numv(j11+2,nchpl,j12,j13,j14)
	kl=ise(j14)+1-isb(j14)
	if(kk.gt.kl)goto 12
	km=isb(j14)+kk-1
	nt1=nt+1
	ktop=min0(km+num-1,ise(j14))
	call copyt(km,ktop,1,1)
	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,ii4)
	call srchnb(j6+1,nchpl,j8)
	call srch(j8,nchpl,' ',j14)
	call getint(j8,j14-1,j15,ii4)
	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(isb(j18).eq.0)call bomb(3,'extrct6 ')
	if(j15.gt.ise(j18)-isb(j18)+1)j15=ise(j18)-isb(j18)+1
	if(j15.lt.j13)call bomb(3,'extrct7 ')
	nt1=nt+1
	call copyt(isb(j18)+j13-1,isb(j18)+j15-1,1,1)
	call defsum(j4,nt1,nt)
	call outqu(j4)
	return
	end
 
 
	subroutine fort(it,ib)
$include: 'picom.for'
	character*1 pa,aa
	character*22 scr
	common/pr1/pa(232),aa(22)
	dimension dum(4)
	if(np.eq.0)return
	if(itd(it).eq.1 .and. itn(it).eq.-1)then
		ib=ib+1
		a(ib)='-'
		goto 42
		endif
	if(itd(it).eq.1 .and. itn(it).eq. 1)then
		ib=ib+1
		a(ib)='+'
		goto 42
		endif
	write(scr,43)itn(it)
	read(scr,9)(aa(kk),kk=1,20)
9	format(20a1)
	j=0
	do 10 i=1,20
	if(aa(i).eq.' ')goto 10
	if(j.eq.0 .and. aa(i).ne.'-')then
		ib=ib+1
		a(ib)='+'
		endif
	j=1
	ib=ib+1
	a(ib)=aa(i)
10	continue
	ib=ib+1
	a(ib)='.'
	if(itd(it).eq.1)goto 42
	ib=ib+1
	a(ib)='/'
	write(scr,43)itd(it)
	read(scr,9)(aa(kk),kk=1,20)
	do 15 i=1,20
	if(aa(i).eq.' ')goto 15
	ib=ib+1
	a(ib)=aa(i)
15	continue
42	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(scr,43)ip
	read(scr,9)(aa(kk),kk=1,20)
43	format(i20)
	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
	do 48 j=75,1,-1
	if(a(j).ne.' ')goto 49
48	continue
49	write(8,50)(a(i),i=1,j)
50	format(80a1)
	do 51 i=1,80
51	a(i)=' '
	a(6)='.'
	ib=11
52	ib=ib+1
	if(a(ib-1).ne.'+' .and. a(ib-1).ne.'-')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
	do 62 j=75,1,-1
	if(a(j).ne.' ')goto 63
62	continue
63	write(8,50)(a(i),i=1,j)
	return
	end
 
 
	subroutine fortrn
$include: 'picom.for'
	call nxt(1,nchpl,' ',j1)
	call numv(j1+1,nchpl,j2,j3,j4)
	if(np.eq.0 .or. isb(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=isb(j4)
	id=ise(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 .or. it.eq.id)goto 4
	a(6)=' '
	call copyc(pv(1,j4),a(7),ia)
	a(7+ia)='='
	call copyc(pv(1,j4),a(8+ia),ia)
	ib=7+ia+ia
	j=0
4	continue
10	return
	end
 
 
	subroutine gcd(inn,idd,ig)
	implicit integer*4 (i-n)
	ig=1
	in=iabs(inn)
	id=idd
	if(in.lt.id)goto 10
5	in=mod(in,id)
	if(in.eq.0)goto 20
10	if(in.eq.1)return
	id=mod(id,in)
	if(id.eq.0)goto 15
	if(id.eq.1)return
	goto 5
15	ig=in
	return
20	ig=id
	return
	end
 
 
	subroutine gcdout(in,id)
	implicit integer*4 (i-n)
	if(id.le.0)call bomb(3,'gcdout1 ')
	if(in.eq.1 .or. in.eq.-1 .or. id.eq.1)return
	if(in.eq.0)goto 1
	call gcd(in,id,ig)
	if(ig.eq.1)return
	in=in/ig
	id=id/ig
	return
1	id=1
	return
	end
 
 
	subroutine getint(j1,j2,ii2,j3)
c  on return j3=integer value of a(j1)-a(j2)
	implicit integer*2 (i-n)
	integer*4 j3
	j3=0
	if(j2.lt.j1)call bomb(3,' getint1')
	do 1 i=j1,j2
	call intval(i,j)
1	j3=10*j3+j
	ii2=j3
	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: 'picom.for'
	ie=1
	call srchnb(j1,nchpl,k)
	if(k.eq.-1)return
	if(a(k).eq.'*' .and. a(k+1).ne.'*')return
	isgn=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(isb(j5).eq.0)call bomb(3,'getexp10')
	ie=itn(isb(j5))
	if(itd(isb(j5)).ne.1)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	isgn=-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*isgn
	call srchnb(j3,nchpl,k)
	return
60	k=j3
	if(j3.eq.nchpl)call bomb(3,' getexp9')
	ie=me*isgn
	return
70	ie=me*isgn
	k=-1
	return
	end
 
 
	function icomp(ia,ib)
$include: 'picom.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(igoflg)
$include: 'picom.for'
	integer*4 i5,i8
	character*1 aj1,aj2
	igoflg=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)=' '
	igoflg=1
	call copyc(a,ap,nchpl)
	return
10	k1=itb(isb(j5))
	k2=ite(ise(j5))
	k3=itb(isb(j8))
	k4=ite(ise(j8))
	i5=itn(isb(j5))*itd(isb(j8))
	i8=itd(isb(j5))*itn(isb(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=isb(j5)
	k2=ise(j5)
	k3=isb(j8)
	k4=ise(j8)
	if(k2-k1.ne.k4-k3)goto 30
	k=k3
	do 12 i=k1,k2
	if(itn(i).ne.itn(k))goto 30
	if(itd(i).ne.itd(k))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(i5.lt.i8)goto 5
	return
24	if(i5.le.i8)goto 5
	return
25	if(aj2.eq.'e')goto 26
	if(i5.gt.i8)goto 5
	return
26	if(i5.ge.i8)goto 5
	return
30	if(aj1.eq.'n')goto 5
	return
	end
 
 
	subroutine incnt
$include: 'picom.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: 'picom.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: 'picom.for'
	call datime
	tprev=itime/100.
	call readin
	ntp=0
10	if(ha(1).ne.'a')call bomb(3,'monitr1 ')
	if(isb(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(isb(is).eq.0 .or. isb(js).eq.0)call bomb(3,'    mon3')
	ia=ise(is)
	ib=isb(js)
	if(ite(ia)+1.ne.itb(ib))call bomb(3,'    mon4')
	do 16 it=ib,ise(js)
	if(itb(it).gt.ite(it))call bomb(3,'    mon5')
	if(itd(it).le.0)call bomb(3,'   mon22')
	if(itn(it).ne.0)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(isb(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(ise(ic).ne.nt)call bomb(3,'   mon14')
	if(ite(ise(ic)).ne.itetop)call bomb(3,'   mon15')
11	call read(ib)
	if(ib.eq.1 .and. np.ge.1)call prntap
1	format(1x,127a1)
	if(ib.eq.1)goto 10
	call datime
	ttot=(itime-itbeg)/100.
	t=ttot-tprev
	tprev=ttot
	ntmntp=nt-ntp
	if(np.ge.3)write(8,3)t,ttot,ntmntp,nt,itetop,maxnt,maxite
3	format(10x,'t=',f7.2,' total t=',f9.1,' ntmntp=',i5,' nt=',i5,
     +       ' itetop=',i6,' maxnt=',i5,' maxite=',i6)
	if(np.ge.2)call prntap
5	ntp=nt
	call copyc(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(igoflg)
	if(igoflg.eq.1)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
	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 mulrat(n1,i1,n2,i2,n3,i3)
	implicit integer*4 (i-n)
	n3=n1*n2
	i3=i1*i2
	call gcdout(n3,i3)
	return
	end
 
 
	subroutine mult(ia,ib,ic,id)
c  multiply terms ia-ib by terms ic-id
c  result goes to nt+1,.....
$include: 'picom.for'
	common/nwt/iitb,iite,jitb,jite,ntb,if
	integer*4 inp,idp
	if(ia.eq.0 .or. ic.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 mulrat(itn(i),itd(i),itn(j),itd(j),inp,idp)
	call sub2(inp,idp)
10	continue
20	continue
	return
	end
 
 
	subroutine multn(i,j,k,l)
$include: 'picom.for'
	nt1=nt+1
	call mult(i,j,k,l)
	if(nt.ge.nt1)return
	call niltrm
	return
	end
 
 
	subroutine niltrm
$include: 'picom.for'
	k=itetop+1
	call incnt
	itn(nt)=0
	itd(nt)=1
	itb(nt)=k
	call nite(k)
	ivn(k)=2
	ivp(k)=0
	return
	end
 
 
	subroutine nite(j)
$include: 'picom.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: 'picom.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 copyc(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: 'picom.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(isb(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 copyc(a(j3),pv(1,nv),l)
	return
	end
 
 
	subroutine nxt(j1,j2,h,j3)
	implicit integer*2 (i-n)
	character*1 h
	call srch(j1,j2,h,j3)
	if(j3.eq.-1)call bomb(3,'     nxt')
	end
 
 
	subroutine nxtch(j1,j2,j3)
	implicit integer*2 (i-n)
	call srchch(j1,j2,j3)
	if(j3.eq.-1)call bomb(3,'   nxtch')
	return
	end
 
 
	subroutine nxtnb(j1,j2,j3)
	implicit integer*2 (i-n)
	call srchnb(j1,j2,j3)
	if(j3.eq.-1)call bomb(3,'   nxtnb')
	end
 
 
	subroutine nxtnnu(j1,j2,j3)
$include: 'picom.for'
	do 10 i=j1,j2
	do 5 j=1,10
	if(a(i).eq.hn(j))goto 10
5	continue
	j3=i
	return
10	continue
	j3=-1
	return
	end
 
 
	subroutine order(is)
$include: 'picom.for'
	ina=isb(is)
	inb=ise(is)
	if(inb.lt.ina .or. ina.eq.0)call bomb(3,'  order1')
	it=ina
2	ia=itb(it)
	if(itn(it).eq.0)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
	itn(ina)=0
	itd(ina)=1
	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(itn(it).eq.0)goto 75
	if(icomp(ina,it))71,73,74
71	call tzap(it,ina)
	goto 68
73	call addrat(itn(ina),itd(ina),itn(it),itd(it),itn(ina),itd(ina))
75	call bkup(it,inb)
	goto 69
74	if(icomp(it-1,it))62,61,68
61	call addrat(itn(it-1),itd(it-1),itn(it),itd(it),
     ,     itn(it-1),itd(it-1))
	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	call addrat(itn(lc),itd(lc),itn(it),itd(it),itn(lc),itd(lc))
	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: 'picom.for'
	call nxt(2,nchpl,' ',j1)
	call srchnb(j1,nchpl,j2)
	if(j2.eq.-1)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
1	format(1x,127a1)
	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 copyc(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: 'picom.for'
1	do 6 it=isb(j5),ise(j5)
	do 5 iv=itb(it),ite(it)
	if(isb(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: 'picom.for'
	if(minns.le.ns)goto 9
	is=iosi(ns)
	nt=ise(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=ise(iosi(js-1))
	kv=ite(kt)
13	do 30 ms=js,ns
	is=iosi(ms)
	if(is.le.0)call bomb(3,'   pack6')
	ia=isb(is)
	ib=ise(is)
	lv=0
	lt=0
	do 14 it=ia,ib
	if(itn(it).eq.0)goto 14
	lt=lt+1
	if(lt.gt.ntermt)call bomb(4,'NTERMT  ')
	jtn(lt)=itn(it)
	jtd(lt)=itd(it)
	jtb(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
	jvn(lv)=2
	jvp(lv)=0
	goto 33
31	lv=lv+1
	jvn(lv)=ivn(i)
	jvp(lv)=ivp(i)
12	continue
33	jte(lt)=lv
	if(lv.gt.ntott)call bomb(4,'NTOTT   ')
14	continue
	if(lv.ne.0)goto 15
	kt=kt+1
	isb(is)=kt
	ise(is)=kt
	kv=kv+1
	itn(kt)=0
	itd(kt)=1
	itb(kt)=kv
	ite(kt)=kv
	ivn(kv)=2
	ivp(kv)=0
	goto 30
15	isb(is)=kt+1
	do 18 it=1,lt
	kt=kt+1
	itn(kt)=jtn(it)
	itd(kt)=jtd(it)
	itb(kt)=kv+1
	do 16 i=jtb(it),jte(it)
	kv=kv+1
	ivn(kv)=jvn(i)
16	ivp(kv)=jvp(i)
18	ite(kt)=kv
	ise(is)=kt
30	continue
	nt=kt
	itetop=kv
	minns=ns+1
	return
	end
 
 
	subroutine print(ii)
$include: 'picom.for'
	character*127 scr
	character*1 pa,aa,bl,sp,sm,cf
	common/pr1/pa(232),aa(22)
	dimension dum(4)
	data bl,sp,sm,cf/' ','+','-','^'/
	if(np.eq.0)return
	ia=isb(ii)
	if(ia.eq.0)write(8,2)
	if(ia.eq.0)return
2	format(' there is nothing to print')
	ib=ise(ii)
	nst=23
	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)
	if(iap.le.99)goto 21
	pa(n+1)='*'
	pa(n+2)='*'
	n=n+2
	goto 30
21	n=n+1
	j=mod(iap,10)
	pa(n)=hn(j+1)
	if(iap.le.9)goto 30
	j=iap/10
	n=n+1
	pa(n)=pa(n-1)
	pa(n-1)=hn(j+1)
30	if(n.gt.232)call bomb(3,'  print3')
	write(scr,31)itn(it)
	read(scr,36)(pa(kk),kk=1,11)
31	format(i11)
	if(itd(it).eq.1)goto 33
	pa(12)='/'
	write(scr,31)itd(it)
	read(scr,36)(aa(kk),kk=1,11)
36	format(22a1)
	i=12
	do 32 j=1,11
	if(aa(j).eq.' ')goto 32
	i=i+1
	if(i.ge.nst)call bomb(3,'print3  ')
	pa(i)=aa(j)
32	continue
33	j=mod(it+1-ia,1000)
	ir=23
	if(pa(ir).ne.' ')goto 41
	do 50 m=ir,4,-1
	if(pa(m).ne.' ')goto 51
50	continue
51	mm=min0(ir-m,8)
	do 55 i=ir,mm+1,-1
	pa(i)=pa(i-mm)
55	pa(i-mm)=' '
41	write(8,34)j,(pa(k),k=1,n)
34	format(1x,i3,126a1,/,26x,106a1)
40	continue
	return
	end
 
 
	subroutine prntap
$include: 'picom.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(1x,127a1)
	return
	end
 
 
	subroutine read(ib)
$include: 'picom.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 copyc(ainpl(ib),a,il)
	call copyc(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: 'picom.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 copyc(a,ainpl(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: 'picom.for'
	integer*4 in,id
	ia=isb(j13)
	iz=ise(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 copyt(ia,ia,1,1)
	in=itn(ja)
	id=itd(ja)
	itn(ja)=id
	if(in.lt.0)itn(ja)=-id
	itd(ja)=iabs(in)
	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 copyt(ib,ic,1,1)
	jc=nt
	if(ic.eq.iz)goto 12
11	jd=nt+1
	call copyt(ic+1,iz,1,1)
	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
	itn(jd)=0
	itd(jd)=1
	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
	itn(jb)=0
	itd(jb)=1
	ic=ib-1
	goto 11
16	jf=nt+1
	call multn(ja,ja,jb,jc)
	jg=nt
	if(maxe.eq.1)goto 18
	jh=nt+1
	call multn(jf,jg,jf,jg)
	ji=nt
	jj=nt+1
	call multn(ja,ja,jd,je)
	jk=nt
	jl=nt+1
	call copyt(jf,jg,-1,1)
	call copyt(jh,ji,1,1)
	call copyt(jj,jk,-1,1)
	jm=nt
	call copyt(ja,ja,1,1)
	call multn(ja,ja,jl,jm)
	nt1=jm+1
	goto 24
18	do 19 i=jf,jg
19	itn(i)=-itn(i)
	call copyt(ja,ja,1,1)
	call multn(ja,ja,jf,jg)
	nt1=jg+1
24	call defsum(j4,nt1,nt)
	call order(j4)
	return
	end
 
 
	subroutine swap(j4)
$include: 'picom.for'
	if=isb(j4)
	ig=ios(j4)
	if(if.ne.0)minns=min0(minns,ios(j4))
	ios(j4)=ios(1)
	iosi(ios(j4))=j4
	isb(j4)=isb(1)
	ise(j4)=ise(1)
	if(if.eq.0)goto 5
	ios(1)=ig
	iosi(ig)=1
	isb(1)=if
	ise(1)=if
	ite(if)=itb(if)
	call pack
	return
5	call niltrm
	isb(1)=0
	call defsum(1,nt,nt)
	return
	end
 
 
	subroutine srch(j1,j2,h,j3)
$include: 'picom.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: 'picom.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: 'picom.for'
	do 1 j3=j1,j2
	if(a(j3).ne.' ')return
1	continue
	j3=-1
	return
	end
 
 
	subroutine sub
$include: 'picom.for'
	integer*4 inc4,idc4
	character*4 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(isb(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(isb(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,inc4,idc4)
	inc=inc4
	idc=idc4
	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(inc,idc,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(in,id,j8,j12)
$include: 'picom.for'
	integer*4 in4,id4
	if(isb(j12).eq.0)call bomb(3,'   subn1')
	if(isb(j8).ne.0)call bomb(3,'   subn2')
	do 8 it=isb(j12),ise(j12)
	do 5 i=itb(it),ite(it)
	if(ivn(i).ne.j8)goto 5
	in4=in
	id4=id
	in4=in4**ivp(i)
	id4=id4**ivp(i)
	call mulrat(itn(it),itd(it),in4,id4,itn(it),itd(it))
	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: 'picom.for'
	integer*4 inj1,idj1
	if(isb(ia).eq.0)call bomb(3,'   subr1')
	if(isb(ib).ne.0)call bomb(3,'   subr2')
	nt1=nt+1
	j1=isb(ia)
	j2=itb(j1)
	j3=ite(j1)
	j4=isb(ic)
	j5=ise(ic)
	do 30 it=j4,j5
	j6=itb(it)
	j7=ite(it)
	minrat=30000
	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
	idj1=itd(j1)
	if(itn(j1).lt.0)idj1=-idj1
	inj1=iabs(itn(j1))
	call mulrat(itn(it),itd(it),idj1,inj1,itn(nt),itd(nt))
	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 copyt(it,it,1,1)
30	continue
	call defsum(ic,nt1,nt)
	call order(ic)
	return
	end
 
 
	subroutine subs(in1,in2,in3)
c  substitute sum=in1 for variable=in2 in sum=in3
c  result goes in nt+1,...
$include: 'picom.for'
	integer*4 inp,idp
	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=isb(in3)
	ic=ise(in3)
	id=isb(in1)
	ie=ise(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(itn(it),itd(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 mulrat(itn(it),itd(it),itn(jt),itd(jt),inp,idp)
	call sub2(inp,idp)
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
	jvn(k)=ivn(i)
71	jvp(k)=ivp(i)*ifp
	jitb=1
	jite=k
	if(ifp.lt.0)goto 72
	call mulrat(itn(it),itd(it),itn(id)**ifp,itd(id)**ifp,inp,idp)
	goto 73
72	ii=-ifp
	call mulrat(itn(it),itd(it),itd(id)**ii,itn(id)**ii,inp,idp)
73	call sub3(inp,idp)
	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
	itn(nt)=jtn(i)
	itd(nt)=jtd(i)
	itb(nt)=k+1
	do 84 j=jtb(i),jte(i)
	k=k+1
	ivn(k)=jvn(j)
84	ivp(k)=jvp(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
	jtn(nts)=itn(i)
	jtd(nts)=itd(i)
	jtb(nts)=ks+1
	do 20 j=itb(i),ite(i)
	ks=ks+1
	jvn(ks)=ivn(j)
20	jvp(ks)=ivp(j)
	jte(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=jtb(jt)
	jite=jte(jt)
	if(jvn(jitb).eq.2)lev2=jvp(jitb)
	if(lev1+lev2.gt.maxe)goto 30
	call mulrat(itn(it),itd(it),jtn(jt),jtd(jt),inp,idp)
	call sub3(inp,idp)
54	continue
30	continue
31	call defsum(in3,ntb,nt)
	return
	end
 
 
	subroutine subv(j4,j8,j12)
$include: 'picom.for'
	if(isb(j12).eq.0)call bomb(3,'   subv1')
	if(isb(j4)+isb(j8).ne.0)call bomb(3,'   subv2')
	do 8 it=isb(j12),ise(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 sub2(in,id)
$include: 'picom.for'
	common/nwt/iitb,iite,jitb,jite,ntb,if
	integer*4 in,id
	if(in.eq.0)return
	call incnt
	k=itetop
	itb(nt)=k+1
	itn(nt)=in
	itd(nt)=id
	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	call addrat(itn(ntb),itd(ntb),itn(nt),itd(nt),itn(ntb),itd(ntb))
	nt=nt-1
	itetop=itesav
	goto 68
34	if(icomp(nt-1,nt))62,61,68
61	call addrat(itn(nt-1),itd(nt-1),itn(nt),itd(nt),
     ,     itn(nt-1),itd(nt-1))
	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	call addrat(itn(lc),itd(lc),itn(nt),itd(nt),itn(lc),itd(lc))
	nt=nt-1
	itetop=itesav
	goto 68
67	la=lc
	goto 63
68	return
	end
 
 
	subroutine sub3(in,id)
$include: 'picom.for'
	integer*4 in,id
	common/nwt/iitb,iite,jitb,jite,ntb,if
	if(in.eq.0)return
	call incnt
	k=itetop
	itb(nt)=k+1
	itn(nt)=in
	itd(nt)=id
	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.jvn(jj))goto 16
	if(ivn(ii).eq.jvn(jj))goto 57
	if(jvp(jj).eq.0)goto 18
	k=k+1
	ivn(k)=jvn(jj)
	ivp(k)=jvp(jj)
18	jj=jj+1
	goto 11
12	ii=ii+1
	goto 11
13	if(jj.gt.jite)goto 60
	if(jvp(jj).eq.0)goto 17
	k=k+1
	ivn(k)=jvn(jj)
	ivp(k)=jvp(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)+jvp(jj).eq.0)goto 58
	k=k+1
	ivn(k)=ivn(ii)
	ivp(k)=ivp(ii)+jvp(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	call addrat(itn(ntb),itd(ntb),itn(nt),itd(nt),itn(ntb),itd(ntb))
	nt=nt-1
	itetop=itesav
	goto 68
34	if(icomp(nt-1,nt))62,61,68
61	call addrat(itn(nt-1),itd(nt-1),itn(nt),itd(nt),
     ,     itn(nt-1),itd(nt-1))
	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	call addrat(itn(lc),itd(lc),itn(nt),itd(nt),itn(lc),itd(lc))
	nt=nt-1
	itetop=itesav
	goto 68
67	la=lc
	goto 63
68	return
	end
 
 
	subroutine switch(j,k,is)
$include: 'picom.for'
	do 2 it=isb(is),ise(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: 'picom.for'
	character*4 star
	integer*4 itf
	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(isb(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=isb(j9),ise(j9)
	call copyt(it,it,1,1)
	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,1,iv,j10)
	itf=1
	do 30 io=1,iorder
	call deriv(j11,iv,j12)
	call outqu(j12)
	if(itn(isb(j12)).eq.0)goto 31
	call copys(j12,j11)
	call subn(0,1,iv,j12)
	itf=itf*io
	call incnt
	k=itetop+1
	itn(nt)=1
	itd(nt)=itf
	itb(nt)=k
	ite(nt)=k
	ivn(k)=iv
	ivp(k)=io
	call nite(k)
	nt1=nt+1
	call mult(nt,nt,isb(j12),ise(j12))
	call defsum(j12,nt1,nt)
	call pack
30	call adds(j12,j10,1,1,j10)
31	if(star.eq.'yes')call switch(nv1,iv,j10)
32	continue
34	call adds(j10,1,1,1,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: 'picom.for'
	integer*4 inc,idc
	call incnt
	call coef(ist,j1,inc,idc)
	itn(nt)=inc
	itd(nt)=idc
	itb(nt)=itetop+1
	if(j1.eq.-1)then
		k=itb(nt)
		ivn(k)=2
		ivp(k)=0
		goto 10
		endif
	k=itetop
5	isgn=1
	if(a(j1).eq.'*')j1=j1+1
	if(a(j1).eq.'/')isgn=-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*isgn
	j1=j5
	if(j5.ne.-1)goto 5
10	call nite(k)
	return
	end
 
 
	subroutine thderv
$include: 'picom.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: 'picom.for'
	integer*4 itnka,itdka
	if(ka.eq.kb)return
	itnka=itn(ka)
	itdka=itd(ka)
	itbka=itb(ka)
	iteka=ite(ka)
	if(ka.gt.kb)goto 5
	do 1 i=ka,kb-1
	itn(i)=itn(i+1)
	itd(i)=itd(i+1)
	itb(i)=itb(i+1)
1	ite(i)=ite(i+1)
	goto 8
5	do 6 i=ka,kb+1,-1
	itn(i)=itn(i-1)
	itd(i)=itd(i-1)
	itb(i)=itb(i-1)
6	ite(i)=ite(i-1)
8	itn(kb)=itnka
	itd(kb)=itdka
	itb(kb)=itbka
	ite(kb)=iteka
	return
	end
 
 
 
	subroutine zsub(j3)
$include: 'picom.for'
c this is available
	return
	end