1 | |
---|
2 | //-*- mode:C++;-*- |
---|
3 | // $Id: lejeune.lib,v 1.6 2005-06-23 15:45:54 cremer Exp $ |
---|
4 | |
---|
5 | |
---|
6 | info=" |
---|
7 | LIBRARY: lejeune.lib Arc space computations |
---|
8 | AUTHOR: Nadine Cremer, nadine.cremer@gmx.de |
---|
9 | [SEE ALSO: <comma-separated words of cross references>] |
---|
10 | [KEYWORDS: <semicolon-separated phrases of index keys>] |
---|
11 | PROCEDURES: |
---|
12 | variables(k,i); creates k*i new var. t,a(1),..,a(i),..,x(1),..,x(i) |
---|
13 | a_z(k); returns kth letter of the alphabet |
---|
14 | tpolys(k,i); creates polyn. a(1)*t+..+a(n)*t^n |
---|
15 | ringchange(i); changes the ring to the one needed in ith step |
---|
16 | plugin_coeffs(i,f) plugs tpolys into f, up to power i |
---|
17 | maxidealstep(i,N); returns ideal needed for contraction in ith step |
---|
18 | N is number of variables of input f |
---|
19 | formaldiff(f,k); computes the formal derivatives D_I with |I|<k |
---|
20 | f_set(f,H); returns the set F corresponding to H as described by |
---|
21 | M. Lejeune |
---|
22 | "; |
---|
23 | |
---|
24 | |
---|
25 | LIB "ring.lib"; // need procedures from these libs |
---|
26 | LIB "general.lib"; |
---|
27 | |
---|
28 | |
---|
29 | |
---|
30 | proc f_set (poly f,intvec H) |
---|
31 | { |
---|
32 | int p; // loop variable |
---|
33 | int m_0=ord(f); |
---|
34 | int b=size(H); |
---|
35 | int c=sum(H,1..b-1); |
---|
36 | if(H[1]!=m_0) // input admissible?! |
---|
37 | { |
---|
38 | "H[1]=ord(f) notwendig!!"; |
---|
39 | return(0); |
---|
40 | } |
---|
41 | for(p=1;p<b;p++) |
---|
42 | { |
---|
43 | if(H[p]<H[p+1]) |
---|
44 | { |
---|
45 | "Unzulaessige Eingabe, H[1]<=...<=H[b] notwendig!"; |
---|
46 | return(0); |
---|
47 | } |
---|
48 | } |
---|
49 | def r=basering; // need that in iteration |
---|
50 | ideal resultf_set,step,I; // save result resp. prelim. result |
---|
51 | for(p=1;p<b-1;p++) // iterating steps |
---|
52 | { |
---|
53 | def tmp=f_set(f,H[1..p+1]); |
---|
54 | setring(tmp); |
---|
55 | def tmp1=resultf_set; |
---|
56 | step=tmp1; |
---|
57 | export(step); |
---|
58 | step;~ |
---|
59 | setring(r); |
---|
60 | } |
---|
61 | def R=formaldiff(f,b-1,c,H[b]); // actual step |
---|
62 | setring R; |
---|
63 | def T=resultdiff; |
---|
64 | ideal resultf_set=T; |
---|
65 | export(resultf_set); |
---|
66 | resultf_set; |
---|
67 | return(R); |
---|
68 | } |
---|
69 | |
---|
70 | |
---|
71 | |
---|
72 | |
---|
73 | |
---|
74 | proc formaldiff (poly f,int i,int a,int k) |
---|
75 | { |
---|
76 | int s,t,v; // loop variables |
---|
77 | int u; |
---|
78 | def R=plugin_coeffs(f,i); // plugs the power series in... |
---|
79 | setring R; // changes the ring |
---|
80 | def Coe=result; |
---|
81 | matrix coe=Coe; // gives the t-coeff. after plugging in |
---|
82 | poly fkv; // need this stuff for the following |
---|
83 | ideal m; // loops... |
---|
84 | ideal m1,m2,J,resultdiff; |
---|
85 | for(v=1;v<=k;v++) // consider the different t-coeff. |
---|
86 | { |
---|
87 | fkv=coe[a+v,1]; |
---|
88 | m=fkv; |
---|
89 | J=fkv; // will save the result in this step |
---|
90 | for(s=1;s<k;s++) |
---|
91 | { |
---|
92 | m1=maxidealstep(i,startvar); // computes the corresponding ideal |
---|
93 | m1=m1^s; |
---|
94 | u=size(m1); |
---|
95 | for(t=1;t<=u;t++) |
---|
96 | { |
---|
97 | m2=contract(m1[t],m); // actual differentiation |
---|
98 | J=J,m2; |
---|
99 | } |
---|
100 | } |
---|
101 | resultdiff=resultdiff,J; |
---|
102 | } |
---|
103 | resultdiff=simplify(resultdiff,2); |
---|
104 | export(resultdiff); // exports the result |
---|
105 | return(R); // return the ring |
---|
106 | |
---|
107 | } |
---|
108 | |
---|
109 | |
---|
110 | |
---|
111 | proc plugin_coeffs (poly f,int i) |
---|
112 | { |
---|
113 | def r=basering; |
---|
114 | def R=ringchange(i); // changes the ring |
---|
115 | setring R; // makes it new basering; |
---|
116 | ideal I=tpolys(i,startvar); |
---|
117 | poly g=imap(r,f); // maps f to new basering |
---|
118 | export(g); // export it |
---|
119 | map h=r,I; // define map according to our purpose |
---|
120 | ideal J=h(f); // gives f with power series plugged in |
---|
121 | export(h); |
---|
122 | matrix result=coeffs(J[1],t); // gives the t-coefficients |
---|
123 | export result; // export it i.o. to use it later on |
---|
124 | return(R); // return ring (ring change!) |
---|
125 | } |
---|
126 | |
---|
127 | |
---|
128 | |
---|
129 | proc ringchange (int i) |
---|
130 | { |
---|
131 | int startvar; |
---|
132 | startvar=nvars(basering); |
---|
133 | export(startvar); |
---|
134 | string str=variables(startvar,i); |
---|
135 | def R=changevar(""+varstr(basering)+",t,"+variables(startvar,i)+"");// change |
---|
136 | return(R); // return the ring, needed in future proc |
---|
137 | } |
---|
138 | |
---|
139 | |
---|
140 | |
---|
141 | proc variables (int k,int i) |
---|
142 | { |
---|
143 | list l; |
---|
144 | int s,u; // loop variables |
---|
145 | string str; |
---|
146 | for (u=1;u<=k;u++) |
---|
147 | { |
---|
148 | for (s=1;s<=i;s++) |
---|
149 | { |
---|
150 | str=""+a_z(u)+"("+string(s)+")"; // creates new variables |
---|
151 | l[(u-1)*i+s]=str; // saves them in a list |
---|
152 | } |
---|
153 | } |
---|
154 | string str1=string(l); // makes the list into a string, |
---|
155 | return(str1); // (needed for ring change) |
---|
156 | } |
---|
157 | |
---|
158 | |
---|
159 | proc a_z (int n) // returns nth letter of the alphabet |
---|
160 | { |
---|
161 | if((n<1)||(n>26)) // input admissible? |
---|
162 | { |
---|
163 | "n must range between 1 and 26!"; |
---|
164 | return(0); |
---|
165 | } |
---|
166 | string s="ring r=0,("+A_Z("a",n)+"),ds;"; |
---|
167 | execute(s); |
---|
168 | return (string(var(n))); |
---|
169 | } |
---|
170 | |
---|
171 | |
---|
172 | |
---|
173 | proc tpolys (int i,int k) // constructs polynomials a(1)*t+... |
---|
174 | { // has to be called from pluin_coeffs |
---|
175 | int s,t; // loop variables |
---|
176 | int v; |
---|
177 | poly sum; |
---|
178 | ideal I; |
---|
179 | for(t=1;t<=k;t++) |
---|
180 | { |
---|
181 | v=(t-1)*i; |
---|
182 | for(s=1;s<=i;s++) |
---|
183 | { |
---|
184 | sum=sum+var(1+k+v+s)*var(k+1)^s; // clumsy: working with "var(1)", |
---|
185 | } // depends on form of basering |
---|
186 | I[t]=sum; |
---|
187 | sum=0; |
---|
188 | } |
---|
189 | return(I); |
---|
190 | } |
---|
191 | |
---|
192 | |
---|
193 | |
---|
194 | |
---|
195 | proc maxidealstep (int i,int N) // returns ideal needed for |
---|
196 | { // differentiation in ith step |
---|
197 | ideal I=var(N+1+i); |
---|
198 | int j; |
---|
199 | for(j=2;j<=N;j++) |
---|
200 | { |
---|
201 | I=I,var(N+1+j*i); |
---|
202 | } |
---|
203 | return(I); |
---|
204 | } |
---|