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
|
subroutine bezous(a,n,c,w,ierr)
c!but
c ce sous programme calcule le second coefficient de bezout du
c couple (a,at) si at=(z**n)*a(1/z) et a est un polynome de degre
c n ayant tous ces poles dans le disque unite.;a et at sont donc
c premiers
c!liste d'appel
c
c subroutine bezous(a,n,c,w,ierr)
c
c a:tableau de taille n+1 contenant les coefficients du polynome a
c ranges par puissances croissantes.
c n:degre de a
c c:tableau de taille n contenant apres execution le facteur de bezout
c les coefficients etant ranges par puissances croissantes
c w:tableau de travail de taille n*n+(n+1)/2
c ierr:indicateur d'erreur:
c si ierr=0 ok
c si ierr=1 a est de degre inferieur a n
c si ierr=2 a et at non premiers ou calcul numeriquement faux.
c
c!methode:
c la methode utilisee ici est de resoudre le systeme lineaire associe
c a la relation de bezout: a*b+at*c=1 c'est a dire:
c [x' y'] [e1] [b]
c [ ].[ ] =[ ]
c [y x ] [ 0] [c]
c ou x et y sont des matrice n*n toeplitz triangulaires superieures
c la premiere ligne de x est formee des n premiers coefficient de a
c la premiere ligne de y des n derniers ranges en ordre inverse
c c est alors solution du systeme : (y'-x'*(y**-1)*x)*c=e1
c!auteur
c serge Steer Inria 1983
c Copyright INRIA
c!sous programmes appeles
c invtpl
c ddot (blas)
c dlslv (linpack.extension)
c!
implicit double precision (a-h,o-z)
dimension a(*),c(*),w(*)
c
c calcul de la premiere ligne de y**-1
call invtpl(a(2),n,n,c,ierr)
if (ierr.ne.0) goto 70
ierr=0
c c contient les coeff de la premiere ligne de y**-1 dans l'ordre invers
c
c calcul de la premiere ligne de la matrice de toeplitz:-(y**-1)*x
j=n+1
do 10 jj=1,n
j=j-1
c(jj)=-ddot(j,a,-1,c(jj),-1)
10 continue
c c contient la premiere ligne du produit ranges dans l'ordre inverse
c
c calcul de x'*(-(y**-1)*x)
i2=0
do 20 i=1,n
w(i2+i)=a(i)*c(n)
20 continue
if(n.eq.1) goto 45
do 40 j=2,n
i1=i2
i2=i2+n
w(i2+1)=a(1)*c(n+1-j)
do 30 i=2,n
w(i2+i)=w(i1+i-1)+a(i)*c(n+1-j)
30 continue
40 continue
c calcul de y'+w
45 i1=-n
do 60 j=1,n
c(j)=0.0d+0
i1=i1+n
do 50 i=j,n
w(i1+i)=w(i1+i)+a(n+1+j-i)
50 continue
60 continue
c(1)=1.0d+0
c w contient la matrice du systeme lineaire et c le second menbre
c
c resolution
call dlslv(w,n,n,c,n,1,w(n*n+1),rcond,ierr,1)
if(ierr.ne.0) goto 80
return
70 ierr=1
return
80 ierr=2
return
end
|
