source: git/Singular/LIB/groups.lib @ cb8103a

spielwiese
Last change on this file since cb8103a was 1e1ec4, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
Updated LIBs according to master add: new LIBs from master fix: updated LIBs due to minpoly/(de)numerator changes fix: -> $Id$ fix: Fixing wrong rebase of SW on master (LIBs)
  • Property mode set to 100644
File size: 28.4 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id$";
3category="Applications";
4info="
5LIBRARY:  groups.lib      Finite Group Theory
6AUTHORS:  Gert-Martin Greuel, email: greuel@mathematik.uni-kl.de @*
7          Gerhard Pfister, email: pfister@mathematik.uni-kl.de
8
9PROCEDURES:
10
11   noSolution(I)         I an ideal in a polynomial ring over Z[x1..xn]
12                         returns a list l of primes <=32003 such that 1
13                         is in IZ/p[x1..xn] for p not in l or an ERROR
14                         if this is wrong or not decided.
15"
16/*
17// ===================== A Problem in Finite Group Theory ===================
18// Posed by Boris Kunyavskii, Bar-Ilan University, Tel Aviv
19//
20// For any word w in X,Y,X^(-1),Y^(-1) consider the sequence U_n
21// of words (depending on w) inductively
22//           U_1 = w
23//         U_n+1 = [X*U_n*X^(-1),Y*U_n*Y^(-1)]
24// with
25//         [X,Y] = X*Y*X^(-1)*Y^(-1)  (the commutator)
26
27
28// Conjecture (1):
29// (by B. Plotkin, slightly modified):
30// A finite group G is solvable <==>
31// there is an n >= 1 such that U_n(x,y) = 1 for all x,y in G
32// and for any of the 4 following words
33//         w1 = X^(-1)*Y*X*Y^(-1)*X
34//         w2 = X^(-2)*Y^(-1)*X
35//         w3 = Y^(-2)*X^(-1)*X
36//         w4 = X*Y^(-2)*X^(-1)*Y*X^(-1)
37//
38// (These words remained as possibly good words by a computer search
39// through about 10 000 candidates, by Fritz Grunewald)
40
41
42//  ==>  is clear
43// The minimal finite non-solvable groups (i.e. every proper subgroup is
44// solvable) have been classified by Thompson in 1968, they are:
45//
46// 1. PSL(2,p)    (p=5 or p=+-2 (mod 5), p!=3)
47// 2. PSL(2,2^p)
48// 3. PSL(2,3^p)  (p odd)
49// 4. Sz(2^p)     (p odd)      (Suzuki group)
50// 5. PSL(3,3)
51//
52// In view of this result Conjecture (1) is equivalent to
53//
54// Conjecture (2):
55// Let G be one of the groups above, then, for at least one of the 4 words w
56// above there are x,y in G such that 1 != U_n(x,y) = U_n+1(x,y) for some n.
57// (then U_n(x,y) != 1 for all n by definitionof U_n)
58
59
60// We give a computer aided proof, using SINGULAR, of Conjecture (2) for
61// the groups 1. - 3. (up to checking for a small, explicit number of primes)
62// Hence only the Suzuki groups have to be checked.
63//
64// We show:  1 != U_1(x,y) = U_2(x,y) for word w1 and some x,y in G
65// We need:
66// Theory
67//   - simple facts from algebraic geometry, singularity theory, finite fields
68//   - theorem of Lang-Weil, estimating the number of rational points on an
69//     absolutely irreducible projective curve C defined over Z (g=genus):
70//     if q >= N(g) ==> #C(F_q) != 0
71//     (N(g) = 4g^2 -2, g arithmetic genus)
72// SINGULAR
73//   - Groebner basis (elimination), multivariate factorization, resolution of
74//     plane curve singularities (Hamburger Noether developement), primary
75//     decomposition
76
77*/
78
79LIB "standard.lib";
80LIB "general.lib";
81LIB "matrix.lib";
82///////////////////////////////////////////////////////////////////////////////
83
84static proc splitS1(ideal I,int s,list #)
85"USAGE: splitS1(I,s[,l]) I ideal, s integer, l list of ideals
86COMPUTE: Factorizes the generators of I. If for one generator f=f1*f2 and
87         fi=ni*gi^ri, ni an integer and gi a polynomial, then
88         compute a standardbasis of I1=I,g1 and I2=I,g2, if this can be
89         done in < s seconds. Then apply splitS to I1 and I2 and continue
90         in the same way. The procedure stops if no generator can be
91         factorized.
92RETURN:  A list L of ideals and prime numbers
93         L[1]: A list of ideals such that the radical of the intersection
94               of these ideals coincides with the radical of I
95               If the optional list l of the input is not empty then
96               the ideals of L[1] which contain an ideal of l are canceled.
97         L[2]: A list of prime numbers appearing as factors of the ni.
98NOTE:   The computation avoids division by integers (by using
99        option(contentSB)) hence the result is correct modulo any prime
100        number which does not appear in the list L[2].
101EXAMPLE: example splitS1; shows an example
102"
103{
104   option(redSB);
105   option(contentSB);
106   int j,k,e;
107   int i=1;
108   int l=attrib(I,"isSB");
109   ideal J;
110   list re,fac,te,pr,qr,w;
111   number n;
112   poly p;
113   re=#;
114
115   if(deg(I[1])==0){return(list(re+list(std(I)),qr));}
116
117   fac=factorize(I[1]);
118
119   while((size(fac[1])==2)&&(i<size(I)))
120   {
121      I[i]=fac[1][2]*fac[1][1];   //not in splitS
122      i++;
123      fac=factorize(I[i]);
124   }
125   if(size(fac[1])==2){I[size(I)]=fac[1][2]*fac[1][1];}  //not in splitS
126   if(size(fac[1])>2)
127   {
128      w=squarefreeP(number(fac[1][1]));
129      n=w[1];
130      qr=insResult(qr,w[2]);
131      for(j=2;j<=size(fac[1]);j++)
132      {
133         I[i]=fac[1][j];
134         attrib(I,"isSB",1);
135         e=1;
136         k=0;
137         while(k<size(re))
138         {
139            k++;
140            if(size(reduce(re[k],I))==0){e=0;break;}
141            attrib(re[k],"isSB",1);
142            if(size(reduce(I,re[k]))==0){re=delete(re,k);k--;}
143         }
144         if(j==2){I[i]=I[i]*n;}
145         if(e)
146         {
147            if(l)
148            {
149               J=I;
150               p=I[i];
151               J[i]=0;
152               J=simplify(J,2);
153               attrib(J,"isSB",1);
154               pr=splitS(std(J,p),s,re);
155               re=pr[1];
156               qr=insResult(qr,pr[2]);
157            }
158            else
159            {
160               J=interred(I);
161               pr=splitS(timeStd(J,s),s,re);
162               re=pr[1];
163               qr=insResult(qr,pr[2]);
164            }
165         }
166      }
167      return(list(re,qr));
168   }
169   J=timeStd(I,s);       //J=std(I) in splitS
170   attrib(I,"isSB",1);
171   if(size(reduce(J,I))==0){return(list(re+list(I),qr));}
172   pr=splitS(J,s,re);
173   return(list(re+pr[1],pr[2]));
174}
175example
176{ "EXAMPLE:"; echo = 2;
177   ring r=0,(b,s,t,u,v,w,x,y,z),dp;
178   ideal i=
179   bv+su,
180   bw+tu,
181   sw+tv,
182   by+sx,
183   bz+tx,
184   sz+ty,
185   uy+vx,
186   uz+wx,
187   vz+wy,
188   bvz;
189   splitS1(i,5);
190}
191
192///////////////////////////////////////////////////////////////////////////////
193static proc splitS(ideal I,int s,list #)
194"USAGE: splitS(I,s[,l]) I ideal, s integer, l list of ideals
195COMPUTE: Factorizes the generators of I. If for one generator f=f1*f2 and
196         fi=ni*gi^ri, ni an integer and gi a polynomial, then
197         compute a standardbasis of I1=I,g1 and I2=I,g2, if this can be
198         done in < s seconds. Then apply splitS to I1 and I2 and continue
199         in the same way. The procedure stops if no generator can be
200         factorized.
201RETURN:  A list L of ideals and prime numbers
202         L[1]: A list of ideals such that the radical of the intersection
203               of these ideals coincides with the radical of I
204               If the optional list l of the input is not empty then
205               the ideals of L[1] which contain an ideal of l are canceled.
206         L[2]: A list of prime numbers appearing as factors of the ni.
207NOTE:   The computation avoids division by integers (by using
208        option(contentSB)) hence the result is correct modulo any prime
209        number which does not appear in the list L[2].
210EXAMPLE: example splitS; shows an example
211"
212{
213   option(redSB);
214   option(contentSB);
215   int j,k,e;
216   int i=1;
217   int l=attrib(I,"isSB");
218   ideal J;
219   list re,fac,te,pr,qr,w;
220   number n;
221   poly p;
222   re=#;
223
224   if(deg(I[1])==0){return(list(re+list(std(I)),qr));}
225
226   fac=factorize(I[1]);
227
228   while((size(fac[1])==2)&&(i<size(I)))
229   {
230      i++;
231      fac=factorize(I[i]);
232   }
233   if(size(fac[1])>2)
234   {
235      w=squarefreeP(number(fac[1][1]));
236      n=w[1];
237      qr=insResult(qr,w[2]);
238      for(j=2;j<=size(fac[1]);j++)
239      {
240         I[i]=fac[1][j];
241         attrib(I,"isSB",1);
242         e=1;
243         k=0;
244         while(k<size(re))
245         {
246            k++;
247            if(size(reduce(re[k],I))==0){e=0;break;}
248            attrib(re[k],"isSB",1);
249            if(size(reduce(I,re[k]))==0){re=delete(re,k);k--;}
250         }
251         if(j==2){I[i]=I[i]*n;}
252         if(e)
253         {
254            if(l)
255            {
256               J=I;
257               p=I[i];
258               J[i]=0;
259               J=simplify(J,2);
260               attrib(J,"isSB",1);
261               pr=splitS(std(J,p),s,re);
262               re=pr[1];
263               qr=insResult(qr,pr[2]);
264            }
265            else
266            {
267               J=interred(I);
268               pr=splitS(timeStd(J,s),s,re);
269               re=pr[1];
270               qr=insResult(qr,pr[2]);
271            }
272         }
273      }
274      return(list(re,qr));
275   }
276   J=std(I);
277   attrib(I,"isSB",1);
278   if(size(reduce(J,I))==0){return(list(re+list(I),qr));}
279   pr=splitS(J,s,re);
280   return(list(re+pr[1],pr[2]));
281}
282example
283{ "EXAMPLE:"; echo = 2;
284   ring r=0,(b,s,t,u,v,w,x,y,z),dp;
285   ideal i=
286   bv+su,
287   bw+tu,
288   sw+tv,
289   by+sx,
290   bz+tx,
291   sz+ty,
292   uy+vx,
293   uz+wx,
294   vz+wy,
295   bvz;
296   splitS(i,5);
297}
298
299///////////////////////////////////////////////////////////////////////////////
300static proc finalSplit(list I,list qr)
301"USAGE: finalSplit(I,qr) I list of ideals, qr list of primes
302RETURN: list l of primes <=32003 such that 1 is in I[j]Z/p[x1..xn]
303        for p not in l and all j or an ERROR if this is wrong or
304        not decided.
305EXAMPLE: example finalSplit; shows an example
306"
307{
308   option(redThrough);
309   option(contentSB);
310   ideal J,K;
311   int i,j,k,n;
312   int count=1;
313   list q,pr,l;
314
315   "trivial splitting";
316   pr=trivialSplit(I,4);
317   l=pr[1];
318   qr=insResult(qr,pr[2]);
319   while(count)
320   {
321      k++;
322      "loop";k;size(l);"======";
323      count=0;list p;
324      for(i=1;i<=size(l);i++)
325      {
326         i;
327         K=changeOrdTest(l[i]);
328         if(deg(K[1])!=0)
329         {
330            "split";
331            pr=splitS1(shortid_L(K,3),2);
332            q=pr[1];
333            qr=insResult(qr,pr[2]);
334            size(q);"out";
335            for(j=1;j<=size(q);j++)
336            {
337               if(deg(q[j][1])==0)
338               {
339                  pr=contentS(q[j]);
340                  q[j]=pr[1];
341                  qr=insResult(qr,pr[2]);
342               }
343               else
344               {
345                  pr=simpliFy(q[j]+K);
346                  q[j]=pr[1];
347                  qr=insResult(qr,pr[2]);
348                  pr=contentS(q[j]);
349                  q[j]=pr[1];
350                  qr=insResult(qr,pr[2]);
351               }
352            }
353            if((size(q)==1)&&(deg(q[1][1])>0))
354            {
355               "split again";
356               K=sort(q[1])[1];
357               J=K[1..20+4*k];
358               pr=splitS(J,5);
359               q=pr[1];
360               qr=insResult(qr,pr[2]);
361               size(q);"out";
362               for(j=1;j<=size(q);j++)
363               {
364                  if(deg(q[j][1])==0)
365                  {
366                     pr=contentS(q[j]);
367                     q[j]=pr[1];
368                     qr=insResult(qr,pr[2]);
369                  }
370                  else
371                  {
372                     pr=simpliFy(q[j]+K);
373                     q[j]=pr[1];
374                     qr=insResult(qr,pr[2]);
375                     pr=contentS(q[j]);
376                     q[j]=pr[1];
377                     qr=insResult(qr,pr[2]);
378                  }
379               }
380
381            }
382            if((size(q)==1)&&(deg(q[1][1])>0))
383            {
384               n++;
385               "split again";
386               K=q[1];
387               J=shortid_L(K,3+n);
388               pr=splitS1(J,5);
389               if(size(pr[1])>1)
390               {
391                  q=pr[1];
392                  qr=insResult(qr,pr[2]);
393                  size(q);"out";
394                  for(j=1;j<=size(q);j++)
395                  {
396                     if(deg(q[j][1])==0)
397                     {
398                       pr=contentS(q[j]);
399                       q[j]=pr[1];
400                       qr=insResult(qr,pr[2]);
401                     }
402                     else
403                     {
404                       pr=simpliFy(q[j]+K);
405                       q[j]=pr[1];
406                       qr=insResult(qr,pr[2]);
407                       pr=contentS(q[j]);
408                       q[j]=pr[1];
409                       qr=insResult(qr,pr[2]);
410                     }
411                  }
412               }
413            }
414            if(size(q)>1)
415            {
416                count++;
417            }
418            else
419            {
420               if(deg(q[1][1])>0)
421               {
422                  q[1]=changeOrdTest(q[1]);
423               }
424            }
425            p=p+q;
426         }
427         else
428         {
429           pr=primefactors(number(K[1]));
430           qr=insResult(qr,pr[1]);
431         }
432      }
433      l=p;
434      kill p;
435   }
436l;
437   for(i=1;i<=size(l);i++)
438   {
439
440      if(deg(l[i][1])>0){l;ERROR("not ready");}
441      pr=primefactors(number(l[i][1]));
442      if(pr[3]!=1){pr;ERROR("not ready");}
443      qr=insResult(qr,pr[1]);
444   }
445   return(qr);
446}
447example
448{ "EXAMPLE:"; echo = 2;
449   ring r=0,x,dp;
450   list qr=2,5,7;
451   ideal i=181x-181,11x+11;
452   list pr=i;
453   finalSplit(pr,qr);
454}
455
456///////////////////////////////////////////////////////////////////////////////
457proc noSolution(ideal I)
458"USAGE: noSolution(I) I ideal
459RETURN: list l of primes <=32003 such that 1 is in IZ/p[x1..xn]
460        for p not in l or an ERROR if this is wrong or not decided.
461EXAMPLE: example noSolution; shows an example
462"
463{
464   int s=1;
465   int t=30;
466   option(redThrough);
467   option(contentSB);
468   ideal J=shortid_L(I,3);
469   ideal K;
470   number n;
471
472   "first splitting";
473   int i,j,k;
474   list l,p,q,re,pr,qr;
475   pr=splitS(J,s);
476   l=pr[1];
477   qr=pr[2];
478   size(l);
479
480   for(i=1;i<=size(l);i++)
481   {
482      if(deg(l[i][1])==0)
483      {
484         K=l[i][1];
485      }
486      else
487      {
488         pr=simpliFy(I+l[i]);
489         J=pr[1];
490         qr=insResult(qr,pr[2]);
491         pr=contentS(J);
492         J=pr[1];
493         qr=insResult(qr,pr[2]);
494         K=timeStd(J,1);
495      }
496      if(deg(K[1])==0)
497      {
498         n=number(K[1]);
499         pr=primefactors(n);
500         if(pr[3]!=1)
501         {
502             J=changeOrdTest(J);
503             p=p+list(J);
504         }
505         else
506         {
507             qr=insResult(qr,pr[1]);
508         }
509      }
510      else
511      {
512         pr=contentS(K);
513         p=p+list(pr[1]);
514         qr=insResult(qr,pr[2]);
515      }
516   }
517   "trivial splitting 1";
518   pr=trivialSplit(p,3);
519   p=pr[1];
520   qr=insResult(qr,pr[2]);
521   size(p);
522   "second splitting";
523
524   for(i=1;i<=size(p);i++)
525   {
526      i;
527      if(size(p[i])<t)
528      {
529         K=timeStd(p[i],10);
530         if(deg(K[1])==0)
531         {
532            n=number(K[1]);
533            pr=primefactors(n);
534            if(pr[3]!=1)
535            {
536                p[i]=changeOrdTest(p[i]);
537                re=re+list(p[i]);
538            }
539            else
540            {
541               qr=insResult(qr,pr[1]);
542            }
543         }
544         else
545         {
546            pr=contentS(K);
547            re=re+list(pr[1]);
548            qr=insResult(qr,pr[2]);
549         }
550      }
551      else
552      {
553         J=p[i];
554         J=J[1..t];
555         "in splitting";
556         pr=splitS(J,s);
557         l=pr[1];
558         qr=insResult(qr,pr[2]);
559         size(l);
560         "out";
561         for(k=1;k<=size(l);k++)
562         {
563            pr=contentS(l[k]);
564            l[k]= pr[1];
565            qr=insResult(qr,pr[2]);
566            pr=simpliFy(l[k]+p[i]);
567            J=pr[1];
568            qr=insResult(qr,pr[2]);
569            K=timeStd(J,10);
570            if(deg(K[1])==0)
571            {
572               n=number(K[1]);
573               pr=primefactors(n);
574               if(pr[3]!=1)
575               {
576                  pr=contentS(J);
577                  qr=insResult(qr,pr[2]);
578                  J=pr[1];
579                  J=changeOrdTest(J);
580                  re=re+list(J);
581               }
582               else
583               {
584                  qr=insResult(qr,pr[1]);
585               }
586            }
587            else
588            {
589               pr=contentS(K);
590               qr=insResult(qr,pr[2]);
591               re=re+list(pr[1]);
592            }
593         }
594      }
595   }
596   "trivial splitting 2";
597   size(re);
598   pr=trivialSplit(re,2);
599   re=pr[1];
600   qr=insResult(qr,pr[2]);
601   "third splitting";
602   size(re);
603   for(i=1;i<=size(re);i++)
604   {
605      if(deg(re[i][1])>0)
606      {
607         i;
608         pr=simpliFy(re[i]);
609         J=pr[1];
610         qr=insResult(qr,pr[2]);
611         "in splitting";
612         J=shortid_L(J,3);
613         pr=splitS(J,s);
614         l=pr[1];
615         if(size(l)==1)
616         {
617            "split again";
618            pr=simpliFy(re[i]);
619            J=pr[1];
620            qr=insResult(qr,pr[2]);
621            J=J[1..t];
622            pr=splitS(J,s);
623            l=pr[1];
624            qr=insResult(qr,pr[2]);
625         }
626         else
627         {
628            qr=insResult(qr,pr[2]);
629         }
630         size(l);
631         "out";
632         for(j=1;j<=size(l);j++)
633         {
634            pr=contentS(l[j]);
635            l[j]=pr[1];
636            qr=insResult(qr,pr[2]);
637            pr=simpliFy(re[i]+l[j]);
638            J=pr[1];
639            qr=insResult(qr,pr[2]);
640            K=timeStd(J,10);
641            if(deg(K[1])==0)
642            {
643               n=number(K[1]);
644               pr=primefactors(n);
645               if(pr[3]!=1)
646               {
647                  pr=contentS(J);
648                  J=pr[1];
649                  qr=insResult(qr,pr[2]);
650                  J=changeOrdTest(J);
651                  q=q+list(J);
652               }
653               else
654               {
655                 qr=insResult(qr,pr[1]);
656               }
657            }
658            else
659            {
660               pr=contentS(K);
661               qr=insResult(qr,pr[2]);
662               q=q+list(pr[1]);
663            }
664         }
665      }
666      else
667      {
668         i;
669         pr=primefactors(number(re[i][1]));
670         qr=insResult(qr,pr[1]);
671      }
672   }
673"jetzt geht es los";
674   return(finalSplit(q,qr));
675}
676example
677{ "EXAMPLE:"; echo = 2;
678   ring r=0,x,dp;
679   ideal i=181x-181,11x+11;
680   noSolution(i);
681}
682
683///////////////////////////////////////////////////////////////////////////////
684
685static proc changeOrdTest(ideal I)
686{
687   def R=basering;
688   if(deg(I[1])==0){return(I);}
689   def RR=changeord(list(list("dp",1:nvars(basering))));
690   setring RR;
691   ideal I=imap(R,I);
692   ideal K=timeStd(I,5);
693   if(deg(K[1])==0)
694   {
695      number n=number(K[1]);
696      if(primefactors(n)[3]==1)
697      {
698         I=K;
699      }
700   }
701   setring R;
702   I=imap(RR,I);
703   kill RR;
704   return(I);
705}
706///////////////////////////////////////////////////////////////////////////////
707
708static proc trivialSplit(list p,int depth, list #)
709"USAGE: trivialSplit(p,s) p list of ideals, s integer the number of iterations
710COMPUTE: Factorizes the monomials among the generators of I.
711         If  one monomial contains the variables x1,..xr, then
712         I1=I(x1=0),...,Ir=I(xr=0) is considered. Then apply
713         trivialSplit to I1 ...Ir and continue in the same way s times.
714RETURN:  A list L of ideals and prime numbers
715         L[1]: A list of ideals such that the radical of the intersection
716               of these ideals at x1=...=xr=0 coincides with the radical of
717               I at x1=...=xr=0.
718         L[2]: A list of prime numbers appearing as factors of the monomials.
719NOTE:   The computation avoids division by integers  hence the result
720        is correct modulo any prime
721        number which does not appear in the list L[2].
722EXAMPLE: example trivialSplit; shows an example
723"
724{
725   list re,l,pr,qr;
726   int i,k;
727   ideal J,K,T,Ke;
728   number n;
729
730   if(size(p)==0){return(list(p,qr));}
731   if(depth<=0){return(list(p,qr));}
732   if(size(#)>0)
733   {
734      T=#[1];
735   }
736   else
737   {
738      T=1;
739   }
740   for(k=1;k<=size(p);k++)
741   {
742      pr=simpliFy(p[k]);
743      p[k]=pr[1];
744      qr=insResult(qr,pr[2]);
745      J=shortid_L(p[k],1);
746      if((size(J)>0)&&(deg(J[1])>=1))
747      {
748         pr=splitS1(J,10);
749         l=pr[1];
750         qr=insResult(qr,pr[2]);
751         for(i=1;i<=size(l);i++)
752         {
753            Ke=l[i];
754            l[i]=trivialSimplify(p[k],l[i]);
755            pr=simpliFy(l[i]);
756            l[i]=pr[1];
757            qr=insResult(qr,pr[2]);
758            K=timeStd(l[i],1);
759            attrib(K,"isSB",1);
760            if(deg(K[1])==0)
761            {
762               n=number(K[1]);
763               if(primefactors(n)[3]!=1)
764               {
765                  l[i]=changeOrdTest(l[i]);
766                  pr=trivialSplit(l[i],depth-1,T);
767                  re=re+pr[1];
768                  qr=insResult(qr,pr[2]);
769               }
770               else
771               {
772                  l[i]=K; //neu
773               }
774            }
775            else
776            {
777               if(size(reduce(trivialSimplify(T,Ke),K))!=0)
778               {
779                  pr=trivialSplit(K,depth-1,trivialSimplify(T,Ke));
780                  re=re+pr[1];
781                  qr=insResult(qr,pr[2]);
782               }
783            }
784            T=intersect(T,Ke);
785         }
786      }
787      else
788      {
789         J=timeStd(p[k],5);
790         if(deg(J[1])==0)
791         {
792            n=number(J[1]);
793            if(primefactors(n)[3]!=1)
794            {
795               p[k]=changeOrdTest(p[k]);
796               re=re+list(p[k]);
797            }
798            else
799            {
800                re=re+list(J);
801            }
802         }
803         else
804         {
805            re=re+list(J);
806         }
807      }
808   }
809   return(list(re,qr));
810}
811example
812{ "EXAMPLE:"; echo = 2;
813   ring r=0,(b,s,t,u,v,w,x,y,z),dp;
814   ideal i=
815   bv+su,
816   bw+tu,
817   sw+tv,
818   by+sx,
819   bz+tx,
820   sz+ty,
821   uy+vx,
822   uz+wx,
823   vz+wy,
824   bvz;
825   trivialSplit(i,2);
826}
827
828///////////////////////////////////////////////////////////////////////////////
829
830static proc trivialSimplify(ideal I, ideal J)
831"USAGE:  trivialSimplify(I,J);  I,J ideals, J generated by variables
832RETURN:  ideal K = eliminate(I,m) m the product of variables of J
833EXAMPLE: example trivialSimplify; shows an example
834"
835{
836   int i;
837   for(i=1;i<=size(J);i++){I=subst(I,J[i],0);}
838   return(simplify(I,2));
839}
840example
841{ "EXAMPLE:"; echo = 2;
842   ring r=0,(x,y,z,w,t,u),dp;
843   ideal i=
844   t,u,
845   x3+y2+2z,
846   x2+3y,
847   x2+y2+z2,
848   w2+z+u;
849   ideal j=t,u;
850   trivialSimplify(i,j);
851}
852
853///////////////////////////////////////////////////////////////////////////////
854
855static proc simpliFy(ideal J,list #)
856"USAGE:  simpliFy(id);  id ideal
857RETURN:  ideal I = eliminate(Id,m) m is a product of variables
858         which are in linear equations involved in the generators of id
859EXAMPLE: example simpliFy; shows an example
860"
861{
862   ideal I=J;
863   if(size(#)!=0){I=#[1];}
864   def R=basering;
865   matrix M=jacob(I);
866   ideal ma=maxideal(1);
867   int i,j,k;
868   map phi;
869   list pr,re;
870
871   for(i=1;i<=nrows(M);i++)
872   {
873      for(j=1;j<=ncols(M);j++)
874      {
875         if((deg(M[i,j])==0)&&(deg(I[i])==1))
876         {
877            ma[j]=(-1/M[i,j])*(I[i]-M[i,j]*var(j));
878            pr=primefactors(number(M[i,j]));
879            if(pr[3]>1){pr[3];ERROR("number too big in simpliFy");}
880            re=insResult(re,pr[1]);
881            phi=R,ma;
882            I=phi(I);
883            J=phi(J);
884            for(k=1;k<=ncols(I);k++)
885            {
886               I[k]=cleardenom(I[k]);
887            }
888            M=jacob(I);
889         }
890      }
891   }
892   J=simplify(J,2);
893   for(i=1;i<=size(J);i++)
894   {
895      if(deg(J[i])==0){J=std(J);break;}
896      J[i]=cleardenom(J[i]);
897   }
898   return(list(J,re));
899}
900example
901{ "EXAMPLE:"; echo = 2;
902   ring r=0,(x,y,z,w),dp;
903   ideal i=
904   x3+y2+2z,
905   x2+3y,
906   x2+y2+z2,
907   w2+z;
908   simpliFy(i);
909}
910///////////////////////////////////////////////////////////////////////////////
911
912static proc squarefreeP(number n)
913"USAGE:  squarefreeP(n);  n number
914RETURN:  list l of numbers. l[2] contains all prime factors of n less
915         then 32003, l[1] is the part of m which has prime factors greater
916         then 32003
917EXAMPLE: example squarefreeP; shows an example
918"
919{
920   list re;
921   if(n<0){n=-n;}
922   if(n==1){return(list(n,re));}
923   list pr=primefactors(n);
924   int i;
925   number m=number(pr[3][1]);
926   re=insResult(re,pr[1]);
927   return(list(m,re));
928}
929example
930{ "EXAMPLE:"; echo = 2;
931   ring r = 0,x,dp;
932   squarefreeP(123456789100);
933}
934
935///////////////////////////////////////////////////////////////////////////////
936
937static proc contentS(ideal I)
938"USAGE:  contentS(I);  I ideal
939RETURN:  list l
940         l[1] the ideal I. Te generators of I are the generators of the
941         input ideal devided by the part of the content which have
942         prime factors less then 32003.
943         l[2] contains the prime numbers which occured in the division
944EXAMPLE: example contentS; shows an example
945"
946
947{
948   option(contentSB);
949   int i,k;
950   number n,m;
951   list pr,re;
952   for(i=1;i<=size(I);i++)
953   {
954      n=leadcoef(I[i]);
955      if(n<0){n=-n;}
956      if(n>1)
957      {
958         pr=primefactors(n);
959         if(n<=32003)
960         {
961            m=n;
962         }
963         else
964         {
965            m=number(pr[1][1])^pr[2][1];
966            for(k=2;k<=size(pr[1]);k++)
967            {
968               m=m*number(pr[1][k])^pr[2][k];
969            }
970         }
971         I[i]=cleardenom(I[i]/m);
972         re=insResult(re,pr[1]);
973      }
974   }
975   return(list(I,re));
976}
977example
978{ "EXAMPLE:"; echo = 2;
979   ring r = 0,(x,y),dp;
980   ideal I=2x+2,3y+3x,1234567891x+1234567891;
981   contentS(I);
982}
983
984///////////////////////////////////////////////////////////////////////////////
985
986static proc insResult(list r,#)
987{
988   if(size(#)==0){return(r);}
989   if(size(r)==0){r[1]=1;}
990   int i,j;
991   for(i=1;i<=size(#);i++)
992   {
993      j=1;
994      while((#[i]>r[j])&&(j<size(r))){j++;}
995      if(#[i]>r[j]){r=insert(r,#[i],j);}
996      if(#[i]<r[j]){r=insert(r,#[i],j-1);}
997   }
998   return(r);
999}
1000example
1001{ "EXAMPLE:"; echo = 2;
1002   list r=list(3,7,13);
1003   intvec v=2,3,5,11,17;
1004   insResult(r,v);
1005}
1006///////////////////////////////////////////////////////////////////////////////
1007
1008static proc shortid_L (ideal id,int n)
1009{
1010  ideal Lc;
1011  intvec v;
1012  int ii;
1013  for(ii= 1; ii<=ncols(id); ii++)
1014  {
1015   if (size(id[ii])<=n and id[ii]!=0)
1016   {
1017       Lc = Lc,id[ii];
1018       v=v,ii;
1019   }
1020  }
1021  Lc = simplify(Lc,2);
1022  list L = Lc,v;
1023  return(Lc);
1024}
1025
1026
1027/*
1028//-------------- the solution of the problem -----------------
1029ring r = 0,(b,c,t),dp;   //global (affine) ring
1030matrix X[2][2] = 0, -1,
1031                 1, t;
1032matrix Y[2][2] = 1, b,
1033                 c, 1+bc;
1034;
1035
1036//------------ create word W1 and ideal I of U1-U2 ------------
1037matrix iX = inverse(X);
1038matrix iY = inverse(Y);
1039matrix U1 = iX*Y*X*iY*X;          //the word w1
1040//matrix U1=iX*iX*iY*X;           //das neue Wort
1041matrix N  = X*U1*iX;
1042matrix M  = Y*U1*iY;
1043matrix iN = inverse(N);
1044matrix iM = inverse(M);
1045matrix U2 = N*M*iN*iM;
1046ideal I = ideal(U1-U2);
1047
1048//list qr=primdecGTZ(I);
1049//I=qr[1][2];
1050
1051ring rh = 0,(b,c,t,h),dp;
1052ideal I = imap(r,I);
1053ideal sI = groebner(I);
1054ideal hI = homog(sI,h);
1055ideal shI =std(hI);
1056ideal J = eliminate(shI,c);       //eliminate c
1057poly f = J[1];
1058
1059ring r1 = 0,(b,t,h),dp;
1060poly hf=b6t6-2b5t7+b4t8+b8t3h-4b7t4h+7b6t5h-6b5t6h+2b4t7h-7b6t4h2+12b5t5h2+b4t6h2-6b3t7h2-3b8th3+12b7t2h3-16b6t3h3+19b4t5h3-12b3t6h3-2b8h4+8b7th4-3b6t2h4+2b5t3h4-45b4t4h4+32b3t5h4+12b2t6h4-12b6th5+50b5t2h5-64b4t3h5+4b3t4h5+26b2t5h5-12b6h6+24b5th6+22b4t2h6+12b3t3h6-73b2t4h6-10bt5h6-8b5h7-6b4th7+88b3t2h7-68b2t3h7-26bt4h7-29b4h8+16b3th8+42b2t2h8+54bt3h8+3t4h8-28b3h9-12b2th9+88bt2h9+10t3h9-38b2h10-8bth10-11t2h10-28bh11-34th11-17h12;
1061//poly hf=b3t4-b2t5+b4t2h-2b3t3h+2b2t4h-bt5h-2b3t2h2+4bt4h2+b2t2h3-bt3h3+t4h3
1062//+2b2th4-6bt2h4-4t3h4+b2h5-2bth5+2t2h5+4th6+h7;
1063
1064int n,m,sA,sB;
1065 n=6;
1066//n=3;
1067//m=7-n;
1068
1069 m = 12-n;
1070 ideal A = maxideal(n);       ideal B = maxideal(m);
1071 sA =size(A);
1072 sB = size(B);
1073
1074ring R = 0,(b(1..sB),a(1..sA),b,t,h),dp;
1075poly hf = imap(r1,hf);
1076ideal A = imap(r1,A);
1077ideal B = imap(r1,B);
1078matrix aa[sA][1]= a(1..sA);
1079matrix bb[sB][1] = b(1..sB);
1080poly f1= (matrix(A)*aa)[1,1];   //Ansatz for degree n
1081poly f2= (matrix(B)*bb)[1,1];   //Ansatz for degree m=12-n
1082poly g = f1*f2-hf;              //assume hf factors
1083matrix C = coef(g,bth);
1084ideal co = submat(C,2,1..ncols(C));//condition for decomposition, size 91
1085
1086ring R1 = 0,(b(1..sB),a(1..sA)),lp;
1087ideal co = imap(R,co);
1088co=subst(co,a(1),1);          //OE a1=1
1089co = subst(co,b(1),-17);      //co1[91]=b(1)+17
1090//co = subst(co,b(1),1);
1091
1092int aa=timer;
1093list pr=noSolution(co);
1094timer-aa;
1095
1096
1097//n=1 0 sec
1098// keine Primzahl
1099//n=2   628 sec
1100// 2,3,5,7,11,13,17,19,23,29,31,37,43,47,61,89,293,347,367,487,491,3463,7498
1101//n=3  1604 sec
1102// 2,3,5,7,11,13,17,19,23,29,31,59,71,79,101,163,197,211,269,281,541,647,863,
1103// 7129,9041,10343,18413,20857
1104//n=4  8296 sec
1105// 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,61,67,71,73,83,89,101,109,127,
1106   131,149,151,157,163,167,173,179,193,199,211,223,229,251,263,283,307,313,331,
1107   347,373,401,601,607,631,643,701,719,727,829,857,887,971,1279,1361,1453,1721,
1108   1847,2003,2069,2213,2671,2693,3299,3373,3391,3517,3593,3701,3779,4111,4409,
1109   4423,4561,4657,4793,5273,5399,5659,5987,6949,7487,8011,8243,8887,9769,10159,
1110   10177,12007,26347
1111//n=5  4715 sec
1112// 2,3,5,7,11,13,17,19,23,31,41,43,61,283,421,631,1609
1113//n=6 18317 sec
1114// 2,3,5,7,11,13,17,19,29,31,37,41,47,61,71,73,79,89,97,127,131,173,181,223,
1115// 269,953,1039,6151,6343,7823
1116
1117
1118//fuer das andere Wort
1119//n=1  keine Primzahl
1120//n=2  2,3,7,109
1121//n=3  2,3,5,7,11,13,19,29,31,41,43,47,89,139,149,167,173,991,1381,1637,27367
1122
1123*/
Note: See TracBrowser for help on using the repository browser.