home *** CD-ROM | disk | FTP | other *** search
/ HomeWare 14 / HOMEWARE14.bin / educate / algtut.arj / PFSAI.FOR < prev    next >
Text File  |  1985-10-03  |  70KB  |  2,951 lines
  1. $include: 'picom.for'
  2.       open(unit=7,file='algin',status='old')
  3.       open(unit=8,file='algout',status='new')
  4.       call datim
  5.       itbeg=itime
  6.       do 4 i=1,nvar
  7.     4 isb(i)=0
  8.       ios(1)=1
  9.       iosi(1)=1
  10.       ipv(1)=1
  11.       pv(1,1)='?'
  12.       isb(1)=1
  13.       ise(1)=1
  14.       itb(1)=1
  15.       ite(1)=1
  16.       ivn(1)=2
  17.       ivp(1)=1
  18.       nd=0
  19.       nv=1
  20.       ns=1
  21.       np=2
  22.       nt=1
  23.       nb=7
  24.       maxe=9999
  25.       maxnt=0
  26.       maxite=0
  27.       itetop=1
  28.       ndump=0
  29.       minns=9999
  30.       lintot=0
  31.       idflg=1
  32.       ha(1)='a'
  33.       call defc(ha,26)
  34.       hc(1)='A'
  35.       call defc(hc,26)
  36.       hn(1)='0'
  37.       call defc(hn,10)
  38.       ma=ntot
  39.       mb=nvar
  40.       mc=nterm
  41.       md=ninch
  42.       me=ninl
  43.       mg=nder
  44.       mh=nlab
  45.       mi=ndind
  46.       mj=mitay
  47.       write(8,1)ma,mb,mc,md,me,mg,mh,mi,mj
  48.     1 format(' version 85/6/20 ',
  49.      +   ' NTOT=',i6,' NVAR=',i4,' NTERM=',i5,' NINCH=',i5
  50.      +   ,/,' NINL=',i6,' NDER=',i3,' NLAB=',i3,' NDIND=',i2
  51.      +   ,' MITAY=',i2)
  52.       call monitr
  53.       call datim
  54.       write(8,90)nd,nv,ns,maxnt,maxite,lintot
  55.    90 format(' nu. derivatives=',i3,'  nu. variables=',i3,
  56.      +    '  nu. sums=',i3,/,' max nu. terms=',i5,
  57.      +  '  max total nu. variables=',i6,'  input lines=',i5)
  58.       t=(itime-itbeg)/100.
  59.       write(8,91)t
  60.    91 format(' time used=',f9.2,' seconds')
  61.       stop
  62.       end
  63.  
  64.  
  65.       subroutine addrat(n1,i1,n2,i2,n3,i3)
  66. c  n3/i3=n1/i1 + n2/i2
  67.       implicit integer*4 (i-n)
  68.       i=i1*i2
  69.       n=n1*i2+n2*i1
  70.       call gcdout(n,i)
  71.       n3=n
  72.       i3=i
  73.       return
  74.       end
  75.  
  76.  
  77.       subroutine adds(is,js,iin,iid,ians)
  78. $include: 'picom.for'
  79.       integer*4 in1,id1,in,id,inn,idd,insav,idsav
  80.       in=iin
  81.       id=iid
  82.       nt1=nt+1
  83.       ib=isb(is)
  84.       ie=ise(is)
  85.       jb=isb(js)
  86.       je=ise(js)
  87.       i=itb(ie)
  88.       if(ivn(i).ne.2 .or. ivp(i).le.maxe)goto 3
  89.       do 1 it=ib,ie
  90.       i=itb(it)
  91.       if(ivn(i).eq.2 .and. ivp(i).gt.maxe)goto 2
  92.     1 continue
  93.     2 ie=it-1
  94.     3 i=itb(je)
  95.       if(ivn(i).ne.2 .or.ivp(i).le.maxe)goto 7
  96.       do 4 it=jb,je
  97.       i=itb(it)
  98.       if(ivn(i).eq.2 .and. ivp(i).gt.maxe)goto 5
  99.     4 continue
  100.     5 je=it-1
  101.     7 if(ib.gt.ie)goto 70
  102.       if(jb.gt.je)goto 12
  103.       if(icomp(ib,jb))8,60,30
  104.     8 jl=jb
  105.       if(jb.eq.je)goto 10
  106.       if(icomp(ib,je))10,16,18
  107.    10 call copyt(jb,je,iin,iid)
  108.    12 call copyt(ib,ie,1,1)
  109.    14 call defsum(ians,nt1,nt)
  110.       return
  111.    16 insav=itn(je)
  112.       idsav=itd(je)
  113.       idd=id
  114.       if(in.lt.0)idd=-idd
  115.       inn=iabs(in)
  116.       call mulrat(itn(ib),itd(ib),idd,inn,in1,id1)
  117.       call addrat(in1,id1,itn(je),itd(je),itn(je),itd(je))
  118.       call copyt(jb,je,iin,iid)
  119.       itn(je)=insav
  120.       itd(je)=idsav
  121.       ib=ib+1
  122.       jb=je+1
  123.       goto 7
  124.    18 jr=je
  125.    20 if(jr-jl.eq.1)goto 28
  126.       jc=(jl+jr)/2
  127.       if(icomp(ib,jc))22,26,24
  128.    22 jl=jc
  129.       goto 20
  130.    24 jr=jc
  131.       goto 20
  132.    26 insav=itn(jc)
  133.       idsav=itd(jc)
  134.       idd=id
  135.       if(in.lt.0)idd=-idd
  136.       inn=iabs(in)
  137.       call mulrat(itn(ib),itd(ib),idd,inn,in1,id1)
  138.       call addrat(in1,id1,itn(jc),itd(jc),itn(jc),itd(jc))
  139.       call copyt(jb,jc,iin,iid)
  140.       itn(jc)=insav
  141.       itd(jc)=idsav
  142.       ib=ib+1
  143.       jb=jc+1
  144.       goto 7
  145.    28 call copyt(jb,jl,iin,iid)
  146.       jb=jr
  147.    30 il=ib
  148.       if(ib.eq.ie)goto 31
  149.       if(icomp(jb,ie))31,36,38
  150.    31 call copyt(ib,ie,1,1)
  151.    32 call copyt(jb,je,iin,iid)
  152.       goto 14
  153.    36 insav=itn(ie)
  154.       idsav=itd(ie)
  155.       call mulrat(in,id,itn(jb),itd(jb),in1,id1)
  156.       call addrat(in1,id1,itn(ie),itd(ie),itn(ie),itd(ie))
  157.       call copyt(ib,ie,1,1)
  158.       itn(ie)=insav
  159.       itd(ie)=idsav
  160.       jb=jb+1
  161.       ib=ie+1
  162.       goto 70
  163.    38 ir=ie
  164.    40 if(ir-il.eq.1)goto 48
  165.       ic=(il+ir)/2
  166.       if(icomp(jb,ic))42,46,44
  167.    42 il=ic
  168.       goto 40
  169.    44 ir=ic
  170.       goto 40
  171.    46 insav=itn(ic)
  172.       idsav=itd(ic)
  173.       call mulrat(in,id,itn(jb),itd(jb),in1,id1)
  174.       call addrat(in1,id1,itn(ic),itd(ic),itn(ic),itd(ic))
  175.       call copyt(ib,ic,1,1)
  176.       itn(ic)=insav
  177.       itd(ic)=idsav
  178.       jb=jb+1
  179.       ib=ic+1
  180.       goto 7
  181.    48 call copyt(ib,il,1,1)
  182.       ib=ir
  183.       goto 8
  184.    60 insav=itn(ib)
  185.       idsav=itd(ib)
  186.       call mulrat(in,id,itn(jb),itd(jb),in1,id1)
  187.       call addrat(in1,id1,itn(ib),itd(ib),itn(ib),itd(ib))
  188.       call copyt(ib,ib,1,1)
  189.       itn(ib)=insav
  190.       itd(ib)=idsav
  191.       ib=ib+1
  192.       jb=jb+1
  193.       goto 7
  194.    70 if(jb.gt.je)goto 14
  195.       goto 32
  196.       end
  197.  
  198.  
  199.       subroutine bomb(i,h)
  200. $include: 'picom.for'
  201.       character*8 h
  202.       goto(100,20,30,40,50),i
  203.    20 write(8,21)
  204.    21 format(' parameter error')
  205.       goto 100
  206.    30 write(8,31)h
  207.    31 format(1x,a8)
  208.       goto 100
  209.    40 write(8,41)h,h
  210.    41 format(1x,a8,' is too small for your problem.',
  211.      +       '  Recompile with bigger ',a8)
  212.    50 continue
  213.   100 write(8,101)line,ap
  214.   101 format(' BOMB: input line nu.',i5,' was:',/,1x,127a1)
  215.       write(8,90)nd,nv,ns,maxnt,maxite,line
  216.    90 format(' nu. derivatives=',i3,
  217.      +   '   nu. variables=',i3,' nu. sums=',i3,/,'  max nu. terms=',
  218.      +      i5,'  max total nu. variables=',i6,' input lines=',i5)
  219.       call datim
  220.       if(ndump.ne.0)call dump
  221.       t=(itime-itbeg)/100.
  222.       write(8,91)t
  223.    91 format(' time used=',f9.2,' seconds')
  224.       stop
  225.       return
  226.       end
  227.  
  228.  
  229.       subroutine bkup(j,k)
  230. $include: 'picom.for'
  231.       if(j.eq.k)goto 2
  232.       do 1 i=j,k-1
  233.       itb(i)=itb(i+1)
  234.       ite(i)=ite(i+1)
  235.       itn(i)=itn(i+1)
  236.     1 itd(i)=itd(i+1)
  237.     2 k=k-1
  238.       return
  239.       end
  240.  
  241.  
  242.       subroutine call
  243. $include: 'picom.for'
  244.       call nxt(3,nchpl,' ',j1)
  245.       call numl(j1,nchpl,j2,j3,j4)
  246.       if(linlab(j4).eq.-1)call bomb(3,'call    ')
  247.       linret=line
  248.       line=linlab(j4)
  249.       return
  250.       end
  251.  
  252.  
  253.       subroutine coef(is,jv,in,id)
  254. c is=start, on return jv points beyond coef
  255. $include: 'picom.for'
  256.       integer*4 in,id,ii,ij
  257.       id=1
  258.       isgn=1
  259.       call srchnb(is,nchpl,j1)
  260.       if(j1.ne.-1)goto 5
  261.     1 jv=-1
  262.     2 in=isgn
  263.       return
  264.     5 isgn=-1
  265.       if(a(j1).eq.'-')goto 30
  266.       isgn=1
  267.       if(a(j1).eq.'+')goto 30
  268.       j2=j1
  269.     7 jv=j2
  270.       call intval(j2,i)
  271.       if(i.eq.-1)goto 2
  272.       call nxtnnu(j2,nchpl,j3)
  273.       call getint(j2,j3-1,in2,in)
  274.       j5=j3
  275.       if(a(j3).eq.'.')goto 10
  276.       if(a(j3).eq.'/')goto 25
  277.       if(a(j3).eq.' ')goto 20
  278.       goto 15
  279.    10 call intval(j3+1,i)
  280.       if(i.ne.-1)goto 16
  281.       j5=j3+1
  282.       if(a(j5).eq.'/')goto 25
  283.       if(a(j5).eq.' ')goto 20
  284.    15 jv=j5
  285.       in=in*isgn
  286.       return
  287.    16 call nxtnnu(j3+1,nchpl,j7)
  288.       call getint(j3+1,j7-1,in2,ii)
  289.       ij=10**(j7-j3-1)
  290.       in=ij*in+ii
  291.       id=ij
  292.       call gcdout(in,id)
  293.       goto 26
  294.    20 call srchnb(j5,nchpl,j1)
  295.       j5=j1
  296.       if(j5.eq.-1)goto 15
  297.       if(a(j5).eq.'/')goto 25
  298.       goto 15
  299.    25 call nxtnb(j5+1,nchpl,j6)
  300.       call intval(j6,i)
  301.       if(i.eq.-1)goto 15
  302.       call nxtnnu(j6,nchpl,j7)
  303.       call getint(j6,j7-1,in2,id)
  304.    26 if(a(j7).eq.'.')j7=j7+1
  305.       call srchnb(j7,nchpl,jv)
  306.       in=in*isgn
  307.       return
  308.    30 call srchnb(j1+1,nchpl,j2)
  309.       if(j2.eq.-1)goto 1
  310.       goto 7
  311.       end
  312.  
  313.  
  314.       subroutine comput
  315. $include: 'picom.for'
  316.   300 call nxt(1,nchpl,' ',j1)
  317.       call numv(j1,nchpl,j2,j3,j4)
  318.       call nxt(j3,nchpl,'=',j5)
  319.       call srch(j5,nchpl,'/',j6)
  320.       if(j6.eq.-1)goto 320
  321. c  /
  322.       call nxtnb(j5+1,j6-1,j7)
  323.       if(a(1).eq.'d')goto 301
  324.       call numv(j7,j6-1,j8,j9,j10)
  325.       call numv(j6+1,nchpl,j11,j12,j13)
  326.       if(isb(j10).eq.0 .or. isb(j13).eq.0)call bomb(3,'not sum ')
  327.       call recip(j13,1)
  328.       nt1=nt+1
  329.       call mult(isb(j10),ise(j10),isb(1),ise(1))
  330.       call defsum(j4,nt1,nt)
  331.       call outqu(j4)
  332.       minns=min0(minns,ios(1))
  333.       it=isb(1)
  334.       ise(1)=it
  335.       ite(it)=itb(it)
  336.       k=itb(it)
  337.       ivn(k)=2
  338.       ivp(k)=0
  339.       call pack
  340.       return
  341. c  compute sj4=dsj10/dfj14
  342.   301 call numv(j7+1,j6-1,j8,j9,j10)
  343.       if(a(j7).ne.'d')call bomb(3,'comput1 ')
  344.       if(idflg.ne.0)call dertab
  345.       call nxt(j6,nchpl,'d',j11)
  346.       call numv(j11+1,nchpl,j12,j13,j14)
  347.       call srch(j13,nchpl,'=',j15)
  348.       iord=1
  349.       if(j15.eq.-1)goto 311
  350.       call nxtnb(j15+1,nchpl,j16)
  351.       call intval(j16,j17)
  352.       if(j17.ne.-1)goto 360
  353.       call numv(j16,nchpl,j19,j20,j18)
  354.       if(isb(j18).eq.0)call bomb(3,'comput14')
  355.       iord=itn(isb(j18))
  356.       if(itd(isb(j18)).ne.1)call bomb(3,'comput15')
  357.       if(iord.gt.0)goto 311
  358.       call bomb(3,'comput16')
  359.   360 call getexp(j15+1,iord,j16,0)
  360.       if(j16.ne.-1 .or. iord.le.0)call bomb(3,'comput2 ')
  361.   311 do 313 ii=1,iord
  362.       if(ii.eq.1)call deriv(j10,j14,1)
  363.       if(ii.ne.1)call deriv(1,j14,1)
  364.       call outqu(1)
  365.       if(itn(isb(1)).eq.0)goto 314
  366.   313 continue
  367.   314 call swap(j4)
  368.       return
  369. c  not /
  370.   320 call srch(j5,nchpl,'+',j6)
  371.       if(a(1).eq.'d')call bomb(3,'comput3 ')
  372.       id12=1
  373.       if(j6.ne.-1)goto 321
  374.       call srch(j5,nchpl,'-',j6)
  375.       in12=-1
  376.       if(j6.ne.-1)goto 322
  377.       call srch(j5,nchpl,'*',j6)
  378.       if(j6.ne.-1)goto 339
  379.       call srch(j5,nchpl,'^',j6)
  380.       if(j6.ne.-1)goto 339
  381. c  compute sj4=sj12
  382.       call numv(j5,nchpl,j6,j11,j12)
  383.       if(isb(j12).eq.0)call bomb(3,'comput4 ')
  384.       call copys(j12,j4)
  385.       call outqu(j4)
  386.       return
  387. c  compute sj4=sj9  + or -  sj12
  388.   321 in12=1
  389.   322 call numv(j5+1,j6-1,j7,j8,j9)
  390.       call numv(j6+1,nchpl,j10,j11,j12)
  391.       call adds(j9,j12,in12,id12,j4)
  392.       call outqu(j4)
  393.       return
  394. c  compute sj4=sj9*sj12
  395.   339 call numv(j5+1,j6-1,j7,j8,j9)
  396.       if(a(j6).eq.'^')goto 350
  397.       if(a(j6+1).eq.'*')goto 349
  398.       call numv(j6+1,nchpl,j10,j11,j12)
  399.       nt1=nt+1
  400.       call mult(isb(j9),ise(j9),isb(j12),ise(j12))
  401.       call defsum(j4,nt1,nt)
  402.       call outqu(j4)
  403.       return
  404.   349 j6=j6+1
  405.   350 call numv(j6+1,nchpl,j10,j11,j12)
  406.       call expnd(j9,j12,j4)
  407.       return
  408.       end
  409.  
  410.  
  411.       subroutine copyc(a,b,n)
  412.       integer*2 i,n
  413.       character*1 a,b
  414.       dimension a(2),b(2)
  415.       do 1 i=1,n
  416.     1 b(i)=a(i)
  417.       return
  418.       end
  419.  
  420.  
  421.       subroutine copys(j1,j2)
  422. $include: 'picom.for'
  423.       nt1=nt+1
  424.       call copyt(isb(j1),ise(j1),1,1)
  425.       call defsum(j2,nt1,nt)
  426.       return
  427.       end
  428.  
  429.  
  430.       subroutine count
  431. $include: 'picom.for'
  432.       call nxt(3,nchpl,' ',j1)
  433.       call numv(j1,nchpl,j3,j4,j5)
  434.       do 12 i=nchpl,j4+1,-1
  435.       if(a(i).ne.' ')goto 14
  436.    12 continue
  437.       call bomb(3,'count1  ')
  438.    14 do 16 j=i,j4+1,-1
  439.       if(a(j).eq.' ')goto 18
  440.    16 continue
  441.    18 call numv(j,i,j6,j7,j8)
  442.       if(isb(j8).eq.0)call bomb(3,'count2  ')
  443.       if(isb(j5).ne.0)goto 20
  444.       k=itetop+1
  445.       call incnt
  446.       itb(nt)=k
  447.       call nite(k)
  448.       itn(nt)=ise(j8)+1-isb(j8)
  449.       itd(nt)=1
  450.       ivn(k)=2
  451.       ivp(k)=0
  452.       call defsum(j5,nt,nt)
  453.       return
  454.    20 ia=isb(j5)
  455.       ib=ise(j5)
  456.       k=itb(ia)
  457.       l=ite(ia)
  458.       ivn(k)=2
  459.       ivp(k)=0
  460.       ise(j5)=ia
  461.       ite(ia)=k
  462.       itn(ia)=ise(j8)+1-isb(j8)
  463.       itd(ia)=1
  464.       if(ia.ne.ib .or. k.ne.l)minns=min0(minns,ios(j5))
  465.       call pack
  466.       return
  467.       end
  468.  
  469.  
  470.       subroutine copyt(ia,ib,in,id)
  471. c  copy terms ia-ib to nt+1,..., coefficients are multiplied by in/id
  472. $include: 'picom.for'
  473.       integer*4 in4,id4
  474.       in4=in
  475.       id4=id
  476.       if(ib.lt.ia)return
  477.       if(ia.lt.1 .or. ib.gt.nt)call bomb(3,'copyt2  ')
  478.       icop=0
  479.       if(id.eq.1 .and. in.eq.-1)icop=-1
  480.       if(id.eq.1 .and. in.eq. 1)icop= 1
  481.       do 15 i=ia,ib
  482.       ja=itb(i)
  483.       if(ivn(ja).eq.2 .and. ivp(ja).gt.maxe)goto 15
  484.       k=itetop
  485.       call incnt
  486.       itn(nt)=itn(i)*icop
  487.       itd(nt)=itd(i)
  488.       if(icop.eq.0)call
  489.      +     mulrat(itn(i),itd(i),in4,id4,itn(nt),itd(nt))
  490.       itb(nt)=k+1
  491.       do 10 j=ja,ite(i)
  492.       k=k+1
  493.       ivn(k)=ivn(j)
  494.    10 ivp(k)=ivp(j)
  495.       call nite(k)
  496.    15 continue
  497.    16 return
  498.       end
  499.  
  500.  
  501.       subroutine datim
  502. $include: 'picom.for'
  503.       integer*4 mult
  504.       call date(iy,mon,id)
  505.       call time(ih,min,is,if)
  506.       mult=60
  507.       itime=100*(is+mult*(min+mult*ih))+if
  508.       write(8,1)iy,mon,id,ih,min,is
  509.     1 format(' job run  ',i4,'/',i2,'/',i2,i9,':',i2,':',i2)
  510.       return
  511.       end
  512.  
  513.  
  514.       subroutine datime
  515. $include: 'picom.for'
  516.       integer*4 mult
  517.       call time(ih,min,is,if)
  518.       mult=60
  519.       itime=100*(is+mult*(min+mult*ih))+if
  520.       return
  521.       end
  522.  
  523.  
  524.       subroutine defc(a,n)
  525.       character*1 a(2)
  526.       integer*2 i,j
  527.       do 1 i=2,n
  528.       j=a(i-1)
  529.     1 a(i)=j+1
  530.       return
  531.       end
  532.  
  533.  
  534.       subroutine defsum(is,ittbv,ittev)
  535. $include: 'picom.for'
  536.       if(ittev.ne.nt .or. ittev+1.ne.ittbv)goto 1
  537.       call niltrm
  538.       ittev=nt
  539.     1 if(ittbv.gt.ittev)call bomb(3,'defsum1 ')
  540.       if(isb(is).eq.0)goto 20
  541.       j=ios(is)
  542.       minns=min0(minns,j)
  543.       if(ittbv.ge.isb(is) .and. ittev.le.ise(is))goto 8
  544.       if(j.eq.ns .and. ittbv.ge.isb(is))goto 8
  545.       if(ittbv.gt.ise(iosi(ns)))goto 4
  546.       call bomb(3,'defsum2 ')
  547.     4 ios(is)=ns
  548.       do 6 i=j,ns-1
  549.       k=iosi(i+1)
  550.       iosi(i)=k
  551.     6 ios(k)=i
  552.       iosi(ns)=is
  553.     8 isb(is)=ittbv
  554.       ise(is)=ittev
  555.       call pack
  556.       return
  557.    20 if(ittbv.le.ise(iosi(ns)))call bomb(3,'defsum3 ')
  558.       ns=ns+1
  559.       ios(is)=ns
  560.       iosi(ns)=is
  561.       isb(is)=ittbv
  562.       ise(is)=ittev
  563.       minns=min0(minns,ns)
  564.       call pack
  565.       return
  566.       end
  567.  
  568.  
  569.       subroutine delet
  570. $include: 'picom.for'
  571. c  delete sj4
  572.       call nxt(4,nchpl,' ',j1)
  573.     1 call numv(j1,nchpl,j2,j3,j4)
  574.       call delsum(j4)
  575.       if(j3.eq.nchpl)goto 10
  576.       call srchnb(j3+1,nchpl,j1)
  577.       if(j1.ne.-1)goto 1
  578.    10 return
  579.       end
  580.  
  581.  
  582.       subroutine delsum(is)
  583. $include: 'picom.for'
  584.       if(isb(is).eq.0)return
  585.       j=ios(is)
  586.       minns=min0(minns,j)
  587.       if(j.eq.0 .or. ns.eq.0)call bomb(3,'delsum1 ')
  588.       ns=ns-1
  589.       if(j.gt.ns)goto 10
  590.       do 2 i=j,ns
  591.       k=iosi(i+1)
  592.       iosi(i)=k
  593.     2 ios(k)=i
  594.    10 iosi(ns+1)=0
  595.       ios(is)=0
  596.       isb(is)=0
  597.       ise(is)=0
  598.       call pack
  599.       return
  600.       end
  601.  
  602.  
  603.       subroutine deriv(j10,ind,jans)
  604. c  d(sum j10)/d(var ind) is put in sum jans
  605. c  ider(1,nd) is dependent variable in derivative
  606. c  ider(2,nd) is independent variable in derivative
  607. c  ider(3,nd) is variable equal to derivative
  608. $include: 'picom.for'
  609.       integer*4 id1,ivpic
  610.       character*1 bl
  611.       common/i/if
  612.       data id1/1/
  613.       nt1=nt+1
  614.       if(isb(j10).lt.1 .or. ise(j10).gt.nt)call bomb(3,'deriv2  ')
  615.       do 90 it=isb(j10),ise(j10)
  616.       ia=itb(it)
  617.       ib=ite(it)
  618.       do 85 ic=ia,ib
  619.       ie=ivn(ic)
  620.       ivpic=ivp(ic)
  621.       if(ie.ne.ind)goto 40
  622.       call incnt
  623.       itb(nt)=itetop+1
  624.       call mulrat(itn(it),itd(it),ivpic,id1,itn(nt),itd(nt))
  625.       if=itetop
  626.       do 36 i=ia,ib
  627.       if(i.eq.ic)goto 32
  628.       if=if+1
  629.       ivn(if)=ivn(i)
  630.       ivp(if)=ivp(i)
  631.       goto 36
  632.    32 if(ivp(i).eq.1 .and. ia.ne.ib)goto 36
  633.       if=if+1
  634.       ivn(if)=ivn(i)
  635.       ivp(if)=ivp(i)-1
  636.    36 continue
  637.       call nite(if)
  638.       if(if.lt.itb(nt))call bomb(3,'deriv3  ')
  639.       goto 85
  640.    40 if(nd.eq.0)goto 85
  641.       i1=1
  642.    41 j1=idtab(i1,ie)
  643.       if(j1.eq.0)goto 85
  644.       k1=ider(2,j1)
  645.       level=1
  646.       if(k1.eq.ind)goto 50
  647.       i2=1
  648.    42 j2=idtab(i2,k1)
  649.       if(j2.eq.0)goto 77
  650.       k2=ider(2,j2)
  651.       level=2
  652.       if(k2.eq.ind)goto 50
  653.       i3=1
  654.    43 j3=idtab(i3,k2)
  655.       if(j3.eq.0)goto 76
  656.       k3=ider(2,j3)
  657.       level=3
  658.       if(k3.eq.ind)goto 50
  659.       i4=1
  660.    44 j4=idtab(i4,k3)
  661.       if(j4.eq.0)goto 75
  662.       k4=ider(2,j4)
  663.       level=4
  664.       if(k4.eq.ind)goto 50
  665.       i5=1
  666.    45 j5=idtab(i5,k4)
  667.       if(j5.eq.0)goto 74
  668.       k5=ider(2,j5)
  669.       level=5
  670.       if(k5.eq.ind)goto 50
  671.       l1=ipv(ie)
  672.       l2=ipv(k1)
  673.       l3=ipv(k2)
  674.       l4=ipv(k3)
  675.       l5=ipv(k4)
  676.       l6=ipv(k5)
  677.       bl=' '
  678.       write(8,49)(pv(i,ie),i=1,l1),bl,(pv(i,k1),i=1,l2),bl
  679.      +  ,(pv(i,k2),i=1,l3),bl,(pv(i,k3),i=1,l4),bl,(pv(i,k4),i=1,l5)
  680.      +  ,bl,(pv(i,k5),i=1,l6)
  681.    49 format(' chain too long: ',90a1)
  682.       call bomb(3,'deriv4  ')
  683.    50 call incnt
  684.       if=itetop
  685.       itb(nt)=if+1
  686.       call mulrat(itn(it),itd(it),ivpic,id1,itn(nt),itd(nt))
  687.       do 55 i=ia,ib
  688.       if(i.eq.ic)goto 52
  689.       if=if+1
  690.       ivn(if)=ivn(i)
  691.       ivp(if)=ivp(i)
  692.       goto 55
  693.    52 if(ivp(i).eq.1)goto 55
  694.       if=if+1
  695.       ivn(if)=ivn(i)
  696.       ivp(if)=ivp(i)-1
  697.    55 continue
  698.       call deriv2(j1)
  699.       if(level.eq.1)goto 77
  700.       call deriv2(j2)
  701.       if(level.eq.2)goto 76
  702.       call deriv2(j3)
  703.       if(level.eq.3)goto 75
  704.       call deriv2(j4)
  705.       if(level.eq.4)goto 74
  706.       call deriv2(j5)
  707.       i5=i5+1
  708.       if(i5.le.ndind)goto 45
  709.    74 i4=i4+1
  710.       if(i4.le.ndind)goto 44
  711.    75 i3=i3+1
  712.       if(i3.le.ndind)goto 43
  713.    76 i2=i2+1
  714.       if(i2.le.ndind)goto 42
  715.    77 i1=i1+1
  716.       if(i1.le.ndind)goto 41
  717.    85 continue
  718.    90 continue
  719.       call defsum(jans,nt1,nt)
  720.       call order(jans)
  721.       return
  722.       end
  723.  
  724.  
  725.       subroutine deriv2(j)
  726. $include: 'picom.for'
  727.       common/i/if
  728.       if=if+1
  729.       call nite(if)
  730.       ivn(if)=ider(3,j)
  731.       ivp(if)=1
  732.       return
  733.       end
  734.  
  735.  
  736.       subroutine dertab
  737. $include: 'picom.for'
  738.       idflg=0
  739.       do 302 i=1,ndind
  740.       do 302 j=1,nvar
  741.   302 idtab(i,j)=0
  742.       do 304 i=1,nd
  743.       idepv=ider(1,i)
  744.       do 303 j=1,ndind
  745.       if(idtab(j,idepv).eq.0)goto 304
  746.   303 continue
  747.       call bomb(4,'NDIND   ')
  748.   304 idtab(j,idepv)=i
  749.       return
  750.       end
  751.  
  752.  
  753.       subroutine dfine
  754. $include: 'picom.for'
  755.       call nxt(4,nchpl,' ',j1)
  756.       call nxtch(j1,nchpl,j2)
  757.       call numv(j2,nchpl,j4,j5,j6)
  758.       ittb6=nt+1
  759.       call srch(j5,nchpl,'=',j7)
  760.       if(j7.eq.-1)goto 410
  761.       call term(j7+1)
  762.       goto 412
  763.   410 call read(ib)
  764.       if(np.ge.2)call prntap
  765.     1 format(1x,127a1)
  766.       if(ib.eq.1)goto 412
  767.       call term(1)
  768.       goto 410
  769.   412 call defsum(1,ittb6,nt)
  770.       call order(1)
  771.       call outqu(1)
  772.       call swap(j6)
  773.       return
  774.       end
  775.  
  776.  
  777.       subroutine dump
  778. $include: 'picom.for'
  779.       character*1 bl
  780.       data bl/' '/
  781.       write(8,1)nd,nv,ns,np,nt,line,maxe,maxnt,maxite
  782.     1 format(/,' nd=',i3,' nv=',i3,' ns=',i3,' np=',i2,' nt=',i5,
  783.      +     ' line=',i3,' maxe=',i4,' maxnt=',i5,' maxite=',i6)
  784.       do 10 i=1,nv
  785.       j=ipv(i)
  786.       if(j.eq.nch)goto 3
  787.       jp=j+1
  788.       do 4 jj=jp,nch
  789.     4 pv(jj,i)=bl
  790.     3 write(8,5)i,(pv(k,i),k=1,nch),isb(i),ise(i),ios(i),iosi(i)
  791.     5 format(' i=',i5,'   pv=',15a1,'   isb=',i5,'   ise=',i5
  792.      +      ,'  ios=',i3,'  iosi=',i3)
  793.    10 continue
  794.       write(8,11)
  795.    11 format(' ')
  796.       write(8,13)
  797.    13 format('        i  itn  itd  itb  ite')
  798.       do 20 i=1,nt,4
  799.       j=i+1
  800.       k=i+2
  801.       l=i+3
  802.       write(8,21)i,itn(i),itd(i),itb(i),ite(i),j,itn(j),itd(j),itb(j),
  803.      +ite(j),k,itn(k),itd(k),itb(k),ite(k),l,itn(l),itd(l),itb(l),ite(l)
  804.    21 format(1x,4(i8,4i5))
  805.    20 continue
  806.       write(8,11)
  807.       do 22 i=1,nd
  808.       j=ider(3,i)
  809.       l=ipv(j)
  810.    22 write(8,24)ider(1,i),ider(2,i),(pv(m,j),m=1,l)
  811.    24 format(' the derivative of',i4,' wrt',i4,' is ', 15a1)
  812.       write(8,11)
  813.       write(8,40)(i,ivn(i),ivp(i),i=1,itetop)
  814.    40 format(10(1x,i5,2i3))
  815.       return
  816.       end
  817.  
  818.  
  819.       subroutine expnd(jb,jp,jans)
  820. $include: 'picom.for'
  821.       integer*4 kn,kd,ij,ik,jprvn,jprvd
  822.       ia=isb(jb)
  823.       if(ia.eq.0)call bomb(3,'expnd1  ')
  824.       if(itn(ia).ne.1 .or. itd(ia).ne.1)call bomb(3,'expnd2  ')
  825.       ic=itb(ia)
  826.       if(ivn(ic).ne.2 .or. ivp(ic).ne.0)call bomb(3,'expnd3  ')
  827.       ic=ia+1
  828.       id=ise(jb)
  829.     5 call copyt(ia,ia,1,1)
  830.       if(ic.le.id)goto 10
  831.     6 call defsum(jans,nt,nt)
  832.       return
  833.    10 ig=nt
  834.       ie=isb(jp)
  835.       kn=itn(ie)
  836.       kd=itd(ie)
  837.       if(kn.eq.0)goto 6
  838.       jprvb=ig
  839.       jprve=ig
  840.       jprvn=1
  841.       jprvd=1
  842.       do 20 kp=1,nb
  843.       if(kn.eq.0)goto 25
  844.       call mulrat(jprvn,jprvd,kn,kd*kp,jprvn,jprvd)
  845.       jtn(kp)=jprvn
  846.       jtd(kp)=jprvd
  847.       kn=kn-kd
  848.       jtb(kp)=nt+1
  849.       call mult(jprvb,jprve,ic,id)
  850.       jte(kp)=nt
  851.       if(jtb(kp).gt.nt)goto 25
  852.       jprvb=jtb(kp)
  853.    20 jprve=nt
  854.       kp=nb+1
  855.    25 kp=kp-1
  856.       do 40 i=1,kp
  857.       ih=jtb(i)
  858.       ii=jte(i)
  859.       ij=jtn(i)
  860.       ik=jtd(i)
  861.       do 30 il=ih,ii
  862.    30 call mulrat(itn(il),itd(il),ij,ik,itn(il),itd(il))
  863.    40 continue
  864.       call defsum(jans,ig,ii)
  865.       call order(jans)
  866.       return
  867.       end
  868.  
  869.  
  870.       subroutine extrct
  871. $include: 'picom.for'
  872.       call nxt(2,nchpl,' ',j1)
  873.       call srch(j1,nchpl,'-',j6)
  874.       if(j6.ne.-1)goto 300
  875.       num=1
  876.       do 1 i=2,10
  877.       if(a(j1-1).eq.hn(i))num=i-1
  878.     1 continue
  879.       call numv(j1,nchpl,j3,j4,j5)
  880.       call nxt(j4,nchpl,'=',j6)
  881.       call nxt(j6,nchpl,'m',j7)
  882.       call numv(j7+1,nchpl,j8,j9,j10)
  883.       k=isb(j10)
  884.       if(k.eq.0)call bomb(3,'extrct1 ')
  885.       kk=itn(k)
  886.       if(kk.gt.0 .and. itd(k).eq.1)goto 15
  887.    12 write(8,14)itn(k),itd(k),(a(i),i=j8,j9)
  888.    14 format(' bad term number  ',2i12,5x, 15a1)
  889.       call bomb(3,'extrct2 ')
  890.    15 call nxt(j9+1,nchpl,'o',j11)
  891.       call numv(j11+2,nchpl,j12,j13,j14)
  892.       kl=ise(j14)+1-isb(j14)
  893.       if(kk.gt.kl)goto 12
  894.       km=isb(j14)+kk-1
  895.       nt1=nt+1
  896.       ktop=min0(km+num-1,ise(j14))
  897.       call copyt(km,ktop,1,1)
  898.       call defsum(j5,nt1,nt)
  899.       call outqu(j5)
  900.       return
  901.   300 call numv(j1,nchpl,j2,j3,j4)
  902.       call nxt(j3,j6,'=',j5)
  903.       call nxt(j5,j6,'t',j9)
  904.       call nxt(j9,j6,' ',j10)
  905.       call nxtnb(j10,j6,j11)
  906.       call srch(j11,j6,' ',j12)
  907.       if(j12.eq.-1)j12=j6
  908.       call getint(j11,j12-1,j13,ii4)
  909.       call srchnb(j6+1,nchpl,j8)
  910.       call srch(j8,nchpl,' ',j14)
  911.       call getint(j8,j14-1,j15,ii4)
  912.       call nxt(j14,nchpl,'o',j19)
  913.       if(a(j19+1).ne.'f' .and. a(j19+2).ne.' ')call bomb(3,'extrct5 ')
  914.   330 call numv(j19+2,nchpl,j16,j17,j18)
  915.       if(isb(j18).eq.0)call bomb(3,'extrct6 ')
  916.       if(j15.gt.ise(j18)-isb(j18)+1)j15=ise(j18)-isb(j18)+1
  917.       if(j15.lt.j13)call bomb(3,'extrct7 ')
  918.       nt1=nt+1
  919.       call copyt(isb(j18)+j13-1,isb(j18)+j15-1,1,1)
  920.       call defsum(j4,nt1,nt)
  921.       call outqu(j4)
  922.       return
  923.       end
  924.  
  925.  
  926.       subroutine fort(it,ib)
  927. $include: 'picom.for'
  928.       character*1 pa,aa
  929.       character*1 scr
  930.       common/pr1/pa(232),aa(22)
  931.       dimension dum(4)
  932.       if(np.eq.0)return
  933.       if(itd(it).eq.1 .and. itn(it).eq.-1)then
  934.          ib=ib+1
  935.          a(ib)='-'
  936.          goto 42
  937.          endif
  938.       if(itd(it).eq.1 .and. itn(it).eq. 1)then
  939.          ib=ib+1
  940.          a(ib)='+'
  941.          goto 42
  942.          endif
  943.       write(scr,43)itn(it)
  944.       read(scr,9)(aa(kk),kk=1,20)
  945.     9 format(20a1)
  946.       j=0
  947.       do 10 i=1,20
  948.       if(aa(i).eq.' ')goto 10
  949.       if(j.eq.0 .and. aa(i).ne.'-')then
  950.          ib=ib+1
  951.          a(ib)='+'
  952.          endif
  953.       j=1
  954.       ib=ib+1
  955.       a(ib)=aa(i)
  956.    10 continue
  957.       ib=ib+1
  958.       a(ib)='.'
  959.       if(itd(it).eq.1)goto 42
  960.       ib=ib+1
  961.       a(ib)='/'
  962.       write(scr,43)itd(it)
  963.       read(scr,9)(aa(kk),kk=1,20)
  964.       do 15 i=1,20
  965.       if(aa(i).eq.' ')goto 15
  966.       ib=ib+1
  967.       a(ib)=aa(i)
  968.    15 continue
  969.    42 ic=itb(it)
  970.       id=ite(it)
  971.       do 60 iv=ic,id
  972.       ip=iabs(ivp(iv))
  973.       in=ivn(iv)
  974.       ie=ipv(in)
  975.       iq=0
  976.       if(ip.eq.1)goto 47
  977.       write(scr,43)ip
  978.       read(scr,9)(aa(kk),kk=1,20)
  979.    43 format(i20)
  980.       il=20
  981.       do 44 if=1,20
  982.       if(aa(if).ne.' ')goto 45
  983.    44 continue
  984.    45 iq=il+1-if+2
  985.    47 if(iq+ie+1.le.72)goto 52
  986.       do 48 j=75,1,-1
  987.       if(a(j).ne.' ')goto 49
  988.    48 continue
  989.    49 write(8,50)(a(i),i=1,j)
  990.    50 format(80a1)
  991.       do 51 i=1,80
  992.    51 a(i)=' '
  993.       a(6)='.'
  994.       ib=11
  995.    52 ib=ib+1
  996.       if(a(ib-1).ne.'+' .and. a(ib-1).ne.'-')a(ib)='*'
  997.       if(ivp(iv).lt.0)a(ib)='/'
  998.       do 53 i=1,ie
  999.       ib=ib+1
  1000.    53 a(ib)=pv(i,in)
  1001.       if(iq.eq.0)goto 60
  1002.       ib=ib+1
  1003.       a(ib)='*'
  1004.       ib=ib+1
  1005.       a(ib)='*'
  1006.       do 54 i=if,il
  1007.       ib=ib+1
  1008.    54 a(ib)=aa(i)
  1009.    60 continue
  1010.       do 62 j=75,1,-1
  1011.       if(a(j).ne.' ')goto 63
  1012.    62 continue
  1013.    63 write(8,50)(a(i),i=1,j)
  1014.       return
  1015.       end
  1016.  
  1017.  
  1018.       subroutine fortrn
  1019. $include: 'picom.for'
  1020.       call nxt(1,nchpl,' ',j1)
  1021.       call numv(j1+1,nchpl,j2,j3,j4)
  1022.       if(np.eq.0 .or. isb(j4).eq.0)goto 10
  1023.       do 1 i=1,nchpl
  1024.     1 a(i)=' '
  1025.       ia=ipv(j4)
  1026.       ib=6
  1027.       do 2 i=1,ia
  1028.       ib=ib+1
  1029.     2 a(ib)=pv(i,j4)
  1030.       ib=ib+1
  1031.       a(ib)='='
  1032.       ic=isb(j4)
  1033.       id=ise(j4)
  1034.       j=0
  1035.       do 4 it=ic,id
  1036.       call fort(it,ib)
  1037.       do 3 i=1,nchpl
  1038.     3 a(i)=' '
  1039.       a(6)='.'
  1040.       ib=7
  1041.       j=j+1
  1042.       if(j.lt.10 .or. it.eq.id)goto 4
  1043.       a(6)=' '
  1044.       call copyc(pv(1,j4),a(7),ia)
  1045.       a(7+ia)='='
  1046.       call copyc(pv(1,j4),a(8+ia),ia)
  1047.       ib=7+ia+ia
  1048.       j=0
  1049.     4 continue
  1050.    10 return
  1051.       end
  1052.  
  1053.  
  1054.       subroutine gcd(inn,idd,ig)
  1055.       implicit integer*4 (i-n)
  1056.       ig=1
  1057.       in=iabs(inn)
  1058.       id=idd
  1059.       if(in.lt.id)goto 10
  1060.     5 in=mod(in,id)
  1061.       if(in.eq.0)goto 20
  1062.    10 if(in.eq.1)return
  1063.       id=mod(id,in)
  1064.       if(id.eq.0)goto 15
  1065.       if(id.eq.1)return
  1066.       goto 5
  1067.    15 ig=in
  1068.       return
  1069.    20 ig=id
  1070.       return
  1071.       end
  1072.  
  1073.  
  1074.       subroutine gcdout(in,id)
  1075.       implicit integer*4 (i-n)
  1076.       if(id.le.0)call bomb(3,'gcdout1 ')
  1077.       if(in.eq.1 .or. in.eq.-1 .or. id.eq.1)return
  1078.       if(in.eq.0)goto 1
  1079.       call gcd(in,id,ig)
  1080.       if(ig.eq.1)return
  1081.       in=in/ig
  1082.       id=id/ig
  1083.       return
  1084.     1 id=1
  1085.       return
  1086.       end
  1087.  
  1088.  
  1089.       subroutine getint(j1,j2,ii2,j3)
  1090. c  on return j3=integer value of a(j1)-a(j2)
  1091.       implicit integer*2 (i-n)
  1092.       integer*4 j3
  1093.       j3=0
  1094.       if(j2.lt.j1)call bomb(3,' getint1')
  1095.       do 1 i=j1,j2
  1096.       call intval(i,j)
  1097.     1 j3=10*j3+j
  1098.       ii2=j3
  1099.       return
  1100.       end
  1101.  
  1102.  
  1103.       subroutine getexp(j1,ie,k,if)
  1104. c  begin looking at a(j1), exponent put in ie, a(k) is next nonblank
  1105. c     (or k=-1 if only blanks are found)
  1106. $include: 'picom.for'
  1107.       ie=1
  1108.       call srchnb(j1,nchpl,k)
  1109.       if(k.eq.-1)return
  1110.       if(a(k).eq.'*' .and. a(k+1).ne.'*')return
  1111.       isgn=1
  1112.       if(a(k).eq.'*' .and. a(k+1).eq.'*')goto 35
  1113.       if(a(k).eq.'/')return
  1114.       call srchch(j1,nchpl,j2)
  1115.       if(j2.eq.k .and. if.ne.0)return
  1116.       if(j2.ne.k)goto 9
  1117.    32 call numv(j2,nchpl,j3,j4,j5)
  1118.       if(isb(j5).eq.0)call bomb(3,'getexp10')
  1119.       ie=itn(isb(j5))
  1120.       if(itd(isb(j5)).ne.1)call bomb(3,'getexp11')
  1121.       k=-1
  1122.       if(j4.ge.nchpl-1)return
  1123.       call srchnb(j4+1,nchpl,k)
  1124.       return
  1125.    9  if(a(k).eq.'+')goto 20
  1126.       if(a(k).eq.'-')goto 10
  1127.       if(a(k).eq.'^')goto 30
  1128.       j3=k
  1129.     1 me=0
  1130.     5 call intval(j3,j4)
  1131.       me=10*me+j4
  1132.       j3=j3+1
  1133.       if(j3.gt.nchpl)goto 70
  1134.       if(a(j3).eq.' ')goto 50
  1135.       if(a(j3).eq.'*')goto 60
  1136.       if(a(j3).eq.'/')goto 60
  1137.       goto 5
  1138.    10 isgn=-1
  1139.    20 if(k+1.gt.nchpl)call bomb(3,' getexp2')
  1140.       call nxtnb(k+1,nchpl,j3)
  1141.       if(j3.eq.-1)call bomb(3,' getexp4')
  1142.       goto 1
  1143.    30 call nxtnb(k+1,nchpl,j3)
  1144.    31 k=j3
  1145.       j2=j3
  1146.       call srchch(j3,nchpl,j4)
  1147.       if(j3.eq.j4)goto 32
  1148.       if(a(k).eq.'-')goto 10
  1149.       if(a(k).eq.'+')goto 20
  1150.       goto 1
  1151.    35 if(k.ge.nchpl-1)call bomb(3,' getexp6')
  1152.       call nxtnb(k+2,nchpl,j3)
  1153.       if(j3.eq.-1)call bomb(3,' getexp8')
  1154.       goto 31
  1155.    40 k=k+1
  1156.       goto 30
  1157.    50 ie=me*isgn
  1158.       call srchnb(j3,nchpl,k)
  1159.       return
  1160.    60 k=j3
  1161.       if(j3.eq.nchpl)call bomb(3,' getexp9')
  1162.       ie=me*isgn
  1163.       return
  1164.    70 ie=me*isgn
  1165.       k=-1
  1166.       return
  1167.       end
  1168.  
  1169.  
  1170.       function icomp(ia,ib)
  1171. $include: 'picom.for'
  1172. c  if ia should precede ib icomp=1
  1173. c  if ib should precede ia icomp=-1
  1174. c  if terms are the same except possibly coefficients, icomp=0
  1175.       ja=itb(ia)
  1176.       jb=itb(ib)
  1177.       ka=ite(ia)
  1178.       kb=ite(ib)
  1179.     1 if(ivn(ja)-ivn(jb))2,8,9
  1180.     2 if(ivp(ja))3,10,4
  1181.     4 icomp=-1
  1182.       return
  1183.     5 if(jb.gt.kb)goto 7
  1184.     3 icomp=1
  1185.       return
  1186.     6 ja=ja+1
  1187.    11 jb=jb+1
  1188.       if(ja.gt.ka)goto 5
  1189.       if(jb.gt.kb)goto 4
  1190.       goto 1
  1191.     7 icomp=0
  1192.       return
  1193.     8 if(ivp(ja)-ivp(jb))3,6,4
  1194.     9 if(ivp(jb))4,11,3
  1195.    10 ja=ja+1
  1196.       if(ja.gt.ka)goto 3
  1197.       goto 1
  1198.       end
  1199.  
  1200.  
  1201.       subroutine iftest(igoflg)
  1202. $include: 'picom.for'
  1203.       integer*4 i5,i8
  1204.       character*1 aj1,aj2
  1205.       igoflg=0
  1206.       call nxt(2,nchpl,' ',j1)
  1207.       call nxt(j1,nchpl,'.',j)
  1208.       call numv(j1,j-1,j3,j4,j5)
  1209.       aj1=a(j+1)
  1210.       aj2=a(j+2)
  1211.       if(aj1.eq.'n')goto 20
  1212.       if(aj1.eq.'l' .or. aj1.eq.'g')goto 22
  1213.       if(aj1.ne.'e')call bomb(3,' iftest1')
  1214.       if(aj2.ne.'q')call bomb(3,' iftest2')
  1215.     4 if(a(j+3).ne.'.')call bomb(3,' iftest3')
  1216.       call numv(j+4,nchpl,j6,j7,j8)
  1217.       if(j5.ne.j8)goto 10
  1218.       if(aj1.eq.'n')return
  1219.     5 call nxtch(j7+1,nchpl,j9)
  1220.       do 6 j=j9,nchpl
  1221.       a(j+1-j9)=a(j)
  1222.     6 a(j)=' '
  1223.       igoflg=1
  1224.       call copyc(a,ap,nchpl)
  1225.       return
  1226.    10 k1=itb(isb(j5))
  1227.       k2=ite(ise(j5))
  1228.       k3=itb(isb(j8))
  1229.       k4=ite(ise(j8))
  1230.       i5=itn(isb(j5))*itd(isb(j8))
  1231.       i8=itd(isb(j5))*itn(isb(j8))
  1232.       if(aj1.eq.'l')goto 23
  1233.       if(aj1.eq.'g')goto 25
  1234.       if(k2-k1.ne.k4-k3)goto 30
  1235.       k=k3
  1236.       do 11 i=k1,k2
  1237.       if(ivn(i).ne.ivn(k))goto 30
  1238.       if(ivp(i).ne.ivp(k))goto 30
  1239.    11 k=k+1
  1240.       k1=isb(j5)
  1241.       k2=ise(j5)
  1242.       k3=isb(j8)
  1243.       k4=ise(j8)
  1244.       if(k2-k1.ne.k4-k3)goto 30
  1245.       k=k3
  1246.       do 12 i=k1,k2
  1247.       if(itn(i).ne.itn(k))goto 30
  1248.       if(itd(i).ne.itd(k))goto 30
  1249.    12 k=k+1
  1250.       if(aj1.eq.'n')return
  1251.       goto 5
  1252.    20 if(a(j+2).ne.'e')call bomb(3,' iftest4')
  1253.       goto 4
  1254.    22 if(aj2.eq.'e' .or. aj2.eq.'t')goto 4
  1255.       call bomb(5,' iftest5')
  1256.    23 if(aj2.eq.'e')goto 24
  1257.       if(i5.lt.i8)goto 5
  1258.       return
  1259.    24 if(i5.le.i8)goto 5
  1260.       return
  1261.    25 if(aj2.eq.'e')goto 26
  1262.       if(i5.gt.i8)goto 5
  1263.       return
  1264.    26 if(i5.ge.i8)goto 5
  1265.       return
  1266.    30 if(aj1.eq.'n')goto 5
  1267.       return
  1268.       end
  1269.  
  1270.  
  1271.       subroutine incnt
  1272. $include: 'picom.for'
  1273.       nt=nt+1
  1274.       maxnt=max0(nt,maxnt)
  1275.       if(maxnt.gt.nterm)call bomb(4,'NTERM   ')
  1276.       return
  1277.       end
  1278.  
  1279.  
  1280.       subroutine intval(j1,i)
  1281. c  on return i=integer value of a(j1)
  1282. $include: 'picom.for'
  1283.       do 1 j=1,10
  1284.       if(a(j1).eq.hn(j))goto 2
  1285.     1 continue
  1286.       i=-1
  1287.       return
  1288.     2 i=j-1
  1289.       return
  1290.       end
  1291.  
  1292.  
  1293.       subroutine monitr
  1294. $include: 'picom.for'
  1295.       call datime
  1296.       tprev=itime/100.
  1297.       call readin
  1298.       ntp=0
  1299.    10 if(ha(1).ne.'a')call bomb(3,'monitr1 ')
  1300.       if(isb(2).ne.0)call bomb(3,'useorder')
  1301.       if(ha(1).eq.'a')goto 11
  1302.       if(ns.lt.2)goto 14
  1303.       do 20 ms=2,ns
  1304.       is=iosi(ms-1)
  1305.       js=iosi(ms)
  1306.       if(is.eq.0 .or. js.eq.0)call bomb(3,'    mon1')
  1307.       if(ios(is).ne.ms-1 .or. ios(js).ne.ms)call bomb(3,'    mon2')
  1308.       if(isb(is).eq.0 .or. isb(js).eq.0)call bomb(3,'    mon3')
  1309.       ia=ise(is)
  1310.       ib=isb(js)
  1311.       if(ite(ia)+1.ne.itb(ib))call bomb(3,'    mon4')
  1312.       do 16 it=ib,ise(js)
  1313.       if(itb(it).gt.ite(it))call bomb(3,'    mon5')
  1314.       if(itd(it).le.0)call bomb(3,'   mon22')
  1315.       if(itn(it).ne.0)goto 15
  1316.       if(itb(it).ne.ite(it))call bomb(3,'    mon6')
  1317.       if(ivn(itb(it)).ne.2)call bomb(3,'    mon7')
  1318.       if(itb(it).gt.ite(it))call bomb(3,'    mon8')
  1319.    15 do 17 i=itb(it),ite(it)
  1320.       if(ivn(i).eq.1)call bomb(3,'    mon9')
  1321.       if(isb(ivn(i)).ne.0)call bomb(3,'   mon10')
  1322.       if(i.eq.itb(it))goto 17
  1323.       if(ivp(i).eq.0)call bomb(3,'   mon11')
  1324.       if(ivn(i-1).ge.ivn(i))call bomb(3,'   mon12')
  1325.    17 continue
  1326.    16 continue
  1327.    20 continue
  1328.    14 if(ns.eq.0)call bomb(3,'   mon13')
  1329.       ic=iosi(ns)
  1330.       if(ise(ic).ne.nt)call bomb(3,'   mon14')
  1331.       if(ite(ise(ic)).ne.itetop)call bomb(3,'   mon15')
  1332.    11 call read(ib)
  1333.       if(ib.eq.1 .and. np.ge.1)call prntap
  1334.     1 format(1x,127a1)
  1335.       if(ib.eq.1)goto 10
  1336.       call datime
  1337.       ttot=itime/100.
  1338.       t=ttot-tprev
  1339.       tprev=ttot
  1340.       ntmntp=nt-ntp
  1341.       if(np.ge.3)write(8,3)t,ttot,ntmntp,nt,itetop,maxnt,maxite
  1342.     3 format(10x,'t=',f7.2,' total t=',f9.1,' ntmntp=',i5,' nt=',i5,
  1343.      +       ' itetop=',i6,' maxnt=',i5,' maxite=',i6)
  1344.       if(np.ge.2)call prntap
  1345.     5 ntp=nt
  1346.       call copyc(a,ap,nchpl)
  1347.       if(a(1).eq.'c')goto 300
  1348.       if(a(1).eq.'d')goto 400
  1349.       if(a(1).eq.'e')goto 500
  1350.       if(a(1).eq.'f')goto 600
  1351.       if(a(1).eq.'g')goto 700
  1352.       if(a(1).eq.'i')goto 900
  1353.       if(a(1).eq.'l')goto 10
  1354.       if(a(1).eq.'m')goto 1300
  1355.       if(a(1).eq.'n')goto 1400
  1356.       if(a(1).eq.'o')goto 1500
  1357.       if(a(1).eq.'p')goto 1600
  1358.       if(a(1).eq.'r')goto 1800
  1359.       if(a(1).eq.'s')goto 1900
  1360.       if(a(1).eq.'t')goto 2000
  1361.       if(a(1).eq.'z')goto 2600
  1362.       call bomb(3,'bad cmnd')
  1363.   300 if(a(2).eq.'o')goto 310
  1364.       if(a(2).ne.'a')call bomb(3,' monitr2')
  1365.       call call
  1366.       goto 10
  1367.   310 if(a(3).eq.'m')goto 320
  1368.       if(a(3).ne.'u')call bomb(3,' monitr3')
  1369.       call count
  1370.       goto 10
  1371.   320 call comput
  1372.       call pack
  1373.       goto 10
  1374.   400 if(a(2).eq.'u')goto 470
  1375.       if(a(2).eq.'i')goto 320
  1376.       if(a(2).ne.'e')call bomb(3,' monitr4')
  1377.       if(a(3).eq.'l')goto 450
  1378.       if(a(3).eq.'f')goto 440
  1379.       call bomb(3,' monitr5')
  1380.   440 call dfine
  1381.       goto 10
  1382.   450 call delet
  1383.       call pack
  1384.       goto 10
  1385.   470 call dump
  1386.       goto 10
  1387.   500 if(a(2).eq.'n')return
  1388.       if(a(2).ne.'x')call bomb(3,' monitr6')
  1389.       call extrct
  1390.       goto 10
  1391.   600 if(a(2).eq.'l')goto 610
  1392.       call fortrn
  1393.       goto 10
  1394.   610 call readin
  1395.       goto 10
  1396.   700 call nxt(1,nchpl,' ',j9)
  1397.       call numl(j9,nchpl,j10,j11,j12)
  1398.       if(linlab(j12).eq.-1)write(8,710)
  1399.       if(linlab(j12).eq.-1)call bomb(3,' monitr7')
  1400.   710 format(' label not found')
  1401.       line=linlab(j12)-1
  1402.       goto 10
  1403.   900 call iftest(igoflg)
  1404.       if(igoflg.eq.1)goto 5
  1405.       goto 10
  1406.  1300 call srch(1,nchpl,'=',j)
  1407.       call srchnb(j+1,nchpl,k)
  1408.       call getexp(k,maxe,j,0)
  1409.       goto 10
  1410.  1400 call nxt(3,nchpl,'=',j1)
  1411.       call getexp(j1+1,j3,j2,0)
  1412.       if(a(2).eq.'d')goto 1405
  1413.       if(a(2).eq.'p')goto 1410
  1414.       if(a(2).eq.'b')goto 1415
  1415.       call bomb(3,' monitr8')
  1416.  1405 if(a(3).eq.'u')goto 1406
  1417.       call bomb(3,' monitr9')
  1418.  1406 ndump=j3
  1419.       goto 10
  1420.  1410 np=j3
  1421.       goto 10
  1422.  1415 nb=j3
  1423.       goto 10
  1424.  1500 if(nv.ne.1)call bomb(3,'putfirst')
  1425.       call ordset
  1426.       goto 10
  1427.  1600 if(a(3).eq.'i')goto 1620
  1428.       if(a(3).ne.'o')call bomb(3,'monitr10')
  1429.  1605 call read(ib)
  1430.       if(np.ge.2)call prntap
  1431.       if(a(1).ne.'r' .or. a(2).ne.'e' .or. a(3).ne.'t')goto 1605
  1432.       goto 10
  1433.  1620 call nxt(4,nchpl,' ',j1)
  1434.       call numv(j1,nchpl,j2,j3,j4)
  1435.       if(np.eq.1)write(8,1)
  1436.       if(np.eq.1)call prntap
  1437.       call print(j4)
  1438.       goto 10
  1439.  1800 if(a(2).ne.'e' .or. a(3).ne.'t')call bomb(3,' return?')
  1440.       line=linret
  1441.       goto 10
  1442.  1900 call sub
  1443.       call pack
  1444.       goto 10
  1445.  2000 if(a(2).eq.'a')goto 2010
  1446.       call thderv
  1447.       goto 10
  1448.  2010 call taylor
  1449.       goto 10
  1450.  2600 call numv(2,nchpl,j1,j2,j3)
  1451.       call zsub(j3)
  1452.       call pack
  1453.       goto 10
  1454.       end
  1455.  
  1456.  
  1457.       subroutine mulrat(n1,i1,n2,i2,n3,i3)
  1458.       implicit integer*4 (i-n)
  1459.       n3=n1*n2
  1460.       i3=i1*i2
  1461.       call gcdout(n3,i3)
  1462.       return
  1463.       end
  1464.  
  1465.  
  1466.       subroutine mult(ia,ib,ic,id)
  1467. c  multiply terms ia-ib by terms ic-id
  1468. c  result goes to nt+1,.....
  1469. $include: 'picom.for'
  1470.       common/nwt/iitb,iite,jitb,jite,ntb,if
  1471.       integer*4 inp,idp
  1472.       if(ia.eq.0 .or. ic.eq.0)call bomb(3,'   mult1')
  1473.       ntb=nt+1
  1474.       if=0
  1475.       do 20 i=ia,ib
  1476.       iitb=itb(i)
  1477.       iite=ite(i)
  1478.       do 10 j=ic,id
  1479.       jitb=itb(j)
  1480.       jite=ite(j)
  1481.       if(ivn(iitb)+ivn(jitb).eq.4 .and. ivp(iitb)+ivp(jitb).gt.maxe)
  1482.      +                                                       goto 20
  1483.       call mulrat(itn(i),itd(i),itn(j),itd(j),inp,idp)
  1484.       call sub2(inp,idp)
  1485.    10 continue
  1486.    20 continue
  1487.       return
  1488.       end
  1489.  
  1490.  
  1491.       subroutine multn(i,j,k,l)
  1492. $include: 'picom.for'
  1493.       nt1=nt+1
  1494.       call mult(i,j,k,l)
  1495.       if(nt.ge.nt1)return
  1496.       call niltrm
  1497.       return
  1498.       end
  1499.  
  1500.  
  1501.       subroutine niltrm
  1502. $include: 'picom.for'
  1503.       k=itetop+1
  1504.       call incnt
  1505.       itn(nt)=0
  1506.       itd(nt)=1
  1507.       itb(nt)=k
  1508.       call nite(k)
  1509.       ivn(k)=2
  1510.       ivp(k)=0
  1511.       return
  1512.       end
  1513.  
  1514.  
  1515.       subroutine nite(j)
  1516. $include: 'picom.for'
  1517.       itetop=j
  1518.       ite(nt)=j
  1519.       maxite=max0(maxite,j)
  1520.       if(maxite.gt.ntot)call bomb(4,'NTOT    ')
  1521.       return
  1522.       end
  1523.  
  1524.  
  1525.  
  1526.       subroutine numl(j1,j2,j3,j4,j5)
  1527. c  search j1-j2  find name in j3-j4, j5 is index of label
  1528. $include: 'picom.for'
  1529.       character*1 c
  1530.       if(j2.lt.j1)call bomb(3,'    numl')
  1531.       call nxtch(j1,j2,j3)
  1532.       do 10 j=j3,j2
  1533.       c=a(j)
  1534.       do 5 k=1,26
  1535.       if(c.eq.ha(k) .or. c.eq.hc(k))goto 10
  1536.     5 continue
  1537.       do 6 k=1,10
  1538.       if(c.eq.hn(k))goto 10
  1539.     6 continue
  1540.       if(c.eq.' ' .or. c.eq.'+' .or. c.eq.'-')goto 12
  1541.       if(c.eq.'*' .or. c.eq.'/' .or. c.eq.'=' .or. c.eq.'^')goto 12
  1542.    10 continue
  1543.       j4=j2
  1544.       goto 14
  1545.    12 j4=j-1
  1546.    14 l=j4-j3+1
  1547.       do 18 i=1,labels
  1548.       if(ipl(i).ne.l)goto 18
  1549.       do 16 j=1,l
  1550.       if(pl(j,i).ne.a(j+j3-1))goto 18
  1551.    16 continue
  1552.       j5=i
  1553.       return
  1554.    18 continue
  1555.       labels=labels+1
  1556.       if(labels.gt.nlab)call bomb(4,'NLAB    ')
  1557.       j5=labels
  1558.       ipl(j5)=l
  1559.       call copyc(a(j3),pl(1,j5),l)
  1560.       return
  1561.       end
  1562.  
  1563.  
  1564.       subroutine numv(j1,j2,j3,j4,j5)
  1565. c  search j1-j2  find name in j3-j4, j5 is index of variable
  1566. $include: 'picom.for'
  1567.       character*1 c
  1568.       if(j2.lt.j1)call bomb(3,'    numv')
  1569.       call nxtch(j1,j2,j3)
  1570.       do 10 j=j3,j2
  1571.       c=a(j)
  1572.       do 5 k=1,26
  1573.       if(c.eq.ha(k) .or. c.eq.hc(k))goto 10
  1574.     5 continue
  1575.       do 6 k=1,10
  1576.       if(c.eq.hn(k))goto 10
  1577.     6 continue
  1578.       if(c.eq.' ' .or. c.eq.'+' .or. c.eq.'-')goto 12
  1579.       if(c.eq.'*' .or. c.eq.'/' .or. c.eq.'=' .or. c.eq.'^')goto 12
  1580.    10 continue
  1581.       j4=j2
  1582.       goto 14
  1583.    12 j4=j-1
  1584.    14 l=j4-j3+1
  1585.       do 28 ii=ns,1,-1
  1586.       i=iosi(ii)
  1587.       if(ipv(i).ne.l)goto 28
  1588.       do 26 j=1,l
  1589.       if(pv(j,i).ne.a(j+j3-1))goto 28
  1590.    26 continue
  1591.       j5=i
  1592.       return
  1593.    28 continue
  1594.       do 18 i=1,nv
  1595.       if(isb(i).ne.0)goto 18
  1596.       if(ipv(i).ne.l)goto 18
  1597.       do 16 j=1,l
  1598.       if(pv(j,i).ne.a(j+j3-1))goto 18
  1599.    16 continue
  1600.       j5=i
  1601.       return
  1602.    18 continue
  1603.       nv=nv+1
  1604.       j5=nv
  1605.       if(nv.gt.nvar)call bomb(4,'NVAR    ')
  1606.       ipv(nv)=l
  1607.       call copyc(a(j3),pv(1,nv),l)
  1608.       return
  1609.       end
  1610.  
  1611.  
  1612.       subroutine nxt(j1,j2,h,j3)
  1613.       implicit integer*2 (i-n)
  1614.       character*1 h
  1615.       call srch(j1,j2,h,j3)
  1616.       if(j3.eq.-1)call bomb(3,'     nxt')
  1617.       end
  1618.  
  1619.  
  1620.       subroutine nxtch(j1,j2,j3)
  1621.       implicit integer*2 (i-n)
  1622.       call srchch(j1,j2,j3)
  1623.       if(j3.eq.-1)call bomb(3,'   nxtch')
  1624.       return
  1625.       end
  1626.  
  1627.  
  1628.       subroutine nxtnb(j1,j2,j3)
  1629.       implicit integer*2 (i-n)
  1630.       call srchnb(j1,j2,j3)
  1631.       if(j3.eq.-1)call bomb(3,'   nxtnb')
  1632.       end
  1633.  
  1634.  
  1635.       subroutine nxtnnu(j1,j2,j3)
  1636. $include: 'picom.for'
  1637.       do 10 i=j1,j2
  1638.       do 5 j=1,10
  1639.       if(a(i).eq.hn(j))goto 10
  1640.     5 continue
  1641.       j3=i
  1642.       return
  1643.    10 continue
  1644.       j3=-1
  1645.       return
  1646.       end
  1647.  
  1648.  
  1649.       subroutine order(is)
  1650. $include: 'picom.for'
  1651.       ina=isb(is)
  1652.       inb=ise(is)
  1653.       if(inb.lt.ina .or. ina.eq.0)call bomb(3,'  order1')
  1654.       it=ina
  1655.     2 ia=itb(it)
  1656.       if(itn(it).eq.0)goto 41
  1657.       ib=ite(it)
  1658.     3 if(ib.eq.ia)goto 25
  1659.       i=ia
  1660.     4 if(ivp(i).eq.0)goto 20
  1661.       i=i+1
  1662.       if(i.le.ib)goto 4
  1663.       if(ia.eq.ib)goto 25
  1664.    23 i=ia+1
  1665.     5 if(ivn(i)-ivn(i-1))10,7,6
  1666.     6 i=i+1
  1667.       if(i.le.ib)goto 5
  1668.       goto 40
  1669.     7 ivp(i-1)=ivp(i-1)+ivp(i)
  1670.       if(i.eq.ib)goto 9
  1671.       ib=ib-1
  1672.       do 8 j=i,ib
  1673.       ivn(j)=ivn(j+1)
  1674.     8 ivp(j)=ivp(j+1)
  1675.       if(ivp(i-1).ne.0)goto 5
  1676.       goto 3
  1677.     9 ib=ib-1
  1678.       if(ivp(ib).ne.0)goto 40
  1679.       goto 3
  1680.    10 if(i.eq.ia+1)goto 19
  1681.       ic=i-2
  1682.    11 if(ivn(i)-ivn(ic))12,13,16
  1683.    12 if(ic.eq.ia)goto 18
  1684.       ic=ic-1
  1685.       goto 11
  1686.    13 ivp(ic)=ivp(ic)+ivp(i)
  1687.       if(i.eq.ib)goto 15
  1688.       ib=ib-1
  1689.       do 14 j=i,ib
  1690.       ivn(j)=ivn(j+1)
  1691.    14 ivp(j)=ivp(j+1)
  1692.       if(ivp(ic).eq.0)goto 3
  1693.       goto 5
  1694.    15 ib=ib-1
  1695.       if(ivp(ic).ne.0)goto 40
  1696.       goto 3
  1697.    16 ip=ivp(i)
  1698.       in=ivn(i)
  1699.       do 17 j=i,ic+2,-1
  1700.       ivn(j)=ivn(j-1)
  1701.    17 ivp(j)=ivp(j-1)
  1702.       ivn(ic+1)=in
  1703.       ivp(ic+1)=ip
  1704.       goto 6
  1705.    18 ic=ia-1
  1706.       goto 16
  1707.    19 ip=ivp(i)
  1708.       in=ivn(i)
  1709.       ivp(i)=ivp(ia)
  1710.       ivn(i)=ivn(ia)
  1711.       ivp(ia)=ip
  1712.       ivn(ia)=in
  1713.       goto 6
  1714.    20 if(i.eq.ib)goto 22
  1715.       ib=ib-1
  1716.       do 21 j=i,ib
  1717.       ivn(j)=ivn(j+1)
  1718.    21 ivp(j)=ivp(j+1)
  1719.       goto 4
  1720.    22 ib=ib-1
  1721.       if(ib.gt.ia)goto 23
  1722.       ib=ia
  1723.    25 if(ivp(ia).eq.0)ivn(ia)=2
  1724.    40 ite(it)=ib
  1725.       if(ivn(ia).ne.2 .or. ivp(ia).le.maxe)goto 52
  1726.    41 call bkup(it,inb)
  1727.       goto 53
  1728.    52 it=it+1
  1729.    53 if(it.le.inb)goto 2
  1730. c  sort terms ina-inb
  1731.       if(ina.lt.inb)goto 70
  1732.       if(ina.eq.inb)goto 82
  1733.    51 inb=ina
  1734.       itn(ina)=0
  1735.       itd(ina)=1
  1736.       k=itb(ina)
  1737.       ite(ina)=k
  1738.       ivn(k)=2
  1739.       ivp(k)=0
  1740.       goto 82
  1741.    70 it=ina
  1742.    68 it=it+1
  1743.    69 if(it.gt.inb)goto 80
  1744.       if(itn(it).eq.0)goto 75
  1745.       if(icomp(ina,it))71,73,74
  1746.    71 call tzap(it,ina)
  1747.       goto 68
  1748.    73 call addrat(itn(ina),itd(ina),itn(it),itd(it),itn(ina),itd(ina))
  1749.    75 call bkup(it,inb)
  1750.       goto 69
  1751.    74 if(icomp(it-1,it))62,61,68
  1752.    61 call addrat(itn(it-1),itd(it-1),itn(it),itd(it),
  1753.      ,     itn(it-1),itd(it-1))
  1754.       call bkup(it,inb)
  1755.       goto 69
  1756.    62 la=ina
  1757.       lb=it-1
  1758.    63 if(lb-la.eq.1)goto 64
  1759.       lc=(la+lb)/2
  1760.       if(icomp(lc,it))65,66,67
  1761.    64 call tzap(it,lb)
  1762.       goto 68
  1763.    65 lb=lc
  1764.       goto 63
  1765.    66 call addrat(itn(lc),itd(lc),itn(it),itd(it),itn(lc),itd(lc))
  1766.       call bkup(it,inb)
  1767.       goto 69
  1768.    67 la=lc
  1769.       goto 63
  1770.    80 if(ina.gt.inb)goto 51
  1771.    82 call defsum(is,ina,inb)
  1772.       return
  1773.       end
  1774.  
  1775.  
  1776.       subroutine ordset
  1777. $include: 'picom.for'
  1778.       call nxt(2,nchpl,' ',j1)
  1779.       call srchnb(j1,nchpl,j2)
  1780.       if(j2.eq.-1)goto 20
  1781.       ntsav=nt
  1782.       itesav=itetop
  1783.       call term(j1)
  1784.       nt=ntsav
  1785.       itetop=itesav
  1786.       return
  1787.    20 if=1
  1788.       do 10 il=1,300
  1789.       call read(ib)
  1790.       if(np.ge.2)call prntap
  1791.     1 format(1x,127a1)
  1792.       if(ib.eq.1)return
  1793.       is=1
  1794.     2 call srchch(is,nchpl,j1)
  1795.       if(j1.eq.-1)goto 10
  1796.       call numv(j1,nchpl,j2,j3,j4)
  1797.       if=if+1
  1798.       if(j4.lt.if)call bomb(3,'ordset1 ')
  1799.       if(j4.eq.if)goto 8
  1800.       do 5 i=nv,if,-1
  1801.       ipv(i)=ipv(i-1)
  1802.       do 4 ic=1,nch
  1803.     4 pv(ic,i)=pv(ic,i-1)
  1804.     5 continue
  1805.       call copyc(a(j2),pv(1,if),j3-j2+1)
  1806.       ipv(if)=j3-j2+1
  1807.     8 is=j3+1
  1808.       if(is.le.nchpl)goto 2
  1809.    10 continue
  1810.       return
  1811.       end
  1812.  
  1813.  
  1814.       subroutine outqu(j5)
  1815. $include: 'picom.for'
  1816.     1 do 6 it=isb(j5),ise(j5)
  1817.       do 5 iv=itb(it),ite(it)
  1818.       if(isb(ivn(iv)).ne.0)goto 9
  1819.     5 continue
  1820.     6 continue
  1821.       return
  1822.     9 jv=ivn(iv)
  1823.       call subs(jv,jv,j5)
  1824.       goto 1
  1825.       end
  1826.  
  1827.  
  1828.       subroutine pack
  1829. $include: 'picom.for'
  1830.       if(minns.le.ns)goto 9
  1831.       is=iosi(ns)
  1832.       nt=ise(is)
  1833.       itetop=ite(nt)
  1834.       minns=ns+1
  1835.       return
  1836.     9 if(minns.le.0)call bomb(3,'   pack4')
  1837.       js=minns
  1838.       kt=0
  1839.       kv=0
  1840.       if(js.eq.1)goto 13
  1841.       kt=ise(iosi(js-1))
  1842.       kv=ite(kt)
  1843.    13 do 30 ms=js,ns
  1844.       is=iosi(ms)
  1845.       if(is.le.0)call bomb(3,'   pack6')
  1846.       ia=isb(is)
  1847.       ib=ise(is)
  1848.       lv=0
  1849.       lt=0
  1850.       do 14 it=ia,ib
  1851.       if(itn(it).eq.0)goto 14
  1852.       lt=lt+1
  1853.       if(lt.gt.ntermt)call bomb(4,'NTERMT  ')
  1854.       jtn(lt)=itn(it)
  1855.       jtd(lt)=itd(it)
  1856.       jtb(lt)=lv+1
  1857.       do 12 i=itb(it),ite(it)
  1858.       if(ivp(i).ne.0)goto 31
  1859.       do 32 j=itb(it),ite(it)
  1860.       if(ivp(j).ne.0)goto 12
  1861.    32 continue
  1862.       lv=lv+1
  1863.       jvn(lv)=2
  1864.       jvp(lv)=0
  1865.       goto 33
  1866.    31 lv=lv+1
  1867.       jvn(lv)=ivn(i)
  1868.       jvp(lv)=ivp(i)
  1869.    12 continue
  1870.    33 jte(lt)=lv
  1871.       if(lv.gt.ntott)call bomb(4,'NTOTT   ')
  1872.    14 continue
  1873.       if(lv.ne.0)goto 15
  1874.       kt=kt+1
  1875.       isb(is)=kt
  1876.       ise(is)=kt
  1877.       kv=kv+1
  1878.       itn(kt)=0
  1879.       itd(kt)=1
  1880.       itb(kt)=kv
  1881.       ite(kt)=kv
  1882.       ivn(kv)=2
  1883.       ivp(kv)=0
  1884.       goto 30
  1885.    15 isb(is)=kt+1
  1886.       do 18 it=1,lt
  1887.       kt=kt+1
  1888.       itn(kt)=jtn(it)
  1889.       itd(kt)=jtd(it)
  1890.       itb(kt)=kv+1
  1891.       do 16 i=jtb(it),jte(it)
  1892.       kv=kv+1
  1893.       ivn(kv)=jvn(i)
  1894.    16 ivp(kv)=jvp(i)
  1895.    18 ite(kt)=kv
  1896.       ise(is)=kt
  1897.    30 continue
  1898.       nt=kt
  1899.       itetop=kv
  1900.       minns=ns+1
  1901.       return
  1902.       end
  1903.  
  1904.  
  1905.       subroutine print(ii)
  1906. $include: 'picom.for'
  1907.       character*127 scr
  1908.       character*1 pa,aa,bl,sp,sm,cf
  1909.       common/pr1/pa(232),aa(22)
  1910.       dimension dum(4)
  1911.       data bl,sp,sm,cf/' ','+','-','^'/
  1912.       if(np.eq.0)return
  1913.       ia=isb(ii)
  1914.       if(ia.eq.0)write(8,2)
  1915.       if(ia.eq.0)return
  1916.     2 format(' there is nothing to print')
  1917.       ib=ise(ii)
  1918.       nst=23
  1919.       do 40 it=ia,ib
  1920.       n=nst
  1921.       do 4 i=1,232
  1922.     4 pa(i)=' '
  1923.       ka=itb(it)
  1924.       kb=ite(it)
  1925.       do 30 k=ka,kb
  1926.       n=n+1
  1927.       pa(n)=bl
  1928.       if=ivn(k)
  1929.       jp=ivp(k)
  1930.       ifl=ipv(if)
  1931.       do 20 j=1,ifl
  1932.       n=n+1
  1933.    20 pa(n)=pv(j,if)
  1934.       if(jp.eq.1)goto 30
  1935.       n=n+1
  1936.       if(jp.ge.0)pa(n)=cf
  1937.       if(jp.lt.0)pa(n)=sm
  1938.       iap=iabs(jp)
  1939.       if(iap.le.99)goto 21
  1940.       pa(n+1)='*'
  1941.       pa(n+2)='*'
  1942.       n=n+2
  1943.       goto 30
  1944.    21 n=n+1
  1945.       j=mod(iap,10)
  1946.       pa(n)=hn(j+1)
  1947.       if(iap.le.9)goto 30
  1948.       j=iap/10
  1949.       n=n+1
  1950.       pa(n)=pa(n-1)
  1951.       pa(n-1)=hn(j+1)
  1952.    30 if(n.gt.232)call bomb(3,'  print3')
  1953.       write(scr,31)itn(it)
  1954.       read(scr,36)(pa(kk),kk=1,11)
  1955.    31 format(i11)
  1956.       if(itd(it).eq.1)goto 33
  1957.       pa(12)='/'
  1958.       write(scr,31)itd(it)
  1959.       read(scr,36)(aa(kk),kk=1,11)
  1960.    36 format(22a1)
  1961.       i=12
  1962.       do 32 j=1,11
  1963.       if(aa(j).eq.' ')goto 32
  1964.       i=i+1
  1965.       if(i.ge.nst)call bomb(3,'print3  ')
  1966.       pa(i)=aa(j)
  1967.    32 continue
  1968.    33 j=mod(it+1-ia,1000)
  1969.       ir=23
  1970.       if(pa(ir).ne.' ')goto 41
  1971.       do 50 m=ir,4,-1
  1972.       if(pa(m).ne.' ')goto 51
  1973.    50 continue
  1974.    51 mm=min0(ir-m,8)
  1975.       do 55 i=ir,mm+1,-1
  1976.       pa(i)=pa(i-mm)
  1977.    55 pa(i-mm)=' '
  1978.    41 write(8,34)j,(pa(k),k=1,n)
  1979.    34 format(1x,i3,126a1,/,26x,106a1)
  1980.    40 continue
  1981.       return
  1982.       end
  1983.  
  1984.  
  1985.       subroutine prntap
  1986. $include: 'picom.for'
  1987.       do 1 n=nchpl,1,-1
  1988.       if(ap(n).ne.' ')goto 2
  1989.     1 continue
  1990.       n=1
  1991.     2 write(8,3)(ap(i),i=1,n)
  1992.     3 format(1x,127a1)
  1993.       return
  1994.       end
  1995.  
  1996.  
  1997.       subroutine read(ib)
  1998. $include: 'picom.for'
  1999. c  ib=0 on return if line has nonblanks
  2000. c  ib=1 on return if line is all blanks
  2001.       line=line+1
  2002.       ib=inpba(line)
  2003.       ie=inpea(line)
  2004.       il=ie-ib+1
  2005.       do 1 i=1,nchpl
  2006.     1 a(i)=' '
  2007.       call copyc(ainpl(ib),a,il)
  2008.       call copyc(a,ap,nchpl)
  2009.       i=1
  2010.       if(a(1).eq.'*')goto 8
  2011.       do 6 i=1,nchpl
  2012.       if(a(i).eq.';')goto 8
  2013.     6 continue
  2014.       goto 14
  2015.     8 do 10 j=i,nchpl
  2016.    10 a(j)=' '
  2017.    14 ib=0
  2018.       do 16 i=1,nchpl
  2019.       if(a(i).ne.' ')goto 18
  2020.    16 continue
  2021.    17 ib=1
  2022.       return
  2023.    18 if(i.eq.1)return
  2024.       do 20 j=i,nchpl
  2025.       a(j+1-i)=a(j)
  2026.    20 a(j)=' '
  2027.       return
  2028.       end
  2029.  
  2030.  
  2031.       subroutine readin
  2032. $include: 'picom.for'
  2033.       do 1 i=1,nlab
  2034.     1 linlab(i)=-1
  2035.       line=0
  2036.       lines=0
  2037.       labels=0
  2038.       iinpl=1
  2039.     5 read(7,6)a
  2040.     6 format(127a1)
  2041.       do 7 k=nchpl,1,-1
  2042.       if(a(k).ne.' ')goto 10
  2043.     7 continue
  2044.       k=1
  2045.    10 lines=lines+1
  2046.       lintot=lintot+1
  2047.       if(lines.gt.ninl)call bomb(4,'NINL    ')
  2048.       if(k.eq.1)goto 15
  2049.       call srchnb(1,k,m)
  2050.       if(a(m).eq.'p' .and. a(m+1).eq.'r' .and. a(m+2).eq.'o')goto 11
  2051.       if(a(m).ne.'l' .or. a(m+1).ne.'a' .or. a(m+2).ne.'b')goto 15
  2052.    11 call nxt(m+2,nchpl,' ',j1)
  2053.       call nxtch(j1,nchpl,j3)
  2054.       call nxt(j3,nchpl,' ',j4)
  2055.       call numl(j1,nchpl,j3,j4,j5)
  2056.       if(linlab(j5).ne.-1)write(8,12)(a(i),i=j3,j4)
  2057.   12  format(' duplicate label: ', 15a1)
  2058.       if(linlab(j5).ne.-1)call bomb(3,' readin1')
  2059.       linlab(labels)=lines
  2060.    15 inpba(lines)=iinpl
  2061.       inpea(lines)=iinpl+k-1
  2062.       if(inpea(lines).gt.ninch)call bomb(4,'NINCH   ')
  2063.       call copyc(a,ainpl(iinpl),k)
  2064.       iinpl=iinpl+k
  2065.       if(k.eq.1)goto 5
  2066.       if(a(m).eq.'f' .and. a(m+1).eq.'l' .and. a(m+2).eq.'u')goto 16
  2067.       if(a(m).ne.'e' .or. a(m+1).ne.'n' .or. a(m+2).ne.'d')goto 5
  2068.       if(a(m+3).ne.' ' .or. a(m+4).ne.'i' .or. a(m+5).ne.'n')goto 5
  2069.    16 write(8,20)lines,labels,iinpl
  2070.    20 format(' number of lines in input file=',i5,' labels='
  2071.      +      ,i2,'  characters=',i6)
  2072.       return
  2073.       end
  2074.  
  2075.  
  2076.       subroutine recip(j13,j4)
  2077. $include: 'picom.for'
  2078.       integer*4 in,id
  2079.       ia=isb(j13)
  2080.       iz=ise(j13)
  2081.       k=itb(ia)
  2082.       if(ivn(k).eq.2 .and. ivp(k).ne.0 .and. ia.ne.iz)then
  2083.          write(8,2)
  2084.          call bomb(3,'  recip1')
  2085.     2 format(' this simple program cannot process the denominator.',
  2086.      +     '  denominators should have 1 term without variable 1')
  2087.       endif
  2088.       ja=nt+1
  2089.       call copyt(ia,ia,1,1)
  2090.       in=itn(ja)
  2091.       id=itd(ja)
  2092.       itn(ja)=id
  2093.       if(in.lt.0)itn(ja)=-id
  2094.       itd(ja)=iabs(in)
  2095.       do 4 i=itb(ja),ite(ja)
  2096.     4 ivp(i)=-ivp(i)
  2097.       if(iz.eq.ia)call defsum(j4,ja,nt)
  2098.       if(iz.eq.ia)return
  2099.       if(maxe.gt.2)write(8,6)
  2100.       if(maxe.gt.2)call bomb(3,'  recip2')
  2101.     6 format(' this simple program cannot divide if maxe is',
  2102.      +      ' greater than 2')
  2103.       ib=ia+1
  2104.       k=itb(ib)
  2105.       if(ivn(k).ne.2 .or. ivp(k).lt.1)write(8,2)
  2106.       if(ivn(k).ne.2 .or. ivp(k).lt.1)call bomb(3,'  recip3')
  2107.       if(ivp(k).eq.2)goto 14
  2108.       ic=ib
  2109.       if(ib.eq.iz)goto 10
  2110.     8 k=itb(ic+1)
  2111.       if(ic.eq.iz .or. ivp(k).eq.2)goto 10
  2112.       ic=ic+1
  2113.       goto 8
  2114.    10 jb=nt+1
  2115.       call copyt(ib,ic,1,1)
  2116.       jc=nt
  2117.       if(ic.eq.iz)goto 12
  2118.    11 jd=nt+1
  2119.       call copyt(ic+1,iz,1,1)
  2120.       je=nt
  2121.       goto 16
  2122.    12 jd=nt+1
  2123.       je=jd
  2124.       itb(jd)=itetop+1
  2125.       call incnt
  2126.       call nite(itb(jd))
  2127.       k=itetop
  2128.       ivn(k)=2
  2129.       ivp(k)=2
  2130.       itn(jd)=0
  2131.       itd(jd)=1
  2132.       goto 16
  2133.    14 jb=nt+1
  2134.       jc=jb
  2135.       itb(jb)=itetop+1
  2136.       call incnt
  2137.       call nite(itb(jb))
  2138.       k=itetop
  2139.       ivn(k)=2
  2140.       ivp(k)=1
  2141.       itn(jb)=0
  2142.       itd(jb)=1
  2143.       ic=ib-1
  2144.       goto 11
  2145.    16 jf=nt+1
  2146.       call multn(ja,ja,jb,jc)
  2147.       jg=nt
  2148.       if(maxe.eq.1)goto 18
  2149.       jh=nt+1
  2150.       call multn(jf,jg,jf,jg)
  2151.       ji=nt
  2152.       jj=nt+1
  2153.       call multn(ja,ja,jd,je)
  2154.       jk=nt
  2155.       jl=nt+1
  2156.       call copyt(jf,jg,-1,1)
  2157.       call copyt(jh,ji,1,1)
  2158.       call copyt(jj,jk,-1,1)
  2159.       jm=nt
  2160.       call copyt(ja,ja,1,1)
  2161.       call multn(ja,ja,jl,jm)
  2162.       nt1=jm+1
  2163.       goto 24
  2164.    18 do 19 i=jf,jg
  2165.    19 itn(i)=-itn(i)
  2166.       call copyt(ja,ja,1,1)
  2167.       call multn(ja,ja,jf,jg)
  2168.       nt1=jg+1
  2169.    24 call defsum(j4,nt1,nt)
  2170.       call order(j4)
  2171.       return
  2172.       end
  2173.  
  2174.  
  2175.       subroutine swap(j4)
  2176. $include: 'picom.for'
  2177.       if=isb(j4)
  2178.       ig=ios(j4)
  2179.       if(if.ne.0)minns=min0(minns,ios(j4))
  2180.       ios(j4)=ios(1)
  2181.       iosi(ios(j4))=j4
  2182.       isb(j4)=isb(1)
  2183.       ise(j4)=ise(1)
  2184.       if(if.eq.0)goto 5
  2185.       ios(1)=ig
  2186.       iosi(ig)=1
  2187.       isb(1)=if
  2188.       ise(1)=if
  2189.       ite(if)=itb(if)
  2190.       call pack
  2191.       return
  2192.     5 call niltrm
  2193.       isb(1)=0
  2194.       call defsum(1,nt,nt)
  2195.       return
  2196.       end
  2197.  
  2198.  
  2199.       subroutine srch(j1,j2,h,j3)
  2200. $include: 'picom.for'
  2201.       character*1 h
  2202.       do 1 j3=j1,j2
  2203.       if(a(j3).eq.h)return
  2204.     1 continue
  2205.       j3=-1
  2206.       return
  2207.       end
  2208.  
  2209.  
  2210.       subroutine srchch(j1,j2,j3)
  2211. $include: 'picom.for'
  2212.       do 2 j3=j1,j2
  2213.       do 1 i=1,26
  2214.       if(a(j3).eq.ha(i))return
  2215.       if(a(j3).eq.hc(i))return
  2216.     1 continue
  2217.     2 continue
  2218.       j3=-1
  2219.       end
  2220.  
  2221.  
  2222.       subroutine srchnb(j1,j2,j3)
  2223. $include: 'picom.for'
  2224.       do 1 j3=j1,j2
  2225.       if(a(j3).ne.' ')return
  2226.     1 continue
  2227.       j3=-1
  2228.       return
  2229.       end
  2230.  
  2231.  
  2232.       subroutine sub
  2233. $include: 'picom.for'
  2234.       integer*4 inc4,idc4
  2235.       character*4 aj4,aj8
  2236.       call nxt(1,nchpl,' ',j1)
  2237.       call nxtnb(j1,nchpl,j2)
  2238.       call nxtch(j1,nchpl,j3)
  2239.       if(j2.ne.j3)goto 10
  2240.       call numv(j1,nchpl,j2,j3,j4)
  2241.       aj4='sum'
  2242.       if(isb(j4).eq.0)aj4='var'
  2243.     1 call nxt(j3+1,nchpl,'f',j5)
  2244.       if(a(j5+1).ne.'o')call bomb(3,' readin2')
  2245.       if(a(j5+2).ne.'r')call bomb(3,' readin4')
  2246.       if(a(j5+3).ne.' ')call bomb(3,' readin6')
  2247.       call numv(j5+3,nchpl,j6,j7,j8)
  2248.       aj8='var'
  2249.       if(isb(j8).gt.0)aj8='sum'
  2250.       call nxt(j7+1,nchpl,'i',j9)
  2251.       if(a(j9+1).ne.'n')call bomb(3,' readin8')
  2252.       if(a(j9+2).ne.' ')call bomb(3,'readin10')
  2253.       call numv(j9+2,nchpl,j10,j11,j12)
  2254.       minns=min0(minns,ios(j12))
  2255.       if(aj4.eq.'num')goto 35
  2256.       if(aj4.eq.'var')goto 30
  2257.       if(aj8.eq.'sum')goto 50
  2258.       call subs(j4,j8,j12)
  2259.       return
  2260.    10 call coef(j1,j4,inc4,idc4)
  2261.       inc=inc4
  2262.       idc=idc4
  2263.       aj4='num'
  2264.       j3=j4-1
  2265.       goto 1
  2266.    30 if(aj8.eq.'sum')goto 40
  2267.       call subv(j4,j8,j12)
  2268.       return
  2269.    35 if(aj8.eq.'sum')call bomb(3,'readin12')
  2270.       call subn(inc,idc,j8,j12)
  2271.       return
  2272.    40 call subr(j8,j4,j12)
  2273.       return
  2274.    50 a(71)='A'
  2275.       a(72)='?'
  2276.       a(73)='0'
  2277.       call numv(71,73,j15,j16,j17)
  2278.       call subr(j8,j17,j12)
  2279.       call subs(j4,j17,j12)
  2280.       return
  2281.       end
  2282.  
  2283.  
  2284.       subroutine subn(in,id,j8,j12)
  2285. $include: 'picom.for'
  2286.       integer*4 in4,id4
  2287.       if(isb(j12).eq.0)call bomb(3,'   subn1')
  2288.       if(isb(j8).ne.0)call bomb(3,'   subn2')
  2289.       do 8 it=isb(j12),ise(j12)
  2290.       do 5 i=itb(it),ite(it)
  2291.       if(ivn(i).ne.j8)goto 5
  2292.       in4=in
  2293.       id4=id
  2294.       in4=in4**ivp(i)
  2295.       id4=id4**ivp(i)
  2296.       call mulrat(itn(it),itd(it),in4,id4,itn(it),itd(it))
  2297.       ivp(i)=0
  2298.       if(ite(it).eq.itb(it))goto 5
  2299.       ite(it)=ite(it)-1
  2300.       if(i.gt.ite(it))goto 5
  2301.       do 1 j=i,ite(it)
  2302.       ivn(j)=ivn(j+1)
  2303.     1 ivp(j)=ivp(j+1)
  2304.     5 continue
  2305.     8 continue
  2306.       call order(j12)
  2307.       return
  2308.       end
  2309.  
  2310.  
  2311.       subroutine subr(ia,ib,ic)
  2312. c  replace pattern=first term in sum ia by variable ib in sum ic
  2313. $include: 'picom.for'
  2314.       integer*4 inj1,idj1
  2315.       if(isb(ia).eq.0)call bomb(3,'   subr1')
  2316.       if(isb(ib).ne.0)call bomb(3,'   subr2')
  2317.       nt1=nt+1
  2318.       j1=isb(ia)
  2319.       j2=itb(j1)
  2320.       j3=ite(j1)
  2321.       j4=isb(ic)
  2322.       j5=ise(ic)
  2323.       do 30 it=j4,j5
  2324.       j6=itb(it)
  2325.       j7=ite(it)
  2326.       minrat=30000
  2327.       do 5 j8=j2,j3
  2328.       jfp=ivn(j8)
  2329.     1 jfm=ivn(j6)
  2330.       if(jfm.gt.jfp)goto 22
  2331.       irat=ivp(j6)/ivp(j8)
  2332.       if(jfm.eq.jfp .and. irat.ge.1)goto 4
  2333.       jt=it+1-j4
  2334.       if(j2.eq.j3 .and. np.ge.2 .and. jfm.eq.jfp .and. irat.le.-1)
  2335.      +  write(8,2)jt,(pv(i,ic),i=1,nch),(pv(i,jfp),i=1,nch),
  2336.      +            (pv(i,ia),i=1,nch)
  2337.     2 format(' term',i4,' of ', 15a1,' has exponent of ', 15a1,
  2338.      +       ' with different sign from ', 15a1)
  2339.       j6=j6+1
  2340.       if(j6.le.j7)goto 1
  2341.       goto 22
  2342.     4 minrat=min0(minrat,ivp(j6)/ivp(j8))
  2343.     5 continue
  2344.       call incnt
  2345.       idj1=itd(j1)
  2346.       if(itn(j1).lt.0)idj1=-idj1
  2347.       inj1=iabs(itn(j1))
  2348.       call mulrat(itn(it),itd(it),idj1,inj1,itn(nt),itd(nt))
  2349.       itb(nt)=itetop+1
  2350.       j6=itb(it)
  2351.       k=itb(nt)
  2352.       do 20 j=j6,j7
  2353.       do 12 j8=j2,j3
  2354.       if(ivn(j8).eq.ivn(j))goto 16
  2355.    12 continue
  2356.       ivn(k)=ivn(j)
  2357.       ivp(k)=ivp(j)
  2358.       goto 20
  2359.    16 ivn(k)=ivn(j)
  2360.       ivp(k)=ivp(j)-minrat*ivp(j8)
  2361.    20 k=k+1
  2362.       ivn(k)=ib
  2363.       ivp(k)=minrat
  2364.       call nite(k)
  2365.       goto 30
  2366.    22 call copyt(it,it,1,1)
  2367.    30 continue
  2368.       call defsum(ic,nt1,nt)
  2369.       call order(ic)
  2370.       return
  2371.       end
  2372.  
  2373.  
  2374.       subroutine subs(in1,in2,in3)
  2375. c  substitute sum=in1 for variable=in2 in sum=in3
  2376. c  result goes in nt+1,...
  2377. $include: 'picom.for'
  2378.       integer*4 inp,idp
  2379.       common/nwt/iitb,iite,jitb,jite,ntb,if
  2380. c  dimension of ittbs ittes is max jp
  2381.       dimension ittbs(30),ittes(30)
  2382.       jp=1
  2383.       ia=in2
  2384.       ib=isb(in3)
  2385.       ic=ise(in3)
  2386.       id=isb(in1)
  2387.       ie=ise(in1)
  2388.       if(ib.eq.0 .or. id.eq.0)call bomb(3,'   subs1')
  2389.       ntb=nt+1
  2390.       lev1=0
  2391.       do 30 it=ib,ic
  2392.       iitb=itb(it)
  2393.       iite=ite(it)
  2394.       if(ivn(iitb).eq.2)lev1=ivp(iitb)
  2395.       if(lev1.gt.maxe)goto 31
  2396.       do 5 if=iitb,iite
  2397.       if(ivn(if).ne.ia)goto 5
  2398.       ifp=ivp(if)
  2399.       lev2=0
  2400.       if(ifp.eq.1)goto 10
  2401.       if(id.eq.ie)goto 70
  2402.       if(ifp.gt.0)goto 17
  2403.       itpr=it-ib+1
  2404.       if(np.ge.2)write(8,1)(pv(i,ia),i=1,nch),
  2405.      +       (pv(i,in3),i=1,nch),itpr
  2406.     1 format(2x, 15a1,' occurs to a negative power in ', 15a1
  2407.      +      '  term ',i4)
  2408.     5 continue
  2409. c  here copy term it to nt+1
  2410.       if=0
  2411.       jitb=0
  2412.       call sub2(itn(it),itd(it))
  2413.       goto 30
  2414. c  here copy term it with substitution
  2415.    10 do 68 jt=id,ie
  2416.       jitb=itb(jt)
  2417.       jite=ite(jt)
  2418.       if(ivn(jitb).eq.2)lev2=ivp(jitb)
  2419.       if(lev1+lev2.gt.maxe)goto 30
  2420.       call mulrat(itn(it),itd(it),itn(jt),itd(jt),inp,idp)
  2421.       call sub2(inp,idp)
  2422.    68 continue
  2423.       goto 30
  2424.    70 if(ivn(itb(id)).eq.2)lev2=ivp(itb(id))*ifp
  2425.       if(lev1+lev2.gt.maxe)goto 30
  2426.       k=0
  2427.       do 71 i=itb(id),ite(id)
  2428.       k=k+1
  2429.       jvn(k)=ivn(i)
  2430.    71 jvp(k)=ivp(i)*ifp
  2431.       jitb=1
  2432.       jite=k
  2433.       if(ifp.lt.0)goto 72
  2434.       call mulrat(itn(it),itd(it),itn(id)**ifp,itd(id)**ifp,inp,idp)
  2435.       goto 73
  2436.    72 ii=-ifp
  2437.       call mulrat(itn(it),itd(it),itd(id)**ii,itn(id)**ii,inp,idp)
  2438.    73 call sub3(inp,idp)
  2439.       goto 30
  2440. c  here get in1**ifp
  2441.    17 jd=id
  2442.       je=ie
  2443.       ifpsav=ifp
  2444.       ntbsav=ntb
  2445.       iitbsv=iitb
  2446.       iitesv=iite
  2447.       ifsav=if
  2448.       ntbmu=nt
  2449.       itesav=itetop
  2450. c  test whether this power is in memory
  2451.       if(jp.eq.1)goto 80
  2452.       if(ifp.lt.jp)goto 27
  2453.       if(ittbs(jp).eq.0)goto 30
  2454.       if(ifp.eq.jp)goto 27
  2455.       jd=nt+1
  2456.       k=itetop
  2457.       do 86 i=ittbs(jp),ittes(jp)
  2458.       call incnt
  2459.       itn(nt)=jtn(i)
  2460.       itd(nt)=jtd(i)
  2461.       itb(nt)=k+1
  2462.       do 84 j=jtb(i),jte(i)
  2463.       k=k+1
  2464.       ivn(k)=jvn(j)
  2465.    84 ivp(k)=jvp(j)
  2466.       call nite(k)
  2467.    86 continue
  2468.       je=nt
  2469.       goto 18
  2470.    80 ks=0
  2471.       nts=0
  2472.    18 jpb=nt+1
  2473.       ntbtmu=nt
  2474.       call mult(id,ie,jd,je)
  2475.       if(nt.eq.ntbtmu)then
  2476.         jp=jp+1
  2477.         ittbs(jp)=0
  2478.         nt=ntbmu
  2479.         ntb=ntbsav
  2480.         itetop=itesav
  2481.         goto 30
  2482.         endif
  2483.       jpe=nt
  2484.       jp=jp+1
  2485.       if(jp.gt.30)call bomb(4,'  max jp')
  2486.       jd=jpb
  2487.       je=jpe
  2488.       ittbs(jp)=nts+1
  2489.       do 21 i=jpb,jpe
  2490.       nts=nts+1
  2491.       jtn(nts)=itn(i)
  2492.       jtd(nts)=itd(i)
  2493.       jtb(nts)=ks+1
  2494.       do 20 j=itb(i),ite(i)
  2495.       ks=ks+1
  2496.       jvn(ks)=ivn(j)
  2497.    20 jvp(ks)=ivp(j)
  2498.       jte(nts)=ks
  2499.       if(ks.gt.ntott)call bomb(4,'NTOTT   ')
  2500.       if(nts.gt.ntermt)call bomb(4,'NTERMT  ')
  2501.    21 continue
  2502.       ittes(jp)=nts
  2503.       if(jp.ne.ifp)goto 18
  2504.       nt=ntbmu
  2505.       itetop=itesav
  2506.       if=ifsav
  2507.       ifp=ifpsav
  2508.       iitb=iitbsv
  2509.       iite=iitesv
  2510.       ntb=ntbsav
  2511.    27 do 54 jt=ittbs(ifp),ittes(ifp)
  2512.       jitb=jtb(jt)
  2513.       jite=jte(jt)
  2514.       if(jvn(jitb).eq.2)lev2=jvp(jitb)
  2515.       if(lev1+lev2.gt.maxe)goto 30
  2516.       call mulrat(itn(it),itd(it),jtn(jt),jtd(jt),inp,idp)
  2517.       call sub3(inp,idp)
  2518.    54 continue
  2519.    30 continue
  2520.    31 call defsum(in3,ntb,nt)
  2521.       return
  2522.       end
  2523.  
  2524.  
  2525.       subroutine subv(j4,j8,j12)
  2526. $include: 'picom.for'
  2527.       if(isb(j12).eq.0)call bomb(3,'   subv1')
  2528.       if(isb(j4)+isb(j8).ne.0)call bomb(3,'   subv2')
  2529.       do 8 it=isb(j12),ise(j12)
  2530.       do 5 i=itb(it),ite(it)
  2531.       if(ivn(i).eq.j8)ivn(i)=j4
  2532.     5 continue
  2533.     8 continue
  2534.       call order(j12)
  2535.       return
  2536.       end
  2537.  
  2538.  
  2539.       subroutine sub2(in,id)
  2540. $include: 'picom.for'
  2541.       common/nwt/iitb,iite,jitb,jite,ntb,if
  2542.       integer*4 in,id
  2543.       if(in.eq.0)return
  2544.       call incnt
  2545.       k=itetop
  2546.       itb(nt)=k+1
  2547.       itn(nt)=in
  2548.       itd(nt)=id
  2549.       ii=iitb
  2550.       if(jitb.eq.0)goto 15
  2551.       jj=jitb
  2552.    11 if(ii.gt.iite)goto 13
  2553.       if(ii.eq.if)goto 12
  2554.       if(jj.gt.jite)goto 15
  2555.       if(ivn(ii).lt.ivn(jj))goto 16
  2556.       if(ivn(ii).eq.ivn(jj))goto 57
  2557.       if(ivp(jj).eq.0)goto 18
  2558.       k=k+1
  2559.       ivn(k)=ivn(jj)
  2560.       ivp(k)=ivp(jj)
  2561.    18 jj=jj+1
  2562.       goto 11
  2563.    12 ii=ii+1
  2564.       goto 11
  2565.    13 if(jj.gt.jite)goto 60
  2566.       if(ivp(jj).eq.0)goto 17
  2567.       k=k+1
  2568.       ivn(k)=ivn(jj)
  2569.       ivp(k)=ivp(jj)
  2570.    17 jj=jj+1
  2571.       goto 13
  2572.    14 if(ii.gt.iite)goto 60
  2573.    15 if(ivp(ii).eq.0)goto 9
  2574.       k=k+1
  2575.       ivn(k)=ivn(ii)
  2576.       ivp(k)=ivp(ii)
  2577.     9 ii=ii+1
  2578.       if(ii.eq.if)goto 9
  2579.       goto 14
  2580.    16 if(ivp(ii).eq.0)goto 12
  2581.       k=k+1
  2582.       ivn(k)=ivn(ii)
  2583.       ivp(k)=ivp(ii)
  2584.       ii=ii+1
  2585.       goto 11
  2586.    57 if(ivp(ii)+ivp(jj).eq.0)goto 58
  2587.       k=k+1
  2588.       ivn(k)=ivn(ii)
  2589.       ivp(k)=ivp(ii)+ivp(jj)
  2590.    58 ii=ii+1
  2591.       jj=jj+1
  2592.       goto 11
  2593.    60 if(k.gt.itetop)goto 56
  2594.       k=itetop+1
  2595.       ivn(k)=2
  2596.       ivp(k)=0
  2597.    56 itesav=itetop
  2598.       call nite(k)
  2599.       if(nt.eq.ntb)goto 68
  2600.       if(icomp(ntb,nt))31,33,34
  2601.    31 call tzap(nt,ntb)
  2602.       goto 68
  2603.    33 call addrat(itn(ntb),itd(ntb),itn(nt),itd(nt),itn(ntb),itd(ntb))
  2604.       nt=nt-1
  2605.       itetop=itesav
  2606.       goto 68
  2607.    34 if(icomp(nt-1,nt))62,61,68
  2608.    61 call addrat(itn(nt-1),itd(nt-1),itn(nt),itd(nt),
  2609.      ,     itn(nt-1),itd(nt-1))
  2610.       nt=nt-1
  2611.       itetop=itesav
  2612.       goto 68
  2613.    62 la=ntb
  2614.       lb=nt-1
  2615.    63 if(lb-la.eq.1)goto 64
  2616.       lc=(la+lb)/2
  2617.       if(icomp(lc,nt))65,66,67
  2618.    64 call tzap(nt,lb)
  2619.       goto 68
  2620.    65 lb=lc
  2621.       goto 63
  2622.    66 call addrat(itn(lc),itd(lc),itn(nt),itd(nt),itn(lc),itd(lc))
  2623.       nt=nt-1
  2624.       itetop=itesav
  2625.       goto 68
  2626.    67 la=lc
  2627.       goto 63
  2628.    68 return
  2629.       end
  2630.  
  2631.  
  2632.       subroutine sub3(in,id)
  2633. $include: 'picom.for'
  2634.       integer*4 in,id
  2635.       common/nwt/iitb,iite,jitb,jite,ntb,if
  2636.       if(in.eq.0)return
  2637.       call incnt
  2638.       k=itetop
  2639.       itb(nt)=k+1
  2640.       itn(nt)=in
  2641.       itd(nt)=id
  2642.       ii=iitb
  2643.       if(jitb.eq.0)goto 15
  2644.       jj=jitb
  2645.    11 if(ii.gt.iite)goto 13
  2646.       if(ii.eq.if)goto 12
  2647.       if(jj.gt.jite)goto 15
  2648.       if(ivn(ii).lt.jvn(jj))goto 16
  2649.       if(ivn(ii).eq.jvn(jj))goto 57
  2650.       if(jvp(jj).eq.0)goto 18
  2651.       k=k+1
  2652.       ivn(k)=jvn(jj)
  2653.       ivp(k)=jvp(jj)
  2654.    18 jj=jj+1
  2655.       goto 11
  2656.    12 ii=ii+1
  2657.       goto 11
  2658.    13 if(jj.gt.jite)goto 60
  2659.       if(jvp(jj).eq.0)goto 17
  2660.       k=k+1
  2661.       ivn(k)=jvn(jj)
  2662.       ivp(k)=jvp(jj)
  2663.    17 jj=jj+1
  2664.       goto 13
  2665.    14 if(ii.gt.iite)goto 60
  2666.    15 if(ivp(ii).eq.0)goto 9
  2667.       k=k+1
  2668.       ivn(k)=ivn(ii)
  2669.       ivp(k)=ivp(ii)
  2670.     9 ii=ii+1
  2671.       if(ii.eq.if)goto 9
  2672.       goto 14
  2673.    16 if(ivp(ii).eq.0)goto 12
  2674.       k=k+1
  2675.       ivn(k)=ivn(ii)
  2676.       ivp(k)=ivp(ii)
  2677.       ii=ii+1
  2678.       goto 11
  2679.    57 if(ivp(ii)+jvp(jj).eq.0)goto 58
  2680.       k=k+1
  2681.       ivn(k)=ivn(ii)
  2682.       ivp(k)=ivp(ii)+jvp(jj)
  2683.    58 ii=ii+1
  2684.       jj=jj+1
  2685.       goto 11
  2686.    60 if(k.gt.itetop)goto 56
  2687.       k=itetop+1
  2688.       ivn(k)=2
  2689.       ivp(k)=0
  2690.    56 itesav=itetop
  2691.       call nite(k)
  2692.       if(nt.eq.ntb)goto 68
  2693.       if(icomp(ntb,nt))31,33,34
  2694.    31 call tzap(nt,ntb)
  2695.       goto 68
  2696.    33 call addrat(itn(ntb),itd(ntb),itn(nt),itd(nt),itn(ntb),itd(ntb))
  2697.       nt=nt-1
  2698.       itetop=itesav
  2699.       goto 68
  2700.    34 if(icomp(nt-1,nt))62,61,68
  2701.    61 call addrat(itn(nt-1),itd(nt-1),itn(nt),itd(nt),
  2702.      ,     itn(nt-1),itd(nt-1))
  2703.       nt=nt-1
  2704.       itetop=itesav
  2705.       goto 68
  2706.    62 la=ntb
  2707.       lb=nt-1
  2708.    63 if(lb-la.eq.1)goto 64
  2709.       lc=(la+lb)/2
  2710.       if(icomp(lc,nt))65,66,67
  2711.    64 call tzap(nt,lb)
  2712.       goto 68
  2713.    65 lb=lc
  2714.       goto 63
  2715.    66 call addrat(itn(lc),itd(lc),itn(nt),itd(nt),itn(lc),itd(lc))
  2716.       nt=nt-1
  2717.       itetop=itesav
  2718.       goto 68
  2719.    67 la=lc
  2720.       goto 63
  2721.    68 return
  2722.       end
  2723.  
  2724.  
  2725.       subroutine switch(j,k,is)
  2726. $include: 'picom.for'
  2727.       do 2 it=isb(is),ise(is)
  2728.       do 1 i=itb(it),ite(it)
  2729.       if(ivn(i).eq.j)ivn(i)=k
  2730.     1 continue
  2731.     2 continue
  2732.       return
  2733.       end
  2734.  
  2735.  
  2736.  
  2737.       subroutine taylor
  2738. $include: 'picom.for'
  2739.       character*4 star
  2740.       integer*4 itf
  2741.       dimension ltay(mitay)
  2742.       data iorder,ifirst/1,0/
  2743.       if(idflg.ne.0)call dertab
  2744.       call nxt(2,nchpl,' ',j1)
  2745.       call srch(j1,nchpl,'=',j2)
  2746.       j4=j1
  2747.       if(j2.eq.-1)goto 5
  2748.       call getexp(j2+1,iorder,j3,0)
  2749.       if(iorder.lt.1)call bomb(3,' taylor1')
  2750.       return
  2751.     5 call nxt(2,nchpl,' ',j13)
  2752.       star='no'
  2753.       if(a(j13-1).eq.'*')star='yes'
  2754.    10 call nxt(j4,nchpl,'i',j3)
  2755.       if(j3.eq.-1)call bomb(3,' taylor2')
  2756.       if(a(j3-1).eq.' ' .and. a(j3+1).eq.'n' .and. a(j3+2).eq.' ')
  2757.      +       goto 12
  2758.       j4=j3+1
  2759.       if(j4.lt.nchpl)goto 10
  2760.       call bomb(3,' taylor3')
  2761.    12 do 14 j5=j3-1,j1,-1
  2762.       if(a(j5).ne.' ')goto 16
  2763.    14 continue
  2764.       call bomb(3,' taylor4')
  2765.    16 do 18 j6=j5,j1,-1
  2766.       if(a(j6).eq.' ')goto 20
  2767.    18 continue
  2768.       call bomb(3,' taylor5')
  2769.    20 call numv(j6,j5,j7,j8,j9)
  2770.       if(isb(j9).eq.0)call bomb(3,' taylor6')
  2771.       itesav=itetop
  2772.       call term(j3+2)
  2773.       itay=itetop-itesav
  2774.       if(itay.gt.mitay)call bomb(4,'MITAY   ')
  2775.       do 21 i=1,itay
  2776.    21 ltay(i)=ivn(itesav+i)
  2777.       nt=nt-1
  2778.       itetop=itesav
  2779.       call niltrm
  2780.       call defsum(1,nt,nt)
  2781.       if(ifirst.ne.0)goto 22
  2782.       a(1)='A'
  2783.       a(2)='?'
  2784.       a(3)='0'
  2785.       call numv(1,3,j14,j15,j10)
  2786.       a(3)='1'
  2787.       call numv(1,3,j14,j15,j11)
  2788.       a(3)='2'
  2789.       call numv(1,3,j14,j15,j12)
  2790.    22 ifirst=1
  2791.       nv1=nv+1
  2792.       do 34 it=isb(j9),ise(j9)
  2793.       call copyt(it,it,1,1)
  2794.       call defsum(j10,nt,nt)
  2795.       do 32 ivv=1,itay
  2796.       iv=ltay(ivv)
  2797.       if(star.eq.'yes')call switch(iv,nv1,j10)
  2798.       call copys(j10,j11)
  2799.       call subn(0,1,iv,j10)
  2800.       itf=1
  2801.       do 30 io=1,iorder
  2802.       call deriv(j11,iv,j12)
  2803.       call outqu(j12)
  2804.       if(itn(isb(j12)).eq.0)goto 31
  2805.       call copys(j12,j11)
  2806.       call subn(0,1,iv,j12)
  2807.       itf=itf*io
  2808.       call incnt
  2809.       k=itetop+1
  2810.       itn(nt)=1
  2811.       itd(nt)=itf
  2812.       itb(nt)=k
  2813.       ite(nt)=k
  2814.       ivn(k)=iv
  2815.       ivp(k)=io
  2816.       call nite(k)
  2817.       nt1=nt+1
  2818.       call mult(nt,nt,isb(j12),ise(j12))
  2819.       call defsum(j12,nt1,nt)
  2820.       call pack
  2821.    30 call adds(j12,j10,1,1,j10)
  2822.    31 if(star.eq.'yes')call switch(nv1,iv,j10)
  2823.    32 continue
  2824.    34 call adds(j10,1,1,1,1)
  2825.       call swap(j9)
  2826.       call order(j9)
  2827.       call delsum(j10)
  2828.       return
  2829.       end
  2830.  
  2831.  
  2832.       subroutine term(ist)
  2833. c  each term is on a separate line, free format with coefficient first
  2834. $include: 'picom.for'
  2835.       integer*4 inc,idc
  2836.       call incnt
  2837.       call coef(ist,j1,inc,idc)
  2838.       itn(nt)=inc
  2839.       itd(nt)=idc
  2840.       itb(nt)=itetop+1
  2841.       if(j1.eq.-1)then
  2842.         k=itb(nt)
  2843.         ivn(k)=2
  2844.         ivp(k)=0
  2845.         goto 10
  2846.         endif
  2847.       k=itetop
  2848.     5 isgn=1
  2849.       if(a(j1).eq.'*')j1=j1+1
  2850.       if(a(j1).eq.'/')isgn=-1
  2851.       if(a(j1).eq.'/')j1=j1+1
  2852.       call numv(j1,nchpl,j2,j3,j4)
  2853.       ie=1
  2854.       j5=-1
  2855.       if(j3.eq.nchpl-1 .and. a(nchpl).ne.' ')call bomb(3,'   term1')
  2856.       if(j3.le.nchpl-2)call getexp(j3+1,ie,j5,1)
  2857.       k=k+1
  2858.       ivn(k)=j4
  2859.       ivp(k)=ie*isgn
  2860.       j1=j5
  2861.       if(j5.ne.-1)goto 5
  2862.    10 call nite(k)
  2863.       return
  2864.       end
  2865.  
  2866.  
  2867.       subroutine thderv
  2868. $include: 'picom.for'
  2869.       idflg=1
  2870.       call nxt(2,nchpl,'o',j1)
  2871.       call nxt(j1,nchpl,' ',j2)
  2872.       call numv(j2,nchpl,j3,j4,j5)
  2873.       call nxt(j4+1,nchpl,'w',j6)
  2874.       call nxt(j6,nchpl,' ',j7)
  2875.       call numv(j7,nchpl,j8,j9,j10)
  2876.       call srch(j9+1,nchpl,'=',j11)
  2877.       if(j11.ne.-1)goto 1
  2878.       call srch(j9+1,nchpl,'i',j12)
  2879.       if(j12.eq.-1)call bomb(3,'thderv1')
  2880.       call nxt(j12,nchpl,' ',j11)
  2881.     1 call srchnb(j11+1,nchpl,j12)
  2882.       if(j12.eq.-1)call bomb(3,'thderv2')
  2883.       ig=0
  2884.       if(a(j12).eq.'0')ig=1
  2885.       if(ig.eq.1)goto 8
  2886.       call numv(j11+1,nchpl,j12,j13,j14)
  2887.     8 if=0
  2888.       do 10 i=1,nd
  2889.       if(ider(1,i).eq.j5 .and. ider(2,i).eq.j10)if=1
  2890.       if(ider(1,i).eq.j5 .and. ider(2,i).eq.j10)in=i
  2891.    10 continue
  2892.       if(if.eq.0 .and. ig.eq.1)goto 50
  2893.       if(if.eq.0)goto 14
  2894.       ja=ider(3,in)
  2895.       jb=ipv(ja)
  2896.       if(np.ge.2)write(8,11)(pv(i,ja),i=1,jb)
  2897.    11 format(' REDEFINITION: previously the  derivative was = ',15a1)
  2898.       if(ig.eq.0)goto 18
  2899.       nd=nd-1
  2900.       if(in.gt.nd)goto 50
  2901.       do 28 ih=in,nd
  2902.       ider(1,ih)=ider(1,ih+1)
  2903.       ider(2,ih)=ider(2,ih+1)
  2904.    28 ider(3,ih)=ider(3,ih+1)
  2905.       goto 50
  2906.    14 nd=nd+1
  2907.       in=nd
  2908.       if(nd.gt.nder)call bomb(4,'NDER    ')
  2909.       ider(1,in)=j5
  2910.       ider(2,in)=j10
  2911.    18 ider(3,in)=j14
  2912.    50 return
  2913.       end
  2914.  
  2915.  
  2916.       subroutine tzap(ka,kb)
  2917. c  zap term at ka to position kb
  2918. $include: 'picom.for'
  2919.       integer*4 itnka,itdka
  2920.       if(ka.eq.kb)return
  2921.       itnka=itn(ka)
  2922.       itdka=itd(ka)
  2923.       itbka=itb(ka)
  2924.       iteka=ite(ka)
  2925.       if(ka.gt.kb)goto 5
  2926.       do 1 i=ka,kb-1
  2927.       itn(i)=itn(i+1)
  2928.       itd(i)=itd(i+1)
  2929.       itb(i)=itb(i+1)
  2930.     1 ite(i)=ite(i+1)
  2931.       goto 8
  2932.     5 do 6 i=ka,kb+1,-1
  2933.       itn(i)=itn(i-1)
  2934.       itd(i)=itd(i-1)
  2935.       itb(i)=itb(i-1)
  2936.     6 ite(i)=ite(i-1)
  2937.     8 itn(kb)=itnka
  2938.       itd(kb)=itdka
  2939.       itb(kb)=itbka
  2940.       ite(kb)=iteka
  2941.       return
  2942.       end
  2943.  
  2944.  
  2945.  
  2946.       subroutine zsub(j3)
  2947. $include: 'picom.for'
  2948. c this is available
  2949.       return
  2950.       end
  2951.