Fix polynomials display.
[scilab.git] / scilab / modules / output_stream / src / fortran / dmpdsp.f
1 c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 c Copyright (C) ????-2011 - INRIA - Serge Steer
3 c
4 c This file must be used under the terms of the CeCILL.
5 c This source file is licensed as described in the file COPYING, which
6 c you should have received as part of this distribution.  The terms
7 c are also available at
8 c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
9       subroutine dmpdsp(mp,d,nl,mm,nn,var,lvar,maxc,mode,ll,lunit,
10      1    cw,iw)
11 c!but
12 c     dmpdsp ecrit une matrice polynomiale (ou un polynome) sous
13 c     la forme d'un tableau de polynomes, avec gestion automatique de
14 c     l'espace disponible.
15 c!liste d'appel
16 c
17 c     subroutine dmpdsp(mp,d,nl,m,n,var,lvar,maxc,mode,ll,lunit,
18 c     1                  cw,iw)
19 c
20 c     double precision mp(*)
21 c     integer d(nl*n+1),nl,m,n,lvar,maxc,mode,iw(*),ll,lunit
22 c     character var*(*),cw*(*)
23 c
24 c     pm : tableau reel contenant les coefficients des polynomes,
25 c     le coefficient de degre k du polynome pm(i,j) est range
26 c     dans pm( d(i + (j-1)*nl + k) )
27 c     pm doit etre de taille au moins d(nl*n+1)-d(1)
28 c     d : tableau entier de taille nl*n+1,  si k=i+(j-1)*nl alors
29 c     d(k)) contient  l'adresse dans pm du coeff de degre 0
30 c     du polynome pm(i,j). Le degre du polynome pm(i,j) vaut:
31 c     d(k+1)-d(k) -1
32 c     nl : entier definissant le rangement dans d
33 c     m : nombre de ligne de la matrice polynomiale
34 c     n : nombre de colonnes de la matrice polynomiale
35 c     var : nom de la variable muette
36 c     lvar : nombre de caracteres de var
37 c     maxc : nombre de caracteres maximum autorise pour
38 c     representer un nombre
39 c     mode : si mode =1 alors representation variable
40 c     si mode = 0 representation d(maxc).(maxc-7)
41 c     ll : longueur de ligne maximum admissible
42 c     lunit : etiquette logique du support d'edition
43 c     cw : chaine de caracteres de travail de longueur au moins 2*ll
44 c     iw : tableau de travail entier de taille au moins egale a
45 c     d(nl*n+1)-d(1) + m*n+1
46 c!
47
48       double precision mp(*),a
49       integer d(*),iw(*),maxc,mode,fl,c1,c2,typ
50       integer sl,sk
51       logical first
52       character var*(*),cw*(*),sgn*1,dl*1
53       character*10 fexp,expo
54       double precision dlamch
55 c
56       m=abs(mm)
57       n=abs(nn)
58
59       cw=' '
60       dl=' '
61       if(m*n.gt.1) dl=' '
62 c
63 c     phase d'analyse: pour chaque coefficient a representer on determine
64 c     format avec lequel on va l'editer, on en deduit la longueur
65 c     de la representation de chacun des polynomes.
66 c     les differents formats sont stockes sous forme codee dans iw
67 c     a partir de lf
68 c     la taille respective des representation des chacun des polynomes
69 c     est contenue dans iw a partir de 1 .
70 c
71       lines=0
72       lbloc=n
73       lf=lbloc+2+n
74       nbloc=1
75       iw(lbloc+nbloc)=n
76       sk=0
77 c
78       ldg=-nl
79       ldef=lf
80       k0=1
81       do 21 k=1,n
82          sl=0
83          iw(k)=0
84          ldg=ldg+nl
85          do 20 l=1,m
86 c
87 c     traitement du polynome (l,k)
88             lp=d(ldg+l)-1
89             np=d(ldg+l+1)-d(ldg+l)
90             lgh=0
91             first=.true.
92             do 10 i=1,np
93                a=abs(mp(lp+i))
94                iw(ldef)=0
95                if(a.eq.0.0d+0) goto 09
96                first=.false.
97 c     determination du format devant representer a
98                typ=1
99                if(mode.eq.1) then
100                   call fmt(a,maxc,typ,n1,n2)
101                else
102                   if (isanan(a).eq.1) then
103                      typ=-2
104                   elseif (a.gt.dlamch('o')) then
105                      typ=-1
106                   endif
107                endif
108
109                if(typ.eq.2) then
110                   fl=n1
111                   iw(ldef)=n2+32*n1
112                elseif(typ.lt.0) then
113                   iw(ldef)=typ
114                   fl=5
115                else
116                   iw(ldef)=1
117                   fl=maxc
118                   n2=maxc-7
119                endif
120
121 c
122 c     determination de la longueur de la representation du monome,
123 c     cette longueur est a priori fl+2 (' '//sgn//rep(a)//var).
124 c     mais peut etre reduite dans des cas particulier
125                lgh=lgh+fl+2
126                if(n2.eq.0) then
127                   lgh=lgh-1
128                   if(i.ne.1.and.int(a+0.1).eq.1) lgh=lgh-1
129                endif
130                if(i.ne.1) lgh=lgh+lvar
131  09            ldef=ldef+1
132  10         continue
133 c
134 c     cas particulier du dernier exposant du polynome
135             nd=ifix(log10(0.5+np))+1
136             lgh=lgh+nd
137 c     cas particulier d'un polynome reduit a 0
138             if(first) lgh=4
139 c
140             iw(k)=max(iw(k),lgh)
141             sl=sl+(lgh/(ll-2))+1
142 c
143  20      continue
144          sk=sk+iw(k)
145          if(sk.gt.ll-2) then
146             if(k.eq.k0) then
147                iw(lbloc+nbloc)=k
148                sk=0
149                k0=k+1
150             else
151                iw(lbloc+nbloc)=k-1
152                sk=iw(k)
153                k0=k
154             endif
155             nbloc=nbloc+1
156             iw(lbloc+nbloc)=n
157             lines=lines+2*sl+m+2
158          endif
159  21   continue
160       nbloc=min(nbloc,n)
161 c
162       l1=1
163       if(mm.lt.0) then
164          write(cw(l1:l1+4),'(''eye *'')')
165          l1=l1+5
166          call basout(io,lunit,cw(1:l1-1))
167          call basout(io,lunit,' ')
168          if(io.eq.-1) goto 99
169       endif
170
171 c
172 c     phase d'edition : les deux chaines de caracteres representant
173 c     la ligne des exposants et la ligne des coefficients,sont
174 c     constituees puis imprimees.
175 c
176       k1=1
177       do 70 ib=1,nbloc
178          k2=iw(lbloc+ib)
179          ll1=0
180          if(nbloc.ne.1) then
181             call blktit(lunit,k1,k2,io)
182             if (io.eq.-1) goto 99
183          endif
184 c
185          cw(1:1)=dl
186          c1=2
187          cw(1+ll:1+ll)=dl
188          c2=2+ll
189 c
190          do 60 l=1,m
191             l1=c1
192             l2=c2
193             if(iw(k1).gt.ll-2) ll1=ll
194             do 50 k=k1,k2
195                ldg=(k-1)*nl+l
196                lp=d(ldg)-1
197                np=d(ldg+1)-d(ldg)
198                ldef=lf-1+d(ldg)-d(1)
199                first=.true.
200 c
201                l0=l1
202                do 45 j=1,np
203                   ifmt=iw(ldef+j)
204                   if(ifmt.eq.0) goto 45
205                   sgn='+'
206                   if(first) sgn=' '
207                   first=.false.
208                   if(mp(lp+j).lt.0.0d+0) sgn='-'
209                   a=abs(mp(lp+j))
210 c
211                   if(ifmt.eq.1) then
212                      fl=maxc
213                      n2=1
214                   elseif(ifmt.ge.0) then
215                      n1=ifmt/32
216                      n2=ifmt-32*n1
217                      fl=n1
218                   elseif(ifmt.lt.0) then
219 c     Inf/Nan
220                      fl=4
221                      n2=1
222                   endif
223 c
224                   nd=0
225                   if(j.gt.2) nd=ifix(log10(0.5+j))+1
226                   if(l2+fl+2+lvar+nd.gt.c2+ll-2) then
227 c     gestion des lignes suites
228                      if(l1.le.ll-1) cw(l1:ll-1)=' '
229                      if(l2.le.c2+ll-3) cw(l2:c2+ll-3)=' '
230                      cw(ll:ll)=dl
231                      call basout(io,lunit,cw(c1-1:ll))
232                      cw(c2+ll-2:c2+ll-2)=dl
233                      call basout(io,lunit,cw(c2-1:c2+ll-2))
234                      if(io.eq.-1) goto 99
235                      cw(c2:c2+9)=' '
236                      l2=c2+10
237                      cw(c1:c1+9)=' '
238                      l1=c1+10
239                   endif
240 c     representation du monome
241                   cw(l2:l2+1)=' '//sgn
242                   l2=l2+1
243                   call formatnumber(a,ifmt,maxc,cw(l2+1:),fl)
244                   l2=l2+fl
245                   if(n2.eq.0) l2=l2-1
246                   if(j.gt.1) then
247                      if(n2.eq.0.and.int(a+0.1).eq.1) l2=l2-1
248                      cw(l2+1:l2+lvar)=var(1:lvar)
249                      l2=l2+lvar
250                   endif
251                   nl1=l2+c1-c2
252                   cw(l1:nl1)=' '
253                   if(j.gt.2) then
254                      write(fexp,110) nd
255                      write(expo,fexp) j-1
256                      cw(nl1+1:nl1+nd)=expo(1:nd)
257                      l1=nl1+nd
258                   endif
259                   l1=l1+1
260                   l2=l2+1
261  45            continue
262                if(first) then
263 c     cas particulier du polynome nul
264                   cw(l1:l1+3)=' '
265                   cw(l2:l2+3)='   0'
266                   l1=l1+4
267                   l2=l2+4
268                   nd=0
269                endif
270                if(nd.ne.0) cw(l2:l2+nd-1)=' '
271                nl1=l0+iw(k)
272                if(ll1.eq.ll) nl1=ll-1
273                cw(l1:nl1)=' '
274                l1=nl1+1
275                cw(l2:c2+nl1-c1)=' '
276                l2=c2+nl1-c1+1
277  50         continue
278             if(cw(c1:l1-1).ne.' ') then
279 c     write(6,'(a)') cw(1:l2)
280 c     write(6,'(''c1,c2,l1,l2 '',4i4)') c1,c2,l1,l2
281                cw(l1:l1)=dl
282                call basout(io,lunit,cw(c1-1:l1))
283             endif
284             cw(l2:l2)=dl
285             call basout(io,lunit,cw(c2-1:l2))
286             if(l.ne.m) then
287                cw(c2:l2-1)=' '
288                call basout(io,lunit,cw(c2-1:l2))
289             endif
290             if(io.eq.-1) goto 99
291  60      continue
292          k1=k2+1
293  70   continue
294 c
295  99   return
296 c
297  110  format('(i',i2,')')
298 c
299       end