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