source: git/Singular/LIB/ringgb.lib @ 1b2216

spielwiese
Last change on this file since 1b2216 was 66d68c, checked in by Hans Schoenemann <hannes@…>, 14 years ago
format git-svn-id: file:///usr/local/Singular/svn/trunk@13499 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 10.4 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id$";
3category="Miscellaneous";
4info="
5LIBRARY:  ringgb.lib     Functions for coefficient rings
6AUTHOR:  Oliver Wienand, email: wienand@mathematik.uni-kl.de
7
8KEYWORDS: vanishing polynomial; zeroreduce; polynomial functions; library, ringgb.lib; ringgb.lib, functions for coefficient rings
9
10PROCEDURES:
11 findZeroPoly(f);        finds a vanishing polynomial for reducing f
12 zeroReduce(f);          normal form of f concerning the ideal of vanishing polynomials
13 testZero(poly f);       tests f defines the constant zero function
14 noElements(def r);      the number of elements of the coefficient ring, if of type (integer, ...)
15";
16
17LIB "general.lib";
18///////////////////////////////////////////////////////////////////////////////
19
20proc findZeroPoly (poly f)
21"USAGE:   findZeroPoly(f); f - a polynomial
22RETURN:  zero polynomial with the same leading term as f if exists, otherwise 0
23EXAMPLE: example findZeroPoly; shows an example
24"
25{
26  list data = getZeroCoef(f);
27  if (data[1] == 0)
28  {
29    return(0);
30  }
31  number q = leadcoef(f) / data[1];
32  if (q == 0)
33  {
34    return(0);
35  }
36  poly g = getZeroPolyRaw(data[2]);
37  g = leadmonom(f) / leadmonom(g) * g;
38  return(q * data[1] * g);
39  //return(system("findZeroPoly", f));
40}
41example
42{ "EXAMPLE:"; echo = 2;
43  ring r = (integer, 65536), (y,x), dp;
44  poly f = 1024*x^8*y^2+11264*x^8*y+28672*x^8+45056*x^7*y^2+36864*x^7*y+16384*x^7+40960*x^6*y^2+57344*x^6*y+32768*x^6+30720*x^5*y^2+10240*x^5*y+8192*x^5+35840*x^4*y^2+1024*x^4*y+20480*x^4+30720*x^3*y^2+10240*x^3*y+8192*x^3+4096*x^2*y^2+45056*x^2*y+49152*x^2+40960*x*y^2+57344*x*y+32768*x;
45  findZeroPoly(f);
46}
47
48proc zeroReduce(poly f, list #)
49"USAGE:   zeroReduce(f, [i = 0]); f - a polynomial, i - noise level (if != 0 prints all steps)
50RETURN:  reduced normal form of f modulo zero polynomials
51EXAMPLE: example zeroReduce; shows an example
52"
53{
54   int i = 0;
55   if (size(#) > 0)
56   {
57     i = #[1];
58   }
59   poly h = f;
60   poly n = 0;
61   poly g = findZeroPoly(h);
62   if (i <> 0) {
63     printf("reducing polyfct  : %s", h);
64   }
65   while ( h <> 0 ) {
66      while ( g <> 0 ) {
67         h = h - g;
68         if (i <> 0) {
69            printf(" reduce with: %s", g);
70            printf(" to: %s", h);
71         }
72         g = findZeroPoly(h);
73      }
74      n = lead(h) + n;
75      if (i <> 0) {
76         printf("head irreducible  : %s", lead(h));
77         printf("irreducible start : %s", n);
78         printf("remains to check  : %s", h - lead(h));
79      }
80      h = h - lead(h);
81      g = findZeroPoly(h);
82   }
83   return(n);
84}
85example
86{ "EXAMPLE:"; echo = 2;
87  ring r = (integer, 65536), (y,x), dp;
88  poly f = 1024*x^8*y^2+11264*x^8*y+28672*x^8+45056*x^7*y^2+36864*x^7*y+16384*x^7+40960*x^6*y^2+57344*x^6*y+32768*x^6+30720*x^5*y^2+10240*x^5*y+8192*x^5+35840*x^4*y^2+1024*x^4*y+20480*x^4+30720*x^3*y^2+10240*x^3*y+8192*x^3+4096*x^2*y^2+45056*x^2*y+49152*x^2+40960*x*y^2+57344*x*y+32768*x;
89  zeroReduce(f);
90  kill r;
91  ring r = (integer, 2, 32), (x,y,z), dp;
92  // Polynomial 1:
93  poly p1 = 3795162112*x^3+587202566*x^2*y+2936012853*x*y*z+2281701376*x+548767119*y^3+16777216*y^2+268435456*y*z+1107296256*y+4244635648*z^3+4244635648*z^2+16777216*z;
94  // Polynomial 2:
95  poly p2 = 1647678464*x^3+587202566*x^2*y+2936012853*x*y*z+134217728*x+548767119*y^3+16777216*y^2+268435456*y*z+1107296256*y+2097152000*z^3+2097152000*z^2+16777216*z;
96  zeroReduce(p1-p2);
97}
98
99proc testZero(poly f)
100"USAGE:   testZero(f); f - a polynomial
101RETURN:  returns 1 if f is zero as a function and otherwise a counterexample as a list [f(x_1, ..., x_n), x_1, ..., x_n]
102EXAMPLE: example testZero; shows an example
103"
104{
105  poly g;
106  int j;
107  bigint i = 0;
108  bigint modul = noElements(basering);
109  printf("Teste %s Belegungen ...", modul^nvars(basering));
110  for (; i < modul^nvars(basering); i = i + 1)
111  {
112    if ((i + 1) % modul^(nvars(basering)/2) == 0)
113    {
114      printf("bisher: %s", i);
115    }
116    g = f;
117    for (j = 1; j <= nvars(basering); j++)
118    {
119      g = subst(g, var(j), number((i / modul^(j-1)) % modul));
120    }
121    if (g != 0)
122    {
123      list counter = g;
124      for (j = 1; j <= nvars(basering); j++)
125      {
126        counter = insert(counter, (i / modul^(j-1)) % modul);
127      }
128      return(counter);
129    }
130  }
131  return(1);
132}
133example
134{ "EXAMPLE:"; echo = 2;
135  ring r = (integer, 12), (y,x), dp;
136  poly f = 1024*x^8*y^2+11264*x^8*y+28672*x^8+45056*x^7*y^2+36864*x^7*y+16384*x^7+40960*x^6*y^2+57344*x^6*y+32768*x^6+30720*x^5*y^2+10240*x^5*y+8192*x^5+35840*x^4*y^2+1024*x^4*y+20480*x^4+30720*x^3*y^2+10240*x^3*y+8192*x^3+4096*x^2*y^2+45056*x^2*y+49152*x^2+40960*x*y^2+57344*x*y+32768*x;
137  zeroReduce(f);
138  testZero(f);
139  poly g = findZeroPoly(x2y3);
140  g;
141  testZero(g);
142}
143
144proc noElements(def r)
145"USAGE:   noElements(r); r - a ring with a finite coefficient ring of type integer
146RETURN:  returns the number of elements of the coefficient ring of r
147EXAMPLE: example noElements; shows an example
148"
149
150{
151  list l = ringlist(basering);
152  return(l[1][2][1]^l[1][2][2]);
153}
154example
155{ "EXAMPLE:"; echo = 2;
156  ring r = (integer, 233,6), (y,x), dp;
157  noElements(r);
158}
159
160static proc getZeroCoef(poly f)
161{
162  if (f == 0)
163  {
164    return(0, leadexp(f))
165  }
166  list data = sort(leadexp(f));
167  intvec exp = data[1];
168  intvec index = data[2];
169  intvec nec = 0:size(exp);
170  int i = 1;
171  int j = 2;
172  bigint g;
173  bigint G = 1;
174  bigint modul = noElements(basering);
175  bigint B = modul;
176
177  for (; exp[i] < 2; i++) {if (i == size(exp)) break;}
178  for (; i <= size(exp); i++)
179  {
180    g = gcd(B, G);
181    G = G * g;
182    B = B / g;
183    if (g != 1)
184    {
185      nec[index[i]] = j - 1;
186    }
187    if (B == 1)
188    {
189      return(B, nec);
190    }
191    for (; j <= exp[i]; j++)
192    {
193      g = gcd(B, bigint(j));
194      G = G * g;
195      B = B / g;
196      if (g != 1)
197      {
198        nec[index[i]] = j;
199      }
200      if (B == 1)
201      {
202        return(B, nec);
203      }
204    }
205  }
206  if (B == modul)
207  {
208    nec = 0;
209    return(0, nec);
210  }
211  return(B, nec);
212}
213
214static proc getZeroPolyRaw(intvec fexp)
215{
216  list data = sort(fexp);
217  intvec exp = data[1];
218  intvec index = data[2];
219  int j = 0;
220  poly res = 1;
221  poly tillnow = 1;
222  int i = 1;
223  for (; exp[i] < 2; i++) {if (i == size(exp)) break;}
224  for (; i <= size(exp); i++)
225  {
226    for (; j < exp[i]; j++)
227    {
228      tillnow = tillnow * (var(1) - j);
229    }
230    res = res * subst(tillnow, var(1), var(index[i]));
231  }
232  return(res);
233}
234
235static proc getZeroPoly(poly f)
236{
237  list data = getZeroCoef(f);
238  poly g = getZeroPolyRaw(data[2]);
239  g = leadmonom(f) / leadmonom(g) * g;
240  return(data[1] * g);
241}
242
243static proc findZeroPolyWrap (poly f)
244"USAGE:   findZeroPolyWrap(f); f - a polynomial
245RETURN:  zero polynomial with the same leading term as f if exists, otherwise 0
246NOTE:    just a wrapper, work only in Z/2^n with n < int_machine_size - 1
247EXAMPLE: example findZeroPoly; shows an example
248"
249{
250  return(system("findZeroPoly", f));
251}
252example
253{ "EXAMPLE:"; echo = 2;
254  ring r = (integer, 2, 16), (y,x), dp;
255  poly f = 1024*x^8*y^2+11264*x^8*y+28672*x^8+45056*x^7*y^2+36864*x^7*y+16384*x^7+40960*x^6*y^2+57344*x^6*y+32768*x^6+30720*x^5*y^2+10240*x^5*y+8192*x^5+35840*x^4*y^2+1024*x^4*y+20480*x^4+30720*x^3*y^2+10240*x^3*y+8192*x^3+4096*x^2*y^2+45056*x^2*y+49152*x^2+40960*x*y^2+57344*x*y+32768*x;
256  findZeroPoly(f);
257}
258
259///////////////////////////////////////////////////////////////////////////////
260
261/*
262                           Examples:
263
264
265// POLYNOMIAL EXAMPLES (Singular ready)
266// ===================
267//
268// For each of the examples below, there are three equivalent polynomials. 'm' indicates the bit-widths of the
269// input/output variables. For some of the polynomials, I have attached the RTL as well.
270//
271//
272// 1) VOLTERRA MODELS:
273//
274//        A) CUBIC FILTER: (m = 32, 3 Vars)
275
276LIB "ringgb.lib";
277ring r = (integer, 2, 32), (x,y,z), dp;
278poly p1 = 3795162112*x^3+587202566*x^2*y+2936012853*x*y*z+2281701376*x+548767119*y^3+16777216*y^2+268435456*y*z \
279         +1107296256*y+4244635648*z^3+4244635648*z^2+16777216*z;
280poly p2 = 1647678464*x^3+587202566*x^2*y+2936012853*x*y*z+134217728*x+548767119*y^3+16777216*y^2+268435456*y*z \
281         +1107296256*y+2097152000*z^3+2097152000*z^2+16777216*z;
282poly p3 = 1647678464*x^3+587202566*x^2*y+2936012853*x*y*z+134217728*x+548767119*y^3+16777216*y^2+268435456*y*z \
283         +1107296256*y+2097152000*z^3+2097152000*z^2+16777216*z;
284zeroReduce(p1-p2);
285zeroReduce(p1-p3);
286zeroReduce(p2-p3);
287
288//        B) DEGREE-4 FILTER: (m=16 , 3 Vars)
289
290LIB "ringgb.lib";
291ring r = (integer, 2, 16), (x,y,z), dp;
292poly p1 = 16384*x^4+y^4+57344*z^4+64767*x*y^3+16127*y^2*z^2+8965*x^3*z+19275*x^2*y*z+51903*x*y*z+32768*x^2*y  \
293         +40960*z^2+32768*x*y^2+49152*x^2+4869*y;
294poly p2 = 8965*x^3*z+19275*x^2*y*z+31999*x*y^3+51903*x*y*z+32768*x*y+y^4+32768*y^3+16127*y^2*z^2+32768*y^2 \
295         +4869*y+57344*z^4+40960*z^2;
296poly p3 = 8965*x^3*z+19275*x^2*y*z+31999*x*y^3+51903*x*y*z+32768*x*y+y^4+16127*y^2*z^2+4869*y+16384*z^3+16384*z;
297zeroReduce(p1-p2);
298zeroReduce(p1-p3);
299zeroReduce(p2-p3);
300
301
302// 2) Savitzsky Golay filter(m=16,5 Vars)
303
304LIB "ringgb.lib";
305ring r = (integer, 2, 16), (v,w,x,y,z), dp;
306poly p1 = 25000*v^2*y+37322*v^2+22142*v*w*z+50356*w^3+58627*w^2+17797*w+17797*x^3+62500*x^2*z+41667*x \
307         +22142*y^3+23870*y^2+59464*y+41667*z+58627;
308poly p2 = 25000*v^2*y+4554*v^2+22142*v*w*z+32768*v+17588*w^3+25859*w^2+17797*w+17797*x^3+29732*x^2*z+32768*x^2 \
309         +32768*x*z+8899*x+22142*y^3+23870*y^2+59464*y+41667*z+58627;
310poly p3 = 25000*v^2*y+4554*v^2+22142*v*w*z+32768*v+17588*w^3+25859*w^2+17797*w+17797*x^3+29732*x^2*z+32768*x*z \
311         +41667*x+22142*y^3+23870*y^2+59464*y+41667*z+58627;
312zeroReduce(p1-p2);
313zeroReduce(p1-p3);
314zeroReduce(p2-p3);
315
316
317// 3) Anti-alias filter:(m=16, 1 Var)
318
319LIB "ringgb.lib";
320ring r = (integer, 2, 16), c, dp;
321poly p1 = 156*c^6+62724*c^5+17968*c^4+18661*c^3+43593*c^2+40224*c+13281;
322poly p2 = 156*c^6+5380*c^5+1584*c^4+43237*c^3+27209*c^2+40224*c+13281;
323poly p3 = 156*c^6+5380*c^5+1584*c^4+10469*c^3+27209*c^2+7456*c+13281;
324zeroReduce(p1-p2);
325zeroReduce(p1-p3);
326zeroReduce(p2-p3);
327
328
329// 4) PSK:(m=16, 2 Var)
330
331LIB "ringgb.lib";
332ring r = (integer, 2, 16), (x,y), dp;
333poly p1 = 4166*x^4+16666*x^3*y+25000*x^2*y^2+15536*x^2+16666*x*y^4+31072*x*y+4166*y^4+15536*y^2+34464;
334poly p2 = 4166*x^4+16666*x^3*y+8616*x^2*y^2+16384*x^2*y+15536*x^2+282*x*y^4+47456*x*y+53318*y^4+31920*y^2+34464;
335poly p3 = 4166*x^4+16666*x^3*y+8616*x^2*y^2+16384*x^2*y+15536*x^2+282*x*y^4+47456*x*y+4166*y^4+15536*y^2+34464;
336zeroReduce(p1-p2);
337zeroReduce(p1-p3);
338zeroReduce(p2-p3);
339
340// Ref: A. Peymandoust G. De Micheli, “Application of Symbolic Computer Algebra in High-Level Data-Flow
341// Synthesis,” IEEE Transactions on CAD/ICAS, Vol. 22, No. 9, September 2003, pp.1154-1165.
342
343*/
Note: See TracBrowser for help on using the repository browser.