1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
|
c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
c Copyright (C) ????-2008 - INRIA - Francois DELBECQUE
c
c This file must be used under the terms of the CeCILL.
c This source file is licensed as described in the file COPYING, which
c you should have received as part of this distribution. The terms
c are also available at
c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
subroutine residu(p,np,a,na,b,nb,v,tol,ierr)
c but! calcul de residus
c calcul de la somme des residus de p/(a.b)
c aux zeros de a
c p=polynome de degre np
c a= na
c b= nb
c
c les zeros de b sont supposes tous differents des
c zeros de a....
c
c a,b et p dimensionnes au moins a leur degre+1 dans le pgm
c appelant.
c rangement par degres croissants.
c v=resultat
c ierr=0 O.K.
c ierr=1 mauvais appel
c si a et b ont une racine commune
c on teste la division par zero par rapport a tol
c principe du calcul:si a (resp b) est une constante on a
c v=p(nb)/b(nb+1)/a(1) (resp v=p(na)/a(na+1)/b(1) )
c sinon on remplace p et a par le reste de la division
c euclidienne de p et a par b,puis on inverse les roles
c de a et b en changeant le signe de v.
c on itere jusqu a trouver degre de a ou degre de b=0.
c routines appelees:dpodiv,idegre (bibli blaise Inria)
c
dimension a(*),b(*),p(*)
double precision a,b,p,v,r,b1,tol
v=0.0d+0
ierr=0
npp=np
call idegre(a,na,na)
call idegre(b,nb,nb)
if(na.eq.0) return
c
c b=constante v=... et return
if(nb.eq.0) then
b1=b(1)
if(b1.eq.0.0d+0) then
ierr=1
return
endif
if(npp.ge.na-1) then
v=p(na)/a(na+1)/b1
return
else
v=0.0d+0
return
endif
endif
c
c degre de b >= 1
c
if(na.le.np) then
c p=p/a (reste de la division euclidienne...)
call dpodiv(p,a,np,na)
call idegre(p,na-1,np)
endif
if(na.le.nb) then
c b=b/a (reste de la div euclidienne...)
call dpodiv(b,a,nb,na)
call idegre(b,na-1,nb)
endif
c ici nb=degre de b < na=degre de a
c et np=degre de p < na=degre de a
c si degre de a=1 (b et p =cstes) v=... et return
if(na.eq.1) then
b1=b(1)
if(abs(b1).le.tol) then
ierr=1
return
endif
v=p(na)/a(na+1)/b1
return
endif
c
c si degre de b=0 v=... et return
call idegre(b,min(na-1,nb),nb)
if(nb.eq.0) then
b1=b(1)
if(abs(b1).le.tol) then
ierr=1
return
endif
if(npp.ge.na-1) then
v=p(na)/a(na+1)/b1
return
else
v=0.0d+0
return
endif
endif
c si degre de b>0 on itere
nit=0
20 continue
if(nit.ge.1) na=nbb
nit=nit+1
nbb=nb
c a=a/b (reste de la division euclidienne...)
c p=p/b (reste de la division euclidienne...)
c
call dpodiv(a,b,na,nb)
call idegre(a,nb-1,na)
call dpodiv(p,b,np,nb)
call idegre(p,nb-1,np)
c b=-a
do 30 k=1,nb+1
r=b(k)
b(k)=-a(k)
a(k)=r
30 continue
c
c si degre de b=0 v=...
c
call idegre(b,na,nb)
if(nb.eq.0) then
b1=b(1)
if(abs(b1).le.tol) then
ierr=1
v=0.0d+0
return
endif
v=p(nbb)/a(nbb+1)/b1
return
endif
c sinon goto 20
goto 20
end
|
