# source:git/Singular/LIB/ncfactor.lib@0a2f7d

spielwiese
Last change on this file since 0a2f7d was 0a2f7d, checked in by Viktor Levandovskyy <levandov@…>, 13 years ago
*levandov: checklib changes of libs form Aachen git-svn-id: file:///usr/local/Singular/svn/trunk@13356 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to 100755
File size: 42.1 KB
Line
1////////////////////////////////////////////////////////////
2version = "$Id: ncfactor.lib, v 1.00 2010/02/15 10:06 heinle Exp$";
3category="Noncommutative";
4info="
5LIBRARY: ncfactor.lib   Tools for factorization in the first Weyl algebra
6AUTHORS: Albert Heinle,     albert.heinle@rwth-aachen.de
7
8PROCEDURES:
9facfirstwa(h);     factorization in the first Weyl algebra
10";
11
12LIB "general.lib";
13LIB "nctools.lib";
14LIB "involut.lib";
15
16//static homogfac(h);         computes one factorization of a homogeneous polynomial h
17//static homogfac_all(h);     computes all factorizations of a homogeneous polynomial h
18//static facwa(h);            computes all factorizations of an inhomogeneous polynomial h
19/* ring R = 0,(x,y),Ws(-1,1); */
20/* def r = nc_algebra(1,1); */
21/* setring(r); */
22
23////////////////////////////////////////////////////////////////////////////////////////////////////
24//==================================================*
25//deletes double-entries in a list of factorization
26//without evaluating the product.
27static proc delete_dublicates_noteval(list l)
28{//proc delete_dublicates_noteval
29     list result= l;
30     int j; int k; int i;
31     int deleted = 0;
32     int is_equal;
33     for (i = 1; i<= size(l); i++)
34     {//Iterate over the different factorizations
35          for (j = i+1; j<= size(l); j++)
36          {//Compare the i'th factorization to the j'th
37               if (size(l[i])!= size(l[j]))
38               {//different sizes => not equal
39                    j++;
40                    continue;
41               }//different sizes => not equal
42               is_equal = 1;
43               for (k = 1; k <= size(l[i]);k++)
44               {//Compare every entry
45                    if (l[i][k]!=l[j][k])
46                    {
47                         is_equal = 0;
48                         break;
49                    }
50               }//Compare every entry
51               if (is_equal == 1)
52               {//Delete this entry, because there is another equal one int the list
53                    result = delete(result, i-deleted);
54                    deleted = deleted+1;
55                    break;
56               }//Delete this entry, because there is another equal one int the list
57          }//Compare the i'th factorization to the j'th
58     }//Iterate over the different factorizations
59     return(result);
60}//proc delete_dublicates_noteval
61
62//==================================================
63//deletes the double-entries in a list with
64//evaluating the products
65static proc delete_dublicates_eval(list l)
66{//proc delete_dublicates_eval
67     list result=l;
68     int j; int k; int i;
69     int deleted = 0;
70     int is_equal;
71     for (i = 1; i<= size(result); i++)
72     {//Iterating over all elements in result
73          for (j = i+1; j<= size(result); j++)
74          {//comparing with the other elements
75               if (product(result[i]) == product(result[j]))
76               {//There are two equal results; throw away that one with the smaller size
77                    if (size(result[i])>=size(result[j]))
78                    {//result[i] has more entries
79                         result = delete(result,j);
80                         continue;
81                    }//result[i] has more entries
82                    else
83                    {//result[j] has more entries
84                         result = delete(result,i);
85                         i--;
86                         break;
87                    }//result[j] has more entries
88               }//There are two equal results; throw away that one with the smaller size
89          }//comparing with the other elements
90     }//Iterating over all elements in result
91     return(result);
92}//proc delete_dublicates_eval
93
94
95//==================================================*
96//given a list of factors g and a desired size nof, the following
97//procedure combines the factors, such that we recieve a
98//list of the length nof.
99static proc combinekfinlf(list g, int nof, intvec limits) //nof stands for "number of factors"
100{//Procedure combinekfinlf
101     list result;
102     int i; int j; int k; //iteration variables
103     list fc; //fc stands for "factors combined"
104     list temp; //a temporary store for factors
105     def nofgl = size(g); //nofgl stands for "number of factors of the given list"
106     if (nofgl == 0)
107     {//g was the empty list
108          return(result);
109     }//g was the empty list
110     if (nof <= 0)
111     {//The user wants to recieve a negative number or no element as a result
112          return(result);
113     }//The user wants to recieve a negative number or no element as a result
114     if (nofgl == nof)
115     {//There are no factors to combine
116          if (limitcheck(g,limits))
117          {
118               result = result + list(g);
119          }
120          return(result);
121     }//There are no factors to combine
122     if (nof == 1)
123     {//User wants to get just one factor
124          if (limitcheck(list(product(g)),limits))
125          {
126               result = result + list(list(product(g)));
127          }
128          return(result);
129     }//User wants to get just one factor
130     for (i = nof; i > 1; i--)
131     {//computing the possibilities that have at least one original factor from g
132          for (j = i; j>=1; j--)
133          {//shifting the window of combinable factors to the left
134               fc = combinekfinlf(list(g[(j)..(j+nofgl - i)]),nof - i + 1,limits); //fc stands for "factors combined"
135               for (k = 1; k<=size(fc); k++)
136               {//iterating over the different solutions of the smaller problem
137                    if (j>1)
138                    {//There are g_i before the combination
139                         if (j+nofgl -i < nofgl)
140                         {//There are g_i after the combination
141                              temp = list(g[1..(j-1)]) + fc[k] + list(g[(j+nofgl-i+1)..nofgl]);
142                         }//There are g_i after the combination
143                         else
144                         {//There are no g_i after the combination
145                              temp = list(g[1..(j-1)]) + fc[k];
146                         }//There are no g_i after the combination
147                    }//There are g_i before the combination
148                    if (j==1)
149                    {//There are no g_i before the combination
150                         if (j+ nofgl -i <nofgl)
151                         {//There are g_i after the combination
152                              temp = fc[k]+ list(g[(j + nofgl - i +1)..nofgl]);
153                         }//There are g_i after the combination
154                    }//There are no g_i before the combination
155                    if (limitcheck(temp,limits))
156                    {
157                         result = result + list(temp);
158                    }
159               }//iterating over the different solutions of the smaller problem
160          }//shifting the window of combinable factors to the left
161     }//computing the possibilities that have at least one original factor from g
162     for (i = 2; i<=nofgl/nof;i++)
163     {//getting the other possible results
164          result = result + combinekfinlf(list(product(list(g[1..i])))+list(g[(i+1)..nofgl]),nof,limits);
165     }//getting the other possible results
166     result = delete_dublicates_noteval(result);
167     return(result);
168}//Procedure combinekfinlf
169
170
171//==================================================*
172//merges two sets of factors ignoring common
173//factors
174static proc merge_icf(list l1, list l2, intvec limits)
175{//proc merge_icf
176     list g;
177     list f;
178     int i; int j;
179     if (size(l1)==0)
180     {
181          return(list());
182     }
183     if (size(l2)==0)
184     {
185          return(list());
186     }
187     if (size(l2)<=size(l1))
188     {//l1 will be our g, l2 our f
189          g = l1;
190          f = l2;
191     }//l1 will be our g, l2 our f
192     else
193     {//l1 will be our f, l2 our g
194          g = l2;
195          f = l1;
196     }//l1 will be our f, l2 our g
197     def result = combinekfinlf(g,size(f),limits);
198     for (i = 1 ; i<= size(result); i++)
199     {//Adding the factors of f to every possibility listed in temp
200          for (j = 1; j<= size(f); j++)
201          {
202               result[i][j] = result[i][j]+f[j];
203          }
204          if(!limitcheck(result[i],limits))
205          {
206               result = delete(result,i);
207               i--;
208          }
209     }//Adding the factors of f to every possibility listed in temp
210     return(result);
211}//proc merge_icf
212
213//==================================================*
214//merges two sets of factors with respect to the occurrence
215//of common factors
216static proc merge_cf(list l1, list l2, intvec limits)
217{//proc merge_cf
218     list g;
219     list f;
220     int i; int j;
221     list pre;
222     list post;
223     list candidate;
224     list temp;
225     int temppos;
226     if (size(l1)==0)
227     {//the first list is empty
228          return(list());
229     }//the first list is empty
230     if(size(l2)==0)
231     {//the second list is empty
232          return(list());
233     }//the second list is empty
234     if (size(l2)<=size(l1))
235     {//l1 will be our g, l2 our f
236          g = l1;
237          f = l2;
238     }//l1 will be our g, l2 our f
239     else
240     {//l1 will be our f, l2 our g
241          g = l2;
242          f = l1;
243     }//l1 will be our f, l2 our g
244     list M;
245     for (i = 2; i<size(f); i++)
246     {//finding common factors of f and g...
247        for (j=2; j<size(g);j++)
248        {//... with g
249             if (f[i] == g[j])
250             {//we have an equal pair
251                  M = M + list(list(i,j));
252             }//we have an equal pair
253        }//... with g
254     }//finding common factors of f and g...
255     if (g[1]==f[1])
256     {//Checking for the first elements to be equal
257          M = M + list(list(1,1));
258     }//Checking for the first elements to be equal
259     if (g[size(g)]==f[size(f)])
260     {//Checking for the last elements to be equal
261          M = M + list(list(size(f),size(g)));
262     }//Checking for the last elements to be equal
263     list result;//= list(list());
264     while(size(M)>0)
265     {//set of equal pairs is not empty
266          temp = M[1];
267          temppos = 1;
268          for (i = 2; i<=size(M); i++)
269          {//finding the minimal element of M
270               if (M[i][1]<=temp[1])
271               {//a possible candidate that is smaller than temp could have been found
272                    if (M[i][1]==temp[1])
273                    {//In this case we must look at the second number
274                         if (M[i][2]< temp[2])
275                         {//the candidate is smaller
276                              temp = M[i];
277                              temppos = i;
278                         }//the candidate is smaller
279                    }//In this case we must look at the second number
280                    else
281                    {//The candidate is definately smaller
282                         temp = M[i];
283                         temppos = i;
284                    }//The candidate is definately smaller
285               }//a possible candidate that is smaller than temp could have been found
286          }//finding the minimal element of M
287          M = delete(M, temppos);
288          if(temp[1]>1)
289          {//There are factors to combine before the equal factor
290               if (temp[1]<size(f))
291               {//The most common case
292                    //first the combinations ignoring common factors
293                    pre = merge_icf(list(f[1..(temp[1]-1)]),list(g[1..(temp[2]-1)]),limits);
294                    post = merge_icf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits);
295                    for (i = 1; i <= size(pre); i++)
296                    {//all possible pre's...
297                         for (j = 1; j<= size(post); j++)
298                         {//...combined with all possible post's
299                              candidate = pre[i]+list(f[temp[1]])+post[j];
300                              if (limitcheck(candidate,limits))
301                              {
302                                   result = result + list(candidate);
303                              }
304                         }//...combined with all possible post's
305                    }//all possible pre's...
306                    //Now the combinations with respect to common factors
307                    post = merge_cf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits);
308                    if (size(post)>0)
309                    {//There are factors to combine
310                         for (i = 1; i <= size(pre); i++)
311                         {//all possible pre's...
312                              for (j = 1; j<= size(post); j++)
313                              {//...combined with all possible post's
314                                   candidate= pre[i]+list(f[temp[1]])+post[j];
315                                   if (limitcheck(candidate,limits))
316                                   {
317                                        result = result + list(candidate);
318                                   }
319                              }//...combined with all possible post's
320                         }//all possible pre's...
321                    }//There are factors to combine
322               }//The most common case
323               else
324               {//the last factor is the common one
325                    pre = merge_icf(list(f[1..(temp[1]-1)]),list(g[1..(temp[2]-1)]),limits);
326                    for (i = 1; i<= size(pre); i++)
327                    {//iterating over the possible pre-factors
328                         candidate = pre[i]+list(f[temp[1]]);
329                         if (limitcheck(candidate,limits))
330                         {
331                              result = result + list(candidate);
332                         }
333                    }//iterating over the possible pre-factors
334               }//the last factor is the common one
335          }//There are factors to combine before the equal factor
336          else
337          {//There are no factors to combine before the equal factor
338               if (temp[1]<size(f))
339               {//Just a check for security
340                    //first without common factors
341                    post=merge_icf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits);
342                    for (i = 1; i<=size(post); i++)
343                    {
344                         candidate = list(f[temp[1]])+post[i];
345                         if (limitcheck(candidate,limits))
346                         {
347                              result = result + list(candidate);
348                         }
349                    }
350                    //Now with common factors
351                    post = merge_cf(list(f[(temp[1]+1)..size(f)]),list(g[(temp[2]+1..size(g))]),limits);
352                    if(size(post)>0)
353                    {//we could find other combinations
354                         for (i = 1; i<=size(post); i++)
355                         {
356                              candidate = list(f[temp[1]])+post[i];
357                              if (limitcheck(candidate,limits))
358                              {
359                                   result = result + list(candidate);
360                              }
361                         }
362                    }//we could find other combinations
363               }//Just a check for security
364          }//There are no factors to combine before the equal factor
365     }//set of equal pairs is not empty
366     return(result);
367}//proc merge_cf
368
369
370//==================================================*
371//merges two sets of factors
372
373static proc mergence(list l1, list l2, intvec limits)
374{//Procedure mergence
375     list g;
376     list f;
377     int l; int k;
378     list F;
379     if (size(l2)<=size(l1))
380     {//l1 will be our g, l2 our f
381          g = l1;
382          f = l2;
383     }//l1 will be our g, l2 our f
384     else
385     {//l1 will be our f, l2 our g
386          g = l2;
387          f = l1;
388     }//l1 will be our f, l2 our g
389     list result;
390     for (l = size(f); l>=1; l--)
391     {//all possibilities to combine the factors of f
392          F = combinekfinlf(f,l,limits);
393          for (k = 1; k<= size(F);k++)
394          {//for all possibilities of combinations of the factors of f
395               result = result + merge_cf(F[k],g,limits);
396               result = result + merge_icf(F[k],g,limits);
397          }//for all possibilities of combinations of the factors of f
398     }//all possibilities to combine the factors of f
399     return(result);
400}//Procedure mergence
401
402
403//==================================================
404//Checks, whether a list of factors doesn't exceed the given limits
405static proc limitcheck(list g, intvec limits)
406{//proc limitcheck
407     int i;
408     if (size(limits)!=3)
409     {//check the input
410          return(0);
411     }//check the input
412     if(size(g)==0)
413     {
414          return(0);
415     }
416     def prod = product(g);
417     def limg = intvec(deg(prod,intvec(1,1)) ,deg(prod,intvec(1,0)),deg(prod,intvec(0,1)));
418     for (i = 1; i<=size(limg);i++)
419     {//the final check
420          if(limg[i]>limits[i])
421          {
422               return(0);
423          }
424     }//the final check
425     return(1);
426}//proc limitcheck
427
428
429//==================================================*
430//one factorization of a homogeneous polynomial
431//in the first Weyl Algebra
432static proc homogfac(poly h)
433"USAGE: homogfac(h); h is a homogeneous polynomial in the first Weyl algebra with respect to the weight vector [-1,1]
434RETURN: list
435PURPOSE: Computes a factorization of a homogeneous polynomial h with respect to the weight vector [-1,1] in the first Weyl algebra
436THEORY: @code{homogfac} returns a list with a factorization of the given, homogeneous polynomial. If the degree of the polynomial is k with k positive, the last k entries in the output list are the second variable. If k is positive, the last k entries will be x. The other entries will be irreducible polynomials of degree zero or 1 resp. -1.
437EXAMPLE: example homogfac; shows examples
439"{//proc homogfac
440     def r = basering;
441     poly hath;
442     int i; int j;
443     if (!homogwithorder(h,intvec(-1,1)))
444     {//The given polynomial is not homogeneous
445          return(list());
446     }//The given polynomial is not homogeneous
447     if (h==0)
448     {
449          return(list(0));
450     }
451     list result;
452     int m = deg(h,intvec(-1,1));
453     if (m!=0)
454     {//The degree is not zero
455          if (m <0)
456          {//There are more x than y
457               hath = lift(var(1)^(-m),h)[1,1];
458               for (i = 1; i<=-m; i++)
459               {
460                    result = result + list(var(1));
461               }
462          }//There are more x than y
463          else
464          {//There are more y than x
465               hath = lift(var(2)^m,h)[1,1];
466               for (i = 1; i<=m;i++)
467               {
468                    result = result + list(var(2));
469               }
470          }//There are more y than x
471     }//The degree is not zero
472     else
473     {//The degree is zero
474          hath = h;
475     }//The degree is zero
476     //beginning to transform x^i*y^i in teta(teta-1)...(teta-i+1)
477     list mons;
478     for(i = 1; i<=size(hath);i++)
479     {//Putting the monomials in a list
480          mons = mons+list(hath[i]);
481     }//Putting the monomials in a list
482     ring tempRing = 0,(x,y,teta),dp;
483     setring tempRing;
484     map tetamap = r,x,y;
485     list mons = tetamap(mons);
486     poly entry;
487     for (i = 1; i<=size(mons);i++)
488     {//transforming the monomials as monomials in theta
490          for (j = 0; j<leadexp(mons[i])[2];j++)
491          {
492               entry = entry * (teta-j);
493          }
494          mons[i] = entry;
495     }//transforming the monomials as monomials in theta
496     list azeroresult = factorize(sum(mons));
497     list azeroresult_return_form;
498     for (i = 1; i<=size(azeroresult[1]);i++)
499     {//rewrite the result of the commutative factorization
500          for (j = 1; j <= azeroresult[2][i];j++)
501          {
502               azeroresult_return_form = azeroresult_return_form + list(azeroresult[1][i]);
503          }
504     }//rewrite the result of the commutative factorization
505     setring(r);
506     map finalmap = tempRing,var(1),var(2),var(1)*var(2);
507     list tempresult = finalmap(azeroresult_return_form);
508     for (i = 1; i<=size(tempresult);i++)
509     {//factorizations of theta resp. theta +1
510          if(tempresult[i]==var(1)*var(2))
511          {
512               tempresult = insert(tempresult,var(1),i-1);
513               i++;
514               tempresult[i]=var(2);
515          }
516          if(tempresult[i]==var(2)*var(1))
517          {
518               tempresult = insert(tempresult,var(2),i-1);
519               i++;
520               tempresult[i]=var(1);
521          }
522     }//factorizations of theta resp. theta +1
523     result = tempresult+result;
524     return(result);
525}//proc homogfac
526/* example */
527/* { */
528/*      "EXAMPLE:";echo=2; */
529/*      ring R = 0,(x,y),Ws(-1,1); */
530/*      def r = nc_algebra(1,1); */
531/*      setring(r); */
532/*      poly h = (x^2*y^2+1)*(x^4); */
533/*      homogfac(h); */
534/* } */
535
536//==================================================
537//Computes all possible homogeneous factorizations
538static proc homogfac_all(poly h)
539"USAGE: homogfac_all(h); h is a homogeneous polynomial in the first Weyl algebra with respect to the weight vector [-1,1]
540RETURN: list
541PURPOSE: Computes all factorizations of a homogeneous polynomial h with respect to the weight vector [-1,1] in the first Weyl algebra
542THEORY: @code{homogfac} returns a list with all factorization of the given, homogeneous polynomial. It uses the output of homogfac and permutes its entries with respect to the commutation rule. Furthermore, if a factor of degree zero is irreducible in K[\theta], but reducible in the first Weyl algebra, the permutations of this element with the other entries will also be computed.
543EXAMPLE: example homogfac; shows examples
545"{//proc HomogFacAll
546     if (deg(h,intvec(1,1)) <= 0 )
547     {//h is a constant
548          return(list(list(h)));
549     }//h is a constant
550     def r = basering;
551     list one_hom_fac; //stands for one homogeneous factorization
552     int i; int j; int k;
553     //Compute again a homogeneous factorization
554     one_hom_fac = homogfac(h);
555     if (size(one_hom_fac) == 0)
556     {//there is no homogeneous factorization or the polynomial was not homogeneous
557          return(list());
558     }//there is no homogeneous factorization or the polynomial was not homogeneous
559     //divide list in A0-Part and a list of x's resp. y's
560     list list_not_azero = list();
561     list list_azero;
562     list k_factor;
563     int is_list_not_azero_empty = 1;
564     int is_list_azero_empty = 1;
565     k_factor = list(one_hom_fac[1]);
566     if (absValue(deg(h,intvec(-1,1)))<size(one_hom_fac)-1)
567     {//There is a nontrivial A_0-part
568          list_azero = one_hom_fac[2..(size(one_hom_fac)-absValue(deg(h,intvec(-1,1))))];
569          is_list_azero_empty = 0;
570     }//There is a nontrivial A_0 part
571     for (i = 1; i<=size(list_azero)-1;i++)
572     {//in homogfac, we factorized theta, and this will be made undone
573          if (list_azero[i] == var(1))
574          {
575               if (list_azero[i+1]==var(2))
576               {
577                    list_azero[i] = var(1)*var(2);
578                    list_azero = delete(list_azero,i+1);
579               }
580          }
581          if (list_azero[i] == var(2))
582          {
583               if (list_azero[i+1]==var(1))
584               {
585                    list_azero[i] = var(2)*var(1);
586                    list_azero = delete(list_azero,i+1);
587               }
588          }
589     }//in homogfac, we factorized theta, and this will be made undone
590     if(deg(h,intvec(-1,1))!=0)
591     {//list_not_azero is not empty
592          list_not_azero = one_hom_fac[(size(one_hom_fac)-absValue(deg(h,intvec(-1,1)))+1)..size(one_hom_fac)];
593          is_list_not_azero_empty = 0;
594     }//list_not_azero is not empty
595     //Map list_azero in K[theta]
596     ring tempRing = 0,(x,y,theta), dp;
597     setring(tempRing);
598     poly entry;
599     map thetamap = r,x,y;
600     if(!is_list_not_azero_empty)
601     {//Mapping in Singular is only possible, if the list before contained at least one element of the other ring
602          list list_not_azero = thetamap(list_not_azero);
603     }//Mapping in Singular is only possible, if the list before contained at least one element of the other ring
604     if(!is_list_azero_empty)
605     {//Mapping in Singular is only possible, if the list before contained at least one element of the other ring
606          list list_azero= thetamap(list_azero);
607     }//Mapping in Singular is only possible, if the list before contained at least one element of the other ring
608     list k_factor = thetamap(k_factor);
609     list tempmons;
610     for(i = 1; i<=size(list_azero);i++)
611     {//rewrite the polynomials in A1 as polynomials in K[theta]
612          tempmons = list();
613          for (j = 1; j<=size(list_azero[i]);j++)
614          {
615               tempmons = tempmons + list(list_azero[i][j]);
616          }
617          for (j = 1 ; j<=size(tempmons);j++)
618          {
620               for (k = 0; k < leadexp(tempmons[j])[2];k++)
621               {
622                    entry = entry*(theta-k);
623               }
624               tempmons[j] = entry;
625          }
626          list_azero[i] = sum(tempmons);
627     }//rewrite the polynomials in A1 as polynomials in K[theta]
628     //Compute all permutations of the A0-part
629     list result;
630     int shift_sign;
631     int shift;
632     poly shiftvar;
633     if (size(list_not_azero)!=0)
634     {//Compute all possibilities to permute the x's resp. the y's in the list
635          if (list_not_azero[1] == x)
636          {//h had a negative weighted degree
637               shift_sign = 1;
638               shiftvar = x;
639          }//h had a negative weighted degree
640          else
641          {//h had a positive weighted degree
642               shift_sign = -1;
643               shiftvar = y;
644          }//h had a positive weighted degree
645          result = permpp(list_azero + list_not_azero);
646          for (i = 1; i<= size(result); i++)
648               shift = 0;
649               for (j=1; j<=size(result[i]);j++)
650               {
651                    if (result[i][j]==shiftvar)
652                    {
653                         shift = shift + shift_sign;
654                    }
655                    else
656                    {
657                         result[i][j] = subst(result[i][j],theta,theta + shift);
658                    }
659               }
661     }//Compute all possibilities to permute the x's resp. the y's in the list
662     else
663     {//The result is just all the permutations of the a_0-part
664          result = permpp(list_azero);
665     }//The result is just all the permutations of the a_0 part
666     if (size(result)==0)
667     {
668          return(result);
669     }
670     //Now we are going deeper and search for theta resp. theta + 1, substitute them by xy resp. yx and go on permuting
671     int found_theta;
672     int thetapos;
673     list leftpart;
674     list rightpart;
675     list lparts;
676     list rparts;
678     for (i = 1; i<=size(result) ; i++)
679     {//checking every entry of result for theta or theta +1
680          found_theta = 0;
681          for(j=1;j<=size(result[i]);j++)
682          {
683               if (result[i][j]==theta)
684               {//the jth entry is theta and can be written as x*y
685                    thetapos = j;
686                    result[i]= insert(result[i],x,j-1);
687                    j++;
688                    result[i][j] = y;
689                    found_theta = 1;
690                    break;
691               }//the jth entry is theta and can be written as x*y
692               if(result[i][j] == theta +1)
693               {
694                    thetapos = j;
695                    result[i] = insert(result[i],y,j-1);
696                    j++;
697                    result[i][j] = x;
698                    found_theta = 1;
699                    break;
700               }
701          }
702          if (found_theta)
703          {//One entry was theta resp. theta +1
704               leftpart = result[i];
705               leftpart = leftpart[1..thetapos];
706               rightpart = result[i];
707               rightpart = rightpart[(thetapos+1)..size(rightpart)];
708               lparts = list(leftpart);
709               rparts = list(rightpart);
710               //first deal with the left part
711               if (leftpart[thetapos] == x)
712               {
713                    shift_sign = 1;
714                    shiftvar = x;
715               }
716               else
717               {
718                    shift_sign = -1;
719                    shiftvar = y;
720               }
721               for (j = size(leftpart); j>1;j--)
722               {//drip x resp. y
723                    if (leftpart[j-1]==shiftvar)
724                    {//commutative
725                         j--;
726                         continue;
727                    }//commutative
728                    if (deg(leftpart[j-1],intvec(-1,1,0))!=0)
729                    {//stop here
730                         break;
731                    }//stop here
732                    //Here, we can only have a a0- part
733                    leftpart[j] = subst(leftpart[j-1],theta, theta + shift_sign);
734                    leftpart[j-1] = shiftvar;
735                    lparts = lparts + list(leftpart);
736                }//drip x resp. y
737                //and now deal with the right part
738                if (rightpart[1] == x)
739                {
740                     shift_sign = 1;
741                     shiftvar = x;
742                }
743                else
744                {
745                     shift_sign = -1;
746                     shiftvar = y;
747                }
748                for (j = 1 ; j < size(rightpart); j++)
749                {
750                     if (rightpart[j+1] == shiftvar)
751                     {
752                          j++;
753                          continue;
754                     }
755                     if (deg(rightpart[j+1],intvec(-1,1,0))!=0)
756                     {
757                          break;
758                     }
759                     rightpart[j] = subst(rightpart[j+1], theta, theta - shift_sign);
760                     rightpart[j+1] = shiftvar;
761                     rparts = rparts + list(rightpart);
762                }
763                //And now, we put all possibilities together
765                for (j = 1; j<=size(lparts); j++)
766                {
767                   for (k = 1; k<=size(rparts);k++)
768                   {
770                   }
771                }
773                result = result + tempadd;
774                continue; //We can may be not be done already with the ith entry
775          }//One entry was theta resp. theta +1
776     }//checking every entry of result for theta or theta +1
777     //map back to the basering
778     setring(r);
779     map finalmap = tempRing, var(1), var(2),var(1)*var(2);
780     list result = finalmap(result);
781     for (i=1; i<=size(result);i++)
783          result[i] = k_factor + result[i];
785     result = delete_dublicates_noteval(result);
786     return(result);
787}//proc HomogFacAll
788/* example */
789/* { */
790/*      "EXAMPLE:";echo=2; */
791/*      ring R = 0,(x,y),Ws(-1,1); */
792/*      def r = nc_algebra(1,1); */
793/*      setring(r); */
794/*      poly h = (x^2*y^2+1)*(x^4); */
795/*      homogfac_all(h); */
796/* } */
797
798//==================================================*
799//Computes all permutations of a given list
800static proc perm(list l)
801{//proc perm
802     int i; int j;
803     list tempresult;
804     list result;
805     if (size(l)==0)
806     {
807          return(list());
808     }
809     if (size(l)==1)
810     {
811          return(list(l));
812     }
813     for (i = 1; i<=size(l); i++ )
814     {
815          tempresult = perm(delete(l,i));
816          for (j = 1; j<=size(tempresult);j++)
817          {
818               tempresult[j] = list(l[i])+tempresult[j];
819          }
820          result = result+tempresult;
821     }
822     return(result);
823}//proc perm
824
825//==================================================
826//computes all permutations of a given list by
827//ignoring equal entries (faster than perm)
828static proc permpp(list l)
829{//proc permpp
830     int i; int j;
831     list tempresult;
832     list l_without_double;
833     list l_without_double_pos;
834     int double_entry;
835     list result;
836     if (size(l)==0)
837     {
838          return(list());
839     }
840     if (size(l)==1)
841     {
842          return(list(l));
843     }
844     for (i = 1; i<=size(l);i++)
845     {//Filling the list with unique entries
846          double_entry = 0;
847          for (j = 1; j<=size(l_without_double);j++)
848          {
849               if (l_without_double[j] == l[i])
850               {
851                    double_entry = 1;
852                    break;
853               }
854          }
855          if (!double_entry)
856          {
857               l_without_double = l_without_double + list(l[i]);
858               l_without_double_pos = l_without_double_pos + list(i);
859          }
860     }//Filling the list with unique entries
861     for (i = 1; i<=size(l_without_double); i++ )
862     {
863          tempresult = permpp(delete(l,l_without_double_pos[i]));
864          for (j = 1; j<=size(tempresult);j++)
865          {
866               tempresult[j] = list(l_without_double[i])+tempresult[j];
867          }
868          result = result+tempresult;
869     }
870     return(result);
871}//proc permpp
872
873//==================================================
874//factorization of the first Weyl Algebra
875
876//The following procedure just serves the purpose to
877//transform the input into an appropriate input for
878//the procedure sfacwa, where the ring must contain the
879//variables in a certain order.
880proc facfirstwa(poly h)
881{//proc facfirstwa
882      //Redefine the ring in my standard form
883     if (!isWeyl())
884     {//Our basering is not the Weyl algebra
885          return(list());
886     }//Our basering is not the Weyl algebra
887     if(nvars(basering)!=2)
888     {//Our basering is the Weyl algebra, but not the first
889          return(list());
890     }//Our basering is the Weyl algebra, but not the first
891     if (ringlist(basering)[6][1,2] == -1) //manual of ringlist will tell you why
892     {
893          def r = basering;
894          ring tempRing = ringlist(r)[1][1],(x,y),Ws(-1,1);
895          def tempRingnc = nc_algebra(1,1);
896          setring(tempRingnc);
897          map transf = r, var(2), var(1);
898          list result = sfacwa(transf(h));
899          setring(r);
900          map transfback = tempRingnc, var(2),var(1);
901          return(transfback(result));
902     }
903     else { return(sfacwa(h));}
904}//proc facfirstwa
905example
906{
907     "EXAMPLE:";echo=2;
908     ring R = 0,(x,y),dp;
909     def r = nc_algebra(1,1);
910     setring(r);
911     poly h = (x^2*y^2+x)*(y-1);
912     facfirstwa(h);
913}
914
915//This is the main program
916static proc sfacwa(poly h)
917"USAGE: facwa(h); h is a polynomial in the first Weyl algebra
918RETURN: list
919PURPOSE: Computes a factorization of a polynomial h in the first Weyl algebra
920THEORY: @code{homogfac} returns a list with some factorizations of the given polynomial. The possibilities of the factorization of the highest homogeneous part and those of the lowest will be merged.
921EXAMPLE: example facwa; shows examples
923"{//proc sfacwa
924     if(homogwithorder(h,intvec(-1,1)))
925     {
926          return(homogfac_all(h));
927     }
928     def r = basering;
929     map invo = basering,-var(1),var(2);
930     int i; int j; int k;
931     intvec limits = deg(h,intvec(1,1)) ,deg(h,intvec(1,0)),deg(h,intvec(0,1));
932     def prod;
933     //end finding the limits
934     poly maxh = jet(h,deg(h,intvec(-1,1)),intvec(-1,1))-jet(h,deg(h,intvec(-1,1))-1,intvec(-1,1));
935     poly minh = jet(h,deg(h,intvec(1,-1)),intvec(1,-1))-jet(h,deg(h,intvec(1,-1))-1,intvec(1,-1));
936     list result;
937     list temp;
938     list homogtemp;
939     list M; list hatM;
940     list f1 = homogfac_all(maxh);
941     list f2 = homogfac_all(minh);
942     int is_equal;
943     poly hath;
944     for (i = 1; i<=size(f1);i++)
945     {//Merge all combinations
946          for (j = 1; j<=size(f2); j++)
947          {
948               M = M + mergence(f1[i],f2[j],limits);
949          }
950     }//Merge all combinations
951     for (i = 1 ; i<= size(M); i++)
952     {//filter valid combinations
953          if (product(M[i]) == h)
954          {//We have one factorization
955               result = result + list(M[i]);
956               M = delete(M,i);
957               continue;
958          }//We have one factorization
959          else
960          {
961               if (deg(h,intvec(-1,1))<=deg(h-product(M[i]),intvec(-1,1)))
962               {
963                    M = delete(M,i);
964                    continue;
965               }
966               if (deg(h,intvec(1,-1))<=deg(h-product(M[i]),intvec(1,-1)))
967               {
968                    M = delete(M,i);
969                    continue;
970               }
971          }
972     }//filter valid combinations
973     M = delete_dublicates_eval(M);
974     while(size(M)>0)
975     {//work on the elements of M
976          hatM = list();
977          for(i = 1; i<=size(M); i++)
978          {//iterate over all elements of M
979               hath = h-product(M[i]);
980               temp = list();
981               //First check for common inhomogeneous factors between hath and h
982               if (involution(NF(involution(hath,invo), std(involution(ideal(M[i][1]),invo))),invo)==0)
983               {//hath and h have a common factor on the left
984                    j = 1;
985                    f1 = M[i];
986                    if (j+1<=size(f1))
987                    {//Checking for more than one common factor
988                         while(involution(NF(involution(hath,invo),std(involution(ideal(product(f1[1..(j+1)])),invo))),invo)==0)
989                         {
990                              if (j+1<size(f1))
991                              {
992                                   j++;
993                              }
994                              else
995                              {
996                                   break;
997                              }
998                         }
999                    }//Checking for more than one common factor
1000                    f2 = list(f1[1..j])+list(involution(lift(involution(product(f1[1..j]),invo),involution(hath,invo))[1,1],invo));
1001                    temp = temp + merge_cf(f2,f1,limits);
1002               }//hath and h have a common factor on the left
1003               if (reduce(hath, std(ideal(M[i][size(M[i])])))==0)
1004               {//hath and h have a common factor on the right
1005                    j = size(M[i]);
1006                    f1 = M[i];
1007                    if (j-1>0)
1008                    {//Checking for more than one factor
1009                         while(reduce(hath,std(ideal(product(f1[(j-1)..size(f1)]))))==0)
1010                         {
1011                              if (j-1>1)
1012                              {
1013                                   j--;
1014                              }
1015                              else
1016                              {
1017                                   break;
1018                              }
1019                         }
1020                    }//Checking for more than one factor
1021                    f2 = list(lift(product(f1[j..size(f1)]),hath)[1,1])+list(f1[j..size(f1)]);
1022                    temp = temp + merge_cf(f2,M[i],limits);
1023               }//hath and h have a common factor on the right
1024               //and now the homogeneous
1025               maxh = jet(hath,deg(hath,intvec(-1,1)),intvec(-1,1))-jet(hath,deg(hath,intvec(-1,1))-1,intvec(-1,1));
1026               minh = jet(hath,deg(hath,intvec(1,-1)),intvec(1,-1))-jet(hath,deg(hath,intvec(1,-1))-1,intvec(1,-1));
1027               f1 = homogfac_all(maxh);
1028               f2 = homogfac_all(minh);
1029               for (j = 1; j<=size(f1);j++)
1030               {
1031                    for (k=1; k<=size(f2);k++)
1032                    {
1033                         homogtemp = mergence(f1[j],f2[k],limits);
1034                    }
1035               }
1036               for (j = 1; j<= size(homogtemp); j++)
1037               {
1038                    temp = temp + mergence(homogtemp[j],M[i],limits);
1039               }
1040               for (j = 1; j<=size(temp); j++)
1041               {//filtering invalid entries in temp
1042                    if(product(temp[j])==h)
1043                    {//This is already a result
1044                         result = result + list(temp[j]);
1045                         temp = delete(temp,j);
1046                         continue;
1047                    }//This is already a result
1048                    if (deg(hath,intvec(-1,1))<=deg(hath-product(temp[j]),intvec(-1,1)))
1049                    {
1050                         temp = delete(temp,j);
1051                         continue;
1052                    }
1053               }//filtering invalid entries in temp
1054               hatM = hatM + temp;
1055          }//iterate over all elements of M
1056          M = hatM;
1057          for (i = 1; i<=size(M);i++)
1058          {//checking for complete factorizations
1059               if (h == product(M[i]))
1060               {
1061                    result = result + list(M[i]);
1062                    M = delete(M,i);
1063                    continue;
1064               }
1065          }//checking for complete factorizations
1066          M = delete_dublicates_eval(M);
1067     }//work on the elements of M
1068     //In the case, that there is none, write a constant factor before the factor of interest.
1069     for (i = 1 ; i<=size(result);i++)
1071          if (deg(result[i][1],intvec(1,1))!=0)
1072          {
1073               result[i] = insert(result[i],1);
1074          }
1076     result = delete_dublicates_noteval(result);
1077     return(result);
1078}//proc sfacwa
1079
1080
1081//==================================================
1082/*Singular has no way implemented to test polynomials
1083 for homogenity with respect to a weight vector.
1084 The following procedure does exactly this*/
1085static proc homogwithorder(poly h, intvec weights)
1086{//proc homogwithorder
1087     if(size(weights) != nvars(basering))
1088     {//The user does not know how many variables the current ring has
1089          return(0);
1090     }//The user does not know how many variables the current ring has
1091     int i;
1092     int dofp = deg(h,weights); //degree of polynomial
1093     for (i = 1; i<=size(h);i++)
1094     {
1095          if (deg(h[i],weights)!=dofp)
1096          {
1097               return(0);
1098          }
1099     }
1100     return(1);
1101}//proc homogwithorder
1102
1103//==================================================
1104//Testfac: Given a list with different factorizations of
1105// one polynomial, the following procedure checks
1106// whether they all refer to the same polynomial.
1107// If they do, the output will be a list, that contains
1108// the product of each factorization. If not, the empty
1109// list will be returned.
1110// If the optional argument # is given (i.e. the polynomial
1111// which is factorized by the elements of the given list),
1112// then we look, if the entries are factorizations of p
1113// and if not, a list with the products subtracted by p
1114// will be returned
1115proc testfac(list l, list #)
1116"USAGE:  - testfac(l); l is a list, that contains lists with entries in the first Weyl algebra or
1117         @*- testfac(l,p); l is a list, that contains lists with entries in the first Weyl algebra and p is a polynomial in the first Weyl algebra
1118RETURN: list@*
1119PURPOSE: Checks whether a list of factorizations contains factorizations of the same element in the first Weyl algebra@*
1120THEORY: @code{testfac} multiplies out each factorization and checks whether each factorization was a factorization of the same element.
1121@* - if there is only a list given, the output will be the empty list, if it does not contain factorizations of the same element. Otherwise it will multiply out all the factorizations and return a list of the length if the given one, in which each entry is the factorized polynomial
1122@* - if there is a polynomial in the second argument, then the procedure checks whether the given list contains factorizations of this polynomial. If it does, then the output will be a list with the length of the given one and each entry contains this polynomial. If it does not, then the procedure returns a list, in which each entry is the given polynomial subtracted by the multiplied candidate for its factorization.
1123EXAMPLE: example testfac; shows examples
1125"{//proc testfac
1126     if (size(l)==0)
1127     {//The empty list is given
1128          return(list());
1129     }//The empty list is given
1130     if (size(#)>1)
1131     {//We want max. one optional argument
1132          return(list());
1133     }//We want max. one optional argument
1134     list result;
1135     int i; int j;
1136     if (size(#)==0)
1137     {//No optional argument is given
1138          int valid = 1;
1139          for (i = size(l);i>=1;i--)
1140          {//iterate over the elements of the given list
1141               if (size(result)>0)
1142               {
1143                    if (product(l[i])!=result[size(l)-i])
1144                    {
1145                         valid = 0;
1146                         break;
1147                    }
1148               }
1149               result = insert(result, product(l[i]));
1150          }//iterate over the elements of the given list
1151          if (valid)
1152          {
1153               return(result);
1154          }
1155          return(list());
1156     }//No optional argument is given
1157     else
1158     {
1159          int valid = 1;
1160          for (i = size(l);i>=1;i--)
1161          {//iterate over the elements of the given list
1162               if (product(l[i])!=#[1])
1163               {
1164                    valid = 0;
1165               }
1166               result = insert(result, product(l[i]));
1167          }//iterate over the elements of the given list
1168          if (valid)
1169          {
1170               return(result);
1171          }
1172          for (i = 1 ; i<= size(result);i++)
1173          {//subtract p from every entry in result
1174               result[i] = result[i] - #[1];
1175          }//subtract p from every entry in result
1176          return(result);
1177     }
1178}//proc testfac
1179example
1180{
1181     "EXAMPLE:";echo=2;
1182     ring R = 0,(x,y),Ws(-1,1);
1183     def r = nc_algebra(1,1);
1184     setring(r);
1185     poly h = (x^2*y^2+1)*(x^4);
1186     def t1 = facfirstwa(h);
1187     //fist a correct list
1188     testfac(t1);
1189     //now a correct list with the factorized polynomial
1190     testfac(t1,h);
1191     //now we put in an incorrect list without a polynomial
1192     t1[5][3] = y;
1193     testfac(t1);
1194     //and last but not least we take h as additional input
1195     testfac(t1,h);
1196}
1197//==================================================
1198//Procedure facsubwa:
1199//This procedure serves the purpose to compute a
1200//factorization of a given polynomial in a ring, whose subring
1201//is the first Weyl algebra. The polynomial must only contain
1202//the two arguments, which are also given by the user.
1203proc facsubwa(poly h, X, D)
1204"USAGE:  facsubwa(h,x,y) - h is a polynomial in the first Weyl algebra, X and D are the variables in h, which have the commutation rule DX = XD + 1@*
1205RETURN: list@*
1206PURPOSE: Given a polynomial in a Ring, which contains the first Weyl algebra as a subring. Furthermore let this polynomial be in the first Weyl algebra with the variables X and D with DX = XD +1.@*
1207THEORY: @code{facsubwa} computes some factorizations of the given polynomial, which is a polynomial in the variables X and D@*
1208EXAMPLE: example facsubwa; shows examples
1210"{//proc facsubwa
1211     if(!isWeyl())
1212     {
1213          return(list());
1214     }
1215     def r = basering;
1216     //We begin to check the input
1217     if (size(X)>1)
1218     {
1219          return(list());
1220     }
1221     if (size(D)>1)
1222     {
1223          return(list());
1224     }
1227     if (sum(lexpofX)!=1)
1228     {
1229          return(list());
1230     }
1231     if (sum(lexpofD)!=1)
1232     {
1233          return(list());
1234     }
1235     int varnumX=1;
1236     int varnumD=1;
1237     while(lexpofX[varnumX] != 1)
1238     {
1239          varnumX++;
1240     }
1241     while(lexpofD[varnumD] != 1)
1242     {
1243          varnumD++;
1244     }
1245     lexpofX = lexpofX +1;
1246     lexpofX[varnumX] = 0;
1247     lexpofX[varnumD] = 0;
1248     if (deg(h,lexpofX)!=0)
1249     {
1250          return(list());
1251     }
1252     //Input successfully checked
1253     ring firstweyl = 0,(var(varnumX),var(varnumD)),dp;
1254     def Firstweyl = nc_algebra(1,1);
1255     setring Firstweyl;
1256     poly h = imap(r,h);
1257     list result = facfirstwa(h);
1258     setring r;
1259     list result = imap(Firstweyl,result);
1260     return(result);
1261}//proc facsubwa
1262example
1263{
1264     "EXAMPLE:";echo=2;
1265     ring R = 0,(x,y),Ws(-1,1);
1266     def r = nc_algebra(1,1);
1267     setring(r);
1268     poly h = (x^2*y^2+x)*x;
1269     facsubwa(h,x,y);
1270}
1271//==================================================
1272/*
1273Example polynomials where one can find factorizations
1274(x^2+y)*(x^2+y);
1275(x^2+x)*(x^2+y);
1276(x^3+x+1)*(x^4+y*x+2);
1277(x^2*y+y)*(y+x*y);
1278y^3+x*y^3+2*y^2+2*(x+1)*y^2+y+(x+2)*y; //Example 5 Gregoriev Schwarz.
1279(y+1)*(y+1)*(y+x*y); //Landau Example projected to the first dimension.
1280 */
Note: See TracBrowser for help on using the repository browser.