← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/poly/dpsimp.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/poly/dpsimp.f
 
1
      subroutine dpsimp(a,na,b,nb,a1,na1,b1,nb1,w,ierr)
 
2
c!    but
 
3
c     Etant donnes une fraction rationnelle donnee par ses polynomes 
 
4
c     numerateur et denominateurs, ce sous programme retourne le numerateur
 
5
c     et le denominateur de sa representation simplifiee.
 
6
c!    liste d'appel
 
7
c     subroutine dpsimp(a,na,b,nb,as,nas,bs,nbs,w,ierr)
 
8
c     
 
9
c     double precision a(na+1),b(nb+1),as(*),bs(*),w,er
 
10
c     integer na,nb,nas,nbs,ierr
 
11
c     
 
12
c     a    :  tableau contenant les coefficients du numerateur range
 
13
c             par puissance croissante.(entree) 
 
14
c     na   :  degre du numerateur (entree) 
 
15
c     b    :  tableau contenant les coefficients du denominateur range
 
16
c             par puissance croissante. (entree) 
 
17
c     nb   :  degre du denominateur (entree) 
 
18
c
 
19
c     a1   :  tableau contenant les coefficients du numerateur range
 
20
c             par puissance croissante.(sortie) 
 
21
c     na1  :  degre+1 du numerateur (sortie) 
 
22
c     b1   :  tableau contenant les coefficients du denominateur range
 
23
c             par puissance croissante. (sortie) 
 
24
c     nb1  :  degre+1 du denominateur (sortie) 
 
25
c
 
26
c     les implantations de a et a1, b et b1 peuvent etre confondues.
 
27
c     Dans les cas ou les zones memoires de a (resp b) et a1 (resp b1) se
 
28
c     chevauchent, l'adresse de a1 (resp b1) doit etre au moins egale a 
 
29
c     l'adresse de a  (resp b)
 
30
c
 
31
c     w    :  tableau de travail de taille:
 
32
c             2*(na+nb)+min(na,nb)+10*max(na,nb)+3*max(na,nb)**2+4
 
33
c     ierr :  
 
34
c             en entree ierr specifie l'espace memoire disponible dans w
 
35
c             en sortie:
 
36
c     ierr=0 : ok
 
37
c     ierr=1 : denominateur nul
 
38
c     ierr=2 : espace memoire insuffisant on retourne les polynomes
 
39
c!    
 
40
c     origine S Steer INRIA 1990
 
41
c     
 
42
c     Copyright INRIA
 
43
      double precision a(na+1),b(nb+1),w(*),a1(*),b1(*),t,er,t1,t2
 
44
      integer na,nb,ierr,ipb(6)
 
45
c     
 
46
      lw=1+2*(na+nb)+min(na,nb)+3
 
47
c     n0=max(na,nb)+1
 
48
c     lfree=lw+10*n0+3*n0*n0
 
49
      
 
50
c     
 
51
      maxw=ierr
 
52
      ierr=0
 
53
c
 
54
c degre reel des polynomes
 
55
c
 
56
 
 
57
c
 
58
      nnb=nb+1
 
59
 08   nnb=nnb-1
 
60
      if(nnb.lt.0) then
 
61
         ierr=1
 
62
         return
 
63
      endif
 
64
      if(b(nnb+1).eq.0.0d+0) goto 08
 
65
 
 
66
      nna=na+1
 
67
 09   nna=nna-1
 
68
      if(nna.lt.0) goto 20
 
69
      if(a(nna+1).eq.0.0d+0) goto 09
 
70
c     
 
71
c     elimination des racines en zero
 
72
      la0=0
 
73
 10   la0=la0+1
 
74
      if(a(la0).eq.0.0d+0) goto 10
 
75
      na1=nna-(la0-1)
 
76
      nz=la0-1
 
77
c     
 
78
      lb0=0
 
79
 11   lb0=lb0+1
 
80
      if(b(lb0).eq.0.0d+0) goto 11
 
81
      nb1=nnb-(lb0-1)
 
82
      nz=nz-(lb0-1)
 
83
c
 
84
      
 
85
      n0=max(na1,nb1)+1
 
86
      lfree=lw+10*n0+3*n0*n0
 
87
      if(lfree.ge.maxw.and.na1.gt.0.and.nb1.gt.0) ierr=2
 
88
      if(lfree.ge.maxw.or.na1.eq.0.or.nb1.eq.0) then
 
89
         if(nz.eq.0) then
 
90
            call dcopy(na1+1,a(la0),1,a1,1)
 
91
            call dcopy(nb1+1,b(lb0),1,b1,1)
 
92
         elseif(nz.gt.0) then
 
93
            call dset(nz,0.0d0,a1,1)
 
