source: git/Singular/LIB/chern.lib @ 20f2239

spielwiese
Last change on this file since 20f2239 was 20f2239, checked in by Hans Schoenemann <hannes@…>, 4 years ago
spelling p7
  • Property mode set to 100644
File size: 126.2 KB
Line 
1////////////////////////////////////////////////////////////////
2version = "version chern.lib 4.1.3.0 Apr_2020 ";  // $Id$
3category = "Chern classes";
4info="
5LIBRARY:  chern.lib    Symbolic Computations with Chern classes,
6                       Computation of Chern classes
7
8AUTHOR:  Oleksandr Iena,      o.g.yena@gmail.com
9
10OVERVIEW:
11  A toolbox for symbolic computations with Chern classes.
12  The Aluffi's algorithms for computation of characteristic classes of algebraic varieties
13  (Segre, Fulton, Chern-Schwartz-MacPherson classes) are implemented as well.
14
15
16REFERENCES:
17  [1] Aluffi, Paolo     Computing characteristic classes of projective schemes.
18                        Journal of Symbolic Computation, 35 (2003), 3-19.
19  [2] Iena, Oleksandr,  On symbolic computations with Chern classes:
20                        remarks on the library chern.lib for Singular,
21                        http://hdl.handle.net/10993/22395, 2015.
22  [3] Lascoux, Alain,   Classes de Chern d'un produit tensoriel.
23                        C. R. Acad. Sci., Paris, Ser. A 286, 385-387 (1978).
24  [4] Manivel, Laurent  Chern classes of tensor products, arXiv 1012.0014, 2010.
25
26PROCEDURES:
27  symm(l [,N]);             symmetric functions in the entries of l
28  symNsym(f, c);            symmetric and non-symmetric parts of a polynomial f
29  CompleteHomog(N, l);      complete homogeneous symmetric functions
30  segre(N, c);              Segre classes in terms of Chern classes
31  chern(N, s);              Chern classes in terms of Segre classes
32  chNum(N, c);              the non-zero Chern numbers in degree N in the entries of c
33  chNumbers(N, c);          the Chern numbers in degree N in the entries of c
34  sum_of_powers(k, l);      the sum of k-th powers of the entries of l
35  powSumSym(c [,N]);        the sums of powers [up to degree N] in terms
36                            of the elementary symmetric polynomials (entries of l)
37  chAll(c [,N]);            Chern character in terms of the Chern classes
38  chAllInv(c);              Chern classes in terms of the Chern character
39  chHE(c);                  the highest term of the Chern character
40  ChernRootsSum(a, b);      the Chern roots of a direct sum
41  chSum(c, C);              the Chern classes of a direct sum
42  ChernRootsDual(l);        the Chern roots of the dual vector bundle
43  chDual(c);                the Chern classes of the dual vector bundle
44  ChernRootsProd(l, L);     the Chern roots of a tensor product of vector bundles
45  chProd(r, c, R, C [,N]);  Chern classes of a tensor product of vector bundles
46  chProdE(c, C);            Chern classes of a tensor product of vector bundles
47  chProdL(r, c, R, C);      Chern classes of a tensor product of vector bundles
48  chProdLP(r, c, R, C);     total Chern class of a tensor product of vector bundles
49  chProdM(r, c, R, C);      Chern classes of a tensor product of vector bundles
50  chProdMP(r, c, R, C);     total Chern class of a tensor product of vector bundles
51  ChernRootsHom(l, L);      the Chern roots of a Hom vector bundle
52  chHom(r, c, R, C [,N]);   Chern classes of the Hom-vector bundle
53  ChernRootsSymm(n, l);     the Chern roots of the n-th symmetric power
54                            of a vector bundle with Chern roots from l
55  ChernRootsWedge(n, l);    the Chern roots of the n-th exterior power
56                            of a vector bundle with Chern roots from l
57  chSymm(k, r, c [,p]);     the rank and the Chern classes of the k-th symmetric power
58                            of a vector bundle of rank r with Chern classes c
59  chSymm2L(r, c);           the rank and the Chern classes of the second symmetric power
60                            of a vector bundle of rank r with Chern classes c
61  chSymm2LP(r, c);          the total Chern class of the second symmetric power
62                            of a vector bundle of rank r with Chern classes c
63  chWedge(k, r, c [,p]);    the rank and the Chern classes of the k-th exterior power
64                            of a vector bundle of rank r with Chern classes c
65  chWedge2L(r, c);          the rank and the Chern classes of the second exterior power
66                            of a vector bundle of rank r with Chern classes c
67  chWedge2LP(r, c);         the total Chern class of the second exterior power
68                            of a vector bundle of rank r with Chern classes c
69  todd(c [,n]);             the Todd class
70  toddE(c);                 the highest term of the Todd class
71  Bern(n);                  the second Bernoulli numbers
72  tdCf(n);                  the coefficients of the Todd class of a line bundle
73  tdTerms(n, f);            the terms of the Todd class of a line bundle
74                            corresponding to the Chern root t
75  tdFactor(n, t);           the Todd class of a line bundle corresponding
76                            to the Chern root t
77  cProj(n);                 the total Chern class of (the tangent bundle on)
78                            the projective space P_n
79  chProj(n);                the Chern character of (the tangent bundle on)
80                            the projective space P_n
81  tdProj(n);                the Todd class of (the tangent bundle on)
82                            the projective space P_n
83  eulerChProj(n, r, c);     Euler characteristic of a vector bundle on
84                            the projective space P_n
85                            via Hirzebruch-Riemann-Roch theorem
86  chNumbersProj(n);         the Chern numbers of the projective space P_n
87  classpoly(l, t);          polynomial in t with coefficients from l
88                            (without constant term)
89  chernPoly(l, t);          Chern polynomial (constant term 1)
90  chernCharPoly(r, l, t);   polynomial in t corresponding to the Chern character
91                            (constant term r)
92  toddPoly(td, t);          polynomial in t corresponding to the Todd class
93                            (constant term 1)
94  rHRR(N, ch, td);          the main ingredient of the right-hand side
95                            of the Hirzebruch-Riemann-Roch formula
96  SchurS(I, S);             the Schur polynomial corresponding to partition I
97                            in terms of the Segre classes S
98  SchurCh(I, C);            the Schur polynomial corresponding to partition I
99                            in terms of the Chern classes C
100  part(m, n);               partitions of integers not exceeding n
101                            into m non-negative summands
102  dualPart(I [,N]);         partition dual to I
103  PartC(I, m);              the complement of a partition with respect to m
104  partOver(n, J);           partitions over a given partition J with summands not exceeding n
105  partUnder(J);             partitions under a given partition J
106  SegreA(I);                Segre class of the projective subscheme defined by I
107  FultonA(I);               Fulton class of the projective subscheme defined by I
108  CSMA(I);                  Chern-Schwartz-MacPherson class of the
109                            projective subscheme defined by I
110  EulerAff(I);              Euler characteristic of the affine subvariety defined by I
111  EulerProj(I);             Euler characteristic of the projective subvariety defined by I
112";
113
114LIB "general.lib";
115LIB "lrcalc.lib"; // needed for chProdM(..) and chProdMP(..)
116//----------------------------------------------------------
117
118proc symm(list l, list #)
119"USAGE:   symm(l [,n]);  l a list of polynomials, n integer
120RETURN:   list of polynomials
121PURPOSE:  computes the list of elementary symmetric functions in the entries of l
122EXAMPLE:  example symm; shows an example
123NOTE:     makes sense only for a list of polynomials
124"
125{
126  int N=size(l);
127  int n=size(l);
128  if(size(#)!=0)
129  {
130    if( is_integer(#[1]) )
131    {
132      N = #[1];
133    }
134  }
135  if(n==0) // if the list is empty, return the empty list
136  {
137    return(list());
138  }
139  else
140  {
141    int i, j;
142    list rez=list(1, l[1]);
143    for(i=2; i<=n; i++)
144    {
145      if( i<=N )
146      {
147        rez=rez+list(0);
148      }
149      for(j = min(i, N); j>=1; j--)
150      {
151        rez[j+1] = rez[j+1] + rez[j]*l[i];
152      }
153    }
154    return(delete(rez, 1));
155  }
156}
157example
158{
159  "EXAMPLE:";echo =2;
160  // elementary symmetric functions in x, y, z:
161  ring r = 0, (x, y, z), dp;
162  list l=(x, y, z);
163  print(symm(l));
164
165  //now let us compute only the first two symmetric polynomials in a(1), ... , a(10)
166  ring q= 0,(a(1..10)), dp;
167  list l=a(1..10);
168  print(symm(l, 2));
169}
170//-----------------------------------------------------------------------
171
172proc symNsym(poly f, list c)
173"USAGE:   symNsym(f, c);  f polynomial; c list of polynomials
174RETURN:   list with 2 poly entries
175PURPOSE:  computes a symmetric and a non-symmetric part of f
176          in terms of the elementary symmetric functions from c
177          as well a non-symmetric remainder
178EXAMPLE:  example symNsym; shows an example
179NOTE:     constants are not considered symmetric
180"
181{
182  ideal V=variables(f); // variables f depends on
183  int nV=size(V); // their number
184  if(nV==0)
185  {
186    return(list(f, 0));
187  }
188  // now f is non-constant and does depend on some variables
189  c=append_by_zeroes(nV, c); // append c by zeroes if it is too short
190  def br@=basering; // remember the base ring
191  // add additional variables to the base ring
192  int ii;
193  list l2 = ringlist(basering)[2];
194  for (ii = 1; ii <= nV; ii++)
195  {
196   l2[size(l2)+1] = "c@("+string(ii)+")";
197  }
198    for (ii = 1; ii <= nV; ii++)
199  {
200   l2[size(l2)+1] = "A@("+string(ii)+")";
201  }
202  ring r@ = create_ring(ringlist(basering)[1], l2, "dp", "no_minpoly");
203  execute( "map F= br@,"+varstr(br@)+";" ); // define the corresponding inclusion of rings
204  ideal V=F(V);
205  poly f=F(f);
206  int i;
207  for(i=1; i<=nV; i++)
208  {
209    f=subst(f, V[i], A@(i) ); // rename the variables of f into A@(1..nV)
210  }
211  int N1=nvars(basering)-nV+1; // the number of variable A@(1)
212  poly rez1=0; // to be the expression in c@(i) of the symmetric part of f
213  poly rez2=0; // to be the remainder
214  poly mon; // monomial in c@(i)
215  poly monc; // the corresponding expression in A@(i)
216  list l=symm(list(A@(1..nV) )); // symmetric functions in A@(i)
217  intvec v=leadexp(f), 0; // the exponent of the leading monomial
218  while(v[N1]!=0)
219  {
220    mon=leadcoef(f); // leading coefficient of f
221    monc=mon;
222    for(i=1; v[N1+i-1]!=0 ;i++ )
223    {
224      mon = mon*c@(i)^( v[N1+i-1]-v[N1+i] );
225      monc = monc*l[i]^( v[N1+i-1]-v[N1+i] ); // has the same leading coefficient as f
226    }
227    rez1=rez1+mon; // add a monomial
228    f=f-monc; // subtract the monomial
229    v=leadexp(f), 0;
230  }
231  while( leadexp(f)!=0 )
232  {
233    rez2=rez2+lead(f);
234    f=f-lead(f);
235  }
236  rez1=rez1+f;
237  setring br@; // come back to the initial base ring
238  // define the specialization homomorphism
239  execute("map FF = r@,"+varstr(br@)+",c[1..nV], V[1..nV];");
240  return( list( FF(rez1), FF(rez2) ) );
241}
242example
243{
244  "EXAMPLE:";echo=2;
245  ring r=0, (x,y,z, c(1..3)), dp;
246  list l=c(1..3);
247  // The symmetric part of f = 3x2 + 3y2 + 3z2 + 7xyz + y
248  // in terms of the elemenatary symmetric functions c(1), c(2), c(3)
249  // and the remainder
250  poly f = 3x2 + 3y2 + 3z2 + 7xyz + y;
251  print( symNsym(f, l) );
252  // Take a symmetrix polynomial in variables x and z
253  f=x2+xz+z2;
254  // Express it in terms of the elementary the symmetric functions
255  print( symNsym(f, l)[1]);
256}
257//-------------------------------------------------------------------------------------------
258
259proc CompleteHomog(int N, list c)
260"USAGE:   CompleteHomog(N, c);  N integer, c list of polynomials
261RETURN:   list of polynomials
262PURPOSE:  computes the list of the complete homogeneous symmetric polynomials
263          in terms of the elementary symmetric polynomials (entries of c)
264EXAMPLE:  example CompleteHomog; shows an example
265NOTE:
266"
267{
268  c=append_by_zeroes(N, c);
269  if(N<0) // if N is negative, return the empty list
270  {
271    return(list());
272  }
273  list rez=list(1); // the result will be computed here
274  int i, j;
275  int sign;
276  poly f;
277  for(i=1; i<=N; i++) // use the recursive formula
278  {
279    f=0;
280    sign=1;
281    for(j=1;j<=i; j++) // compute the next complete homogeneous symmetric polynomial
282    {
283      f=f+sign*c[j]*rez[i-j+1];
284      sign=-sign;
285    }
286    rez=rez+( list(f)  );
287  }
288  return(rez);
289}
290example
291{
292  "EXAMPLE:";echo =2;
293  ring r = 0, (x(1..3)), dp;
294  list l=x(1..3);
295  //Complete homogeneous symmetric polynomials up to degree 3 in variables x(1), x(2), x(3)
296  print( CompleteHomog(3, l) );
297}
298//-----------------------------------------------------------------------
299
300proc segre(list c, list #)
301"USAGE:   segre(c[, N]); c list of polynomials, N integer
302RETURN:   list of polynomials
303PURPOSE:  computes the list of the Segre classes up to degree N
304          in terms of the Chern classes from c
305EXAMPLE:  example segre; shows an example
306NOTE:
307"
308{
309  int N;
310  if(size(#)>0)
311  {
312    if( is_integer(#[1]) )
313    {
314      N=#[1];
315    }
316  }
317  else
318  {
319    N=size(c);
320  }
321  c=append_by_zeroes(N, c);
322  if(N<0) // if N is negative, return the empty list
323  {
324    return(list());
325  }
326  list rez=list(1); // the result will be computed here
327  int i, j;
328  poly f;
329  for(i=1; i<=N; i++) // use the recursive formula
330  {
331    f=0;
332    for(j=1;j<=i; j++) // compute the next Segre class
333    {
334      f=f-c[j]*rez[i-j+1];
335    }
336    rez=rez+( list(f)  );
337  }
338  return(delete(rez,1));
339}
340example
341{
342  "EXAMPLE:";echo =2;
343  ring r = 0, (c(1..3)), dp;
344  list l=c(1..3);
345  //Segre classes up to degree 5 in Chern classes c(1), c(2), c(3)
346  print( segre(l, 5) );
347}
348//-----------------------------------------------------------------------
349
350proc chern(list s, list #)
351"USAGE:   chern(s);  s list of polynomials
352RETURN:   list of polynomials
353PURPOSE:  computes the list of the Chern classes up to degree N
354          in terms of the Segre classes from s
355EXAMPLE:  example chern; shows an example
356NOTE:
357"
358{
359  return( segre(s, #) );
360}
361example
362{
363  "EXAMPLE:"; echo =2;
364  ring r = 0, (s(1..3)), dp;
365  list l=s(1..3);
366  // Chern classes in Segre classes s(1), s(2), s(3)
367  print( chern(l) );
368  // This procedure is inverse to segre(...). Indeed:
369  print( segre(chern(l), 3) );
370}
371//-----------------------------------------------------------------------
372
373proc chNum(int N, list c)
374"USAGE:   chNum(N, c);  N integer, c list
375RETURN:   list
376PURPOSE:  computes the Chern numbers of a vector bundle with Chern classes c
377          on a complex manifold (variety) of dimension N,
378          the zeroes corresponding to the higher zero Chern classes are ignored
379EXAMPLE:  example chNumbers; shows an example
380NOTE:     computes basically the partitions of N
381          in summands not greater than the length of c
382"
383{
384  int n=size(c);
385  if(N<0)
386  {
387    print("");
388    return( list() );
389  }
390  if( (n==0) || (N==0)  )
391  {
392    return(list(1));
393  }
394  if(n==1) // if there is only one entry in the list
395  {
396    return(list(c[1]^N));
397  }
398  int i;
399  int j;
400  poly f; // the powers of the last variable will be stored here
401  list l=delete(c, n); // delete the last variable
402  list L;
403  list rez=chNum(N, l); // monomials not involving the last variable
404  for(i=1;i<=(N div n); i++) // add the monomials involving the last variable
405  {
406    f=c[n]^i; // the power of the last variable
407    // monomials without the last variable that,
408    // multiplied by the i-th power of the last variable,
409    // give a monomial of the required type
410    L=chNum(N-n*i, l);
411    for(j=1; j<=size(L) ;j++) // multiply every such monomial
412    {
413      L[j]=L[j]*f; // by the i-th power of the last variable
414    }
415    rez=rez+L; // add the monomials involving the i-th power of the last variable
416  }
417  return(rez);
418}
419example
420{
421  "EXAMPLE:";echo=2;
422  ring r = 0, (c(1..2)), dp;
423  list l=c(1..2);
424  // Let c(1) be a variable of degree 1, let c(2) be a variable of degree 2.
425  // The monomials in c(1) and c(2) of weighted degree 5 are:
426  print( chNum( 5, l ) );
427
428  // Compare the result to the output of chNumbers(...):
429  print( chNumbers(5, l) );
430}
431//----------------------------------------------------------------------------------------
432
433proc chNumbers(int r, list c)
434"USAGE:   chNumbers(r, c);  r integer,  c list
435RETURN:   list
436PURPOSE:  computes the Chern numbers of a vector bundle with Chern classes c
437          on a complex manifold (variety) of dimension r
438EXAMPLE:  example chNumbers; shows an example
439NOTE:     computes basically the partitions of r
440"
441{
442  if(r<0)
443  {
444    print("The dimension of a manifold must be a non-negative integer!");
445    return(list()); // return the empty list in this case
446  }
447  if(r==0)
448  {
449    return(list(1));
450  }
451  //-----------------
452  // from now on r>0
453  //----------------
454  int n=size(c);
455  c=append_by_zeroes(r, c);
456  c=c[1..r]; // throw away redundant data
457  return(chNum(r, c));
458}
459example
460{
461  "EXAMPLE:";echo=2;
462  ring r = 0, (c(1..3)), dp;
463  list l=c(1..3);
464  // The Chern numbers of a vector bundle with Chern classes c(1), c(2), c(3)
465  // on a 3-fold:
466  print( chNumbers( 3, l ) );
467
468  // If the highest Chern class is zero, the Chern numbers are:
469  l=c(1..2);
470  print( chNumbers( 3, l ) );
471
472  // Compare this to the output of chNum(...):
473  print( chNum( 3, l ) );
474}
475//---------------------------------------------------------------------------------------
476
477proc sum_of_powers(int k, list l)
478"USAGE:   sum_of_powers(k, l);  k non-negative integer, l list of polynomials
479RETURN:   polynomial
480PURPOSE:  computes the sum of k-th powers of the entries of l
481EXAMPLE:  example sum_of_powers; shows an example
482NOTE:     returns 0 if k is negative
483"
484{
485  if(k<0) // return 0 if k is negative
486  {
487    print("The exponent must be non-negative; 0 has been returned");
488    return(0);
489  }
490  int i;
491  int n=size(l);
492  poly rez; // the result will be computed here
493  for(i=1;i<=n;i++) // compute the sum of powers
494  {
495    rez=rez+l[i]^k;
496  }
497  return(rez);
498}
499example
500{
501  "EXAMPLE:";echo =2;
502  ring r = 0, (x, y, z), dp;
503  list l=x, y, z;
504  //sum of 7-th powers of x, y, z
505  print( sum_of_powers(7, l) );
506}
507//-----------------------------------------------------------------------
508
509proc powSumSym(list c, list #)
510"USAGE:   powSumSym(l [,N]);  l a list of polynomials, N integer
511RETURN:   list of polynomials
512PURPOSE:  computes the expressions for the sums of powers [up to degree N]
513          in terms of the elementary symmetric polynomials (entries of l),
514EXAMPLE:  example powSumSym; shows an example
515NOTE:     returns the terms of the Chern character
516          multiplied by the corresponding factorials
517"
518{
519  int n;
520  if( size(#) == 0 ) // if there are no optional parameters
521  {
522    n = size(c); // set n to be the length of c
523  }
524  else // if there are optional parameters
525  {
526    if( is_integer(#[1])) // if the first optional parameter is an integer
527    {
528      n = max( #[1], 0 ); // if the parameter is negative, reset it to be zero
529      c = append_by_zeroes(n, c); // if n is greater than the length of c, append c by zeroes
530      if( n != 0 ) // if n is non-zero
531      {
532        c = c[1..n]; // take into account only the first n entries of c
533      }
534    }
535    else // if the optional parameter is not an integer, then
536    {
537      n=size(c); // ignore it and set n to be the length of c
538    }
539  }
540  list rez; // the result will be computed here
541  if(n==0) // return the empty list
542  {
543    return(rez)
544  }
545  else // otherwise proceed as follows:
546  {
547    // first compute the sums of powers of the Chern roots
548    // in terms of the Chern classes using the Newton's identities
549    int i, j, sign;
550    poly f;
551    // the first term of the Chern character coincides with the first Chern class,
552    // or equivalently with the sum of Chern roots
553    rez = rez + list(c[1]);
554    // compute the sums of powers of Chern roots recursively using the Newton's identities
555    for(j=2; j<=n; j++)
556    {
557      sign=1;
558      f=0;
559      for(i=1; i<j; i++)
560      {
561        f=f + c[i]*sign*rez[j-i];
562        sign = -sign;
563      }
564      f=f+sign*j*c[j];
565      rez=rez+list(f); // add the newly computed sum of powers of Chern roots to the list
566    }
567    return(rez); // return the result
568  }
569}
570example
571{
572  "EXAMPLE:";echo =2;
573  // the expressions of the first 3 sums of powers of 3 variables a(1), a(2), a(3)
574  // in terms of the elementary symmetric polynomials c(1), c(2), c(3):
575  ring r = 0, (c(1..3)), dp;
576  list l=(c(1..3));
577  print(powSumSym(l));
578  // The first 5 sums in the same situation
579  print(powSumSym(l, 5));
580}
581//---------------------------------------------------------------------------------
582
583proc chAll(list c, list #)
584"USAGE:   chAll(l [,N]);  l a list of polynomials, N integer
585RETURN:   list of polynomials
586PURPOSE:  computes the list of terms of positive degree [up to degree N] of
587          the Chern character, where the entries of l are considered as the Chern classes
588EXAMPLE:  example chAll; shows an example
589NOTE:     makes sense only for a list of polynomials
590"
591{
592    list rez; // to be the result
593    rez = powSumSym(c, #); // get the sums of powers of the Chern roots
594    int n = size(rez);
595    bigint fct=1;
596    int i;
597    for(i=1;i<=n;i++) // get the terms of the Chern character
598    {
599      fct=fct*i;
600      rez[i]=rez[i]/fct;
601    }
602    return(rez); // return the result
603}
604example
605{
606  "EXAMPLE:";echo =2;
607  // Chern character (terms of degree 1, 2, 3)
608  // corresponding to the Chern classes c(1), c(2), c(3):
609  ring r = 0, (c(1..3)), dp;
610  list l=(c(1..3));
611  print(chAll(l));
612  // terms up to degree 5 in the same situation
613  print(chAll(l, 5));
614}
615//---------------------------------------------------------------------------------
616
617proc chAllInv(list c)
618"USAGE:   chAllInv(l);  l a list of polynomials
619RETURN:   list of polynomials
620PURPOSE:  procedure inverse to chAll(), computes the list of Chern classes
621          from the list of terms of positive degree of the Chern character
622EXAMPLE:  example chAllInv; shows an example
623NOTE:     makes sense only for a list of polynomials
624"
625{
626  int n = size(c);
627  list rez;
628  if(n==0) // if the list of terms of Chern character is empty, return the empty list
629  {
630    return(rez);
631  }
632  else // otherwise compute the result using the Newton's identities
633  {
634    int j, i, sign;
635    poly f;
636    // transform the list of terms of the Chern character
637    // to the list of sums of powers of Chern roots
638    //bigint fct=1;
639    for(i=1; i<=n; i++)
640    {
641      //fct=fct*i;
642      //c[i]=fct*c[i];
643      c[i]=factorial(i)*c[i];
644    }
645    // the first Chern class coincides with the first degree term of the Chern character
646    rez=rez+list(c[1]);
647    // compute the higher Chern classes recursively using the Newton's identities
648    for(j=2;j<=n;j++)
649    {
650      sign=1;f=0;
651      for(i=1;i<j;i++)
652      {
653        f=f+ c[i]*sign*rez[j-i];
654        sign=-sign;
655      }
656      f=f+sign*c[j];
657      rez=rez+list(f/j);
658    }
659    return(rez); // return the result
660  }
661}
662example
663{
664  "EXAMPLE:";echo=2;
665  // first 3 Chern classes in terms of the first 3 terms
666  // of the Chern character Chern  ch(1), ch(2), ch(3):
667  ring r = 0, (ch(1..3)), dp;
668  list l=(ch(1..3));
669  print(chAllInv(l));
670  // let's see that chAllInv() is inverse to chAll()
671  print( chAll( chAllInv(l) ) );
672}
673//---------------------------------------------------------------------------------
674
675proc chHE(list c)
676"USAGE:   chHE(c);  c list of polynomials
677RETURN:   polynomial
678PURPOSE:  computes the highest relevant term of the Chern character
679EXAMPLE:  example chHE; shows an example
680NOTE:     uses the elimination and is extremely inefficient,
681          is included just for comparison with chAll(c)
682"
683{
684  int i;
685  // insure that the entries of c are polynomials
686  // in order to be able to apply maps
687  for(i=1;i<=size(c);i++)
688  {
689    c[i]=poly(c[i]);
690  }
691  int n=size(c);
692  if(n==0) // in this case return the empty list
693  {
694    return(list());
695  }
696  else // otherwise proceed as follows
697  {
698    def br@=basering; // remember the base ring
699    // add additional variables c@, a@(1..n) to the base ring
700    list l3 = "c@";
701    l3 = l3+ ringlist(basering)[2];
702    for (int ii = 1; ii <= n; ii++)
703    {
704     l3[size(l3)+1] = "a@("+string(ii)+")";
705    }
706    ring r@ = create_ring(ringlist(basering)[1], l3, "lp", "no_minpoly");
707    execute( "map F= br@,"+varstr(br@)+";" ); // define the corresponding inclusion of rings
708    list c=F(c); // embed c in the bigger ring
709    poly rez;
710    list A=a@(1..n);
711    list sym=symm(A);
712    ideal I;
713    poly E=1; // to be the product of variables which should be eliminated
714    for(i=1;i<=n;i++)
715    {
716      E=E*a@(i); // compute the product of the variables that must be eliminated
717      I=I, c[i]-sym[i];
718    }
719    I=I, c@-sum_of_powers(n, A);
720    I=simplify(elim(I, E),1);
721    rez = -subst(I[1], c@, 0);
722    setring br@; // come back to the initial base ring
723    execute( "map FF= r@,0,"+varstr(br@)+";" ); // define the specialization homomorphism
724    poly rez=FF(rez); // bring the result to the base ring
725    return( (1/factorial(n))*rez);
726  }
727}
728example
729{
730  "EXAMPLE:";echo =2;
731  ring r = 0, (c(1..3)), dp;
732  list l=c(1..3);
733  //the third degree term of the Chern character
734  print( chHE(l) );
735}
736//----------------------------------------------------------------------
737
738proc ChernRootsSum(list a, list b)
739"USAGE:   ChernRootsSum(a, b); a, b lists of polynomials
740RETURN:   list of polynomials
741PURPOSE:  computes the Chern roots of the direct (Whitney) sum
742          of a vector bundle with Chern roots a and a vector bundle with Chern roots b
743EXAMPLE:  example ChernRootsSum; shows an example
744NOTE:
745"
746{
747  return(a+b);
748}
749example
750{
751  "EXAMPLE:";echo =2;
752  ring r = 0, (a(1..3), b(1..2)), dp;
753  // assume a(1), a(2), a(3) are the Chern roots of a vector bundle E
754  // assume b(1), b(2) are the Chern roots of a vector bundle F
755  list l=a(1..3);
756  list L=b(1..2);
757  // the Chern roots of their direct sum is
758  print( ChernRootsSum(l, L) );
759}
760//----------------------------------------------------------------------
761
762proc chSum(list c, list C)
763"USAGE:   chSum(c, C);  c, C lists of polynomials
764RETURN:   list of polynomials
765PURPOSE:  computes the Chern classes of a direct sum of two vector bundles
766EXAMPLE:  example chSum; shows an example
767NOTE:
768"
769{
770  int N=size(c)+size(C);
771  c=append_by_zeroes(N, c); // append by zeroes if necessary
772  C=append_by_zeroes(N, C); // append by zeroes if necessary
773  list rez; // to be the result
774  int i;
775  int j;
776  poly f;
777  for(i=1;i<=N;i++)
778  {
779    f=c[i]+C[i];
780    for(j=1;j<i;j++)
781    {
782      f=f+c[j]*C[i-j];
783    }
784    rez=rez+list(f);
785  }
786  return(rez);
787}
788example
789{
790  "EXAMPLE:";echo =2;
791  ring r = 0, (c(1..3), C(1..2)), dp;
792  // Let E be a vector bundle with Chern classes c(1), c(2), c(3).
793  // Let F be a vector bundle with Chern classes C(1), C(2).
794  list l=c(1..3);
795  list L=C(1..2);
796  // Then the Chern classes of their direct sum are
797  print( chSum(l, L) );
798}
799//----------------------------------------------------------------------
800
801proc ChernRootsDual(list l)
802"USAGE:   ChernRootsDual(l); l a list of polynomials
803RETURN:   list of polynomials
804PURPOSE:  computes the Chern roots of the dual vector bundle
805          of a vector bundle with Chern roots from l
806EXAMPLE:  example ChernRootsDual; shows an example
807NOTE:
808"
809{
810  int n=size(l);
811  int i;
812  for(i=1;i<=n;i++) // change the sign of the entries of a
813  {
814    l[i]=-l[i];
815  }
816  return(l);
817}
818example
819{
820  "EXAMPLE:";echo =2;
821  ring r = 0, (a(1..3)), dp;
822  // assume a(1), a(2), a(3) are the Chern roots of a vector bundle
823  list l=a(1..3);
824  // the Chern roots of the dual vector bundle
825  print( ChernRootsDual(l) );
826}
827//----------------------------------------------------------------------
828
829proc chDual(list c)
830"USAGE:   chDual(c);   c list of polynomials
831RETURN:   list of polynomials
832PURPOSE:  computes the list of Chern classes of the dual vector bundle
833EXAMPLE:  example chDual; shows an example
834NOTE:
835"
836{
837  int n=size(c);
838  int i;
839  for(i=1;i<=n;i=i+2)
840  {
841    c[i]=-c[i];
842  }
843  return(c);
844}
845example
846{
847  "EXAMPLE:"; echo=2;
848  // Chern classes of a vector bundle that is dual to a vector bundle
849  // with Chern classes c(1), c(2), c(3)
850  ring r=0, (c(1..3)), dp;
851  list l=c(1..3);
852  print(chDual(l));
853}
854//-----------------------------------------------------------------------------------------
855
856proc ChernRootsProd(list a, list b)
857"USAGE:   ChernRootsProd(a, b); a, b lists of polynomials
858RETURN:   list of polynomials
859PURPOSE:  computes the Chern roots of the tensor product of a vector bundle with Chern roots a
860          and a vector bundles with Chern roots b
861EXAMPLE:  example ChernRootsProd; shows an example
862NOTE:
863"
864{
865  int na=size(a);
866  int nb=size(b);
867  int i;
868  int j;
869  list rez; // the result will be computed here
870  for(i=1;i<=na;i++) // compute the result
871  {
872    for(j=1;j<=nb;j++)
873    {
874      rez=rez+list(a[i]+b[j]);
875    }
876  }
877  return(rez);
878}
879example
880{
881  "EXAMPLE:"; echo=2;
882  ring r=0, (a(1..2), b(1..3)), dp;
883  list l=a(1..2);
884  list L=b(1..3);
885  // Chern roots of the tensor product of a vector bundle with Chern roots a(1), a(2)
886  // and a vector bundle with Chern roots b(1), b(2), b(3)
887  print(ChernRootsProd(l, L));
888}
889//-----------------------------------------------------------------------------------------
890
891proc chProd(def r,  list c, def R, list C, list #)
892"USAGE:   chProd(r, c, R, C [, N]);  r, R polynomials (integers);
893          c, C lists of polynomials, N integer
894RETURN:   list of polynomials
895PURPOSE:  computes the list of Chern classes of the product of two vector bundles
896          in terms of their ranks and Chern clases [up to degree N]
897EXAMPLE:  example chProd; shows an example
898NOTE:
899"
900{
901  // check the input data
902  if( is_integer(r) ) // if r is an integer
903  {
904    if(r<=0) // if r is negative or zero return the empty list
905    {
906      return( list() );
907    }
908    //----------------------------
909    //now r is a positive integer
910    //----------------------------
911    c=append_by_zeroes(r, c); // append c by zeroes if r is greater than the length of c
912    c=c[1..r]; // make c shorter (of length r) if r is smaller than the length of c
913  }
914  if( is_integer(R) ) // if R is an integer
915  {
916    if(R<=0) // if R is negative or zero return the empty list
917    {
918      return( list() );
919    }
920    //----------------------------
921    //now R is a positive integer
922    //----------------------------
923    C=append_by_zeroes(R, C); // append C by zeroes if R is greater than the length of C
924    C=C[1..R]; // make C shorter (of length R) if R is smaller than the length of C
925  }
926  //----------------------------------------------------------
927  //now r > 0 if it is an integer; R > 0 if it is an integer
928  //----------------------------------------------------------
929  int n;
930  if( is_integer(r) && is_integer(R) ) // if both r and R are integers
931  {
932    n=r*R; // set n to be the rank of the product bundle
933  }
934  else // otherwise define the rank of the product vector bundle by
935  {
936    n=size(c)*size(C); // looking at the lengths of c and C
937  }
938  if( size(#) != 0 ) // if there is an optional parameter
939  {
940    if( is_integer( #[1] )  ) // if this parameter is an integer
941    {
942      if( #[1]<=0 ) // if it is negative or zero, return the empty list
943      {
944        return( list() );
945      }
946      // now #[1] is positive
947      // the product bundle can only have non-zero Chern classes up to degree n
948      // so ignore the optional parameter if it is greater than n
949      n = min(#[1], n);
950    }
951  }
952  if(n==0) // if n is zero, return the empty list
953  {
954    return( list() );
955  }
956  //-----------------------------------------------------------
957  //now n is positive, we can perform the relevant computations
958  //-----------------------------------------------------------
959  int i, j;
960  c=append_by_zeroes(n, c); // append c by zeroes up to degree n
961  C=append_by_zeroes(n, C); // append C by zeroes up to degree n
962  c=c[1..n]; // throw away the redundant data if needed
963  C=C[1..n]; // throw away the redundant data if needed
964  // build the list of all terms of the Chern characters: for rank r, and Chern classes c
965  list ch = list(r) + chAll(c);
966  list CH = list(R) + chAll(C); // do the same for rank R and Chern classes C
967  poly f;
968  list chP;
969  // compute the list of the non-zero degree terms of the Chern character
970  // of the tensor product of two vector bundles
971  for(i=1;i<=n;i++) // using the multiplicativity of the Chern character
972  {
973    f=0;
974    for(j=0;j<=i;j++)
975    {
976      f=f+ch[j+1]*CH[i-j+1];
977    }
978    chP=chP+list(f);
979  }
980  return( chAllInv(chP) ); // return the corresponding Chern classes
981}
982example
983{
984  "EXAMPLE:"; echo =2;
985  ring H = 0, ( r, R, c(1..3), C(1..2) ), dp;
986  list l=c(1..3);
987  list L=C(1..2);
988  // the Chern classes of the tensor product of a vector bundle E of rank 3
989  // with Chern classes c(1), c(2), c(3)
990  // and a vector bundle F of rank 2 with Chern classes C(1) and C(2):
991  print( chProd(3, l, 2, L) );
992  // the first two Chern classes of the tensor product
993  // of a vector bundle E of rank r with Chern classes c(1) and c(2)
994  // and a vector bundle G of rank R with Chern classes C(1) and C(2)
995  // this gives the Chern classes of a tensor product on a complex surface
996  l=c(1..2);
997  L=C(1..2);
998  print( chProd(r, l, R, L, 2 ) );
999}
1000//---------------------------------------------------------------------------------
1001
1002proc chProdE(list c, list C)
1003"USAGE:   chProdE(c, C);   c, C lists of polynomials
1004RETURN:   list of polynomials
1005PURPOSE:  computes the list of Chern classes of the product
1006          of two vector bundles in terms of their Chern clases
1007EXAMPLE:  example chProdE; shows an example
1008NOTE:     makes sense only for (lists of) polynomials;
1009          uses elimination, hence very inefficient;
1010          included only for comparison with chProd(...)
1011"
1012{
1013  int r=size(c);
1014  int R=size(C);
1015  // insure that the entries of c and C are polynomials
1016  // in order to be able to apply maps
1017  int i,j;
1018  for(i=1;i<=r;i++)
1019  {
1020    c[i]=poly(c[i]);
1021  }
1022  for(i=1;i<=R;i++)
1023  {
1024    C[i]=poly(C[i]);
1025  }
1026  if( (r==0) && (R==0) ) // if one of the ranks is 0,
1027  {
1028    return( list() ); // return the empty list (zero bundles have no Chern classes)
1029  }
1030  //------------------------------------
1031  //now both r and R are greater than 0
1032  //------------------------------------
1033  int n=r*R; // the rank of the product of two vector bundles
1034  def br@=basering; // remember the base ring
1035  // add additional variables a@(1..r), b@(1..R), x@ to the base ring
1036  int ii;
1037  list l4 = "x@";
1038  l4 = l4+ ringlist(basering)[2];
1039  for (ii = 1; ii <= r; ii++)
1040  {
1041   l4[size(l4)+1] = "a@("+string(ii)+")";
1042  }
1043  for (ii = 1; ii <= R; ii++)
1044  {
1045   l4[size(l4)+1] = "b@("+string(ii)+")";
1046  }
1047  ring r@ = create_ring(ringlist(basering)[1], l4, "lp", "no_minpoly");
1048  execute( "map F= br@,"+varstr(br@)+";" ); // define the corresponding inclusion of rings
1049  list c=F(c); // embed c in the bigger ring
1050  list C=F(C); // embed C in the bigger ring
1051  list A=a@(1..r); // list of Chern roots of the first vector bundle
1052  list syma = symm(A); // symmetric functions in the Chern roots of the first vector bundles
1053  list B=b@(1..R); // list of Chern roots of the second vector bundle
1054  list symb=symm(B); // symmetric functions in the Chern roots of the second vector bundles
1055  ideal I;
1056  // the product of variables (all Chern roots) which should be eliminated
1057  poly E=product(A)*product(B);
1058  for(i=1; i<=r; i++)
1059  {
1060    for(j=1; j<=R; j++)
1061    {
1062      I=I, c[i]-syma[i], C[j]-symb[j]; // add the relations
1063    }
1064  }
1065  // the Chern roots of the tensor product in terms of the Chern roots of the factors
1066  list crt=ChernRootsProd(A, B);
1067  list Cf=symm(crt); // Chern classes of the product in terms of the Chern roots of the factors
1068  list rez; // the result will be computed here
1069  ideal J;
1070  for(i=1;i<=n;i++)
1071  {
1072    J = I, x@-Cf[i]; // add the equation for the i-th Chern class to the ideal of relations
1073    J = simplify(elim(J, E), 1); // eliminate the Chern roots
1074    // get the expression for the i-th Chern class of the product
1075    // in terms of the Chern classes of the factors
1076    rez = rez + list( -subst(J[1], x@, 0) );
1077  }
1078  setring br@; // come back to the initial base ring
1079  execute( "map FF= r@, 0,"+varstr(br@)+";" ); // define the specialization homomorphism t@=0
1080  list rez=FF(rez); // bring the result to the base ring
1081  return(rez); // return the corresponding Chern classes
1082}
1083example
1084{
1085  "EXAMPLE:"; echo =2;
1086  ring H = 0, ( c(1..3), C(1..2) ), dp;
1087  list l=c(1..3);
1088  list L=C(1..2);
1089  // the Chern classes of the tensor product of a vector bundle E of rank 3
1090  // with Chern classes c(1), c(2), c(3)
1091  // and a vector bundle F of rank 2 with Chern classes C(1) and C(2):
1092  print( chProdE(l,  L) );
1093}
1094//------------------------------------------------------------------------------------
1095
1096proc chProdL(int r, list c, int R, list C)
1097"USAGE:   chProdL(r, c, R, C);  r, R integers; c, C lists of polynomials
1098RETURN:   list
1099PURPOSE:  computes the list of Chern classes of the product of two vector bundles
1100          in terms of their Chern clases
1101EXAMPLE:  example chProdL; shows an example
1102NOTE:     Implementation of the formula of Lascoux, the Schur polynomials are computed
1103          using the second Jacobi-Trudi formula (in terms of the Chern classes)
1104"
1105{
1106  // check the input data
1107  if(r<=0) // if r is negative or zero return the empty list
1108  {
1109    return( list() );
1110  }
1111  //----------------------------
1112  //now r is a positive integer
1113  //----------------------------
1114  c=append_by_zeroes(r, c); // append c by zeroes if r is greater than the length of c
1115  c=c[1..r]; // make c shorter (of length r) if r is smaller than the length of c
1116  if(R<=0) // if R is negative or zero return the empty list
1117  {
1118    return( list() );
1119  }
1120  //----------------------------
1121  //now R is a positive integer
1122  //----------------------------
1123  C=append_by_zeroes(R, C); // append C by zeroes if R is greater than the length of C
1124  C=C[1..R]; // make C shorter (of length R) if R is smaller than the length of C
1125  //----------------------------------------------------------
1126  // now r > 0 and R > 0
1127  //----------------------------------------------------------
1128  def br@=basering; // remember the base ring
1129  // add  additional variables to the base ring
1130  int ii;
1131  list l5 = ringlist(basering)[2]; 
1132  l5[size(l5)+1] = "t@";
1133  for (ii = 1; ii <= r; ii++)
1134  {
1135   l5[size(l5)+1] = "c@("+string(ii)+")";
1136  }
1137  for (ii = 1; ii <= R; ii++)
1138  {
1139   l5[size(l5)+1] = "C@("+string(ii)+")";
1140  }
1141  ring r@ = create_ring(ringlist(basering)[1], l5, "dp", "no_minpoly");
1142  execute( "map F= br@,"+varstr(br@)+";" ); // define the corresponding inclusion of rings
1143  list c, C;
1144  int i;
1145  for(i=1;i<=r;i++)
1146  {
1147    c[i]=c@(i)*t@^i;
1148  }
1149  for(i=1;i<=R;i++)
1150  {
1151    C[i]=C@(i)*t@^i;
1152  }
1153  poly f = chProdLP(r,c,R,C); // get the total Chern class using the Lascoux formula
1154  matrix CF = coeffs(f, t@); // get its coefficients in front of the powers of t@
1155  int N=r*R;
1156  list rez; // write them in a list
1157  for(i=1;i<=N;i++)
1158  {
1159    rez=rez+list(CF[i+1,1]);
1160  }
1161  setring br@; // come back to the initial base ring
1162  // define the specialization homomorphism
1163  execute("map FF = r@,"+varstr(br@)+",0, c[1..r], C[1..R];");
1164  return( FF( rez ) ); // bring the result to the initial ring
1165}
1166example
1167{
1168  "EXAMPLE:"; echo =2;
1169  // The Chern classes of the tensor product of a vector bundle of rank 3
1170  // with Chern classes c(1), c(2), c(3) and a vector bundle of rank 1 with
1171  // Chern class C(1)
1172  ring r = 0, ( c(1..3), C(1)), dp;
1173  list c=c(1..3);
1174  list C=C(1);
1175  print( chProdL(3,c,1,C) );
1176}
1177//---------------------------------------------------------------------------------------
1178
1179proc chProdLP(int r, list c, int R, list C)
1180"USAGE:   chProdLP(r, c, R, C);  r, R integers; c, C lists of polynomials
1181RETURN:   polynomial
1182PURPOSE:  computes the total Chern class of the product of two vector bundles
1183          in terms of their ranks and Chern clases
1184EXAMPLE:  example chProdLP; shows an example
1185NOTE:     Implementation of the formula of Lascoux, the Schur polynomials are computed
1186          using the second Jacobi-Trudi formula (in terms of the Chern classes)
1187"
1188{
1189  if(r<=0) // if r is negative or zero, return 1
1190  {
1191    return( 1 );
1192  }
1193  if(R<=0) // if R is negative or zero, return 1
1194  {
1195    return( 1 );
1196  }
1197  //-------------------------------------------
1198  // now r and R are positive
1199  //-------------------------------------------
1200  c=append_by_zeroes(r, c);
1201  C=append_by_zeroes(R, C);
1202  c=c[1..r];
1203  C=C[1..R];
1204  list P;
1205  P=part(r, R); // compute the partitions of numbers up to R into r summands
1206  int sz=size(P); // number of such partitions
1207  int szu;
1208  int i, j;
1209  list T;
1210  list PU;
1211  list TU;
1212  poly rez; // the result will be computed here
1213  poly ST;
1214  // implement the formula of Lascoux:
1215  for(i=1;i<=sz;i++) // run through all the partitions from P
1216  {
1217    T=P[i]; // the current partition
1218    ST= SchurS( PartC(T, R) , C ); // compute the corresponding Schur polynomial
1219    PU=partUnder(T); // compute the partitions under T
1220    szu=size(PU); // number of such partitions
1221    for(j=1;j<=szu;j++) // run through all the partitions lying under T
1222    {
1223      TU=PU[j]; // for each of them
1224      rez=rez+IJcoef(T, TU)* SchurCh(TU, c) *ST; // add the corresponding term to the result
1225     }
1226  }
1227  return(rez); // return the result
1228}
1229example
1230{
1231  "EXAMPLE:"; echo =2;
1232  // The total Chern class of the tensor product of a vector bundle of rank 3
1233  // with Chern classes c(1), c(2), c(3) and a vector bundle of rank 1 with
1234  // Chern class C(1)
1235  ring r = 0, ( c(1..3), C(1)), ws(1,2,3, 1);
1236  list c=c(1..3);
1237  list C=C(1);
1238  print( chProdLP(3,c,1,C) );
1239}
1240//---------------------------------------------------------------------------------------
1241
1242proc chProdM(int r, list c, int R, list C)
1243"USAGE:   chProdM(r, c, R, C);  r, R integers; c, C lists of polynomials
1244RETURN:   list
1245PURPOSE:  computes the list of Chern classes of the product of two vector bundles
1246          in terms of their Chern clases
1247EXAMPLE:  example chProdM; shows an example
1248NOTE:     Implementation of the formula of Manivel
1249"
1250{
1251  // check the input data
1252  if(r<=0) // if r is negative or zero return the empty list
1253  {
1254    return( list() );
1255  }
1256  //----------------------------
1257  //now r is a positive integer
1258  //----------------------------
1259  c=append_by_zeroes(r, c); // append c by zeroes if r is greater than the length of c
1260  c=c[1..r]; // make c shorter (of length r) if r is smaller than the length of c
1261  if(R<=0) // if R is negative or zero return the empty list
1262  {
1263    return( list() );
1264  }
1265  //----------------------------
1266  //now R is a positive integer
1267  //----------------------------
1268  C=append_by_zeroes(R, C); // append C by zeroes if R is greater than the length of C
1269  C=C[1..R]; // make C shorter (of length R) if R is smaller than the length of C
1270  //----------------------------------------------------------
1271  // now r > 0 and R > 0
1272  //----------------------------------------------------------
1273  def br@=basering; // remember the base ring
1274  // add  additional variables to the base ring
1275  int ii;
1276  list l6 = ringlist(basering)[2]; 
1277  l6[size(l6)+1] = "t@";
1278  for (ii = 1; ii <= r; ii++)
1279  {
1280   l6[size(l6)+1] = "c@("+string(ii)+")";
1281  }
1282  for (ii = 1; ii <= R; ii++)
1283  {
1284   l6[size(l6)+1] = "C@("+string(ii)+")";
1285  }
1286  ring r@ = create_ring(ringlist(basering)[1], l6, "dp", "no_minpoly");
1287  execute( "map F= br@,"+varstr(br@)+";" ); // define the corresponding inclusion of rings
1288  list c, C;
1289  int i;
1290  for(i=1;i<=r;i++)
1291  {
1292    c[i]=c@(i)*t@^i;
1293  }
1294  for(i=1;i<=R;i++)
1295  {
1296    C[i]=C@(i)*t@^i;
1297  }
1298  poly f = chProdMP(r,c,R,C); // get the total Chern class using the Manivel formula
1299  matrix CF = coeffs(f, t@); // get its coefficients in front of the powers of t@
1300  int N=r*R;
1301  list rez; // write them in a list
1302  for(i=1;i<=N;i++)
1303  {
1304    rez=rez+list(CF[i+1,1]);
1305  }
1306  setring br@; // come back to the initial base ring
1307  // define the specialization homomorphism
1308  execute("map FF = r@,"+varstr(br@)+",0, c[1..r], C[1..R];");
1309  return( FF( rez ) ); // bring the result to the initial ring
1310}
1311example
1312{
1313  "EXAMPLE:"; echo = 2;
1314  // The Chern classes of the tensor product of a vector bundle of rank 3
1315  // with Chern classes c(1), c(2), c(3) and a vector bundle of rank 1 with
1316  // Chern class C(1)
1317  ring r = 0, ( c(1..3), C(1)), dp;
1318  list c=c(1..3);
1319  list C=C(1);
1320  print( chProdM(3,c,1,C) );
1321}
1322//---------------------------------------------------------------------------------------
1323
1324proc chProdMP(int r, list c, int R, list C)
1325"USAGE:   chProdMP(r, c, R, C);  r, R integers; c, C lists of polynomials
1326RETURN:   polynomial
1327PURPOSE:  computes the total Chern class of the product of two vector bundles
1328          in terms of their ranks and Chern clases
1329EXAMPLE:  example chProdMP; shows an example
1330NOTE:     Implementation of the formula of Lascoux, the Schur polynomials are computed
1331          using the second Jacobi-Trudi formula (in terms of the Chern classes)
1332"
1333{
1334  if(r<=0) // if r is negative or zero, return 1
1335  {
1336    return( 1 );
1337  }
1338  if(R<=0) // if R is negative or zero, return 1
1339  {
1340    return( 1 );
1341  }
1342  //-------------------------------------------
1343  // now r and R are positive
1344  //-------------------------------------------
1345  c=append_by_zeroes(r, c);
1346  C=append_by_zeroes(R, C);
1347  c=c[1..r];
1348  C=C[1..R];
1349  list P;
1350  P=part(r, R); // compute the partitions of numbers up to R into r summands
1351  int sz=size(P); // number of such partitions
1352  int szu;
1353  int i, j;
1354  list T;
1355  list PU;
1356  list TU;
1357  poly rez; // the result will be computed here
1358  poly ST;
1359  // implement the formula of Manivel:
1360  for(i=1;i<=sz;i++) // run through all the partitions from P
1361  {
1362    T=P[i]; // the current partition
1363    ST= SchurS( PartC(T, R) , C ); // compute the corresponding Schur polynomial
1364    PU=partUnder(T); // compute the partitions under T
1365    szu=size(PU); // number of such partitions
1366    for(j=1;j<=szu;j++) // run through all the partitions lying under T
1367    {
1368      TU=PU[j]; // for each of them
1369      // add the corresponding term to the result
1370      rez=rez+Pcoef( TU, PartC(dualPart(T, R), r), r, R )* SchurCh(TU, c) *ST;
1371    }
1372  }
1373  return(rez); // return the result
1374}
1375example
1376{
1377  "EXAMPLE:"; echo =2;
1378  // The total Chern class of the tensor product of a vector bundle of rank 3
1379  // with Chern classes c(1), c(2), c(3) and a vector bundle of rank 1 with
1380  // Chern class C(1)
1381  ring r = 0, ( c(1..3), C(1)), ws(1,2,3, 1);
1382  list c=c(1..3);
1383  list C=C(1);
1384  print( chProdMP(3,c,1,C) );
1385}
1386//---------------------------------------------------------------------------------------
1387
1388proc ChernRootsHom(list a, list b)
1389"USAGE:   ChernRootsHom(a, b); a, b lists of polynomials
1390RETURN:   list of polynomials
1391PURPOSE:  for a vector bundle E with Chern roots a and a vector bundle F
1392          with Chern roots b, computes the Chern roots of Hom(E, F)
1393EXAMPLE:  example ChernRootsHom; shows an example
1394NOTE:
1395"
1396{
1397  int na=size(a);
1398  int nb=size(b);
1399  int i;
1400  int j;
1401  list rez; // the result will be computed here
1402  for(i=1;i<=na;i++) // compute the result
1403  {
1404    for(j=1;j<=nb;j++)
1405    {
1406      rez=rez+list(-a[i]+b[j]);
1407    }
1408  }
1409  return(rez);
1410}
1411example
1412{
1413  "EXAMPLE:"; echo=2;
1414  ring r=0, (a(1..2), b(1..3)), dp;
1415  list l=a(1..2);
1416  list L=b(1..3);
1417  // Let E be a vector bundle with Chern roots a(1). a(2),
1418  // let F be a vector bundle with CHern roots b(1), b(2), b(3).
1419  // Then the Chern roots of Hom(E, F) are
1420  print(ChernRootsHom(l, L));
1421}
1422//-----------------------------------------------------------------------------------------
1423
1424proc chHom(def r,  list c, def R, list C, list #)
1425"USAGE:   chHom(r, c, R, C [, N]);  r, R polynomials (integers);
1426          c, C lists of polynomials, N integer
1427RETURN:   list of polynomials
1428PURPOSE:  computes [up to degree N] the list of Chern classes of the vector bundle Hom(E, F)
1429          in terms of the ranks and the Chern classes of E and F
1430EXAMPLE:  example chHom; shows an example
1431NOTE:
1432"
1433{
1434  return( chProd(r, chDual(c), R, C, # ) );
1435}
1436example
1437{
1438  "EXAMPLE:"; echo=2;
1439  ring H = 0, ( r, R, c(1..3), C(1..2) ), dp;
1440  list l=c(1..3);
1441  list L=C(1..2);
1442  // the Chern classes of Hom(E, F) for a vector bundle E of rank 3
1443  // with Chern classes c(1), c(2), c(3)
1444  // and a vector bundle F of rank 2 with Chern classes C(1) and C(2):
1445  print( chHom(3, l, 2, L) );
1446  // the first two Chern classes of Hom(E, F) for a vector bundle E of rank r
1447  // with Chern classes c(1) and c(2)
1448  // and a vector bundle G of rank R with Chern classes C(1) and C(2)
1449  // this gives the Chern classes of a tensor product on a complex surface
1450  l=c(1..2);
1451  L=C(1..2);
1452  print( chHom(r, l, R, L, 2 ) );
1453}
1454//---------------------------------------------------------------------------------
1455
1456proc ChernRootsSymm(int n, list l)
1457"USAGE:   ChernRootsSymm(m, l); m integer, l a list of polynomials
1458RETURN:   list of polynomials
1459PURPOSE:  computes the Chern roots of m-th symmetric power
1460          of a vector bundle with Chern roots from l
1461EXAMPLE:  example ChernRootsSymm; shows an example
1462NOTE:
1463"
1464{
1465  if(n<0) // return the empty list if n is negative
1466  {
1467    return(list(0));
1468  }
1469  int r=size(l);
1470  def br@=basering; // remember the base ring
1471  ring r@=0, (a@(1..r)), dp;
1472  ideal mon = a@(1..r);
1473  mon=mon^n; // all monomials of degree n
1474  list rez;
1475  int i, j;
1476  int N = size(mon);
1477  intvec v;
1478  for(i=1; i<=N; i++) // collect in rez the exponents of the monomials of degree n
1479  {
1480    v = leadexp(mon[i]);
1481    rez = rez + list(v);
1482  }
1483  setring br@;
1484  poly f;
1485  list rez1;
1486  // run over all exponents and construct the corresponding sums of the Chern roots
1487  for(i=1; i<=N; i++)
1488  {
1489    f=0;
1490    for(j=1;j<=r;j++)
1491    {
1492      f=f+rez[i][j]*l[j];
1493    }
1494    rez1=rez1+list(f);
1495  }
1496  return(rez1);
1497}
1498example
1499{
1500  "EXAMPLE:";echo =2;
1501  ring r=0, (a(1..3)), dp;
1502  list l=a(1..3);
1503  // the Chern roots of the second symmetric power of a vector bundle
1504  // with Chern  roots a(1), a(2), a(3)
1505  print( ChernRootsSymm(2, l) );
1506}
1507//------------------------------------------------------------
1508
1509proc ChernRootsWedge( int m, list l)
1510"USAGE:   ChernRootsWedge(m, l); m integer, l a list of polynomials
1511RETURN:   list of polynomials
1512PURPOSE:  computes the Chern roots of m-th exterior power
1513          of a vector bundle with Chern roots from l
1514EXAMPLE:  example ChernRootsWedge; shows an example
1515NOTE:     makes sense only for list of polynomials
1516"
1517{
1518  int n=size(l);
1519  if((m>n)||(m<=0) ) // if m is bigger that n or non-positive
1520  {
1521    return( list(0) ); // return the list with one zero entry
1522  }
1523  else
1524  {
1525    if(m==n) // if m equals n, the only Chern root of the exterior power will be
1526    {
1527      return( list(sum(l)) ); // the sum of the initial Chern roots
1528    }
1529    else // otherwise proceed recursively
1530    {
1531      int i;
1532      list rez;
1533      list rez1;
1534      list l1 = delete(l, 1); // throw away the first element from the list
1535      poly f = l[1]; // remember the first entry of l
1536      // compute the Chern roots of the (m-1)-th exterior power of the smaller list
1537      rez1 = ChernRootsWedge(m-1, l1 );
1538      int s = size( rez1 );
1539      // add the first entry of the bigger list to every entry in the result,
1540      // this will give all Chern roots involving f
1541      for(i=1; i<=s; i++)
1542      {
1543        rez1[i] = f+rez1[i];
1544      }
1545      // return the union of those Chern roots with f and those without f
1546      rez = ChernRootsWedge(m, l1) + rez1;
1547      return( rez );
1548    }
1549  }
1550}
1551example
1552{
1553  "EXAMPLE:";echo =2;
1554  ring r=0, (a(1..3)), dp;
1555  list l=a(1..3);
1556  // the Chern roots of the second exterior power of a vector bundle
1557  // with Chern  roots a(1), a(2), a(3)
1558  print( ChernRootsWedge(2, l) );
1559}
1560//---------------------------------------------------------------------------------
1561
1562proc chSymm(int k, int r, list c, list #)
1563"USAGE:   chSymm(k, r, c[, pos]);  k, r integers, c list of polynomials, pos list of integers
1564RETURN:   list with entries: int N, list of polynomials l
1565PURPOSE:  computes the rank and the Chern classes of the symmetric power of a vector bundle
1566EXAMPLE:  example chSymm; shows an example
1567NOTE:     for the second symmetric power chSymm2L(...) could be faster
1568"
1569{
1570  // insure that the entries of c are polynomials
1571  // in order to be able to apply maps
1572  int i;
1573  for(i=1;i<=size(c);i++)
1574  {
1575    c[i]=poly(c[i]);
1576  }
1577  if(r<0) // if the rank is negative
1578  {
1579    print("The rank of a vector bundle can non be negative");
1580    return(list()); // return the empty list in this case
1581  }
1582  if(r==0) // if we deal with the zero bundle
1583  {
1584    return( list( 0, list() ) ); // return the data corresponding to the zero bundle
1585  }
1586  //-----------------------------------
1587  // from now on we are in the case r>0
1588  //-----------------------------------
1589  // if the length  n of the list of Chern classes is smaller
1590  // than the rank of the vector bundle,
1591  // the higher classes are assumed to be zero and the list is appended by zeroes up to length r
1592  c=append_by_zeroes(r, c);
1593  // if the length of the list of the Chern classes is greater than the rank
1594    c=c[1..r]; // throw away the redundant data
1595  //-----------------------------------
1596  // from now on the length of c is r>0
1597  //-----------------------------------
1598  if(k<0)
1599  {
1600    print("You are trying to compute a negative symmetric power of a vector bundle");
1601    return( list(0, list() ) ); // assume such a power to be just a zero bundle
1602  }
1603  if(k==0) // the zeroth symmetric power is the trivial line bundle
1604  {
1605    return( list(1, list(0)) );
1606  }
1607  if(k==1) // the first symmetric power is equal to the vector bundle itself
1608  {
1609    return(list(r, c));
1610  }
1611  //-----------------------------------
1612  // from now on we are in the case k>2
1613  //-----------------------------------
1614  list LM = integer_list(#);
1615  int M = LM[2]; // maximum among the optional parameters
1616  # = LM[1]; // take into account only the first integer optional parameters that are positive
1617  //-------------------------------
1618  // Perform the computations now
1619  //-------------------------------
1620  def br@=basering; // remember the base ring
1621  // add additional variables to the base ring
1622  int ii;
1623  list l7 = "x@";
1624  l7 = l7+ ringlist(basering)[2];
1625  for (ii = 1; ii <= r; ii++)
1626  {
1627   l7[size(l7)+1] = "a@("+string(ii)+")";
1628  }
1629  ring r@ = create_ring(ringlist(basering)[1], l7, "lp", "no_minpoly");
1630  execute( "map F= br@,"+varstr(br@)+";" ); // define the corresponding inclusion of rings
1631  list c=F(c); // embed c into the bigger ring
1632  list rez; // the Chern classes of the symmetric power are going to be written here
1633  poly E = product( list( a@(1..r ) )  ); // product of the Chern roots
1634  list ss=ChernRootsSymm(k, list( a@(1..r) ) ); // list of the Chern roots of the symmetric power
1635  int N=size(ss); // number of such roots, it equals the rank of the symmetric power
1636  // the entries in C will be the Chern classes of the symmetric power
1637  // expressed in terms of the Chern roots of the initial vector bundle
1638  list C;
1639  ideal I, J;
1640  // list of the Chern classes of the initial vector bundle expressed in its Chern roots
1641  list sym=symm(list(a@(1..r)));
1642  if(size(#)==0) // if there are no optional parameters, compute all Chern classes
1643  {
1644    // the entries here are the Chern classes of the symmetric power
1645    // expressed in terms of Chern roots of the initial vector bundle
1646    C=symm(ss);
1647    for(i=1;i<=N;i++) // eliminate the Chern roots
1648    {
1649      if(i<= r) // first add all relevant  formulas for the Chern classes in terms of Chern roots
1650      {
1651        I=I, c[i]-sym[i];
1652      }
1653      J = I,  x@-C[i];
1654      // Notice that elim(...) is from the library "elim.lib",
1655      // it is loaded as a result of loading "general.lib"
1656      J=simplify(elim(J, E), 1);
1657      // get the expression of the next Chern class
1658      // in terms of the Chern classes of the initial vector bundle
1659      rez=rez+list( -subst(  J[1], x@, 0)  );
1660    }
1661  }
1662  else // otherwise compute only the needed Chern classes
1663  {
1664    C=symm(ss, M); // only the needed Chern classes
1665    int j;
1666    i=1;
1667    // the maximal number of optional parameters to be considered does not exceed N,
1668    // i.e., the rank of the symmetric power
1669    int NN = min( size(#), N);
1670    for(j=1; j <= NN; j++) // process the optional parameters
1671    {
1672      // process the optional parameters only until they are not bigger than N;
1673      // notice they are positive anyway after integer_list(...)
1674      if( #[j]<=N  )
1675      {
1676        for( ; i<=#[j];i++)
1677        {
1678          if(i<=r)
1679          {
1680            // add the relevant formulas for the Chern classes in terms of the Chern roots
1681            I=I, c[i]-sym[i];
1682          }
1683        }
1684        J= I, x@-C[ #[j]];
1685        // Notice that elim(...) is from the library "elim.lib",
1686        // it is loaded as a result of loading "general.lib"
1687        J=simplify(elim(J, E), 1);
1688        // get the expression of the next Chern class
1689        // in terms of the Chern classes of the initial vector bundle
1690        rez=rez+list( -subst(  J[1], x@, 0)  );
1691      }
1692      else // get out from the loop
1693      {
1694        break;
1695      }
1696    }
1697  }
1698  // used because Singular seems not to be able to apply maps to empty lists (see below)
1699  if(size(rez)==0)
1700  {
1701    return(list(N, list()));
1702  }
1703  setring br@; // come back to the initial base ring
1704  // define the specialization homomorphism,
1705  // evaluate the formulas for the Chern classes on their given values
1706  execute( "map FF = r@,0,"+varstr(br@)+";" );
1707  list rez=FF( rez ); // bring the result back to the initial ring
1708  return( list( N, rez  ) ); // return the result together with the rank of the symmetric power
1709}
1710example
1711{
1712  "EXAMPLE:";echo =2;
1713  ring r=0, (c(1..5)), dp;
1714  list l=c(1..5);
1715  // the rank and the Chern classes of the second symmetric power of a vector bundle of rank 3
1716  print( chSymm(2, 3, l) );
1717  // the rank and the first 3 Chern classes
1718  // of the second symmetric power of a vector bundle of rank 5
1719  print( chSymm(2, 5, l, 1, 2, 3) );
1720}
1721//----------------------------------------------------------------------------------
1722
1723proc chSymm2L(int r, list c)
1724"USAGE:   chSymm2L(r, c); r integer, c list of polynomials
1725RETURN:   list of polynomials
1726PURPOSE:  computes the Chern classes of the second symmetric power of a vector bundle
1727EXAMPLE:  example chSymm2L; shows an example
1728NOTE:     Implementation of the formula of Lascoux, the Schur polynomials are computed
1729          using the second Jacobi-Trudi formula (in terms of the Chern classes)
1730"
1731{
1732  // insure that the entries of c are polynomials
1733  // in order to be able to apply maps
1734  int i;
1735  for(i=1;i<=size(c);i++)
1736  {
1737    c[i]=poly(c[i]);
1738  }
1739  if(r<0) // if the rank is negative
1740  {
1741    print("The rank of a vector bundle can non be negative");
1742    return(list()); // return the empty list in this case
1743  }
1744  if(r==0) // if we deal with the zero bundle
1745  {
1746    return( list( 0, list() ) ); // return the data corresponding to the zero bundle
1747  }
1748  //-----------------------------------
1749  // from now on we are in the case r>0
1750  //-----------------------------------
1751  c=append_by_zeroes(r, c);
1752  c=c[1..r];
1753  def br@=basering; // remember the base ring
1754  // add  additional variables to the base ring
1755  list l8 = ringlist(basering)[2]; 
1756  l8[size(l8)+1] = "t@";
1757  for (int ii = 1; ii <= r; ii++)
1758  {
1759   l8[size(l8)+1] = "c@("+string(ii)+")";
1760  }
1761  ring r@ = create_ring(ringlist(basering)[1], l8, "dp", "no_minpoly");
1762  execute( "map F= br@,"+varstr(br@)+";" ); // define the corresponding inclusion of rings
1763  list c;
1764  for(i=1;i<=r;i++)
1765  {
1766    c[i]=c@(i)*t@^i;
1767  }
1768  poly f = chSymm2LP(r,c); // get the total Chern class using the formula of Lascoux
1769  matrix CF = coeffs(f, t@);
1770  int N=r*(r+1) div 2;
1771  list rez; // write the coefficients in front of the powers of t@ into a list
1772  for(i=1;i<=N;i++)
1773  {
1774    rez=rez+list(CF[i+1,1]);
1775  }
1776  setring br@; // come back to the initial base ring
1777  execute("map FF = r@,"+varstr(br@)+",0, c[1..r];"); // define the specialization homomorphism
1778  return( list(N, FF( rez )) ); // bring the result to the initial ring
1779}
1780example
1781{
1782  "EXAMPLE:";echo =2;
1783  ring r=0, (c(1..2)), dp;
1784  list l=c(1..2);
1785  // the Chern classes of the second symmetric power of a vector bundle of rank 2
1786  print( chSymm2L(2, l));
1787}
1788//---------------------------------------------------------------------------------------
1789
1790proc chSymm2LP(int r, list c)
1791"USAGE:   chSymm2LP(r, c); r integer, c list of polynomials
1792RETURN:   poly
1793PURPOSE:  computes the total Chern class of the second symmetric power of a vector bundle
1794EXAMPLE:  example chSymm2LP; shows an example
1795NOTE:     Implementation of the formula of Lascoux, the Schur polynomials are computed
1796          using the second Jacobi-Trudi formula (in terms of the Chern classes)
1797"
1798{
1799  if(r<0) // if the rank is negative
1800  {
1801    print("The rank of a vector bundle can non be negative");
1802    return(1); // return 1 in this case
1803  }
1804  if(r==0) // if we deal with the zero bundle
1805  {
1806    return( 1 ); // return 1 in this case
1807  }
1808  //-------------------------------------------
1809  // from now on we are in the case r > 0
1810  //-------------------------------------------
1811  c=append_by_zeroes(r, c);
1812  c=c[1..r];
1813  list I; // the partition (1,2,...,r) will be stored here
1814  int i;
1815  for(i=1;i<=r;i++)
1816  {
1817    I=I+list(i);
1818  }
1819  list PU = partUnder(I); // compute the partitions under I
1820  int sz=size(PU); // get their number
1821  poly rez; // the result will be computed here
1822  list J;
1823  poly cf;
1824  int ex;
1825  // implement the formula of Lascoux
1826  for(i=1;i<=sz;i++)
1827  {
1828    J=PU[i];
1829    ex=sum(J)- r*(r-1) div 2;
1830    if(ex>=0)
1831    {
1832      cf=bigint(2)^ex*IJcoef(I, J);
1833    }
1834    else
1835    {
1836      cf=IJcoef(I, J)/bigint(2)^(-ex);
1837    }
1838    rez = rez + cf * SchurCh(J, c );
1839  }
1840  return(rez);
1841}
1842example
1843{
1844  "EXAMPLE:";echo =2;
1845  ring r=0, (c(1..2)), ws(1, 2);
1846  list l=c(1..2);
1847  // the total Chern class of the second symmetric power of a vector bundle of rank 2
1848  print( chSymm2LP(2, l));
1849}
1850//---------------------------------------------------------------------------------------
1851
1852proc chWedge(int k, int r, list c, list #)
1853"USAGE:   chWedge(k, r, c [,pos]);  k, r integers, c list of polynomials, pos list of integers
1854RETURN:   list with entries: int N, list of polynomials l
1855PURPOSE:  computes the rank and the Chern classes of the exterior power of a vector bundle
1856EXAMPLE:  example chWedge; shows an example
1857NOTE:     for the second exterior power chWedge2L(...) could be faster
1858"
1859{
1860  // insure that the entries of c are polynomials
1861  // in order to be able to apply maps
1862  int i;
1863  for(i=1;i<=size(c);i++)
1864  {
1865    c[i]=poly(c[i]);
1866  }
1867  if(r<0) // if the rank is negative
1868  {
1869    print("The rank of a vector bundle can non be negative");
1870    return(list()); // return the empty list in this case
1871  }
1872  if(r==0) // if we deal with the zero bundle
1873  {
1874    return( list( 0, list() ) ); // return the data corresponding to the zero bundle
1875  }
1876  //-------------------------------------------
1877  // from now on we are in the case r > 0
1878  //-------------------------------------------
1879  if(k<0)
1880  {
1881    print("You are trying to compute a negative exterior power of a vector bundle");
1882    return( list(0, list() ) ); // assume such a power to be just a zero bundle
1883  }
1884  if(k==0) // the zeroth exterior power is the trivial line bundle
1885  {
1886    return( list(1, list(0)) );
1887  }
1888  if(k==1) // the first exterior power is equal to the vector bundle itself
1889  {
1890    c=append_by_zeroes(r, c);
1891    c=c[1..r];
1892    return(list(r, c));
1893  }
1894  //---------------------------------------
1895  // from now on we are in the case k > 2
1896  //---------------------------------------
1897  // if the length of the list of Chern classes is smaller than the rank of the vector bundle,
1898  // the higher classes are assumed to be zero and the list is appended by zeroes up to length r
1899  c=append_by_zeroes(r, c);
1900  // if the length of the list of the Chern classes is greater than the rank
1901    c=c[1..r]; // throw away the redundant data
1902  //------------------------------------------
1903  // from now on the length of c is r > 0
1904  //------------------------------------------
1905  if( k>r ) // if k>r, the exterior power is zero
1906  {
1907    return( list( int(0), list()  ) );
1908  }
1909  //-----------------------------------------------
1910  // from now on we are in the case 0 < k <= r = n
1911  //-----------------------------------------------
1912  if(k==r)
1913  {
1914    return(list( int(1), list( c(1) ) ) );
1915  }
1916  //-----------------------------------------------
1917  // from now on we are in the case 0 < k < r = n
1918  //-----------------------------------------------
1919  list LM = integer_list(#);
1920  int M=LM[2]; // maximum among the optional parameters if there are any, zero otherwise
1921  # = LM[1]; // take into account only the first integer optional parameters that are positive
1922  //-----------------------------
1923  // Let us compute now
1924  //-----------------------------
1925  def br@=basering; // remember the base ring
1926  // add additional variables a@(1..r), x@ to the base ring
1927  list l9 = "x@";
1928  l9 = l9+ ringlist(basering)[2];
1929  for (int ii = 1; ii <= r; ii++)
1930  {
1931   l9[size(l9)+1] = "a@("+string(ii)+")";
1932  }
1933  ring r@ = create_ring(ringlist(basering)[1], l9, "lp", "no_minpoly");
1934  execute( "map F= br@,"+varstr(br@)+";" ); // define the corresponding inclusion of rings
1935  list c = F(c); // embed c into the bigger ring
1936  list rez; // the result should be computed here
1937  poly E = product( list( a@(1..r ) )  ); // product of the Chern roots to be eliminaned
1938  list ss=ChernRootsWedge(k,  list( a@(1..r) )); // list of the Chern roots of the exterior product
1939  int N=size(ss); // length of ss, equals the rank of the exterior product
1940  // list of the Chern classes of the initial vector bundle in terms of their Chern roots
1941  list sym=symm(list(a@(1..r)));
1942  // the entries here will be the Chern classes we need
1943  // expressed in  terms of the Chern roots of the initial vector bundle
1944  list C;
1945  ideal I, J;
1946  if( size(#) == 0 ) // if there are no optional parameters, compute all Chern classes
1947  {
1948    // the entries here are the Chern classes we need
1949    // expressed in  terms of the Chern roots of the initial vector bundle
1950    C=symm(ss);
1951    for(i=1;i<=N;i++) // eliminate the Chern roots
1952    {
1953      if(i<= r) // first add all relevant  formulas for the Chern classes in terms of Chern roots
1954      {
1955        I=I, c[i]-sym[i];
1956      }
1957      J = I,  x@-C[i];
1958      // Notice that elim(...) is from the library "elim.lib",
1959      // it is loaded as a result of loading "general.lib"
1960      J=simplify(elim(J, E), 1);
1961      // get the expression of the next Chern class
1962      // in terms of the Chern classes of the initial vector bundle
1963      rez=rez+list( -subst(  J[1], x@, 0)  );
1964    }
1965  }
1966  else // otherwise compute only the needed Chern classes
1967  {
1968    // the entries here are the Chern classes we need
1969    // expressed in  terms of the Chern roots of the initial vector bundle
1970    C=symm(ss, M);
1971    int j;
1972    i=1;
1973    // the maximal number of optional parameters to be considered
1974    // does not exceed N, the rank of the exterior power
1975    int NN = min( size(#), N);
1976    for(j=1; j <= NN; j++) // process the optional parameters
1977    {
1978      // process the optional parameters only until they are not bigger than N;
1979      // notice they are positive anyway after integer_list(...)
1980      if( #[j]<=N  )
1981      {
1982        for( ; i<=#[j]; i++)
1983        {
1984          if( i<=r )
1985          {
1986            // add the relevant formulas for the Chern classes in terms of the Chern roots
1987            I=I, c[i]-sym[i];
1988          }
1989        }
1990        J= I, x@-C[ #[j]];
1991        // Notice that elim(...) is from the library "elim.lib",
1992        // it is loaded as a result of loading "general.lib"
1993        J=simplify(elim(J, E), 1);
1994        // get the expression of the next Chern class
1995        // in terms of the Chern classes of the initial vector bundle
1996        rez=rez+list( -subst(  J[1], x@, 0)  );
1997      }
1998      else // get out from the loop
1999      {
2000        break;
2001      }
2002    }
2003  }
2004  // used because Singular seems not to be able to apply maps to empty lists (see below)
2005  if(size(rez)==0)
2006  {
2007    return(list(N, list()));
2008  }
2009  setring br@; // come back to the initial base ring
2010  // define the specialization homomorphism,
2011  // evaluate the formulas for the Chern classes on their given values
2012  execute( "map FF = r@,0,"+varstr(br@)+";" );
2013  list rez=FF( rez ); // bring the result back to the initial ring
2014  return( list( N, rez  ) ); //return the rank and the Chern classes of the exterior product
2015}
2016example
2017{
2018  "EXAMPLE:";echo =2;
2019  ring r=0, (c(1..5)), dp;
2020  list l=c(1..5);
2021  // the rank and the Chern classes of the second exterior power of a vector bundle of rank 3
2022  print( chWedge(2, 3, l) );
2023  // the rank and the first 3 Chern classes
2024  // of the fourth exterior power of a vector bundle of rank 5
2025  print( chWedge(4, 5, l, 1, 2, 3) );
2026}
2027//---------------------------------------------------------------------------------
2028
2029proc chWedge2L(int r, list c)
2030"USAGE:   chWedge2L(r, c );  r integer, c list of polynomials
2031RETURN:   list of polynomials
2032PURPOSE:  computes the Chern classes of the second exterior power of a vector bundle
2033EXAMPLE:  example chWedge2L; shows an example
2034NOTE:     Implementation of the formula of Lascoux, the Schur polynomials are computed
2035          using the second Jacobi-Trudi formula (in terms of the Chern classes)
2036"
2037{
2038  // insure that the entries of c are polynomials
2039  // in order to be able to apply maps
2040  int i;
2041  for(i=1;i<=size(c);i++)
2042  {
2043    c[i]=poly(c[i]);
2044  }
2045  if(r<0) // if the rank is negative
2046  {
2047    print("The rank of a vector bundle can non be negative");
2048    return(list()); // return the empty list in this case
2049  }
2050  if(r==0) // if we deal with the zero bundle
2051  {
2052    return( list( 0, list() ) ); // return the data corresponding to the zero bundle
2053  }
2054  //-------------------------------------------
2055  // from now on we are in the case r > 0
2056  //-------------------------------------------
2057  c=append_by_zeroes(r, c);
2058  c=c[1..r];
2059  def br@=basering; // remember the base ring
2060  // add  additional variables to the base ring
2061  list l10 = ringlist(basering)[2]; 
2062  l10[size(l10)+1] = "t@";
2063  for (int ii = 1; ii <= r; ii++)
2064  {
2065   l10[size(l10)+1] = "c@("+string(ii)+")";
2066  }
2067  ring r@ = create_ring(ringlist(basering)[1], l10, "dp", "no_minpoly");
2068  execute( "map F= br@,"+varstr(br@)+";" ); // define the corresponding inclusion of rings
2069  list c;
2070  for(i=1;i<=r;i++)
2071  {
2072    c[i]=c@(i)*t@^i;
2073  }
2074  poly f = chWedge2LP(r,c); // get the total Chern class using the formula of Lascoux
2075  matrix CF = coeffs(f, t@);
2076  int N=r*(r-1) div 2;
2077  list rez; // write its coefficients in front of the powers of t@ to a list
2078  for(i=1;i<=N;i++)
2079  {
2080    rez=rez+list(CF[i+1,1]);
2081  }
2082  setring br@; // come back to the initial base ring
2083  execute("map FF = r@,"+varstr(br@)+",0, c[1..r];"); // define the specialization homomorphism
2084  return( list(N, FF( rez )) ); // bring the result to the initial ring
2085}
2086example
2087{
2088  "EXAMPLE:";echo =2;
2089  ring r=0, (c(1..3)), dp;
2090  list l=c(1..3);
2091  // the Chern classes of the second exterior power of a vector bundle of rank 3
2092  print(chWedge2L(3, l));
2093}
2094//---------------------------------------------------------------------------------------
2095
2096proc chWedge2LP(int r, list c)
2097"USAGE:   chWedge2LP(r, c );  r integer, c list of polynomials
2098RETURN:   poly
2099PURPOSE:  computes the total Chern class of the second exterior power of a vector bundle
2100EXAMPLE:  example chWedge2LP; shows an example
2101NOTE:     Implementation of the formula of Lascoux, the Schur polynomials are computed
2102          using the second Jacobi-Trudi formula (in terms of the Chern classes)
2103"
2104{
2105  if(r<0) // if the rank is negative
2106  {
2107    print("The rank of a vector bundle can non be negative");
2108    return(1); // return 1 in this case
2109  }
2110  if(r==0) // if we deal with the zero bundle
2111  {
2112    return( 1 ); // return 1 in this case
2113  }
2114  //-------------------------------------------
2115  // from now on we are in the case r > 0
2116  //-------------------------------------------
2117  c=append_by_zeroes(r, c);
2118  c=c[1..r];
2119  list I; // the partition (0,1,...,r-1) will be stored here
2120  int i;
2121  for(i=0;i<=r-1;i++)
2122  {
2123    I=I+list(i);
2124  }
2125  list PU = partUnder(I); // compute the partitions under I
2126  int sz=size(PU); // get their number
2127  poly rez; // the result will be computed here
2128  list J;
2129  poly cf;
2130  // implement the Lascoux formula
2131  for(i=1;i<=sz;i++)
2132  {
2133    J=PU[i];
2134    cf = IJcoef(I,J)/bigint(2)^( r*(r-1) div 2-sum(J)  );
2135    rez = rez + cf * SchurCh(J, c );
2136  }
2137  return(rez);
2138}
2139example
2140{
2141  "EXAMPLE:";echo =2;
2142  ring r=0, (c(1..3)), ws(1,2,3);
2143  list l=c(1..3);
2144  // the total Chern class of the second exterior power of a vector bundle of rank 3
2145  print(chWedge2LP(3, l));
2146}
2147//---------------------------------------------------------------------------------------
2148
2149proc todd(list c, list #)
2150"USAGE:   todd(l [, n] );  l a list of polynomials, n integer
2151RETURN:   list of polynomials
2152PURPOSE:  computes [the first n] terms of the Todd class
2153EXAMPLE:  example todd; shows an example
2154NOTE:     returns an empty list if l is empty
2155"
2156{
2157  int i, j, k;
2158  // insure that the entries of c are polynomials
2159  // in order to be able to apply maps
2160  for(i=1;i<=size(c); i++)
2161  {
2162    c[i]=poly(c[i]);
2163  }
2164  int n;
2165  # = integer_list(#)[1]; // take into account only the first integer entries that are positive
2166  if( size(#) == 0 ) // if there are no  optional parameters
2167  {
2168    n = size(c);
2169  }
2170  else
2171  {
2172    // set n to be 0, if the parameter is non-positive,
2173    // set n to the value of the parameter otherwise
2174    n = max( #[1], 0 );
2175    c = append_by_zeroes(n, c); // append c by zeroes if the length of c  is smaller than n
2176    if(n!=0) // throw away the redundant data if n is positive and smaller than the length of c
2177    {
2178      c = c[1..n];
2179    }
2180  }
2181  if(n==0) // return the empty list
2182  {
2183    return(list());
2184  }
2185  else // otherwise proceed as follows
2186  {
2187    def br@=basering; // remember the base ring
2188    // add  additional variables to the base ring
2189    list l11 = ringlist(basering)[2]; 
2190    l11[size(l11)+1] = "a@";
2191    for (int ii = 1; ii <= n; ii++)
2192    {
2193     l11[size(l11)+1] = "c@("+string(ii)+")";
2194    }
2195    ring r@ = create_ring(ringlist(basering)[1], l11, "dp", "no_minpoly");
2196    execute( "map F= br@,"+varstr(br@)+";" ); // define the corresponding inclusion of rings
2197    list c=F(c); // embed c into the bigger ring
2198    list prev;
2199    list next;
2200    next=tdTerms(n, c@(1)); // the Todd class terms of a line budle
2201    list step = tdTerms(n, a@);
2202    poly f;
2203    list hC=c@(1)-a@; // "old" first Chern class
2204    for(k=2;k<=n;k++) // do n-1 iterations
2205    {
2206      prev=next;
2207      next=list();
2208      hC=hC+list( c@(k)-a@*hC[k-1] ); // "old" k-th Chern class
2209      for(i=0;i<k;i++) // these terms have already been computed in the previous iterations
2210      {
2211        next = next + list(prev[i+1]);
2212      }
2213      for(i=k;i<=n;i++) // new values in terms of "old" Chern classes and the Chern root a
2214      {
2215        f=0;
2216        for(j=0; j<=i; j++)
2217        {
2218          f=f + step[j+1]*prev[i-j+1];
2219        }
2220        // substitute the old values of Chern classes
2221        // by their expressions in the new ones and the Chern root a
2222        for(j=1;j<k;j++)
2223        {
2224          f=subst(f, c@(j), hC[j] );
2225        }
2226        f=reduce(f, std(hC[k]) ); // eliminate the Chern root
2227        next = next + list(f);
2228      }
2229    }
2230    next = delete(next, 1); // throw away the zeroth term which is always equal to 1
2231    setring br@; // come back to the initial base ring
2232    execute("map FF = r@,"+varstr(br@)+",0, c[1..n];"); // define the specialization homomorphism
2233    return( FF( next ) ); // bring the result to the initial ring
2234  }
2235}
2236example
2237{
2238  "EXAMPLE:";echo =2;
2239  // the terms of the Todd class up to degree 5
2240  // in terms of the Chern classes c(1), c(2), c(3), c(4), c(5)
2241  ring r=0, (c(1..5)), dp;
2242  list l=c(1..5);
2243  print( todd( l ) );
2244
2245  // in the same situation compute only first two terms
2246  print( todd(l, 2) );
2247
2248  // compute the first 5 terms corresponding to the Chern classes c(1), c(2)
2249  l=c(1..2);
2250  print( todd(l, 5) );
2251}
2252//------------------------------------------------------------------------------------------
2253
2254proc toddE(list c)
2255"USAGE:   toddE(l);  l a list of polynomials
2256RETURN:   polynomial
2257PURPOSE:  computes the highest relevant term of the Todd class
2258EXAMPLE:  example toddE; shows an example
2259NOTE:     returns an empty list if l is empty,
2260          very inefficient because the elimination is used, included for comparison with todd(c)
2261"
2262{
2263  int i;
2264  for(i=1;i<=size(c);i++)
2265  {
2266    c[i]=poly( c[i] );
2267  }
2268  int n=size(c);
2269  if(n==0) // return the empty list if c is empty
2270  {
2271    return(list());
2272  }
2273  else
2274  {
2275    def br@=basering; // remember the base ring
2276    // add additional variables a@(1..n), x@ to the base ring
2277    list l12 = "x@";
2278    l12 = l12+ ringlist(basering)[2];
2279    for (int ii = 1; ii <= n; ii++)
2280    {
2281     l12[size(l12)+1] = "a@("+string(ii)+")";
2282    }
2283    ring r@ = create_ring(ringlist(basering)[1], l12, "lp", "no_minpoly");
2284    execute( "map F= br@,"+varstr(br@)+";" ); // define the corresponding inclusion of rings
2285    list c=F(c); // embed c into the bigger ring
2286    int j;
2287    int k;
2288    poly E = a@(1); // to be the product of the Chern roots that will be eliminated later
2289    list ss=tdTerms( n, a@(1) );
2290    list next; // to be the terms of the Todd class corresponding to the next Chern root
2291    for(i=2;i<=n;i++) // compute the terms of the Todd class in terms of the Chern roots
2292    {
2293      E=E*a@(i); // to compute the product of variables to be eliminated
2294      next=tdTerms( n, a@(i) );
2295      for(j=n;j>=1;j--)
2296      {
2297        for(k=0;k<j;k++)
2298        {
2299          ss[j+1]=ss[j+1]+ss[k+1]*next[j-k+1];
2300        }
2301      }
2302    }
2303    ideal I=x@ - ss[n+1]; // formula for the highest degree term of the Todd class
2304    list sym=symm(list(a@(1..n))); // expressions for the Chern classes in terms of the Chern roots
2305    for(i=1;i<=n;i++)
2306    {
2307      I=I, c[i]-sym[i]; // add the relations
2308    }
2309    I=simplify(elim(I, E), 1); // eliminate the Chern roots
2310    poly rez=-subst(I[1],x@, 0); // get the required formula
2311    setring br@; // come back to the initial base ring
2312    // define the specialization homomorphism (all added variables are set to zero)
2313    execute( "map FF = r@,0, "+varstr(br@)+";" );
2314    poly rez=FF( rez ); // bring the result back to the initial base ring
2315    return(rez);
2316  }
2317}
2318example
2319{
2320  "EXAMPLE:";echo =2;
2321  // first  3 terms of the Todd class in terms of the Chern classes c(1), c(2), c(3)
2322  ring r=0, (c(1..3)), dp;
2323  list l;
2324  //first term
2325  l=c(1);
2326  print( toddE( l ) );
2327  // second term
2328  l=c(1..2);
2329  print( toddE( l ) );
2330  // third term
2331  l=c(1..3);
2332  print( toddE( l ) );
2333}
2334//---------------------------------------------------------------------------------
2335
2336proc  Bern(int n)
2337"USAGE:   Bern(n);  n non-negative integer
2338RETURN:   list of numbers
2339PURPOSE:  computes the list of (second) Bernoulli numbers from B(0) to B(n)
2340EXAMPLE:  example Bern; shows an example
2341NOTE:     needs a base ring to be defined, returns an empty list if n is negative,
2342          uses the Akiyama-Tanigawa algorithm
2343"
2344{
2345  // the Akiyama-Tanigawa algorithm
2346  //could be replaced by a more efficient one
2347  list rez, steprez;
2348  int i, j;
2349  if(n<0) // if n is negative, return the empty list
2350  {
2351    return(list());
2352  }
2353  for(i=0;i<=n;i++)
2354  {
2355    steprez=steprez+list( 1/number(i+1) );
2356    for(j=i;j>=1;j--)
2357    {
2358      steprez[j]=j*(steprez[j]-steprez[j+1]);
2359    }
2360    rez=rez+list(steprez[1]);
2361  }
2362  return(rez);
2363}
2364example
2365{
2366  "EXAMPLE:";echo =2;
2367  // first 10 Bernoulli numbers: B(0), ..., B(9)
2368  ring r=0,(t), dp;
2369  print( Bern(9) );
2370}
2371//---------------------------------------------------------------------------------
2372
2373proc tdCf(int n)
2374"USAGE:   tdCf(n);  n integer
2375RETURN:   list of rational numbers
2376PURPOSE:  computes up to degree n the coefficients of the Todd class of a line bundle
2377EXAMPLE:  example tdCf; shows an example
2378NOTE:
2379"
2380{
2381  list rez=Bern(n); // notice that Bern(n) is able to take care of negative n
2382  int i;
2383  for(i=1;i<=n+1;i++)
2384  {
2385  rez[i]=rez[i]/factorial(i-1);
2386  }
2387  return(rez);
2388}
2389example
2390{
2391  "EXAMPLE:";echo =2;
2392  // first 5 coefficients
2393  ring r=0,(t), dp;
2394  print( tdCf(4) );
2395}
2396//---------------------------------------------------------------------------------
2397
2398proc tdTerms(int n, poly f)
2399"USAGE:   tdTerms(n, f);  n integer, f polynomial
2400RETURN:   list of polynomials
2401PURPOSE:  computes the terms of the Todd class of the line bundle with the Chern root f
2402EXAMPLE:  example tdTerms; shows an example
2403NOTE:
2404"
2405{
2406  list rez=Bern(n); // notice that Bern(n) takes care of negative n
2407  int i;
2408  for(i=1;i<=n+1;i++)
2409  {
2410  rez[i]=( rez[i]/factorial(i-1) )* f^(i-1);
2411  }
2412  return(rez);
2413}
2414example
2415{
2416  "EXAMPLE:";echo =2;
2417  ring r=0, (t), ls;;
2418  // the terms of the Todd class of a line bundle with Chern root t up to degree 4
2419  print( tdTerms(4, t) );
2420}
2421//---------------------------------------------------------------------------------
2422
2423proc tdFactor(int n, poly t)
2424"USAGE:   tdFactor(n, a);  n integer, a polynomial
2425RETURN:   polynomial
2426PURPOSE:  computes up to degree n the Todd class
2427          of the line bundle corresponding to the Chern root t
2428EXAMPLE:  example tdFactor; shows an example
2429NOTE:     returns 0 if n is negative
2430"
2431{
2432  int i;
2433  poly rez=0;
2434  list l=Bern(n); // get the coefficients
2435  for(i=0; i<=n; i++) // form the polynomial
2436  {
2437    rez=rez+(l[i+1]/factorial(i))*t^i;
2438  }
2439  return(rez);
2440}
2441example
2442{
2443  "EXAMPLE:";echo =2;
2444  // the Todd class up do degree 4
2445  ring r=0,(t), ls;
2446  print( tdFactor(4, t) );
2447}
2448//---------------------------------------------------------------------------------
2449
2450proc cProj(int n)
2451"USAGE:   cProj(n);  n  integer
2452RETURN:   list of integers
2453PURPOSE:  computes the terms of positive degree of the total Chern class
2454          of the tangent bundle on the complex projective space
2455EXAMPLE:  example cProj; shows an example
2456NOTE:
2457"
2458{
2459  if(n<0)
2460  {
2461    print("The dimension of the projective space must be non-negative!");
2462    return(list()); // return the empty list in this case
2463  }
2464  else
2465  {
2466    list rez;
2467    int i;
2468    for(i=1;i<=n;i++)
2469    {
2470      rez=rez+list( binomial(n+1, i) );
2471    }
2472    return(rez);
2473  }
2474}
2475example
2476{
2477  "EXAMPLE:";echo =2;
2478  ring r=0, (t), dp;
2479  // the coefficients of the total Chern class of the complex projective line
2480  print( cProj(1) );
2481
2482  // the coefficients of the total Chern class of the complex projective line
2483  print( cProj(2) );
2484
2485  // the coefficients of the total Chern class of the complex projective line
2486  print( cProj(3) );
2487}
2488//------------------------------------------------------------------------------------------
2489
2490proc chProj(int n)
2491"USAGE:   chProj(n);  n  integer
2492RETURN:   list of (rational) numbers
2493PURPOSE:  computes the terms of the Chern character of the tangent bundle
2494          on the complex projective space
2495EXAMPLE:  example chProj; shows an example
2496NOTE:
2497"
2498{
2499  if(n<0)
2500  {
2501    print("The dimension of the projective space must be non-negative!");
2502    return( list() ); // return the empty list in this case
2503  }
2504  else
2505  {
2506    list rez=list(number(n));
2507    int i;
2508    for(i=1;i<=n;i++)
2509    {
2510      rez=rez+list( (n+1)/factorial(i) );
2511    }
2512    return(rez);
2513  }
2514}
2515example
2516{
2517  "EXAMPLE:";echo =2;
2518  ring r=0, (t), dp;
2519  // the coefficients of the Chern character of the complex projective line
2520  print( chProj(1) );
2521
2522  // the coefficients of the Chern character of the complex projective plane
2523  print( chProj(2) );
2524
2525  // the coefficients of the Chern character of the complex 3-dimentional projectice space
2526  print( chProj(3) );
2527}
2528//------------------------------------------------------------------------------------------
2529
2530proc tdProj(int n)
2531"USAGE:   tdProj(n);  n  integer
2532RETURN:   list of (rational) numbers
2533PURPOSE:  computes the terms of the Todd class
2534          of the (tangent bundle of the) complex projective space
2535EXAMPLE:  example tdProj; shows an example
2536NOTE:
2537"
2538{
2539  if(n<0)
2540  {
2541    print("The dimension of the projective space must be non-negative!");
2542    return( list() ); // return the empty list in this case
2543  }
2544  else
2545  {
2546    def br@=basering; // remember the base ring
2547    ring r@= 0, t@, lp; // ring with one variable t@
2548    ideal T=std( t@^(n+1) );
2549    poly f= tdFactor(n, t@);
2550    f=reduce( f^(n+1), T);
2551    matrix C = coeffs(f, t@);
2552    list rez;
2553    int i;
2554    for(i=0;i<=n;i++)
2555    {
2556     rez=rez+list(C[i+1, 1]);
2557    }
2558    setring br@; // come back to the initial base ring
2559    map FF= r@, 0 ; // define the specialization homomorphism t@=0
2560    return(FF(rez)); // bring the result to the base ring
2561  }
2562}
2563example
2564{
2565  "EXAMPLE:";echo =2;
2566  ring r=0, (t), dp;
2567
2568  // the coefficients of the Todd class of the complex projective line
2569  print( tdProj(1) );
2570
2571  // the coefficients of the Todd class of the complex projective line
2572  print( tdProj(2) );
2573
2574  // the coefficients of the Todd class of the complex projective line
2575  print( tdProj(3) );
2576}
2577//------------------------------------------------------------------------------------------
2578
2579proc eulerChProj(int n, def r, list c)
2580"USAGE:   eulerChProj(n, r, c);  n  integer, r polynomial (or integer), c list of polynomials
2581RETURN:   polynomial
2582PURPOSE:  computes the Euler characteristic of a vector bundle on P_n
2583          in terms of its rank and Chern classes
2584EXAMPLE:  example eulerChProj; shows an example
2585NOTE:
2586"
2587{
2588  if(n<0)
2589  {
2590    print("The dimension of the projective space must be non-negative!");
2591    return(0); // return zero in this case
2592  }
2593  else
2594  {
2595    if(n==0)
2596    {
2597      return(r);
2598    }
2599    // now n is at least 1
2600    c=append_by_zeroes(n, c); // append c by zeroes if its size is smaller than n
2601    c=c[1..n]; // throw away the redundant data
2602    // now the size of c is n
2603    list td = tdProj(n); // terms of the Todd class of P_n
2604    list ch = list(r) + chAll(c); // terms of the Chern character of the vector bundle
2605    return( rHRR(n, ch, td) );
2606  }
2607}
2608example
2609{
2610  "EXAMPLE:";echo =2;
2611  ring h=0, (r, c(1..3)),  ws(0,1,2,3);
2612  list l=c(1..3);
2613  // the Euler characteristic of a vector bundle on the projective line
2614  print( eulerChProj(1, r, l) );
2615
2616  // the Euler characteristic of a vector bundle on the projective plane
2617  print( eulerChProj(2, r, l) );
2618
2619  // the Euler characteristic of a vector bundle on P_3
2620  print( eulerChProj(3, r, l) );
2621
2622  // assume now that we have a bundle framed at a subplane of P_3
2623  // this implies c(1)=c(2)=0
2624  l= 0, 0, c(3);
2625
2626  // the Euler characteristic is
2627  print( eulerChProj(3, r, l) );
2628  // which implies that c(3) must be even in this case
2629}
2630//-------------------------------------------------------
2631
2632proc chNumbersProj(int n)
2633"USAGE:   chNumbersProj(n);  n  integer
2634RETURN:   list of integers
2635PURPOSE:  computes the Chern numbers of the projective space P_n
2636EXAMPLE:  example chNumbersProj; shows an example
2637NOTE:
2638"
2639{
2640  return( chNumbers( n, cProj(n) ) );
2641}
2642example
2643{
2644  "EXAMPLE:";echo =2;
2645  ring h=0, (t), dp;
2646  // The Chern numbers of the projective plane P_2:
2647  print( chNumbersProj(2) );
2648
2649  // The Chern numbers of P_3:
2650  print( chNumbersProj(3) );
2651}
2652//-------------------------------------------------------
2653
2654proc classpoly(list l, poly t)
2655"USAGE:   classpoly(l, t);   l list of polynomials, t polynomial
2656RETURN:   polynomial
2657PURPOSE:  computes the polynomial in t with coefficients being the entries of l
2658EXAMPLE:  example classpoly; shows an example
2659NOTE:
2660"
2661{
2662  int n=size(l);
2663  poly pow=1; // powers of t will be compured here
2664  poly rez=0; // result will be computed here
2665  int i;
2666  for(i=1; i<=n; i++)
2667  {
2668    pow=pow*t; // compute the required power of t
2669    // add the i-th entry of l multiplied by the corresponding power of t to the result
2670    rez=rez + l[i]*pow;
2671  }
2672  return( rez );
2673}
2674example
2675{
2676  "EXAMPLE:";echo=2;
2677  ring r=0, (c(1..5), t), ds;
2678  list l=c(1..5);
2679  // get the polynomial c(1)*t + c(2)*t^2 + ... + c(5)*t^5
2680  print( classpoly(l, t) );
2681}
2682//----------------------------------------------------------------------------------------
2683
2684proc chernPoly(list c, poly t)
2685"USAGE:   chernPoly(c, t);   c list of polynomials, t polynomial
2686RETURN:   polynomial
2687PURPOSE:  computes the Chern polynomial in t
2688EXAMPLE:  example chernPoly; shows an example
2689NOTE: does the same as toddPoly(...)
2690"
2691{
2692  return( 1+classpoly(c, t) );
2693}
2694example
2695{
2696  "EXAMPLE:";echo=2;
2697  ring r=0, (c(1..5), t), ds;
2698  list l=c(1..5);
2699  // get the Chern polynomial 1 + c(1)*t + c(2)*t^2 + ... + c(5)*t^5
2700  print( chernPoly(l, t) );
2701}
2702//----------------------------------------------------------------------------------------
2703
2704proc chernCharPoly(poly r, list ch, poly t)
2705"USAGE:   chernCharPoly(r, ch, t);   r polynomial, ch list of polynomials, t polynomial
2706RETURN:   polynomial
2707PURPOSE:  computes the polynomial in t corresponding to the Chern character
2708EXAMPLE:  example chernpoly; shows an example
2709NOTE:
2710"
2711{
2712  return( r+classpoly(ch, t) );
2713}
2714example
2715{
2716  "EXAMPLE:";echo=2;
2717  ring h=0, (r, ch(1..5), t), ds;
2718  list l=ch(1..5);
2719  // get the polynomial r + ch(1)*t + ch(2)*t^2 + ... + ch(5)*t^5
2720  print( chernCharPoly(r, l, t) );
2721}
2722//----------------------------------------------------------------------------------------
2723
2724proc toddPoly(list td, poly t)
2725"USAGE:   toddPoly(td, t);   td list of polynomials, t polynomial
2726RETURN:   polynomial
2727PURPOSE:  computes the polynomial in t corresponding to the Todd class
2728EXAMPLE:  example toddPoly; shows an example
2729NOTE:     does the same as chernPoly(...)
2730"
2731{
2732  return( 1+classpoly(td, t) );
2733}
2734example
2735{
2736  "EXAMPLE:"; echo=2;
2737  ring r=0, (td(1..5), c(1..5), t), ds;
2738  list l=td(1..5);
2739  // get the polynomial 1 + td(1)*t + td(2)*t^2 + ... + td(5)*t^5
2740  print( toddPoly(l, t) );
2741}
2742//---------------------------------------------------------------------------------------
2743
2744proc rHRR(int N, list ch, list td)
2745"USAGE:   rHRR( N, ch, td);  N integer, ch, td  lists of polynomials
2746RETURN:   polynomial
2747PURPOSE:  computes the the main ingredient of the right-hand side
2748          of the Hirzebruch-Riemann-Roch formula
2749EXAMPLE:  example rHRR; shows an example
2750NOTE:     in order to get the right-hand side of the HRR formula
2751          one needs to be able to compute the degree of the output of this procedure
2752"
2753{
2754  poly rez; // to be the result
2755  int i;
2756  int nch=size(ch); // length of ch
2757  int ntd=size(td); // length of td
2758  for(i=1; i<=N+1; i++) // compute the highest degree term of ch.td
2759  {
2760    if( (i<=nch) && (N-i+2 <= ntd) )
2761    {
2762      rez = rez + ch[i]*td[N-i+2];
2763    }
2764  }
2765  return(rez);
2766}
2767example
2768{
2769  "EXAMPLE:"; echo=2;
2770  ring r=0, (td(0..3), ch(0..3)), dp;
2771  // Let ch(0), ch(1), ch(2), ch(3) be the terms of the Chern character
2772  // of a vector bundle E on a 3-fold X.
2773  list c = ch(0..3);
2774  // Let td(0), td(1), td(2), td(3) be the terms of the Todd class of X.
2775  list t = td(0..3);
2776  // Then the highest term of the product ch(E).td(X) is:
2777  print( rHRR(3, c, t) );
2778}
2779//---------------------------------------------------------------------------------------
2780
2781proc SchurS(list I, list S)
2782"USAGE:   SchurS(I, S);  I list of integers representing a partition, S list of polynomials
2783RETURN:   poly
2784PURPOSE:  computes the Schur polynomial in the Segre classes S (of the dual vector bundle),
2785          i.e., in the  complete homogeneous symmetric polynomials, with respect to the partition I
2786EXAMPLE:  example SchurS; shows an example
2787NOTE:     if S are the Segre classes of the tautological bundle on a grassmanian,
2788          this gives the cohomology class of a Schubert cycle
2789"
2790{
2791  int m=size(I); // size of I
2792  S=list(1)+S; // add the zeroth Segre class
2793  int szS=size(S); // size of S
2794  int h,k;
2795  int in; // variable for the index of the required Segre class
2796  // construct the required m x m matrix from the first determinantal (Jacobi-Trudi) formula
2797  matrix M[m][m];
2798  for(h=1;h<=m;h++)
2799  {
2800    for(k=1;k<=m;k++)
2801    {
2802      in=I[k]+k-h; // compute the index
2803      if(in<0) // if it is negative, assume the corresponding Segre class to be zero
2804      {
2805        M[h,k]=0;
2806      }
2807      else
2808      {
2809        if(in>=szS) // if it is bigger than the number of the highest available Segre class in S
2810        {
2811          M[h, k]=0; // assume the corresponding Segre class is zero
2812        }
2813        else // otherwise
2814        {
2815          M[h, k]= S[in+1]; // use a value from S for the corresponding Segre class
2816        }
2817      }
2818    }
2819  }
2820  return(det(M)); // return the determinant of the computed matrix
2821}
2822example
2823{
2824  "EXAMPLE:"; echo=2;
2825  // The Schur polynomial corresponding to the partition 1,2,4
2826  // and the Segre classes 1, s(1), s(2),..., s(6)
2827  ring r=0,(s(1..6)), dp;
2828  list I=1,2,4;
2829  list S=s(1..6);
2830  print( SchurS(I, S) );
2831  // compare this with the Schur polynomial computed using Chern classes
2832  list C=chDual(chern(S));
2833  print( SchurCh(I, C) );
2834}
2835//---------------------------------------------------------------------------------------
2836
2837proc SchurCh(list I, list C)
2838"USAGE:   SchurCh(I, C);  I list of integers representing a partition, C list of polynomials
2839RETURN:   poly
2840PURPOSE:  computes the Schur polynomial in the Chern classes C,
2841          i.e., in the elementary symmetric polynomials, with respect to the partition I
2842EXAMPLE:  example SchurCh; shows an example
2843NOTE:     if C are the Chern classes of the tautological bundle on a grassmanian,
2844          this gives the cohomology class of a Schubert cycle
2845"
2846{
2847  I=dualPart(I); // dual partition to I
2848  int m=size(I); // size of I
2849  C=list(1)+C; // add the zeroth Chern class
2850  int szC=size(C); // size of C
2851  int h,k;
2852  int in; // variable for the index of the required Chern class
2853  // construct the required m x m matrix from the second determinantal (Jacobi-Trudi) formula
2854  matrix M[m][m];
2855  for(h=1;h<=m;h++)
2856  {
2857    for(k=1;k<=m;k++)
2858    {
2859      in=I[k]+k-h; // compute the index
2860      if(in<0) // if it is negative, assume the corresponding Chern class to be zero
2861      {
2862        M[h,k]=0;
2863      }
2864      else
2865      {
2866        if(in>=szC) // if it is bigger than the number of the highest available Chern class in C
2867        {
2868          M[h, k]=0; // assume the corresponding Chern class is zero
2869        }
2870        else // otherwise
2871        {
2872          M[h, k]= C[in+1]; // use a value from C for the corresponding Chern class
2873        }
2874      }
2875    }
2876  }
2877  return(det(M)); // return the determinant of the computed matrix
2878}
2879example
2880{
2881  "EXAMPLE:"; echo=2;
2882  // The Schur polynomial corresponding to the partition 1,2,4
2883  // and the Chern classes c(1), c(2), c(3)
2884  ring r=0,(c(1..3)), dp;
2885  list I=1,2,4;
2886  list C=c(1..3);
2887  print( SchurCh(I, C) );
2888  // Compare this with the Schur polynomial computed using Segre classes
2889  list S=segre( chDual( list(c(1..3)) ), 6 );
2890  print(SchurS(I,S));
2891}
2892//---------------------------------------------------------------------------------------
2893
2894proc part(int m, int n)
2895"USAGE:   part( m, n );  m positive integer, n non-negative integer
2896RETURN:   list of lists
2897PURPOSE:  computes all partitions of integers not exceeding n into m non-negative summands
2898EXAMPLE:  example part; shows an example
2899NOTE:     if n is negative or m is non-positive, the list with one empty entry is returned
2900"
2901{
2902  if( n<0 ) // if n is negative
2903  {
2904    return(list(list())); // return the list with one empty entry
2905  }
2906  if(n==0) // if n equals 0, there is only one partition of 0 into m non-negative summands
2907  {
2908    return(list(listSame(0,m))); // return the list with one entry consistion of m zeroes
2909  }
2910  // otherwise proceed recursively
2911  list rez=part(m, n-1); // get all partitions for n-1
2912  int i;
2913  for(i=1;i<=m;i++) // for every i between 1 and m, add the partitions with exactly
2914  {
2915    rez=rez + appendToAll( part(m-i, n-1), listSame(n, i) ); // i summands equal to n
2916  }
2917  return(rez); // return the result
2918}
2919example
2920{
2921  "EXAMPLE:"; echo=2;
2922  // partitions into 3 summands of numbers not exceeding 1
2923  print( part(3, 1) );
2924}
2925//---------------------------------------------------------------------------------------
2926
2927proc dualPart(list I, list #)
2928"USAGE:   dualPart( I [,N] );  I list of integers, N integer
2929RETURN:   list of integers
2930PURPOSE:  computes the partition dual (conjugate) to I
2931EXAMPLE:  example dualPart; shows an example
2932NOTE:     the result is extended by zeroes to length N if an optional integer
2933          parameter N is given and the length of the computed  dual partition
2934          is smaller than N
2935"
2936{
2937  int m= size(I); // size of I
2938  if(m==0) // if I is the empty list
2939  {
2940    print("You are trying to compute the dual of the empty partition!");
2941    print("The partition with one zero is returned.");
2942    return(list(0));
2943  }
2944  // compute the dual partition
2945  list J; // the result will be computed here
2946  int i;
2947  int j=I[1];
2948  int k;
2949  for(k=1;k<=j;k++)
2950  {
2951    J=list(m)+J;
2952  }
2953  for(i=2;i<=m;i++)
2954  {
2955    j=I[i]-I[i-1];
2956    for(k=1;k<=j;k++)
2957    {
2958      J=list(m-i+1)+J;
2959    }
2960  }
2961  if(size(J)==0) // if the dual partition J is empty (if I consists of zeroes)
2962  {
2963    J = list(0); // add zero to the result
2964  }
2965  if(size(#)>0) // if there is an optional parameter N
2966  {
2967    if( is_integer( #[1] ) ) // if the parameter is an integer,
2968    {
2969      if( size(J) < #[1] ) // if N  is bigger than the length of J,
2970      {
2971        J=listSame(0, #[1]-size(J))+J; // extend J by zeroes to length N
2972      }
2973    }
2974  }
2975  return(J); // return the result
2976}
2977example
2978{
2979  "EXAMPLE:"; echo =2;
2980  // dual partition to (1, 3, 4):
2981  list I = 1, 3, 4;
2982  print( dualPart(I) );
2983}
2984//---------------------------------------------------------------------------------------
2985
2986proc PartC(list I, int m)
2987"USAGE:   PartC( I, m);  I list of integers, m integer
2988RETURN:   list of integers
2989PURPOSE:  computes the complement of a partition with respect to m
2990EXAMPLE:  example PartC; shows an example
2991NOTE:     returns the zero partition if the maximal element of the partition is smaller than m
2992"
2993{
2994  int n=size(I); // size of I
2995  if( m<I[n] ) // if m is smaller than the last term of I,
2996  {
2997    // give a warning
2998    print("You are trying to compute a complement of a partition with respect");
2999    print("to a number that is smaller than the maximal summand of the partition!");
3000    print("The zero partition is returned.");
3001    return(list(0)); // and return the zero partition
3002  }
3003  list J; // the result will be computed here
3004  int i;
3005  for(i=n;i>=1;i--) // invert the order of numbers
3006  {
3007    J=J+list(m-I[i]); // and substitute them by their complemenst to m
3008  }
3009  return(J); // return the result
3010}
3011example
3012{
3013  "EXAMPLE:"; echo =2;
3014  // Complement of the partition (1, 3, 4) with respect to 5
3015  list I = 1, 3, 4;
3016  print( PartC(I, 5) );
3017}
3018//---------------------------------------------------------------------------------------
3019
3020proc partOver(int n, list J)
3021"USAGE:   partOver( n, J);  n integer, J list of integers (partition)
3022RETURN:   list of lists
3023PURPOSE:  computes the partitions over a given one  with summands not exceeding n
3024EXAMPLE:  example partOver; shows an example
3025NOTE:
3026"
3027{
3028  int m=size(J); // size of J
3029  if( m==0 ) // if J is an empty list
3030  {
3031    // give a warning
3032    print("You are trying to compute partitions over an empty partition!");
3033    return( list() ); // and return the empty list
3034  }
3035  if( J[m] > n ) // if the biggest summand of the partition is bigger than n
3036  {
3037    return( list( ) ); // return the empty list
3038  }
3039  if( J[m] == 0 ) // if J consists of zeroes
3040  {
3041    return( part(m,n) ); // return all partitions of n into m summands
3042  }
3043  // now J is non-empty, contains con-zero summands, has partitions over it
3044  list rez1; // the result will be computed here
3045  int i,j;
3046  if(m==1) // if J has only one element
3047  {
3048    for(j=J[1]; j<=n; j++) // run through the integers from J[1] to n
3049    {
3050      rez1=rez1 + list(j); // add the corresponding one element lists to the result
3051    }
3052    return(rez1); // return the result
3053  }
3054  // now J has at least two elements
3055  // get the partitions over the partition without the last summand
3056  list rez = partOver(n, delete(J, m));
3057  int sz=size(rez); // number of such partitions
3058  list P;
3059  int last;
3060  for(i=1; i<=sz; i++) // run trough all such partitions
3061  {
3062    P=rez[i]; // for each partition P of this type
3063    last = max( P[size(P)], J[m] );
3064    for(j = last;j<= n;j++) // run through the integers exceeding the last summands of P and J
3065    {
3066      // append them to P at the end and add the resulting partition to the result
3067      rez1=rez1 + list(P+list(j));
3068    }
3069  }
3070  return(rez1); // return the result
3071}
3072example
3073{
3074  "EXAMPLE:"; echo =2;
3075  // Partitions over the partition (3, 3, 4) with summands not exceeding 4
3076  list I = 3, 3, 4;
3077  print( partOver(4, I) );
3078}
3079//---------------------------------------------------------------------------------------
3080
3081proc partUnder(list J)
3082"USAGE:   partUnder(J);  J list of integers (partition)
3083RETURN:   list of lists
3084PURPOSE:  computes the partitions under a given one
3085EXAMPLE:  example partUnder; shows an example
3086NOTE:
3087"
3088{
3089  int m=size(J); // size of J
3090  if(m==0) // if J is empty
3091  {
3092    return(list()); // return an empty list
3093  }
3094  list rez1; // the result will be computed here
3095  int i;
3096  if(m==1) // if J contains only one element
3097  {
3098    for(i=0; i<=J[1]; i++)
3099    {
3100      rez1=rez1+list(list(i));
3101    }
3102  }
3103  // now J contains at least two elements
3104  list rez;
3105  int Jlast=J[m]; // last element of J
3106  rez = partUnder(delete(J, m)); // partitions under J without the last element
3107  int j;
3108  int sz=size(rez); // their number
3109  list P;
3110  int last;
3111  for(i=1; i<=sz; i++) // for every such partition
3112  {
3113    P = rez[i];
3114    last = P[size(P)];
3115    for(j = last;j<=Jlast ;j++) // for every number between its last entry and the last entry of J
3116    {
3117      // append that number to the end of the partition
3118      // and append the resulting partition to the final result
3119      rez1 = rez1 + list(P+list(j));
3120    }
3121  }
3122  return(rez1);
3123}
3124example
3125{
3126  "EXAMPLE:"; echo =2;
3127  // Partitions under the partition (0, 1, 1)
3128  list I = 0, 1, 1;
3129  print( partUnder(I) );
3130}
3131//-------------------------------------------------------------------------------------------
3132
3133proc SegreA(ideal I)
3134"USAGE:   SegreA(I);  I an ideal
3135RETURN:   list of integers
3136PURPOSE:  computes the Segre classes of the subscheme defined by I
3137EXAMPLE:  example SegreA; shows an example
3138NOTE:
3139"
3140{
3141  if( !homog(I) ) // if the ideal is not homogeneous
3142  {
3143    print("You are trying to compute the Segre class of a non-homogeneous ideal!");
3144    print("The ideal must be homogeneous, an empty list is returned.");
3145    return( list() );
3146  }
3147  // modify the generators of the ideal so that all of them are of the same degree
3148  I=equal_deg(I);
3149  int degGen=deg(I[1]); // this degree d is stored in this variable
3150  int i;
3151  list rez; // define the variable for the result
3152  int n=nvars(basering)-1; // the dimension of the projective space
3153  if(degGen==-1) // if the ideal is zero
3154  {
3155    for(i=0;i<=n;i++)
3156    {
3157      rez=rez+ list( int((-1)^i*binomial(n+i, i)) );
3158    }
3159    return(rez);
3160  }
3161  int sz=ncols(I); // the number of new generators is stored here
3162  def br@=basering; // remember the base ring
3163  // add additional variables t@(1), ... , t@(sz) and u@ to the base ring
3164  list l13 = ringlist(basering)[2];
3165  for (int ii = 1; ii <= sz; ii++)
3166  {
3167   l13[size(l13)+1] = "t@("+string(ii)+")";
3168  }
3169  l13[size(l13)+1] = "u@";
3170  ring r@ = create_ring(ringlist(basering)[1], l13, "dp", "no_minpoly");
3171  execute( "map F= br@,"+varstr(br@)+";" ); // define the corresponding inclusion of rings
3172  ideal I=F(I); // the ideal generated by I in the new ring
3173  ideal J(0..n); // define n+1 ideals J(0), ... , J(n)
3174  // compute the ideal of the Rees algebra of the ideal I:
3175  for(i=1; i<=sz; i++) // consider the ideal generated by t@(i)-u@*I[i]
3176  {
3177    J(0)=J(0) + ideal( t@(i)-u@*I[i] );
3178  }
3179  J(0)=eliminate(J(0), u@); // and eliminate the variable u@
3180  ideal T=t@(1..sz); // define the ideal generated by the additional variables t@(1), ... , t@(sz)
3181  for(i=1;i<=n;i++)// for all i=1, ... n define J(j) as in 3.6 of the Aluffi's paper
3182  {
3183    // add a random general linear form in variables t@(1), ... , t@(n) to J(i-1)
3184    J(i)=sat( random_hypersurf(J(i-1), T)   , T)[1]; // and saturate with respect to T
3185  }
3186  poly prd=product(T); // compute the product of t@(i)
3187  poly cl;
3188  poly mlt;
3189  for(i=0;i<=n;i++)
3190  {
3191    // eliminate all t@(i) from J(i)
3192    // compute the degree of the scheme defined by this ideal
3193    // and use it to compute the class corresponding to c(O(d))^n * G\otimes O(d)
3194    cl=cl+mult(std(eliminate(J(i), prd)))*u@^i*(1+degGen*u@)^(n-i);
3195    // the (n+1)-st power of the inverse of the Chern class of O(d)
3196    mlt=mlt+binomial(n+i, i)*(-degGen*u@)^i;
3197  }
3198  poly resPoly;
3199  // compute the Segre class by the Aluffi's formula from Proposition 3.1 as polynomial in u@
3200  resPoly= NF( 1-cl*mlt, u@^(n+1));
3201  matrix cf=coeffs(resPoly, u@); // coefficients of the Segre class
3202  rez=list(); // empty the list
3203  for(i=0;i<=n;i++) // fill the list with the the Segre classes in positive degrees
3204  {
3205    if( i < nrows(cf) ) // if i is not bigger than the maximal degree of non-zero Segre classes
3206    {
3207      rez=rez+list(int(cf[i+1,1]));
3208    }
3209    else // otherwise fill the list with zeroes
3210    {
3211      rez=rez+list( int(0) );
3212    }
3213  }
3214  return(rez);
3215}
3216example
3217{
3218  "EXAMPLE:";echo =2;
3219  // Consider a 3-dimensional projective space
3220  ring r = 0, (x, y, z, w), dp;
3221  // Consider 3 non-coplanar lines trough one point and compute the Segre class
3222  ideal I=xy, xz, yz;
3223  I;
3224  SegreA(I);
3225  // Now consider 3 coplanar lines trough one point and its Segre class
3226  ideal J=w, x*y*(x+y);
3227  J;
3228  SegreA(J);
3229}
3230//-------------------------------------------------------------------------------------------
3231
3232proc FultonA(ideal I)
3233"USAGE:   FultonA(I);  I an ideal
3234RETURN:   list of integers
3235PURPOSE:  computes the Fulton classes of the subscheme defined by I
3236EXAMPLE:  example FultonA; shows an example
3237NOTE:
3238"
3239{
3240  if( !homog(I) ) // if the ideal is not homogeneous
3241  {
3242    print("You are trying to compute the Segre class of a non-homogeneous ideal!");
3243    print("The ideal must be homogeneous, an empty list is returned.");
3244    return( list() );
3245  }
3246  // modify the generators of the ideal so that all of them are of the same degree
3247  I=equal_deg(I);
3248  int degGen=deg(I[1]); // this degree d is stored in this variable
3249  int i;
3250  list rez; // define the variable for the result
3251  int n=nvars(basering)-1; // the dimension of the projective space
3252  if(degGen==-1) // if the ideal is zero
3253  {
3254    rez=rez+list(int(1));
3255    for(i=1;i<=n;i++)
3256    {
3257      rez=rez+ list( int(0) );
3258    }
3259    return(rez);
3260  }
3261  int sz=ncols(I); // the number of new generators is stored here
3262  def br@=basering; // remember the base ring
3263  // add additional variables t@(1), ... , t@(sz) and u@ to the base ring
3264  list l14 = ringlist(basering)[2];
3265  for (int ii = 1; ii <= sz; ii++)
3266  {
3267   l14[size(l14)+1] = "t@("+string(ii)+")";
3268  }
3269  l14[size(l14)+1] = "u@";
3270  ring r@ = create_ring(ringlist(basering)[1], l14, "dp", "no_minpoly");
3271  execute( "map F= br@,"+varstr(br@)+";" ); // define the corresponding inclusion of rings
3272  ideal I=F(I); // the ideal generated by I in the new ring
3273  ideal J(0..n); // define n+1 ideals J(0), ... , J(n)
3274  // compute the ideal of the Rees algebra of the ideal I:
3275  for(i=1; i<=sz; i++) // consider the ideal generated by t@(i)-u@*I[i]
3276  {
3277    J(0)=J(0) + ideal( t@(i)-u@*I[i] );
3278  }
3279  J(0)=eliminate(J(0), u@); // and eliminate the variable u@
3280  ideal T=t@(1..sz); // define the ideal generated by the additional variables t@(1), ... , t@(sz)
3281  for(i=1;i<=n;i++)// for all i=1, ... n define J(j) as in 3.6 of the Aluffi's paper
3282  {
3283    // add a random general linear form in variables t@(1), ... , t@(n) to J(i-1)
3284    J(i)=sat(   random_hypersurf(J(i-1), T)  , T)[1]; // and saturate with respect to T
3285  }
3286  poly prd=product(T); // compute the product of t@(i)
3287  poly cl;
3288  poly mlt;
3289  for(i=0;i<=n;i++)
3290  {
3291    // eliminate all t@(i) from J(i)
3292    // compute the degree of the scheme defined by this ideal
3293    // the class corresponding to c(O(d))^n * G\otimes O(d)
3294    cl=cl+ mult(std( eliminate( J(i),  prd) ))*u@^i*(1+degGen*u@)^(n-i);
3295    // the (n+1)-st power of the inverse of the Chern class of O(d)
3296    mlt=mlt+binomial(n+i, i)*(-degGen*u@)^i;
3297  }
3298  poly resPoly;
3299  // compute the Fulton class using the Aluffi's formula for the Segre class and
3300  // multiplying it by the Chern class of the projective space (1+u@)^(n+1)
3301  resPoly= NF( (1+u@)^(n+1)*(1-cl*mlt), u@^(n+1) );
3302  matrix cf=coeffs(resPoly, u@); // coefficients of the Fulton class
3303  rez=list(); // empty the list
3304  for(i=0;i<=n;i++) // fill the list with the the Fulton classes in positive degrees
3305  {
3306    if( i < nrows(cf) ) // if i is not bigger than the maximal degree of non-zero Fulton classes
3307    {
3308      rez=rez+list(int(cf[i+1,1]));
3309    }
3310    else // otherwise fill the list with zeroes
3311    {
3312      rez=rez+list( int(0) );
3313    }
3314  }
3315  return(rez);
3316}
3317example
3318{
3319  "EXAMPLE:";echo =2;
3320  // Consider a 3-dimensional projective space
3321  ring r = 0, (x, y, z, w), dp;
3322  // Consider 3 non-coplanar lines trough one point and compute the Fulton class
3323  ideal I=xy, xz, yz;
3324  I;
3325  FultonA(I);
3326  // Now consider 3 coplanar lines trough one point and its Fulton class
3327  ideal J=w, x*y*(x+y);
3328  J;
3329  FultonA(J);
3330}
3331//-------------------------------------------------------------------------------------------
3332
3333proc CSMA(ideal I)
3334"USAGE:   CSMA(I);  I an ideal
3335RETURN:   list of integers
3336PURPOSE:  computes the Chern-Schwartz-MacPherson classes of the variety defined by I
3337EXAMPLE:  example CSMA; shows an example
3338NOTE:
3339"
3340{
3341  if( !homog(I) ) // if the ideal is not homogeneous
3342  {
3343    print("You are trying to compute the Chern-Schwartz-MacPherson class.");
3344    print("However the input ideal is not homogeneous!");
3345    print("The ideal must be homogeneous, an empty list is returned.");
3346    return( list() );
3347  }
3348  int sz=ncols(I);
3349  int i;
3350  int n=nvars(basering)-1;
3351  bigintmat REZ[1][n+1];
3352  int j;
3353  int szpr;
3354  list pr;
3355  int sgn=-1;
3356  for(i=1;i<=sz;i++)
3357  {
3358    sgn=-sgn;
3359    pr=prds(i,I);
3360    szpr=size(pr);
3361    for(j=1;j<=szpr;j++)
3362    {
3363      REZ=REZ+sgn*CSM_hypersurf(pr[j]);
3364    }
3365  }
3366  list rez;
3367  for(i=0;i<=n;i++)
3368  {
3369    rez=rez+list(REZ[1,i+1]);
3370  }
3371  return(rez);
3372}
3373example
3374{
3375  "EXAMPLE:";echo =2;
3376  // consider the projective plane with homogeneous coordinates x, y, z
3377  ring r = 0, (x, y, z), dp;
3378  // the Chern-Schwartz-MacPherson class of a smooth cubic:
3379  ideal I=x3+y3+z3;
3380  I;
3381  CSMA(I);
3382  // the Chern-Schwartz-MacPherson class of singular cubic
3383  // that is a union of 3 non-collinear lines:
3384  ideal J=x*y*z;
3385  J;
3386  CSMA(J);
3387  // the Chern-Schwartz-MacPherson class of singular cubic
3388  // that is a union of 3 lines passing through one point
3389  ideal K=x*y*(x+y);
3390  K;
3391  CSMA(K);
3392}
3393//-------------------------------------------------------------------------------------------
3394
3395proc EulerAff(ideal I)
3396"USAGE:   EulerAff(I);  I an ideal
3397RETURN:   integer
3398PURPOSE:  computes the Euler characteristic of the affine variety defined by I
3399EXAMPLE:  example EulerAff; shows an example
3400NOTE:
3401"
3402{
3403  int n=nvars(basering);
3404  def br@=basering; // remember the base ring
3405  ring r@ = create_ring(ringlist(basering)[1], "("+varstr(basering)+",homvar@)", "dp", "no_minpoly");
3406  execute( "map F= br@,"+varstr(br@)+";" ); // define the corresponding inclusion of rings
3407  ideal I=F(I);
3408  ideal J=homog(I, homvar@);
3409  ideal JJ=J, homvar@;
3410  return( CSMA(J)[n+1]-CSMA(JJ)[n+1] );
3411}
3412example
3413{
3414  "EXAMPLE:";echo =2;
3415  ring r = 0, (x, y), dp;
3416  // compute the Euler characteristic of the affine elliptic curve y^2=x^3+x+1;
3417  ideal I=y2-x3-x-1;
3418  EulerAff(I);
3419}
3420//-------------------------------------------------------------------------------------------
3421
3422proc EulerProj(ideal I)
3423"USAGE:   EulerProj(I);  I an ideal
3424RETURN:   integer
3425PURPOSE:  computes the highest degree term of the Chern-Schwartz-MacPherson class
3426          of the variety defined by I, which equals the Euler characteristic
3427EXAMPLE:  example EulerProj; shows an example
3428NOTE:     uses CSMA(...)
3429"
3430{
3431  if( !homog(I) ) // if the ideal is not homogeneous
3432  {
3433    print("The ideal must be homogeneous, zero is returned.");
3434    return( int(0) );
3435  }
3436  int n=nvars(basering);
3437  return( CSMA(I)[n] );
3438}
3439example
3440{
3441  "EXAMPLE:";echo =2;
3442  // consider the projective plane with homogeneous coordinates x, y, z
3443  ring r = 0, (x, y, z), dp;
3444  // Euler characteristic of a smooth cubic:
3445  ideal I=x3+y3+z3;
3446  I;
3447  EulerProj(I);
3448  // Euler characteritic of 3 non-collinear lines:
3449  ideal J=x*y*z;
3450  J;
3451  EulerProj(J);
3452  // Euler characteristic of 3 lines passing through one point
3453  ideal K=x*y*(x+y);
3454  K;
3455  EulerProj(K);
3456}
3457//----------------------------------------------------------------------------------------
3458// The procedures below are for the internal usage only
3459//----------------------------------------------------------------------------------------
3460
3461static proc append_by_zeroes(int N, list c)
3462"USAGE:   append_by_zeroes( N, c);  N integer, c a list
3463RETURN:   list
3464PURPOSE:  appends by zeroes up to the length N
3465EXAMPLE:  example append_by_zeroes; shows an example
3466NOTE:
3467"
3468{
3469  int n=size(c);
3470  if(N>n) // if N is greater than the length of c, append c by zeroes up to the length N
3471  {
3472    int i;
3473    for(i=n+1;i<=N;i++)
3474    {
3475      c=c+list( poly(0) );
3476    }
3477  }
3478  return(c);
3479}
3480example
3481{
3482  "EXAMPLE:";echo =2;
3483  ring r = 0, (x, y, z), dp;
3484  list l=(x, y, z);
3485  //append the list by two zeroes and get a list of length 5
3486  print( append_by_zeroes(5, l) );
3487}
3488//-----------------------------------------------------------------------
3489
3490static proc is_integer(def r)
3491"USAGE:   is_integer(r);  r any type
3492RETURN:   1 or 0
3493PURPOSE:  checks whether r is of type int or bigint
3494EXAMPLE:  example is_integer; shows an example
3495NOTE: returns 1 if r is of type int or bigint, otherwise returns 0
3496"
3497{
3498  if( (typeof(r)=="int") || (typeof(r)=="bigint") )
3499  {
3500    return(1);
3501  }
3502  else
3503  {
3504    return(0);
3505  }
3506}
3507example
3508{
3509  "EXAMPLE:";echo =2;
3510  // test on int, bigint, poly
3511  ring r;
3512  int i=12;
3513  bigint j=16;
3514  poly f=x;
3515  print( is_integer(i) );
3516  print( is_integer(j) );
3517  print( is_integer(f) );
3518}
3519//------------------------------------------------------------------------------------
3520
3521static proc integer_list(list l)
3522"USAGE:   integer_list(l);  l list
3523RETURN:   list
3524PURPOSE:  gets the first positive ingerer entries of l, computes their maximum;
3525          used for adjusting the lists of optional parameters that are supposed to be integers
3526EXAMPLE:  example integer_list; shows an example
3527NOTE:     used in chWedge(...) and chSymm(...)
3528"
3529{
3530  int M=0;
3531  int n=size(l);
3532  if(n==0)
3533  {
3534    return(list(l, M));
3535  }
3536  // now n>0
3537  list rez; // the result will be computed here
3538  int i=1;
3539  while( is_integer( l[i] ) ) // take only the first integer entries of l
3540  {
3541    if(l[i]>0) // if they are positive
3542    {
3543      rez=rez+list( l[i] );
3544      if(l[i]>M) // adjust the maximum if necessary
3545      {
3546        M=l[i];
3547      }
3548      i++;
3549    }
3550    else // otherwise get out from the loop
3551    {
3552      break;
3553    }
3554  }
3555  return( list( rez, M)   );
3556}
3557example
3558{
3559  "EXAMPLE:";echo =2;
3560  // the first integer entries of 1,2,3,t are 1,2,3
3561  ring r=0,(t), ls;
3562  list l=1,2,3, t;
3563  print( integer_list(l) );
3564}
3565//---------------------------------------------------------------------------------------
3566
3567static proc appendToAll(list L, list A)
3568"USAGE:   appendToAll( L, A );  L list of lists, A list
3569RETURN:   list
3570PURPOSE:  appends A to every entry of L
3571EXAMPLE:  example appendToAll; shows an example
3572NOTE:
3573"
3574{
3575  int n=size(L);
3576  int i;
3577  for(i=1;i<=n;i++) // run through all elements of L
3578  {
3579    L[i]=L[i]+A; // and append A to each of them
3580  }
3581  return(L);
3582}
3583example
3584{
3585  "EXAMPLE:"; echo=2;
3586  // Consider two lists
3587  list l1, l2;
3588  l1=1,2;
3589  l2=3,4;
3590  // The first one is
3591  print(l1);
3592  // The second one is
3593  print(l2);
3594  // Now consider the list with entries l1 and l2
3595  list L= l1, l2;
3596  print(L);
3597  // and consider a list A
3598  list A = 7,9;
3599  print(A);
3600  // append A to all entries of L
3601  print( appendToAll(L, A) );
3602}
3603//---------------------------------------------------------------------------------------
3604
3605static proc listSame(int n, int k)
3606"USAGE:   listSame( n, k );  n integer, k non-negative integer
3607RETURN:   list
3608PURPOSE:  list with k entries each equal to n
3609EXAMPLE:  example listSame; shows an example
3610NOTE:     if k is negative or zero, the empty list is returned
3611"
3612{
3613  list rez;
3614  int i;
3615  for(i=1;i<=k;i++) // create a list with k entries, each equal to n
3616  {
3617    rez=rez+list(n);
3618  }
3619  return(rez);
3620}
3621example
3622{
3623  "EXAMPLE:"; echo=2;
3624  // list of 5 zeroes
3625  print( listSame(0, 5) );
3626}
3627//---------------------------------------------------------------------------------------
3628
3629static proc IJcoef(list I, list J)
3630"USAGE:   IJcoef( I, J);  J, J lists of integers
3631RETURN:   bigint
3632PURPOSE:  computes the coefficient used in the formula of Lascoux
3633EXAMPLE:  example IJcoef; shows an example
3634NOTE:     these coefficients are denoted (I, J) in the paper of Lascoux
3635"
3636{
3637  int m = size(I);
3638  if(m != size(J)) // if the sizes of I and J are different
3639  {
3640    // give a warning
3641    print("The sizes of the partitions are different!");
3642    print("Zero is returned.");
3643    return( bigint(0) ); // and return zero
3644  }
3645  // now the sizes of I and J are equal m
3646  int h, k;
3647  bigintmat M[m][m]; // construct the required matrix
3648  for(h=1; h<=m; h++)
3649  {
3650    for(k=1; k<=m; k++)
3651    {
3652      M[h,k] = binomial( I[k]+k-1, J[h]+h-1 );
3653    }
3654  }
3655  return( det(M) ); // and return its determinant
3656}
3657example
3658{
3659  "EXAMPLE:"; echo =2;
3660  // The coefficient corresponding to the partitions (1, 3, 4) and (0, 3, 3)
3661  list I = 1, 3, 4;
3662  list J = 1, 3, 3;
3663  print( IJcoef(I, J) );
3664}
3665//---------------------------------------------------------------------------------------
3666
3667static proc invertPart(list l)
3668"USAGE:   invertPart(I);  I list of integers (partition),
3669RETURN:   list of integers
3670PURPOSE:  inverts the ordering of the elements in the list
3671EXAMPLE:  example invertPart; shows an example
3672NOTE:
3673"
3674{
3675  list L;
3676  int sz=size(l);
3677  int i;
3678  for(i=sz;i>=1;i--)
3679  {
3680    L=L+list(l[i]);
3681  }
3682  return(L);
3683}
3684example
3685{
3686  "EXAMPLE:"; echo = 2;
3687  // Invert the ordering of elements in (3, 2, 1)
3688  list l = 3, 2, 1;
3689  print( invertPart(l) );
3690}
3691//---------------------------------------------------------------------------------------
3692
3693static proc LRmul(list I, list J)
3694"USAGE:   LRmul(x, y);  x, y lists of integers (partitions)
3695RETURN:   list of lists
3696PURPOSE:  computes the partitions z for which the Littlewood-Richardson
3697          coefficient c^z_{x,y} is non-zero together with that coefficient;
3698          partitions up to length r
3699EXAMPLE:  example LRmul; shows an example
3700NOTE:     uses LRmult(..) from lrcalc.lib, does the same,
3701          only uses the inverted ordering of the elements in the partition
3702"
3703{
3704  list rez=LRmult(invertPart(I), invertPart(J));
3705  int sz=size(rez);
3706  int i;
3707  for(i=1;i<=sz;i++)
3708  {
3709    rez[i][2]=invertPart(rez[i][2]);
3710  }
3711  return(rez);
3712}
3713example
3714{
3715  "EXAMPLE:"; echo = 2;
3716  // Compute the partitions z for which the Littlewood-Richardson coefficient
3717  // c^z_{x,y} is non-zero together with that coefficient
3718  // for x= (1, 2), y=(1, 2)
3719  list x = 1, 2;
3720  list y = 1, 2;
3721  print( LRmul(x, y) );
3722}
3723//---------------------------------------------------------------------------------------
3724
3725static proc hook(list I)
3726"USAGE:   hook(I);  I list of integers (partition),
3727RETURN:   bigint
3728PURPOSE:  computes the product of the hook lenhths of the partition I
3729EXAMPLE:  example hook; shows an example
3730NOTE:
3731"
3732{
3733  bigint rez=1;
3734  list dI= invertPart( dualPart(I) );
3735  I=invertPart( I );
3736  int szI=size(I);
3737  int szdI=size(dI);
3738  int i, j;
3739  for(i=1;i<=szI;i++)
3740  {
3741    for(j=1;j<=I[i];j++)
3742    {
3743      rez=rez*(I[i]+dI[j]-i-j+1);
3744    }
3745  }
3746  return(rez);
3747}
3748example
3749{
3750  "EXAMPLE:"; echo = 2;
3751  // compute the product of all hook lengths of the partition (1, 1, 3)
3752  list I = 1, 1, 3;
3753  print( hook(I) );
3754}
3755//---------------------------------------------------------------------------------------
3756
3757static proc apn0_int(int N, list c)
3758"USAGE:   apn0_int( N, c);  N integer, c list of integers (partition)
3759RETURN:   list of integers
3760PURPOSE:  appends by integer zeroes up to the length N
3761EXAMPLE:  example apn0_int; shows an example
3762NOTE:
3763"
3764{
3765  int n=size(c);
3766  if(N>n) // if N is greater than the length of c, append c by zeroes up to the length N
3767  {
3768    int i;
3769    for(i=n+1;i<=N;i++)
3770    {
3771      c=c+list( int(0) );
3772    }
3773  }
3774  return(c);
3775}
3776example
3777{
3778  "EXAMPLE:";echo =2;
3779  ring r = 0, (x, y, z), dp;
3780  list l=(1, 2, 3);
3781  //append the list by two zeroes and get a list of length 5
3782  print( apn0_int(5, l) );
3783}
3784//---------------------------------------------------------------------------------------
3785
3786static proc contentPoly(list I, list J, poly e)
3787"USAGE:   contentPoly(I, J, e);  L, M lists of integers (partitions),
3788          e polynomial
3789RETURN:   poly
3790PURPOSE:  computes the content polynomial of the skew partition corresponding to I>J
3791EXAMPLE:  example contentPoly; shows an example
3792NOTE:
3793"
3794{
3795  int i, j;
3796  int szI=size(I);
3797  I=invertPart(I);
3798  J=invertPart(J);
3799  J=apn0_int(szI, J);
3800  poly rez=1;
3801  for(i=1; i<=szI; i++)
3802  {
3803    for(j=J[i]+1; j<=I[i]; j++)
3804    {
3805      rez = rez*(e- i+j);
3806    }
3807  }
3808  return(rez);
3809}
3810example
3811{
3812  "EXAMPLE:"; echo = 2;
3813  // compute the content Polynomial of the skew partition
3814  // corresponding to  (1,2) > (0, 1) with respect to the variable x
3815  ring r = 0, (x), dp;
3816  list L=1,2;
3817  list M=0,1;
3818  print( contentPoly(L, M, x) );
3819}
3820//---------------------------------------------------------------------------------------
3821
3822static proc Pcoef(list L, list M, poly e, poly f)
3823"USAGE:   Pcoef(L, M, e, f);  L, M lists of integers (partitions),
3824          e, f polynomials
3825RETURN:   poly
3826PURPOSE:  computes the polynomial P_{L, M}(e, f) from the paper of Manivel
3827EXAMPLE:  example Pcoef; shows an example
3828NOTE:
3829"
3830{
3831  list P = LRmul(dualPart(L), M);
3832  int  sz=size(P);
3833  poly rez;
3834  //poly h;
3835  list T, DT;
3836  //bigint lrc;
3837  int i;
3838  for(i=1;i<=sz;i++)
3839  {
3840    T=P[i][2];
3841    DT=dualPart(T);
3842    //lrc=P[i][1];
3843    //h=contentPolyM(DT, L, e)*contentPoly(T, M, f);
3844    rez=rez+P[i][1]* contentPoly(DT, L, e)*contentPoly(T, M, f)/hook(DT);
3845  }
3846 return(rez);
3847}
3848example
3849{
3850  "EXAMPLE:"; echo = 2;
3851  // compute P_{L, M}(e, f) from the paper of Manivel
3852  // for L = (0,1) and M = (1, 1)
3853  ring r = 0, (e, f), dp;
3854  list L=1,2,3;
3855  list M=1,2;
3856  print( Pcoef(L, M, e, f) );
3857}
3858//---------------------------------------------------------------------------------------
3859
3860static proc max_deg(ideal I)
3861"USAGE:   max_deg(I);  I an ideal
3862RETURN:   integer
3863PURPOSE:  computes the maximal degree of the generators of I
3864EXAMPLE:  example max_deg; shows an example
3865NOTE:
3866"
3867{
3868  int rez=0;
3869  int i;
3870  int sz = ncols(I);
3871  for(i=1;i<=sz;i++)
3872  {
3873    rez=max(rez, deg(I[i]) );
3874  }
3875  return(rez);
3876}
3877example
3878{
3879  "EXAMPLE:";echo =2;
3880  // the maximal degree of the ideal in k[x, y]:
3881  ring r = 0, (x, y), dp;
3882  ideal I= x4, y7, x2y3;
3883  print(max_deg(I));
3884}
3885//-------------------------------------------------------------------------------------------
3886
3887static proc var_pow(int n)
3888"USAGE:   var_pow(n);  n an integer
3889RETURN:   ideal
3890PURPOSE:  computes the ideal generated by the n-th powers of the variables of the base ring
3891EXAMPLE:  example var_pow; shows an example
3892NOTE:
3893"
3894{
3895  ideal I=maxideal(1);
3896  int sz=ncols(I);
3897  int i;
3898  ideal J;
3899  for(i=1; i<=sz; i++)
3900  {
3901    J=J+ideal(I[i]^n);
3902  }
3903  return(J);
3904}
3905example
3906{
3907  "EXAMPLE:";echo =2;
3908  // the 3-rd powers of the variables in k[x, y]:
3909  ring r = 0, (x, y), dp;
3910  print(var_pow(3));
3911}
3912//-------------------------------------------------------------------------------------------
3913
3914static proc equal_deg(ideal I)
3915"USAGE:   equal_deg(I);  I an ideal
3916RETURN:   ideal
3917PURPOSE:  computes an ideal generated by elements of the same degree
3918          that defines the same projective subscheme as I
3919EXAMPLE:  example equal_deg; shows an example
3920NOTE:
3921"
3922{
3923  I=simplify(I, 8+2);
3924  int sz=ncols(I);
3925  int mxd=max_deg(I);
3926  int i;
3927  ideal J;
3928  for(i=1;i<=sz;i++)
3929  {
3930    J=J+I[i]*var_pow( mxd-deg(I[i]) );
3931  }
3932
3933  return(sort( simplify(J, 8+2) )[1]);
3934}
3935example
3936{
3937  "EXAMPLE:"; echo=2;
3938  // change the ideal (x, y^2) in k[x, y, z]:
3939  ring r = 0, (x, y, z), dp;
3940  ideal I=x, y*z;
3941  // the ideal defines a two points subscheme in the projective plane
3942  // and is generated by elements of different degrees
3943  print(I);
3944  ideal J=equal_deg(I);
3945  // now the ideal is generated by elements of degree 2
3946  // and defines the same subscheme in the projective plane
3947  J;
3948  // notice that both ideals have the same saturation
3949  // with respect to the irrelevant ideal (x, y, z)
3950  // the saturation of the initial ideal coincides with the ideal itself
3951  sat(I, maxideal(1))[1];
3952  // the saturation of the modified ideal
3953  sat(J, maxideal(1))[1];
3954}
3955//-------------------------------------------------------------------------------------------
3956
3957static proc CSM_hypersurf(poly f)
3958"USAGE:   CSM_hypersurf(f);  f a polynomial
3959RETURN:   list of integers
3960PURPOSE:  computes the Chern-Schwartz-MacPherson classes of the hypersurface defined by f
3961EXAMPLE:  example CSM_hypersurf; shows an example
3962NOTE:
3963"
3964{
3965  ideal I=jacob(f);
3966  I=simplify(I, 8+2); // ignore repetitions and zero generators
3967  I=sort(I)[1]; // sort the generators, it speeds up the computations
3968  int degGen=deg(I[1]); // the degree of the generators of the Jacobian ideal
3969  int i;
3970  int n=nvars(basering)-1; // the dimension of the projective space
3971  bigintmat REZ[1][n+1];
3972  if(degGen==-1) // if the Jacobian ideal is zero
3973  {
3974    for(i=0;i<=n;i++)
3975    {
3976      REZ[1,i+1]=binomial(n+1,i);
3977    }
3978    return(REZ);
3979  }
3980  //TODO need to check the zero ideal and I==1;
3981  int sz=ncols(I);
3982  def br@=basering; // remember the base ring
3983  list l15 = ringlist(basering)[2];
3984  for (int ii = 1; ii <= sz; ii++)
3985  {
3986   l15[size(l15)+1] = "t@("+string(ii)+")";
3987  }
3988  l15[size(l15)+1] = "u@";
3989  ring r@ = create_ring(ringlist(basering)[1], l15, "dp", "no_minpoly");
3990  execute( "map F= br@,"+varstr(br@)+";" ); // define the corresponding inclusion of rings
3991  ideal I=F(I);
3992  ideal J(0..n);
3993  for(i=1; i<=sz; i++)
3994  {
3995    J(0)=J(0) + ideal( t@(i)-u@*I[i] );
3996  }
3997  J(0)=eliminate(J(0), u@);
3998  ideal T=t@(1..sz);
3999  for(i=1;i<=n;i++)
4000  {
4001    J(i) = sat( random_hypersurf(J(i-1), T), T )[1];
4002  }
4003  poly prd=product(T);
4004  poly cl;
4005  poly mlt;
4006  for(i=0;i<=n;i++)
4007  {
4008    cl=cl+  mult( std(eliminate( J(i),  prd)) ) *(-u@)^i*(1+u@)^(n-i);
4009  }
4010  cl=(1+u@)^(n+1)-cl;
4011  poly resPoly=NF( cl, u@^(n+1) );
4012  matrix cf=coeffs(resPoly, u@);
4013  for(i=0;i<=n;i++)
4014  {
4015    if( i < nrows(cf) )
4016    {
4017      REZ[1,i+1]=int(cf[i+1,1]);
4018    }
4019    else
4020    {
4021      REZ[1,i+1]=int(0);
4022    }
4023  }
4024  return(REZ);
4025}
4026example
4027{
4028  "EXAMPLE:";echo =2;
4029  // consider the projective plane with homogeneous coordinates x, y, z
4030  ring r = 0, (x, y, z), dp;
4031  // the Chern-Schwartz-MacPherson class of a smooth cubic:
4032  poly f=x3+y3+z3;
4033  f;
4034  CSM_hypersurf(f);
4035  // the Chern-Schwartz-MacPherson class of singular cubic
4036  // that is a union of 3 non-collinear lines:
4037  poly g=x*y*z;
4038  g;
4039  CSM_hypersurf(g);
4040  // the Chern-Schwartz-MacPherson class of singular cubic
4041  // that is a union of 3 lines passing through one point
4042  poly h=x*y*(x+y);
4043  h;
4044  CSM_hypersurf(h);
4045}
4046//-------------------------------------------------------------------------------------------
4047
4048static proc random_hypersurf(ideal I, ideal V)
4049"USAGE:   random_hypersurf(I, V);  I, V ideals
4050RETURN:   ideal
4051PURPOSE:  computes the sum of I with the ideal generated by a random
4052          linear combination of the generators of V such that the dimension decreases
4053EXAMPLE:  example random_hypersurf; shows an example
4054NOTE:     if the ideal I=1 (the whole ring), then I is returned
4055"
4056{
4057  ideal H;
4058  ideal J;
4059 // if(isSubModule(ideal(1), I)) // if I equals 1 (the whole ring);
4060  if( is_zero(I) )
4061  {
4062    return(I);
4063  }
4064  // otherwise
4065  int ok=0;
4066  int ntries; // number of tries
4067  while( !ok ) // give two tries for every b in randomid(V, 1, b)
4068  {
4069    H=randomid(V, 1, ntries div 2 +1);
4070    ntries++; // increase by 1
4071    if( isSubModule( quotient(I, ideal(H) ), I) )
4072    {
4073      J=I + H;
4074      ok=1;
4075    }
4076  }
4077  return(J);
4078}
4079example
4080{
4081  "EXAMPLE:";echo =2;
4082  // Consider an ideal in  k[x, y, z, s, t] and find its intersection with a general hyperplane
4083  // given by as+bt=0
4084  ring r = 0, (x, y, z, s, t), dp;
4085  ideal I=x2, yz, s+t;
4086  I;
4087  ideal V= s, t;
4088  V;
4089  // the ideal of the intersection with the random general hyperplane
4090  random_hypersurf(I, V);
4091}
4092//-------------------------------------------------------------------------------------------
4093
4094static proc prds(int n, def l)
4095"USAGE:   prds(n, l);  n an integer, l list of polynomials or ideal
4096RETURN:   list of polynomials
4097PURPOSE:  computes all possible products of length n (without repetitions) of the entries of l
4098EXAMPLE:  example prds; shows an example
4099NOTE:
4100"
4101{
4102  int sz;
4103  if(typeof(l)=="ideal")
4104  {
4105    sz=ncols(l);
4106  }
4107  else
4108  {
4109    sz=size(l);
4110  }
4111  if( (n>sz)||(sz==0)||(n<0) )
4112  {
4113    return( list() );
4114  }
4115  // otherwise
4116  if(n==0)
4117  {
4118    return( list(int(1)) );
4119  }
4120  if(sz==n)
4121  {
4122    return(product(l));
4123  }
4124  list L, LL, ll;
4125  ll=l[2..sz];
4126  poly f=l[1];
4127  L=prds(n, ll );
4128  LL=prds(n-1,ll);
4129  int i;
4130  sz=size(LL);
4131  for(i=1;i<=sz;i++)
4132  {
4133    LL[i]=f*LL[i];
4134  }
4135  return(L+LL);
4136}
4137example
4138{
4139  "EXAMPLE:";echo =2;
4140  ring r = 0, (x, y, z, w), dp;
4141  // compute all possible 2-products between the variables x,y,z,w
4142  list l=x,y,z,w;
4143  prds(2, l);
4144  // compute all possible 3-products between the variables x,y,z,w
4145  ideal I=x,y,z,w;
4146  prds(3, l);
4147}
4148//-------------------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.