# source:git/Singular/LIB/lejeune.lib@2b2f12

spielwiese
Last change on this file since 2b2f12 was 2b2f12, checked in by Nadine Cremer <cremer@…>, 18 years ago
*cremer: it's running, now save the results from each step. git-svn-id: file:///usr/local/Singular/svn/trunk@8394 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 5.6 KB
Line
1
2//-*- mode:C++;-*-
3// \$Id: lejeune.lib,v 1.7 2005-06-24 07:28:19 cremer Exp \$
4
5
6info="
7LIBRARY: lejeune.lib  Arc space computations
10[KEYWORDS: <semicolon-separated phrases of index keys>]
11PROCEDURES:
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
25LIB "ring.lib";                       // need procedures from these libs
26LIB "general.lib";
27
28
29
30proc 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);
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,I;                 // save result resp. prelim. result
51  for(p=b-1;p>1;p--)                   // iterating steps
52    {
53      def tmp=f_set(f,intvec(H[1..p]));
54      setring(tmp);
55      def tmp1=resultf_set;
56      export(tmp1);
57      setring(r);
58    }
59  def R=formaldiff(f,b-1,c,H[b]);      // actual step
60  setring R;
61  def T=resultdiff;
62  ideal resultf_set=T;
63  export(resultf_set);
64  resultf_set;
65  return(R);
66}
67
68
69
70
71
72proc formaldiff (poly f,int i,int a,int k)
73{
74  int s,t,v;                          // loop variables
75  int u;
76  def R=plugin_coeffs(f,i);           // plugs the power series in...
77  setring R;                          // changes the ring
78  def Coe=result;
79  matrix coe=Coe;                     // gives the t-coeff. after plugging in
80  poly fkv;                           // need this stuff for the following
81  ideal m;                            // loops...
82  ideal m1,m2,J,resultdiff;
83  for(v=1;v<=k;v++)                   // consider the different t-coeff.
84   {
85     fkv=coe[a+v,1];
86     m=fkv;
87     J=fkv;                           // will save the result in this step
88     for(s=1;s<k;s++)
89       {
90         m1=maxidealstep(i,startvar); // computes the corresponding ideal
91         m1=m1^s;
92         u=size(m1);
93         for(t=1;t<=u;t++)
94          {
95            m2=contract(m1[t],m);     // actual differentiation
96            J=J,m2;
97          }
98       }
99     resultdiff=resultdiff,J;
100   }
101  resultdiff=simplify(resultdiff,2);
102  export(resultdiff);                // exports the result
103  return(R);                         // return the ring
104
105}
106
107
108
109proc plugin_coeffs (poly f,int i)
110{
111  def r=basering;
112  def R=ringchange(i);              // changes the ring
113  setring R;                        // makes it new basering;
114  ideal I=tpolys(i,startvar);
115  poly g=imap(r,f);                 // maps f to new basering
116  export(g);                        // export it
117  map h=r,I;                        // define map according to our purpose
118  ideal J=h(f);                     // gives f with power series plugged in
119  export(h);
120  matrix result=coeffs(J[1],t);     // gives the t-coefficients
121  export result;                    // export it i.o. to use it later on
122  return(R);                        // return ring (ring change!)
123}
124
125
126
127proc ringchange (int i)
128{
129  int startvar;
130  startvar=nvars(basering);
131  export(startvar);
132  variables(startvar,i);
133  def R=changevar(""+varstr(basering)+",t,"+variables(startvar,i)+"");// change
134  return(R);                 // return the ring, needed in future proc
135}
136
137
138
139proc variables (int k,int i)
140{
141  list l;
142  int s,u;                              // loop variables
143  string str;
144  for (u=1;u<=k;u++)
145   {
146     for (s=1;s<=i;s++)
147     {
148       str=""+a_z(u)+"("+string(s)+")"; // creates new variables
149       l[(u-1)*i+s]=str;                // saves them in a list
150     }
151   }
152  string str1=string(l);                // makes the list into a string,
153  return(str1);                         // (needed for ring change)
154}
155
156
157proc a_z (int n)                        // returns nth letter of the alphabet
158{
160   {
161     "n must range between 1 and 26!";
162      return(0);
163   }
164  string s="ring r=0,("+A_Z("a",n)+"),ds;";
165  execute(s);
166  return (string(var(n)));
167}
168
169
170
171proc tpolys (int i,int k)             // constructs polynomials a(1)*t+...
172{                                     // has to be called from pluin_coeffs
173  int s,t;                            // loop variables
174  int v;
175  poly sum;
176  ideal I;
177  for(t=1;t<=k;t++)
178   {
179     v=(t-1)*i;
180     for(s=1;s<=i;s++)
181      {
182        sum=sum+var(1+k+v+s)*var(k+1)^s;    // clumsy: working with "var(1)",
183      }                                     // depends on form of basering
184     I[t]=sum;
185     sum=0;
186   }
187  return(I);
188}
189
190
191
192
193proc maxidealstep (int i,int N)       // returns ideal needed for
194{                                     // differentiation in ith step
195  ideal I=var(N+1+i);
196  int j;
197  for(j=2;j<=N;j++)
198   {
199     I=I,var(N+1+j*i);
200   }
201return(I);
202}
Note: See TracBrowser for help on using the repository browser.