# source:git/Singular/LIB/decomp.lib@0df59c8

jengelh-datetimespielwiese
Last change on this file since 0df59c8 was 0df59c8, checked in by Hans Schoenemann <hannes@…>, 11 years ago
• Property mode set to 100644
File size: 45.3 KB
Line
1///////////////////////////////////////////////////////////////////////
2version = "$Id: decomp.lib,v 0.99 2012/05/21 q$";
3
4// last changed  21.5.12 C.G. reversal wieder eingefuegt (standalone)
5category = "general";
6info =
7"
8LIBRARY: decomp.lib  Functional Decomposition of Polynomials
9AUTHOR:  Christian Gorzel, University of Muenster
10email: gorzelc@math.uni-muenster.de
11
12
13OVERVIEW:
14@texinfo
15 This library implements functional uni-multivariate decomposition
16 of multivariate polynomials.
17
18    A (multivariate) polynomial f is a composite if it can be written as
19    @math{g \\circ h} where g is univariate and h is multivariate,
20    where @math{\\deg(g), \\deg(h)>1}.
21
22    Uniqueness for monic polynomials is up to linear coordinate change
23@tex
24     $g\\circ h = g(x/c -d) \\circ c(h(x)+d)$.
25@end tex
26
27    If f is a composite, then @code{decompose(f);} returns an ideal (g,h);
28    such that @math{\\deg(g) < \\deg(f)} is maximal, (@math{\\deg(h)\\geq 2}).
29    The polynomial h is, by the maximality of @math{\\deg(g)}, not a composite.
30
31    The polynomial g is univariate in the (first) variable vvar of f,
32    such that deg_vvar(f) is maximal.
33
34    @code{decompose(f,1);} computes a full decomposition, i.e. if f is a
35    composite, then an ideal @math{(g_1,\\dots ,g_m,h)} is returned, where
36    @math{g_i} are univariate and each entry is primitive such that
37    @math{f=g_1\\circ \\dots \\circ g_m\\circ h}.
38
39    If f is not a composite, for instance if @math{\\deg(f)} is prime,
40    then @code{decompose(f);} returns f.
41
42    The command @code{decompose} is the inverse: @code{compose(decompose(f,1))==f}.
43
44    Recall, that Chebyshev polynomials of the first kind commute by composition. @*
45
46    The decomposition algorithms work in the tame case, that is if
47    char(basering)=0 or p:=char(basering) > 0 but deg(g) is not divisible by
48    p.
49    Additionally, it works for monic polynomials over @math{Z} and in some
50    cases for monic polyomials over coefficient rings. @* See
51    @code{is_composite} for examples.  (It also works over the reals but
52    there it seems not be numerical stable.) @*
53
55
56    Univariate decomposition is created, with the additional assumption
57    @math{\\deg(g), \\deg(h)>1}. @*
58
59    A multivariate polynomial f is a composite, if f can be written as
60    @math{g \\circ h}, where @math{g} is a univariate polynomial and @math{h}
61    is multivariate.  Note, that unlike in the univariate case, the polynomial
62    @math{h} may be of degree @math{1}. @*
63    E.g. @math{f = (x+y)^2+ 2(x+y) +1} is the composite of
64    @math{g = x^2+2x+1} and @math{h = x+y}. @*
65
66    If @code{nvars(basering)>1}, then, by default, a single-variable
67    multivariate polynomial is not considered to be the same as in the
68    one-variable polynomial ring; it will always be decomposed. That is: @*
69   @code{> ring r1=0,x,dp;} @*
70   @code{> decompose(x3+2x+1);} @*
71   @code{x3+2x+1} @*
72    but: @*
73   @code{> ring r2=0,(x,y),dp;} @*
74   @code{> decompose(x3+2x+1);} @*
75   @code{_[1]=x3+2x+1} @*
76   @code{_[2]=x} @*
77
78    In particular: @*
79    @code{is_composite(x3+2x+1)==1;}  in @code{ring r1} but  @*
80    @code{is_composite(x3+2x+1)==0;}  in @code{ring r2}. @*
81
82   This is justified by interpreting the polynomial decomposition as an
83   affine Stein factorization of the mapping @math{f:k^n \\to k, n\\geq 2}.
84
85      The behaviour can changed by the some global variables.
86
87    @code{int DECMETH;} choose von zur Gathen's or Kozen-Landau's method.
88@*    @code{int MINS;} compute f = g o h, such that h(0) = 0. @*
89    @code{int IMPROVE;} simplify the coefficients of g and h if f is
90    not monic. @*
91    @code{int DEGONE;} single-variable multivariate are
92    considered uni-variate. @*
93
95
96    Additional information is displayed if @code{printlevel > 0}.
97@end texinfo
98REFERENCES:
99@texinfo
100@tex
101D. Kozen, S. Landau: Polynomial Decomposition Algorithms, \\par
103J. von zu Gathen: Functional Decomposition of Polynomials: the Tame Case,\\par
105 J. von zur Gathen, J. Gerhard:  Modern computer algebra,  \\par
107@end tex
108@end texinfo
109PROCEDURES:
110// decompunivmonic(f,r);
111// decompmultivmonic(f,var,s);
112 decompopts([\"reset\"]);    displays resp. resets global options
113 decompose(f[,1]);           [complete] functional decomposition of poly f
114 is_composite(f);            predicate, is f a composite polynomial?
115 chebyshev(n[,1]);           the nth Chebyshev polynomial of the first kind
116 compose(f1,..,fn);          compose f1 (f2 (...(fn))), f_i polys of ideal
117
118AUXILIARY PROCEDURES:
119 makedistinguished(f,var);   transforms f to a var-distinguished polynomial
120// divisors(n[,1]);            intvec [increasing] of the divisors d of n
121// gcdv(v);                    the gcd of the entries in intvec v
122// maxdegs(f);                 maximal degree for each variable of the poly f
123// randomintvec(n,a,b[,1]);    random intvec size n, [non-zero] entries in {a,b}
124KEYWORDS: Functional decomposition
125";
126
127/*
128 decompunivpoly(poly f,list #)  // f = goh; r = deg g, s = deg h;
129
130 Ablauf ist:
131
132decompose(f)
133| check whether f is the composite by a monomial
134| check whether f is univariate
135| transformation to a distinguished polynomial
136     decompmultivmonic(f,vvar,r)
137         decompunivmonic(f,r)   // detect vvar by maxdegs
138     |lift univariate decomposition
139| back-transformation
140| fulldecompose, iterate
141   | decompuniv for g
142
143*/
144///////////////////////////////////////////////////////////////////////////////
145
146
147proc decompopts(list #)
148"USAGE:  decompopts(); or decompopts(\"reset\");
149RETURN: nothing
150NOTE:
151@texinfo
152        in the first case, it shows the setting of the control parameters;@*
153        in the second case, it kills the user-defined control parameters and@*
154                            resets to the default setting which will then
155                            be diplayed. @* @*
156        int DECMETH; Method for computing the univariate decomposition@*
157                 0 : (default) Kozen-Landau @*
158                 1 : von zur Gathen @*
159
160        int IMPROVE Choice of coefficients for the decomposition @*
161         @math{(g_1,\ldots,g_l,h)} of a non-monic polynomials f. @*
163              and @math{g_2,\ldots,g_l,h} are monic @*
164         1 : (default), content(@math{g_i}) = 1 @*
165
166        int MINS @*
167         @math{f=g\circ h, (g_1,\ldots,g_m,h)} of a non-monic polynomials f.@*
168               0 : g(0) = f(0), h(0) = 0   [ueberlegen fuer complete] @*
169               1 : (default),  g(0)=0, h(0) = f(0) @*
170               2 : Tschirnhaus @*
171
172        int DECORD; The order in which the decomposition will be computed@*
173               0 : minfirst @*
174               1 : (default) maxfirst @*
175
176      int DEGONE; decompose also polynomials built on linear ones @*
177               0 : (default) @*
178               1 :
179@end texinfo
180EXAMPLE: example decompopts; shows an example
181"
182{
183/*
184  siehe Erlaeuterungen, globale Variablen wie im Header angegeben,
185  suchen mit CTRL-S Top::
186  diese eintragen
187*/
188  if (size(#))
189  {
190    if (string(#[1]) == "reset")
191    {
192      if (defined(DECMETH)) {kill DECMETH;}
193      //   if (defined(DECORD)) {kill DECORD;}
194      if (defined(MINS)) {kill MINS;}
195      if (defined(IMPROVE)) {kill IMPROVE;}
196    }
197  }
198
199 if (voice==2)
200 {
201  "";
202  "  === Global variables for decomp.lib === ";
203  "";
204
205  if (!defined(DECMETH)) {" -- DECMETH (int) not defined, implicitly 1";}
206  else
207  {
208    if (DECMETH!=0 and DECMETH!=1) { DECMETH=1; }
209   " -- DECMETH =", DECMETH;
210  }
211/*
212  if (!defined(DECORD)) {" -- DECORD (int) not defined, implicitly 1";}
213  else
214  {
215    if (DECORD!=0 and DECORD!=1) { DECORD=1; }
216   " -- (int) DECORD =", DECORD;
217  }
218*/
219  if (!defined(MINS)) {" -- MINS (int) not defined, implicitly 0";}
220  else
221  {
222    if (MINS!=0 and MINS!=1) { MINS = 0; }
223   " -- (int) MINS =", MINS;
224  }
225
226  if (!defined(IMPROVE)) {" -- IMPROVE (int) not defined, implicitly 1";}
227  else
228  {
229   if (IMPROVE!=0 and IMPROVE!=1) { IMPROVE=1; }
230   " -- (int) IMPROVE =", IMPROVE;
231  }
232 }
233}
234example;
235{ "EXAMPLE:"; echo =2;
236   decompopts();
237}
238///////////////////////////////////////////////////////////////////////////////
239
240//static
241proc decompmonom(poly f, list #)
242"USAGE: decompmonom(f[,vvar]); f poly, vvar poly
243PURPOSE: compute a maximal decomposition in case that
244         f = g o h, where g is univariate and h is a single monomial
245RETURN: ideal, (g,h); g univariate, h monomial if such a decomposition exist,
246        poly, the input, otherwise
247ASSUME: f is non-constant
248EXAMPLE: example decompmonom; shows an example
249"
250{
251  int i,k;
252  poly g;
253
254  poly vvar = var(1);
255  if (size(#)) { vvar = var(rvar(#[1])); }
256
257  //poly vvar = maxdeg(f);
258  poly zeropart = jet(f,0);
259  poly ff = f - zeropart;
260  int mindeg = -deg(ff,-1:nvars(basering));
261  poly minff = jet(ff,mindeg);
262  if (size(minff)>1) { return(f); }
264  minv = minv/gcdv(minv);
265  for (i=1;i<=size(ff);i++)
266  {
268    if (k==0)  { return(f); }
269    else { g = g + leadcoef(ff[i])*vvar^k; }
270  }
271  g = g + zeropart;
272  dbprint("* Sucessfully multivariate decomposed by a monomial"+newline);
273  return(ideal(g,monomial(minv)));
274}
275example
276{ "EXAMPLE:"; echo =2;
277   ring r = 0,(x,y),dp;
278   poly f = subst((x2+x3)^150,x,x2y3);
279   decompmonom(f);
280
281   ring rxyz = 0,(x,y,z),dp;
282   poly g = 1+x2+x3+x5;
283   poly G = subst(g,x,x7y5z3);
284   ideal I = decompmonom(G^50);
285   I[2];
286}
287///////////////////////////////////////////////////////////////////////////////
288
289static proc divintvecs(intvec v,intvec w)
290"USAGE: divintvecs(v,w); v,w intvec, w!=0
291RETURN: int, k if  v = k*w,
292             0 otherwise
293NOTE: if w==0, then an Error message occurs
294EXAMPLE: example divintevcs; shows an example
295"
296{
297  if (w==0) {
298    ERROR("// Error: proc divintvecs: the second argument has to be non-zero.");
299    return(0);
300  }
301  int i=1;
302  while (w[i]==0) { i++; }
303  int k = v[i] div w[i];
304  if (v == k*w) { return(k); }
305  else { return(0); }
306}
307example
308{ "EXAMPLE:"; echo =2;
309  intvec v = 1,2,3;
310  intvec w = 2,4,6;
311  divintvecs(w,v);
312  divintvecs(intvec(3,2,9),v);
313}
314///////////////////////////////////////////////////////////////////////////////
315
316static proc gcdv(intvec v)
317"USAGE: gcdv(v); intvec v
318RETURN: int, the gcd of the entries in v
319NOTE:   if v=0, then gcdv(v)=1 @*
320        this is different from Singular's builtin gcd, where gcd(0,0)==0
321EXAMPLE: example gcdv; shows an example
322"
323{
324  int ggt;
325  int i,n;
326
327  ggt = v[1];
328  for (i=2;i<=size(v);i++)
329  {
330    ggt = gcd(ggt,v[i]);
331  }
332  if (ggt==0)
333  {
334    ggt = 1;
335  }
336  return(ggt);
337}
338example
339{ "EXAMPLE:"; echo =2;
340  intvec v = 6,15,21;
341  gcdv(v);
342  gcdv(0:3);
343}
344///////////////////////////////////////////////////////////////////////////////
345
346static proc divisors(int n,list #)
347"USAGE:  divisors(n); n int
348         divisors(n,1); n int
349RETURN: intvec, the positive divisors of n  @*
350           in decreasing order (default) @*
351           in increasing order in the second case
352EXAMPLE: example divisors; shows an example
353"
354{
355  int i,j;
356  intvec v = 1;
357
358  list l = primefactors(n);
359  list primesl = l[1];
360  list multl = l[2];
361
362  for (i=1;i<=size(primesl);i++)
363  {
364   for (j=1;j<=multl[i];j++)
365     { v = v,primesl[i]*v;}
366  }
367
368  ring rhelp =0,x,dp;        // sort the intvec
369  poly h;
370  for(i=1;i<=size(v);i++)
371  {
372    h = h+x^v[i];
373  }
374  v=0;
375  for(i=1;i<=size(h);i++)
376  {
378  }
379  if (size(#)) {
380    return(intvec(v[size(v)..1]));
381  }
382
383  return(v);
384}
385example
386{ "EXAMPLE:"; echo = 2;
387  divisors(30);
388  divisors(-24,1);
389}
390///////////////////////////////////////////////////////////////////////////////
391//
392// Dies wirkt sich nur aus wenn Brueche vorhanden sind?!
393// Laeuft dann so statt cleardenom usw. problemlos ueber Z,Z_m
394// ansehen.
395//
396static proc improvecoef(poly g0,poly h0,number lc)
397"USAGE: improvecoef(g0,h0,lc); g0, h0 poly; lc number
398RETURN: poly, poly, number
399ASSUME: global ordering
400EXAMPLE: example improvecoef; shows an example
401"
402{
403  int Zcoefs = find(charstr(basering),"integer");
404  poly vvar = var(univariate(g0));
406  number denom;
407
408  if (Zcoefs and lch0<0)  // da cleardenom fuer integer buggy ist.
409  {
410    h0 = h0/(-1);
411    denom = -1;
412  }
413  else
414  {
415   h0 = cleardenom(h0);
417  }
418  g0 = subst(g0,vvar,1/denom*vvar);
419  g0 = lc*g0;
421  g0= 1/lc*g0;
422  return(g0,h0,lc);
423}
424example
425{ "EXAMPLE:"; echo = 2;
426   ring r = 0,x,dp;
427
428   poly g = 3x2+5x;
429   poly h = 4x3+2/3x;
430   number lc = 7;
431
432   improvecoef(g,h,lc);
433}
434///////////////////////////////////////////////////////////////////////////////
435
436proc compose(list #)
437"USAGE: compose(f1,...,fn); f1,...,fn poly
438        compose(I); I ideal,  @*
439ASSUME: the ideal consists of n=ncols(I) >= 1 entries, @*
440          where I[1],...,I[n-1] are univariate in the same variable @*
441          but I[n] may be multivariate.
442RETURN: poly, the composition  I[1](I[2](...I[n]))
443NOTE:   this procedure is the inverse of decompose
444EXAMPLE: example compose; shows some examples
445SEE: decompose
446"
447{
448  def d = basering;   // Ohne dies kommt es zu Fehler, wenn auf Toplevel
449                      // ring r definiert ist.
450
451  ideal I = ideal(#[1..size(#)]);
452  int n=ncols(I);
453  poly f=I[1];
454  map phisubst;
455  ideal phiid = maxideal(1);
456
457  int varnum = univariate(f);
458
459  if (varnum<0) {
460    " // the first polynomial is a constant";
461    return(f);
462  }
463  if (varnum==0 and n>1) {
464    " // the first polynomial is not univariate";
465    return(f);
466  }
467  // Hier noch einen Test ergaenzen
468
469  poly vvar = var(varnum);
470
471  for(int i=2;i<=n;i++)
472  {
473    phiid[varnum]=I[i];
474    // phisubst=d,phiid;
475    phisubst=basering,phiid;
476    f = phisubst(f);
477  }
478  return(f);
479}
480example
481{ "EXAMPLE:"; echo =2;
482   ring r = 0,(x,y),dp;
483   compose(x3+1,x2,y3+x);
484   // or the input as one ideal
485   compose(ideal(x3+1,x2,x3+y));
486}
487///////////////////////////////////////////////////////////////////////////////
488
489proc is_composite(poly f)
490"USAGE: is_composite(f); f poly
491RETURN: int @*
492     1, if f is decomposable @*
493     0, if f is not decomposable @*
494 -1, if char(basering)>0 and deg(f) is  divisible by char(basering) but no
495      decomposition has been found.
496NOTE:  The last case means that it could exist a decomposition f=g o h  with
497char(basering)|deg(g), but this wild case cannot be decided by the algorithm.@*
498       Some additional information will be displayed when called by the user.
499EXAMPLE: example is_composite; shows some examples
500"
501{
502  int d = deg(f,nvars(basering));
503  int cb = char(basering);
504
505  if (d<1)
506  {
507    " The polynomial is constant ";
508    return(0);
509  }
510  if (d==1)
511  {
512    " The polynomial is linear ";
513    return(0);
514  }
515
516  if (nvars(basering)==1 and d==prime(d))
517  {
518    " The degree is prime.";
519    return(0);
520  }
521
522  if (nvars(basering)>1 and univariate(f))  // and not(defined(DEGONE))
523  {
524    return(1);
525  }
526
527  // else try to decompose
528  int nc = ncols(ideal(decompose(f)));
529
530  if (cb > 0)   // check the not covered wild case
531  {
532    if ((d mod cb == 0) and (nc == 1))
533    {
534      if (voice==2)
535      {
536        "// -- Warning: wild case, cannot decide whether the polynomial has a";
537        "// -- decomposition goh with deg(g) divisible by char(basering) = "
538          + string(cb) + ".";
539       }
540      return(-1);
541    }
542  }
543  // in the tame case, decompose gives the correct result
544  return(nc>1);
545}
546example
547{ "EXAMPLE:"; echo =2;
548
549   ring r0 = 0,x,dp;
551
552   is_composite(2x2+x+1);     // prime degree
553   // -----------------------------------------------------------------------
554   // polynomial ring with several variables
555   ring R = 0,(x,y),dp;
556   // -----------------------------------------------------------------------
557   // single-variable multivariate polynomials
558   is_composite(2x+1);
559   is_composite(2x2+x+1);
560   // -----------------------------------------------------------------------
561   // prime characteristic
562   ring r7 = 7,x,dp;
563   is_composite(compose(ideal(x2+x,x14)));     // is_composite(x14+x7);
564   is_composite(compose(ideal(x14+x,x2)));     // is_composite(x14+x2);
565
566}
567///////////////////////////////////////////////////////////////////////////////
568
569proc decompose(poly f,list #)
570"USAGE:  decompose(f); f poly
571        decompose(f,1); f poly
572RETURN: poly, the input, if f is not a composite
573        ideal, if the input is a composite
574NOTE: computes a full decomposition if called by the second variant
575EXAMPLE: example decompose; shows some examples
576SEE: compose
577"
578{
579 if (!defined(IMPROVE)){ int IMPROVE = 1; }
580 if (!defined(MINFIRST)){ int MINFIRST = 0; }
581 int fulldecompose;
582
583 if (size(#)) {        // cf. ERROR-msg in randomintvec
584   if (typeof(#[1])=="int") {
585       fulldecompose = (#[1]==1);
586   }
587 }
588
589 int m,iscomposed;
590 int globalord = 1;
591 ideal I;
592
593// --- preparatory stuff ----------------------------------------------------
594  // The degree is not independent of the term order
595 int n = deg(f,1:nvars(basering));
596 int varnum = univariate(f);  // to avoid transformation if f is univariate
597
598 // if (deg(f)<=1) {return(f);} //steigt automatisch bei der for-schleife aus m = 2
599 if (n==prime(n) and nvars(basering)==1
600     //  or (varnum>0 and nvars(basering))
601  ) {return(f);}
602
603 if (varnum<0)
604 {
605   ERROR("// -- Error proc decompoly: the polynomial is constant.");
606 }
607 //--------------------------------------------------------------------------
608
609 int minfirst = MINFIRST!=0;
610 list mdeg;
611 intvec maxdegv,degcand;
612
613 // -- switch to global order, necessary for division -- // Weiter nach oben
614 if (typeof(attrib(basering,"global"))!="int") {
615   globalord = 0;
616 }
617 else {
618   globalord = attrib(basering,"global");
619 }
620
621 if (!globalord) {
622   def d = basering;
623   list ll = ringlist(basering);
624   ll[3] = list(list("dp",1:nvars(basering)),list("C",0));
625   def rneu = ring(ll);
626   setring rneu;
627   poly f = fetch(d,f);
628   ideal I;
629 }
630 // -----------------------------------------------------------------------
631
632 map phiback;
633 poly f0,g0,h0,vvar;
634 number lc;
635 ideal J;   // wird erst in fulldecompose benoetigt
636
637 // --- Determine the candidates for deg(g) a decreasing sequence of divisors
638 poly lf = jet(f,n)-jet(f,n-1);
639 //"lf = ",lf;
640 if (size(lf)==1)         // the leading homogeneous part is a monomial
641 {
643 }
644 else
645 {
646   degcand = divisors(n);            // Das ist absteigend
647 }
648
649 if(printlevel>0) {degcand;}
650
651 // --- preparatory steps for the multivariate case -------------------------
652
653 if (varnum>0)            // -- univariate polynomial
654 {
655  vvar = var(varnum);
656  f0 = f;   // save f
657 }
658 else  // i.e. multivariate (varnum==0),the case varnum < 0 is excluded above
659 {
660  // -- find variable with  maximal degree
661  mdeg = maxdegs(f);
662  maxdegv = mdeg[2];
663  varnum = maxdegv[2];
664  vvar = var(varnum);
665  phiback = maxideal(1);
666
667// special case, the polynomial is  a composite of a single monomial //20.6.10
668   if (qhweight(f)!=0) { I = decompmonom(f,vvar); }
669   iscomposed = size(I)>1;
670   if (iscomposed)          // 3.6.11 - dies decompmonom
671   { //I;
672      ideal J = decompunivmonic(I[1],deg(I[1]));
673      I[2]= subst(J[2],vvar,I[2]);
674      I[1] = J[1];
675      //I;
676   }
677
678   if (!iscomposed) // -- transform into a distinguished polynomial
679   {
680     f0,phiback = makedistinguished(f,vvar);
681   }
682 }
683 // ------ Start computation ------------------------------------------------
684 // -- normalize and  save the leading coefficient
685 lc = 1;
686 //f0;
687 //"vvar = ",vvar;
688
689 // --- 11.4.11 hier auch noch gewichteten Grad beruecksichtigen ? --
690
691  if (!iscomposed) { lc = leadcoef(coeffs(f0,vvar)[deg(f0)+1,1]); } // 20.6.10
692
693  // if Z, Z_m, and f is not monic (and  content !=1) // if (f0/lc*lc!=f0)
694  if (find(charstr(basering),"integer") and not(lc==1 or lc==-1)) // 6.4.11
695  {
696    ERROR("// -- Error proc decompose: Can not decompose non-monic polynomial over Z!");
697  }
698
699  if (lc!=1){ f0 = 1/number(lc)*f0;}      // --- normalize the polynomial
700
701  // -- Now the input is prepared to be monic and vvar-distinguished
702  //----------------------------------------------------------------
703  m = 1;
704
705  // --- Special case: a multivariate can be composite of a linear polynom
706  if (univariate(f) and nvars(basering)==1) // 11.8.09 d.h.
707  {    // --- if univariate ----------------------------------------
708    if(minfirst) {degcand = divisors(n,1);} // dies ist aufsteigend
709    m = 2;                                  // skip first entry
710  }
711  // if decomposed as the decomposition with a monomial
712  // then skip the multivariate process // 20.6.10 detected as decompmonomial
713  if (iscomposed) { degcand = 1; }
714
715  if (printlevel>0 and !iscomposed) { "* Degree candidates are", degcand; }
716
717  // -- check succesively for each candidate
718  // whether f is decomposable with deg g = r
719
720  for(;m<size(degcand);m++)   // decreasing
721  { //r = degcand[m];
722    I =  decompmultivmonic(f0,vvar,degcand[m]);
723    if (size(I)>1)
724    {
725     iscomposed = 1;
726     break;
727    }
728  }
729 // -- all candidates have be checked but f is primitive
730 if(!iscomposed) {
731   if (!globalord) { setring d; }   // restore old ring
732   dbprint("** not decomposable: linear / not tame / prime degree --");
733   return(f);
734 }
735
736 // -- the monic vvar-distinguished polynomial f0 is decomposed -------
737 // -- retransformation for the multivariate case ---------------------
738  g0,h0 = I;
739
740  if (!univariate(f)) { h0 = phiback(h0);}
741
742  if (IMPROVE) { g0,h0,lc=improvecoef(g0,h0,lc);}  //  ueber switch
743  I = h0;
744
745 // -- Full decomposition: try to decompose g further ------------------
746  if (fulldecompose) {
747    dbprint(newline+"** Compute a complete decomposition");
748    while (iscomposed) {
749     iscomposed=0;
750     degcand=divisors(deg(g0,1:nvars(basering)));  // absteigend
751     if (printlevel> 0) { "** Degree candidates are now: ", degcand; }
752     for (m=2;m<size(degcand);m++) //OK, ergibt lexicographically ..
753     {
754       J =decompunivmonic(g0,degcand[m]); /*      J =decompuniv(g0);*/
755       g0 = J[1];
756       h0=J[2];
757       iscomposed = deg(h0,1:nvars(basering))>1;
758       if (iscomposed) {
759         if (IMPROVE) { g0,h0,lc=improvecoef(g0,h0,lc); }  //  ueber switch
760         I = h0,I;
761         break;
762       }
763      }
764    }
765  dbprint("** completely decomposed"+newline);
766  }
767  I = lc*g0,I;
768  if (!globalord) {
769   setring d;
770   I = fetch(rneu,I);
771  }
772  return(I);
773}
774example
775{ "EXAMPLE:"; echo =2;
776   ring r2 = 0,(x,y),dp;
777
778   decompose(((x3+2y)^6+x3+2y)^4);
779
780   // complete decomposition
781   decompose(((x3+2y)^6+x3+2y)^4,1);
782   // -----------------------------------------------------------------------
783   // decompose over the integers
784   ring rZ = integer,x,dp;
785   decompose(compose(ideal(x3,x2+2x,x3+2)),1);
786   // -----------------------------------------------------------------------
787   // prime characteristic
788   ring r7 = 7,x,dp;
789   decompose(compose(ideal(x2+x,x7)));   // tame case
790   // -----------------------------------------------------------------------
791   decompose(compose(ideal(x7+x,x2)));   // wild case
792   // -----------------------------------------------------------------------
793   ring ry = (0,y),x,dp;     // y is now a parameter
794   compose(x2+yx+5,x5-2yx3+x);
795   decompose(_);
796
797  // Usage of variable IMPROVE
798  ideal J = x2+10x, 64x7-112x5+56x3-7x, 4x3-3x;
799  decompose(compose(J),1);
800  int IMPROVE=0;
801  exportto(Decomp,IMPROVE);
802  decompose(compose(J),1);
803}
804///////////////////////////////////////////////////////////////////////////////
805/*   ring rt =(0,t),x,dp;
806   poly f = 36*x6+12*x4+15*x3+x2+5/2*x+(-t);
807   decompose(f);
808*/
809
810
811// Dies gibt stets ein ideal zurueck, wenn f composite ist
812// gibt das polynom zurueck, wenn es primitiv ist
813// static
814proc decompmultivmonic(poly f,poly vvar,int r)
815"USAGE:  decompmultivmonic(f,vvar,r); f,vvar poly; r int
816RETURN:  ideal, I = ideal(g,h) if f = g o h with deg(g) = r@*
817           poly f, if  f is not a composite or char(basering) divides r
818ASSUME:  f is monic and distinguished w.r.t. vvar,
819         1<=r<=deg(f) is a divisor of deg(f)
820         and char(basering) does not divide r.
821EXAMPLE: example decompmultivmonic; shows an example
822"
823{
824  def d = basering;
825 int i,isprimitive;
826 int m = nvars(basering);
827 int n = deg(f);
828 int varnum = rvar(vvar);
829 intvec v = 1:m;   // weight-vector for jet
830   v[varnum]=0;
831 int s = n div r;
832 // r = deg g; s = deg h;
833
834 poly f0 = f;
835 poly h,h0,g,gp,fgp,k,t,u;
836 ideal I,rem,phiid;
837 list l;
838 map phisubst;
839
840// -- entscheidet intern, abhaengig von der Anzahl der Ringvariablen,
841// -- ob f0 primitive ist.
842// " r = ",r;
843
844 if (s*r!=n)
845 {
846   ERROR("// -- Error proc decompmultivmonic: r = "+string(r)+
847             " does not divide deg(f) = "+string(n)+".");
848 }
849
850 int cb = char(basering);  // oder dies in decompunivmonic
851 if (cb>0)
852 {
853     if (r mod cb == 0)
854     {
855       if (voice == 2)
856       {
857        "// Warning: wild case in characteristic " + string(cb) +
858         ". We cannot decide";
859         "// whether a decomposition goh with deg(g) = " + string(r)+
860         " exists.";"";
861       }
862       return(f);
863     }
864 }
865//---------------------------------------------------------------------------
866
867 for (i=1;i<=m;i++)
868 {
869  if (i!=varnum) {f0 = subst(f0,var(i),0);}
870 }
871 //" f0 = ",f0;
872 // f0 ist nun das univariate
873
874 // 24.3.09   // 11.8.09  nochmals ansehen
875 if (r==deg(f0)) // the case of a linear multivarcomposite
876 {
877   dbprint("** try to decompose in linear h, deg g = "+string(r));
878  I = f0,vvar;   // Das ist hier wichtig
879 }
880 else   // find decomposition of the univariate f0
881 {
882  I = decompunivmonic(f0,r);
883  // dbprint(" ** monic decomposed");//" I = ";I;
884
885  isprimitive=(deg(I[2])==1);
886  if (isprimitive) {return(f);}
887 }
888
889//---- proceed in the multivariate case
890//---- lift the univariate decomposition
891 if (!univariate(f))
892 {
893   dbprint("* Lift the univariate decomposition");
894   g,h0 = I;
895   k = h0;
896   gp = diff(g,vvar);
897
898   // -- This is substitution ----
899   // t = substitute(gp,vvar,h0);
900   phiid = maxideal(1);
901   phiid[varnum]=h0;
902   phisubst=basering,phiid;
903    t = phisubst(gp);
904   // -- substitution ende
905   fgp = 1;
906   i = 0;
907   while(fgp!=0)
908   {
909     i++;
910     // -- This  is substitution ----
911     //gp = substitute(g,vvar,k);
912     phiid[varnum]=k;
913     phisubst=basering,phiid;
914     gp = phisubst(g);
915     // --  substitution ende
916
917     fgp = f - gp;
918     u = jet(fgp,i,v) - jet(fgp,i-1,v);    // oder mit reduce(maxideal(x))
919     l = division(u,t);                    // die kleineren Terme abschneiden
920     rem = l[2];
921     u = l[1][1,1];    // the factor
922     if (rem!=0)
923     {
924       isprimitive = 1;
925       break;
926     }
927     k = k + u;
928   }
929   h = k;
930   I = g,h;
931   //"decomposed as =";
932   //I;
933 }
934 if (isprimitive) {
935   dbprint(">>> not multivariate decomposed"+newline);
936  return(f);
937 }
938 else {
939   dbprint("* Sucessfully multivariate decomposed"+newline);
940  return(I);
941 }
942}
943example
944{ "EXAMPLE:"; echo = 2;
945   ring r = 0,(x,y),lp;
946   poly f = 3xy4 + 2xy2 + x5y3 + x + y6;
947   decompmultivmonic(f,y,2);
948
949   ring rx = 0,x,lp;
950   decompmultivmonic(x8,x,4);
951}
952///////////////////////////////////////////////////////////////////////////////
953//static
954proc decompunivmonic(poly f,int r)
955"USAGE: decompunivmonic(f,r); f poly, r int
956RETURN:  ideal,  (g,h) such that f = goh and deg(g) = r
957         poly f, if such a decomposition does not exist.
958ASSUME:  f is univariate, r is a divisor of deg(f)      @*
959         and char(basering) does not divide r in case that char(basering) > 0.
960         global order of the basering is assumed.
961EXAMPLE: example decompunivmonic; shows an example
962"
963{
964 int d = deg(f);
965 int s;    // r = deg g; s = deg h;
966 int minf,mins;
967 int iscomposed = 1;
968
969 if (!defined(MINS)) { int MINS = 0; }
970 if (!defined(DECMETH)) { int DECMETH = 1; }
971 int savedecmeth = DECMETH;
972 int Zcoefs =charstr(basering)=="integer";//find(charstr(basering),"integer");
973
974 number cf;
975 poly h,g;
976 ideal I;
977 matrix cc;
978
979 // --- Check input and create the results for the simple cases
980
981 if (deg(f)<1){return(ideal(f,var(1)));}  // wird dies aufgerufen?
982 //-------------------------
983
984 int varnum = univariate(f);
985
986 if (varnum==0)
987 {
988   "// -- The polynomial is not univariate";
989   return(f);
990 }
991
992 poly vvar = var(varnum);
993 I = f,vvar;
994
996 {
997     "// -- Error proc decompunivmonic: the polynomial is not monic.";
998     return(f);
999 }
1000 /* Dies einklammern, wenn (x+1)^2 zerlegt werden sollte
1001  // aus decompose heraus, wird dies gar nicht aufgerufen!
1002 if (deg(f)==1 or deg(f)==prime(deg(f)))
1003 {
1004   "// -- The polynomial is not a composite.";
1005  return(I);
1006 }
1007 */
1008 /* ---------------------------------------------------- */
1009 s = d div r;
1010
1011 if (d!=s*r)
1012 {
1013   ERROR("// -- Error proc decompunivmonic: the second argument does not divide deg f.");
1014 }
1015 int cb = char(basering);
1016 if (cb>0)
1017 {
1018     if (r mod cb ==0)
1019     {
1020         "wild case: cannot determine a decomposition";
1021         return(I);
1022     }
1023 }
1024// -------------------------------------------------------------------------
1025 // The Newton iteration only works over coefficient *fields*
1026 // Therefore use in this case the Kozen-Landau method i.e. set DECMETH = 1;
1027 if (savedecmeth==0 and Zcoefs) { DECMETH=1; }
1028
1029// -- Start the computation ----------------------------------------------
1030
1031  dbprint("* STEP 1: Determine h");
1032  dbprint(" d = deg f = " +string(n) + " f = goh"," r = deg g = "+string(r),
1033          " s = deg h = " +string(s));
1034  int tt = timer;
1035
1036  if(DECMETH==1) {  // Kozen-Landau
1037    dbprint("* Kozen-Landau method");
1038
1039   // Determine ord(f);
1040   //cc  = coef(f,vvar);  // extract coefficents of f
1042
1043   // dbprint("time: "+string(timer-tt)); tt = timer;
1044   // minf = deg(cc[1,ncols(cc)]);   // 11.8.09 Doch OK.
1045   minf = -deg(f,-1:nvars(basering));  // this is local ord 15.3.10
1046
1047   // oder: mins = 1;  if (minf) { .. dies .. }
1048   mins = (minf div r) + (minf mod r) > 0;  // i.e. ceil(minf/r)
1049
1050   if (mins==0 and MINS) { mins=1; } // omit the constant term i.e. h(0) = 0
1051
1052   dbprint("** min f = "+string(minf) + " | min s = "+string(mins) +
1053           " | s-mins = "+ string(s-mins));
1054
1055   // Dies wird wohl nicht benoetigt.
1056   //  int minr=  (minf div s) + ((minf mod s)>0); // ceil
1057   dbprint("** extract the coeffs ");
1058   cc  = coeffs(f,vvar);
1059
1060   dbprint("time: "+ string(timer -tt));
1061
1062   h = vvar^s;
1063   for (int j=1;j<=s-mins;j++)
1064   {
1065/*
1066     timer = 1;H = Power(h,r);  "Power H"; timer;
1067     timer = 1;G = h^r;  "h^r"; timer;
1068*/
1069    cf = (number(cc[d-j+1,1])-number(coeffs(h^r,vvar)[d-j+1,1]));
1070
1071//    d-j+1,"cf =",cf, " r= ",r;
1072// dbprint("*** "+ string(d-j+1) + " cf = "+string(cf) + " r= "+string(r));
1073
1074    if (Zcoefs) { if (bigint(cf) mod r != 0) { iscomposed = 0; break; }}
1075    cf = cf/r;
1076
1077    //else { cf = cf/r; }
1078    h = h + cf*vvar^(s-j);
1079//    " h = ",h;
1080   }
1081  } else {
1082   dbprint("* von zur Gathen-method");
1083   //  "f=",f;
1084   h = reversal(newtonrroot(reversal(f,d),r,s+!MINS),s,vvar);  // verdreht OK
1085   //   " h = ",h;
1086   dbprint("* END STEP 1: time: "+string(timer -tt));
1087  }
1088  DECMETH=savedecmeth;   // restore the original method
1089
1090  if (iscomposed == 0) {
1091    dbprint("** Failed in STEP 1: not decomposed with deg h = "+string(s)+newline);
1092    return(I);
1093  }
1094
1095  // -- Step 2: try to rewrite f as a sum of powers of h ---
1096  dbprint("* STEP 2: Determine g");
1097  poly H = h^r;
1098  int dalt = r;
1099  int ds;
1100  number c;
1101  while (d >= 0)   // i.e. f!=0
1102  {
1103    //dbprint("d = ",d);
1104    ds = d div s;
1105    if (ds * s !=d)   // d mod s != 0, i.e. remaining f is a power of h
1106    {
1107      iscomposed = 0;
1108      break;
1109    }
1111    g = g + c*vvar^ds;
1112    H = division(H,h^(dalt - ds))[1][1,1];   // 10.3.10
1113    // H = H / h^(dalt - ds);
1114    f = f - c*H;
1115    //"f = ",f;
1116
1117    dalt = ds;
1118    d = deg(f);
1119  }
1120  dbprint("* END STEP 2: time: "+string(timer -tt));
1121  if (iscomposed)
1122  {
1123   dbprint("** Sucessfully univariate decomposed with deg g = "+string(r)+newline);
1124    I = g,h;
1125  } else {
1126   dbprint("** Failed in STEP 2: not decomposed with deg g = "+string(r)+newline);
1127  }
1128
1129 return(I);
1130}
1131example
1132{ "EXAMPLE:"; echo = 2;
1133  ring r=0,(x,y),dp;
1134  decompunivmonic((x2+x+1)^3,3);
1135  decompunivmonic((x2+x)^3,3);
1136
1137  decompunivmonic((y2+y+1)^3,3);
1138}
1139///////////////////////////////////////////////////////////////////////////////
1140// aus polyaux.lib
1141proc reversal(poly f,list #)
1142"USAGE: reversal(f); f poly
1143         reversal(f,k); f poly, k int
1144         reversal(f,k,vvar); f poly, k int, vvar poly (a ring variable)
1145RETURN: poly, the reversal x^k*f(1/x) of the input f
1146ASSUME: f is univariate and that  k>=deg(f)
1147@*      since no negative exponents are possible in Singular
1148@*      if k<deg(f) then k = deg(f) is used
1149NOTE: reversal(f); is by default reversal(f,deg(f));
1150      the third variant is needed if f is a non-zero constant and k>0 @*
1151@*    reversal is only idempotent,
1152@*    if called twice with the deg(f) as second argument
1153EXAMPLE: example reversal; shows an example
1154"
1155{
1156  int k = 0;
1157  poly vvar = var(1);
1158
1159  if (size(#)) {
1160    k  = #[1] - deg(f) ;
1161    if (k<0) { k=0; }
1162    if (size(#)==2){            // check whether second optional argument
1163      vvar = var(univariate(#[2]));     // is a ring variable
1164    }
1165  }
1166
1167  int varnum = univariate(f);
1168
1169  if (varnum==0) {
1170    ERROR("// -- the input is not univariate.");
1171  }
1172  if (varnum<0) {  // the polynomial is  constant
1173    return(f*vvar^k);
1174  }
1175
1176  def d = basering;
1177  list l = ringlist(d);
1178  list varl = l[2];
1179  varl = insert(varl,"@z",size(varl));
1180  l[2] = varl;
1181  def rnew = ring(l);
1182  setring rnew;
1183  poly f = fetch(d,f);
1184  f = subst(homog(f,@z),var(varnum),1,@z,var(varnum))*var(varnum)^k;
1185
1186  setring d;
1187  f = fetch(rnew,f);
1188  return(f);
1189}
1190example
1191{ "EXAMPLE:"; echo = 2;
1192   ring r = 0,x,dp;
1193   poly f = x3+2x+5;
1194   reversal(f);
1195   // the same as
1196   reversal(f,3);
1197   reversal(f,5);
1198
1199   poly g  = x3+2x;
1200   reversal(g);
1201
1202   // Not idempotent
1203   reversal(reversal(g));
1204
1205   // idempotent
1206   reversal(reversal(g,deg(g)),deg(g));
1207   // or for short
1208   // reversal(reversal(g),deg(g));
1209}
1210///////////////////////////////////////////////////////////////////////////////
1211// aus polyaux.lib
1212proc newtonrroot(poly f,int r,int l)
1213"USAGE: newtonrroot(f,r,l); f poly; r int; l int
1214RETURN:  poly h, the solution of h^r = f modulo vvar^l
1215ASSUME: f(0) = 1
1216NOTE: this uses p-adic Newton iteration. It is the adaption of Algorithm 9.22@*
1217        of von zur Gathen & Gerhard p. 264 for the special case: phi = Y^r - f
1218EXAMPLE: example newtonrroot; shows some examples
1219"
1220{
1221 // phi = Y^r - f
1222
1223 poly g = 1;  // start polynomial
1224
1225 poly s =  1/number(r); // initial solution
1226 int i = 2;
1227 //"s initial",s;
1228
1229 while(i<l) {
1230   // "iteration i",i;
1231
1232  //  g = (g -(g^r-f)*s) mod x^i;
1233  g = jet((g -(g^r-f)*s), i-1);
1234  //  s = 2*s - (r*g^(r-1)*s^2) mod x^i;
1235  s = jet(2*s - (r*g^(r-1)*s^2),i-1);
1236  // "s is now ",s;
1237
1238  i = 2*i;
1239 }
1240 //"return newtonrroot";
1241 //jet((g -(g^r-f)*s),l-1);
1242
1243 return(jet((g -(g^r-f)*s),l-1));
1244}
1245example
1246{ "EXAMPLE:"; echo = 2;
1247   ring r = 0,x,dp;
1248
1249   ring r3 = 3,x,dp;
1250   poly f = x+1;
1251   // determine square root of f modulo x^4
1252   poly g = newtonrroot(f,2,4);
1253   g;
1254   g^2;
1255   ring R = (0,b,c,d),x,ds;
1256//   poly f = 1 + bx +cx2+dx3;
1257   poly f = 1 + 5bx +5cx2+5dx3;
1258   poly g2 = newtonrroot(f,2,4);
1259   g2;
1260   f-g2^2;
1261   poly f5 = 1 +5*(bx+cx2+dx3);
1262   poly g5 = newtonrroot(f5,5,4);
1263   g5;
1264   f5-g5^5;
1265   // Multivariate polynomials
1266   ring r = 0,(x,y,z),ds;
1267   ring r2 =(0,a,b,c,d,e),(x,y),ds;
1268//   poly f = 1 +ax+by+cx2+dxy+ey2;
1269   poly f3 = 1 +9*(ax+by+cx2+dxy+ey2);
1270   poly g3 = newtonrroot(f3,3,4);
1271   jet(g3^3-f3,5);
1272}
1273///////////////////////////////////////////////////////////////////////////////
1274
1275static proc randomintvec(int n,int a,int b,list #)
1276"USAGE: randomintvec(n,a,b); n,a,b int;
1277        randomintvec(n,a,b,1); n,a,b int;
1278RETURN: intvec, say v, of length n
1279        with entries a<=v[i]<=b, in the first case, resp.
1280        with entries a<=v[i]<=b, where v[i]!=0, in the second case
1281NOTE:   a<=b should be satisfied, otherwise always v[i]=b (due to random).
1282EXAMPLE: example randomintvec; shows some examples
1283"
1284{
1285  int i,randint,nozeroes;
1286  intvec v;
1287
1288  if (size(#)) {
1289   if (typeof(#[1])!="int") {
1290     ERROR("4th argument can only be an integer, assumed 1.");
1291   }
1292    nozeroes = #[1]==1;
1293  }
1294
1295  for (i=1;i<=n;i++)
1296  {
1297    randint = random(a,b);
1298    while (nozeroes and randint==0) { randint = random(a,b); }
1299    v[i] = randint;
1300  }
1301  return(v);
1302}
1303example
1304{ "EXAMPLE:"; echo = 1;
1305  int randval = system("--random");  // store initial value
1306  system("--random",0815);
1307  echo = 2;
1308  randomintvec(7,-1,1);   // 7 entries in {-1,0,1}
1309  randomintvec(7,-1,1,1); // 7 entries either -1 or 1
1310  randomintvec(3,-10,10);
1311  echo = 1;
1312  system("--random",randval);      // reset random generator
1313}
1314///////////////////////////////////////////////////////////////////////////////
1315
1316proc makedistinguished(poly f,poly vvar)
1317"USAGE:  makedistinguished(f,vvar); f, vvar poly; where vvar is a ring variable
1318RETURN:  (poly, ideal): the transformed polynomial and an ideal defining
1319                       the map which reverses the transformation.
1320PURPOSE: let vvar = var(1). Then  f is transformed by a random linear
1321         coordinate change
1322         phi = (var(1), var(2)+c_2*vvar,...,var(n)+c_n*vvar)  @*
1323         such that phi(f) = f o phi becomes distinguished with respect
1324         to vvar. That is, the new polynomial contains the monomial vvar^d,
1325         where d is the degree of f. @*
1326         If already f is distinguished w.r.t. vvar, then f is left unchanged
1327         and the re-transformation is the identity.
1328NOTE 1:  (this proc correctly works independent of the term ordering.)
1329         to apply the reverse transformation, either define a map
1330         or use substitute (to be loaded from poly.lib).
1331NOTE 2:  If p=char(basering) > 0, then there exist polynomials of degree d>=p,
1332         e.g. @math{(p-1)x^p y + xy^p}, that cannot be transformed to a
1333         vvar-distinguished polynomial. @*
1334         In this case, *p random trials will be made and the proc
1335         may leave with an ERROR message.
1336EXAMPLE: example makedistinguished; shows some examples
1337"
1338{
1339  def d = basering;   // eigentlich ueberfluessig // wg Bug mit example part
1340  map phi;            // erforderlich
1341  ideal Db= maxideal(1);
1342  int n,b = nvars(basering),1;
1343  intvec v= 0:n;
1344  intvec w =v;
1345  int varnum = rvar(vvar);
1346  w[varnum]=1;   // weight vector for deg
1347
1348  poly g = f;
1349  int degg = deg(g);
1350
1351  int count = 1; // limit the number of trials in char(p) > 0
1352  //int count =2*char(basering);
1353
1354  while(deg(g,w)!=degg and (count-2*char(basering)))  // do a transformation
1355  {
1356    v = randomintvec(n,-b,b,1);  // n non-zero entries
1357    v[varnum] = 0;
1358    phi = d,ideal(matrix(maxideal(1),n,1) + var(varnum)*v);  // transformation;
1359    g = phi(f);
1360    b++;    // increase the range for the random values
1361    // count--;
1362    count++;
1363  }
1364  if (deg(g,w)!=degg) {
1365   ERROR("it could not be transform to a "+string(vvar)+"-distinguished polynomial.");
1366  }
1367  Db = ideal(matrix(maxideal(1),n,1) - var(varnum)*v); // back transformation
1368  return(g,Db);
1369}
1370example
1371{ "EXAMPLE:";
1372  int randval = system("--random");  // store initial value
1373  system("--random",0815);
1374  echo = 2;
1375
1376   ring r = 0,(x,y),dp;
1377   poly g;
1378   map phi;
1379   // -----------------------------------------------------------------------
1380   // Example 1:
1381   poly f = 3xy4 + 2xy2 + x5y3 + x + y6;    // degree 8
1382   // make the polynomial y-distinguished
1383   g, phi = makedistinguished(f,y);
1384   g;
1385   phi;
1386
1387   // to reverse the transformation apply the map
1388   f == phi(g);
1389
1390   //  -----------------------------------------------------------------------
1391   // Example 2:
1392   // The following polynomial is already x-distinguished
1393   f = x6+y4+xy;
1394   g,phi = makedistinguished(f,x);
1395   g;                         // f is left unchanged
1396   phi;                       // the transformation is the identity.
1397   echo = 1;
1398
1399   system("--random",randval);      // reset random generator
1400   // -----------------------------------------------------------------------
1401   echo = 2;
1402   // Example 3:    // polynomials which cannot be transformed
1403   // If p=char(basering)>0, then (p-1)*x^p*y + x*y^p factorizes completely
1404   // in linear factors, since (p-1)*x^p+x equiv 0 on F_p. Hence,
1405   // such polynomials cannot be transformed to a distinguished polynomial.
1406
1407   ring r3 = 3,(x,y),dp;
1408   makedistinguished(2x3y+xy3,y);
1409}
1410///////////////////////////////////////////////////////////////////////////////
1411
1412static proc maxdegs(poly f)
1413"USAGE: maxdegs(f); f poly
1414RETURN:  list of two intvecs
1415         _[1] intvec: degree for variable i, 1<=i<=nvars(basering) @*
1416         _[2] intvec: max of _[1], index of first variable with this max degree
1417EXAMPLE: example maxdegs; shows an example
1418"
1419{
1420  int i,n;
1421  intvec degs,maxdeg;
1422  list l;
1423
1424  n = nvars(basering);
1425
1426  for (i=1;i<=n;i++)
1427  {
1428   degs[i] = nrows(coeffs(f,var(i)))-1;
1429   if (degs[i] > maxdeg)
1430   {
1431    maxdeg[1] = degs[i];
1432    maxdeg[2] = i;
1433   }
1434  }
1435  return(list(degs,maxdeg));
1436}
1437example
1438{ "EXAMPLE:"; echo =2;
1439   ring r = 0,(x,y,z),lp;
1440   poly f = 3xy4 + 2xy2 + x5y3 + xz6 + y6;
1441   maxdegs(f);
1442}
1443///////////////////////////////////////////////////////////////////////////////
1444
1445proc chebyshev(int n,list #)
1446"USAGE:   chebyshev(n); n int, n >= 0
1447         chebyshev(n,c); n int, n >= 0, c number, c!=0
1448RETURN:  poly, the [monic] nth Chebyshev polynomial of the first kind. @*
1449         The polynomials are defined in the first variable, say x, of the
1450         basering.
1451NOTE:   @texinfo
1452 The (generalized) Chebyshev polynomials of the first kind  can be
1453         defined by the recursion:
1454@tex
1455$C_0 = c,\ C_1 = x,\ C_n = 2/c\cdot x\cdot C_{n-1}-C_{n-2},\ n \geq 2,c\neq 0$.
1456@end tex
1457@end texinfo
1458        These polynomials commute by composition:
1459        @math{C_m \circ C_n = C_n\circ C_m}. @*
1460        For c=1, we obtain the standard (non monic) Chebyshev polynomials
1461        @math{T_n} which satisfy @math{T_n(x)  = \cos(n \cdot \arccos(x))}. @*
1462        For c=2 (default), we obtain the monic Chebyshev polynomials @math{P_n}
1463        which satisfy the relation @math{P_n(x+ 1/x) = x^n+ 1/x^n}. @*
1464        By default the monic Chebyshev polynomials are returned:
1465        @math{P_n =}@code{chebyshev(n)} and @math{T_n=}@code{chebyshev(n,1)}.@*
1466        It holds @math{P_n(x) = 2\cdot T_n(x/2)} and more generally
1467        @math{C_n(c\cdot x) = c\cdot T_n(x)} @*
1468        That is @code{subst(chebyshev(n,c),var(1),c*var(1))= c*chebyshev(n,1)}.
1469
1470        If @code{char(basering) = 2}, then
1471        @math{C_0 = 1, C_1 = x, C_2 = 1, C_3 = x}, and so on.
1472EXAMPLE: example chebyshev; shows some examples
1473"
1474{
1475 number startv = 2;
1476
1477 if (size(#)){ startv = #[1]; }
1478 if (startv == 0) { startv = 1; }
1479
1480 poly f0,f1 = startv,var(1);
1481 poly fneu,falt = f1,f0;
1482 poly fh;
1483
1484 if (n<=0) {return(f0);}
1485 if (n==1) {return(f1);}
1486
1487 for(int i=2;i<=n;i++)
1488 {
1489   fh = 2/startv*var(1)*fneu - falt;
1490  // fh = 2*var(1)*fneu - falt;
1491  falt = fneu;
1492  fneu = fh;
1493 }
1494 return(fh);
1495}
1496example
1497{ "EXAMPLE:"; echo = 2;
1498   ring r = 0,x,lp;
1499
1500   // The monic Chebyshev polynomials
1501   chebyshev(0);
1502   chebyshev(1);
1503   chebyshev(2);
1504   chebyshev(3);
1505
1506   // These polynomials commute
1507   compose(chebyshev(2),chebyshev(6)) ==
1508   compose(chebyshev(6),chebyshev(2));
1509
1510   // The standard Chebyshev polynomials
1511   chebyshev(0,1);
1512   chebyshev(1,1);
1513   chebyshev(2,1);
1514   chebyshev(3,1);
1515   // -----------------------------------------------------------------------
1516   // The relation for the various Chebyshev polynomials
1517   5*chebyshev(3,1)==subst(chebyshev(3,5),x,5x);
1518   // -----------------------------------------------------------------------
1519   // char 2 case
1520   ring r2 = 2,x,dp;
1521   chebyshev(2);
1522   chebyshev(3);
1523}
1524///////////////////////////////////////////////////////////////////////////////
1525
1526/*
1527
1528// Examples for decomp.lib
1529
1530ring r02 = 0,(x,y),dp;
1531
1532decompose(compose(x6,chebyshev(4),x2+y3+x5y7),1);
1533
1534int  MINS = 0;
1535decompose((xy+1)^7);
1536//_[1]=x7
1537//_[2]=xy+1
1538
1539decompose((x2y3+1)^7);
1540//_[1]=y7
1541//_[2]=x2y3+1
1542
1543MINS = 1;
1544ring r01 = 0,x,dp;
1545decompose((x+1)^7);
1546//x7+7x6+21x5+35x4+35x3+21x2+7x+1
1547
1548decompunivmonic((x+1)^7,7);
1549//_[1]=x7
1550//_[2]=x+1
1551
1552int MINS =1;
1553 decompunivmonic((x+1)^7,7);
1554//_[1]=x7+7x6+21x5+35x4+35x3+21x2+7x+1
1555//_[2]=x
1556
1557 // --  Example -------------
1558
1559//  Comparision Kozen-Landau vs. von zur Gathen
1560
1561 ring r02 = 0,(x,y),dp;
1562
1563  // printlevel = 5;
1564
1565 decompopts("reset");
1566
1567 poly F = compose(x6,chebyshev(4)+3,8x2+y3+7x5y7+2);
1568 deg(F);
1569
1570 timer = 1;decompose(F,1);timer;
1571
1572 int MINS = 1;
1573 timer = 1;decompose(F,1);timer;
1574 int IMPROVE  =0;
1575 timer = 1;decompose(F,1);timer;
1576
1577 decompopts("reset");
1578 int DECMETH = 0;  // von zur Gathen
1579
1580 timer = 1;decompose(F,1);timer;
1581
1582decompopts("reset");
1583
1584 // -- Example -------------
1585
1586ring rZ10 = (integer,10),x,dp;
1587chebyshev(2);
1588//x2+8
1589chebyshev(3);
1590//x3+7x
1591
1592compose(chebyshev(2),chebyshev(3));
1593//x6+4x4+9x2+8
1594decompose(_);
1595int MINS =1;
1596decompose(compose(chebyshev(2),chebyshev(3)));
1597compose(_);
1598
1599decompopts("reset");
1600
1601// --  Example -------------
1602
1603ring rT =(0,y),x,dp;
1604compose(x2,x3+y,(y+1)*x2);
1605//(y6+6y5+15y4+20y3+15y2+6y+1)*x12+(2y4+6y3+6y2+2y)*x6+(y2)
1606
1607 decompose(_,1);
1608//_[1]=(y6+6y5+15y4+20y3+15y2+6y+1)*x2
1609//_[2]=x3+(y)/(y3+3y2+3y+1)
1610//_[3]=x2
1611
1612int MINS =1;
1613compose(x2,x3+y,(y+1)*x2);
1614//(y6+6y5+15y4+20y3+15y2+6y+1)*x12+(2y4+6y3+6y2+2y)*x6+(y2)
1615
1616decompose(_,1);
1617//_[1]=(y6+6y5+15y4+20y3+15y2+6y+1)*x2+(2y4+6y3+6y2+2y)*x+(y2)
1618//_[2]=x3
1619//_[3]=x2
1620
1621//ring rt =(0,t),x,dp;
1622//compose(x2+tx+5,x5-2tx3+x);
1623//x10+(-4t)*x8+(4t2+2)*x6+(t)*x5+(-4t)*x4+(-2t2)*x3+x2+(t)*x+5
1624
1625decompose(_);
1626//_[1]=x2+(-1/4t2+5)
1627//_[2]=x5+(-2t)*x3+x+(1/2t)
1628
1629int IMPROVE = 1;
1630compose(x2+tx+5,x5-2tx3+x);
1631//x10+(-4t)*x8+(4t2+2)*x6+(t)*x5+(-4t)*x4+(-2t2)*x3+x2+(t)*x+5
1632
1633decompose(_);
1634//_[1]=x2+(-1/4t2+5)
1635//_[2]=x5+(-2t)*x3+x+(1/2t)
1636
1637int IMPROVE = 0;
1638compose(x2+tx+5,x5-2tx3+x);
1639//x10+(-4t)*x8+(4t2+2)*x6+(t)*x5+(-4t)*x4+(-2t2)*x3+x2+(t)*x+5
1640decompose(_);
1641//_[1]=x2+(-1/4t2+5)
1642//_[2]=x5+(-2t)*x3+x+(1/2t)
1643
1644int MINS = 1;
1645compose(x2+tx+5,x5-2tx3+x);
1646//x10+(-4t)*x8+(4t2+2)*x6+(t)*x5+(-4t)*x4+(-2t2)*x3+x2+(t)*x+5
1647
1648decompose(_);
1649//_[1]=x2+(t)*x+5
1650//_[2]=x5+(-2t)*x3+x
1651
1652*/
1653///////////////////////////////////////////////////////////////////////////////
1654// --- End of decomp.lib --------------------------------------------------- //
1655///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.