94
            call dcopy(na1+1,a(la0),1,a1(nz+1),1)
 
95
            call dcopy(nb1+1,b(lb0),1,b1,1)
 
96
            na1=na1+nz
 
97
         else
 
98
            call dcopy(na1+1,a(la0),1,a1,1)
 
99
            call dset(-nz,0.0d0,b1,1)
 
100
            call dcopy(nb1+1,b(lb0),1,b1(-nz+1),1)
 
101
            nb1=nb1-nz
 
102
         endif
 
103
         na1=na1+1
 
104
         nb1=nb1+1
 
105
         return
 
106
      endif
 
107
c     normalize highest degree coefficients of num and den
 
108
      t1=a(nna+1)
 
109
      t2=b(nnb+1)
 
110
      call dscal(na1+1,1.0d0/t1,a(la0),1)
 
111
      call dscal(nb1+1,1.0d0/t2,b(lb0),1)
 
112
c     
 
113
      call recbez(a(la0),na1,b(lb0),nb1,w,ipb,w(lw),er)
 
114
      if(er.gt.1d-3) goto 30 
 
115
      nden=ipb(5)-ipb(4)
 
116
      nnum=ipb(6)-ipb(5)
 
117
      if(na1.ne.nnum-1) then
 
118
         t=w(ipb(5)-1)
 
119
         t=(1.0d0)/t
 
120
         if(nz.eq.0) then
 
121
            call dcopy(nnum,w(ipb(5)),1,a1,1)
 
122
            call dcopy(nden,w(ipb(4)),1,b1,1)
 
123
            call dscal(nden,t,b1,1)
 
124
         else if(nz.gt.0) then
 
125
            call dcopy(nnum,w(ipb(5)),1,a1(1+nz),1)
 
126
            call dset(nz,0.0d0,a1,1)
 
127
            nnum=nnum+nz
 
128
            call dcopy(nden,w(ipb(4)),1,b1,1)
 
129
            call dscal(nden,t,b1,1)
 
130
         else if(nz.lt.0) then
 
131
            nz=-nz
 
132
            call dcopy(nnum,w(ipb(5)),1,a1,1)
 
133
            call dcopy(nden,w(ipb(4)),1,b1(1+nz),1)
 
134
            call dset(nz,0.0d+0,b1,1)
 
135
            call dscal(nden,t,b1(1+nz),1)
 
136
            nden=nden+nz
 
137
         endif
 
138
         call dscal(nnum,-t*t1/t2,a1,1)
 
139
      else
 
140
         if(nz.eq.0) then
 
141
            call dcopy(nnum,a(la0),1,a1,1)
 
142
            call dcopy(nden,b(lb0),1,b1,1)
 
143
         else if(nz.gt.0) then
 
144
            call dcopy(nnum,a(la0),1,a1(1+nz),1)
 
145
            call dset(nz,0.0d+0,a1,1)
 
146
            nnum=nnum+nz
 
147
            call dcopy(nden,b(lb0),1,b1,1)
 
148
         else 
 
149
            nz=-nz
 
150
            call dcopy(nnum,a(la0),1,a1,1)
 
151
            call dcopy(nden,b(lb0),1,b1(1+nz),1)
 
152
            call dset(nz,0.0d+0,b1,1)
 
153
            nden=nden+nz
 
154
         endif
 
155
         call dscal(nnum,t1,a1,1)
 
156
         call dscal(nden,t2,b1,1)
 
157
 
 
158
      endif
 
159
      na1=nnum
 
160
      nb1=nden
 
161
      return
 
162
 20   a1(1)=0.0d+0
 
163
      b1(1)=1.0d+0
 
164
      na1=1
 
165
      nb1=1
 
166
      return
 
167
 30   if(nz.eq.0) then
 
168
         call dcopy(na1+1,a(la0),1,a1,1)
 
169
         call dcopy(nb1+1,b(lb0),1,b1,1)
 
170
      elseif(nz.gt.0) then
 
171
         call dset(nz,0.0d0,a1,1)
 
172
         call dcopy(na1+1,a(la0),1,a1(nz+1),1)
 
173
         call dcopy(nb1+1,b(lb0),1,b1,1)
 
174
         na1=na1+nz
 
175
      else
 
176
         call dcopy(na1+1,a(la0),1,a1,1)
 
177
         call dset(-nz,0.0d0,b1,1)
 
178
         call dcopy(nb1+1,b(lb0),1,b1(-nz+1),1)
 
179
         nb1=nb1-nz
 
180
      endif
 
181
      na1=na1+1
 
182
      nb1=nb1+1
 
183
      call dscal(na1,t1,a1,1)
 
184
      call dscal(nb1,t2,b1,1)
 
185
      return
 
186
 
 
187
      end