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