Changeset 2e7410 in git


Ignore:
Timestamp:
Sep 24, 2010, 4:52:16 PM (13 years ago)
Author:
Stefan Steidel <steidel@…>
Branches:
(u'spielwiese', '828514cf6e480e4bafc26df99217bf2a1ed1ef45')
Children:
f0e080284bc1a4ba286ba66df791c785e808feb8
Parents:
1492bbce4e131a22de468dd41b9ced71d0fe3f93
Message:
new libs

git-svn-id: file:///usr/local/Singular/svn/trunk@13273 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
Singular/LIB
Files:
1 added
1 edited
1 moved

Legend:

Unmodified
Added
Removed

  • Singular/LIB/assprimeszerodim.lib

    • Property mode changed from 100644 to 100755
    r1492bb r2e7410  
    11////////////////////////////////////////////////////////////////////////////////
    22version="$Id$";
    3 category = "Commutative algebra";
     3category = "Commutative Algebra";
    44info="
    5 LIBRARY:  assPrimes.lib    associated primes of a zero-dimensional ideal
     5LIBRARY:  assprimeszerodim.lib   associated primes of a zero-dimensional ideal
    66
    77AUTHORS:  N. Idrees       nazeranjawwad@gmail.com
    … …  
    1111OVERVIEW:
    1212
    13   A library for computing the associated primes  and the radical of a
     13  A library for computing the associated primes and the radical of a
    1414  zero-dimensional ideal in the polynomial ring over the rational numbers,
    15   Q[x_1,...,x_n] using modular computations.
     15  Q[x_1,...,x_n], using modular computations.
    1616
    1717SEE ALSO: primdec_lib
    1818
    1919PROCEDURES:
     20 zeroR(I);             computes the radical of I
    2021 assPrimes(I);         computes the associated primes of I
    21  zeroR(I);             compute the radical of I
    2222";
    2323
    … …  
    3030"USAGE:  zeroR(I,[n]); I ideal, optional: n number of processors (for parallel computing)
    3131ASSUME:  I is zero-dimensional in Q[variables]
     32NOTE:    Parallelization is just applicable using 32-bit Singular version since
     33         MP-links are not compatible with 64-bit Singular version.
    3234RETURN:  the radical of I
    3335EXAMPLE: example zeroR; shows an example
    … …  
    6466
    6567//--------------------  Initialize the list of primes  -------------------------
    66 // nach Umstellung von modstd.lib auf unsere neue Version: primeList(I,10)
    67    intvec L = primeList(10);
     68   intvec L = primeList(I,10);
    6869   L[5] = prime(random(100000000,536870912));
    6970
    … …  
    141142      }
    142143
    143       //hier deleteUnlucky einbauen, streichen, wenn der Grad nicht stimmt
     144      // insert deleteUnluckyPrimes
    144145      G = chinrem(CO1,CO2);
    145146      N = CO2[1];
    … …  
    160161
    161162         j = size(L) + 1;
    162 
    163 // nach Umstellung von modstd.lib auf unsere neue Version: primeList(I,10,L)
    164          L = primeList(10,L);
     163         L = primeList(I,10,L);
    165164
    166165         if(n > 1)
    … …  
    198197
    199198proc assPrimes(list #)
    200 "USAGE:  assPrimes(I,[n],[a]); I ideal or module,
     199"USAGE:  assPrimes(I,[a],[n]); I ideal or module,
    201200         optional: n number of processors (for parallel computing), a
    202          - a=1: method of Eisenbud/Hunecke/Vasconcelos
    203          - a=2: method of Gianni/Trager/Zacharias
    204          - a=3: mathod of Monico
    205          assPrimes(I)  chooses n=a=1
     201         - a = 1: method of Eisenbud/Hunecke/Vasconcelos
     202         - a = 2: method of Gianni/Trager/Zacharias
     203         - a = 3: mathod of Monico
     204         assPrimes(I) chooses n = a = 1
    206205ASSUME:  I is zero-dimensional over Q[variables]
     206NOTE:    Parallelization is just applicable using 32-bit Singular version since
     207         MP-links are not compatible with 64-bit Singular version.
    207208RETURN:  a list pr of associated primes of I:
    208209EXAMPLE: example assPrimes; shows an example
    … …  
    214215
    215216   ideal I;
    216    if(typeof(#[1])=="ideal")
     217   if(typeof(#[1]) == "ideal")
    217218   {
    218219      I = #[1];
    … …  
    223224      I = Ann(M);
    224225   }
     226   
     227   //---------------------  Initialize optional parameter  ------------------------
     228   if(size(#) > 1)
     229   {
     230      if(size(#) == 2)
     231      {
     232         int alg = #[2];
     233         int n = 1;
     234      }
     235      if(size(#) == 3)
     236      {
     237         int alg = #[2];
     238         int n = #[3];
     239         if((n > 1) && (1 - system("with","MP")))
     240         {
     241            "========================================================================";
     242            "There is no MP available on your system. Since this is necessary to     ";
     243            "parallelize the algorithm, the computation will be done without forking.";
     244            "========================================================================";
     245            n = 1;
     246         }
     247      }
     248   }
     249   else
     250   {
     251      int alg = 1;
     252      int n = 1;
     253   }
    225254
    226255   def SPR = basering;
    227    I = std(I);
     256   I = modStd(I,n);
    228257   int d = vdim(I);
    229258   if(d == -1) { ERROR("Ideal is not zero-dimensional."); }
    … …  
    263292      "CPU-time for radical-check is "+string(timer - T)+" seconds.";
    264293   }
    265 //---------------------  Initialize optional parameter  ------------------------
    266    if(size(#) > 1)
    267    {
    268       if(size(#) == 2)
    269       {
    270          int alg = #[2];
    271          int n = 1;
    272       }
    273       if(size(#) == 3)
    274       {
    275          int alg = #[2];
    276          int n = #[3];
    277          if((n > 1) && (1 - system("with","MP")))
    278          {
    279             "========================================================================";
    280             "There is no MP available on your system. Since this is necessary to     ";
    281             "parallelize the algorithm, the computation will be done without forking.";
    282             "========================================================================";
    283             n = 1;
    284          }
    285       }
    286    }
    287    else
    288    {
    289       int alg = 1;
    290       int n = 1;
    291    }
    292294
    293295   export(SPR);
    … …  
    297299   export(I_for_fork);           // I available for each link
    298300
     301//--------------------  Initialize the list of primes  -------------------------
     302   intvec L = primeList(I,10);
     303   L[5] = prime(random(1000000000,2134567879));
     304
    299305   list Re;
    300306
    … …  
    308314   int j = 1;
    309315   int index = 1;
    310 
    311 
    312 //--------------------  Initialize the list of primes  -------------------------
    313 // nach Umstellung von modstd.lib auf unsere neue Version: primeList(I,10)
    314    intvec L = primeList(10);
    315    L[5] = prime(random(1000000000,2134567879));
    316316
    317317//-----  If there is more than one processor available, we parallelize the  ----
    … …  
    359359            for(i = 1; i <= n; i++)
    360360            {
    361                if(status(l(i), "read", "ready"))                         // ask if link l(i) is ready otherwise sleep for t seconds
     361               if(status(l(i), "read", "ready"))           // ask if link l(i) is ready otherwise sleep for t seconds
    362362               {
    363                   P = read(l(i));                                        // read the result from l(i)
     363                  P = read(l(i));                          // read the result from l(i)
    364364                  CO1[index] = P[1];
    365365                  CO2[index] = bigint(P[2]);
    … …  
    378378               }
    379379            }
    380             if(k == n)                                                  // k describes the number of closed links
     380            if(k == n)                                     // k describes the number of closed links
    381381            {
    382382               j++;
    … …  
    401401         "CPU-time for computing list in assPrimes is "+string(timer - tt)+" seconds.";
    402402      }
    403       //hier deleteUnlucky einbauen, streichen, wenn der Grad nicht stimmt
     403     
     404      // insert deleteUnluckyPrimes
    404405      G = chinrem(CO1,CO2);
    405406      N = CO2[1];
    … …  
    461462                  {
    462463                     H[1][i] = quickSubst(H[1][i],f,I);
    463 
    464                      if(typeof(#[1])=="ideal")
    465                      {
    466                         Re[i-1] = #[1] + ideal(H[1][i]);
    467                      }
    468                      else
    469                      {
    470                         Re[i-1] = M + H[1][i]*freemodule(nrows(M));
    471                      }
     464                     Re[i-1] = I + ideal(H[1][i]);
    472465                  }
    473466
    … …  
    514507      j = size(L) + 1;
    515508
    516 // nach Umstellung von modstd.lib auf unsere neue Version: primeList(I,10,L)
    517       L = primeList(10,L);
     509      setring(SPR);
     510      L = primeList(I,10,L);
     511      setring rHelp;
    518512
    519513      if(n > 1)
    … …  
    746740////////////////////////////////////////////////////////////////////////////////
    747741
    748 static proc primeList(int n, list #)
    749 "USAGE:  primeList(n); (resp. primeList(n,L); )
    750 RETURN: the intvec of n greatest primes  <= 536870912 (resp. n greatest
    751         primes < L[size(L)] union with L)
    752 EXAMPLE:example primList; shows an example
    753 "
    754 {
    755    intvec L;
    756    int i,p;
    757    if(size(#) == 0)
    758    {
    759       p = prime(536870912);
    760       L[1] = p;
    761    }
    762    else
    763    {
    764       L = #[1];
    765       p = prime(L[size(L)]-1);
    766       L[size(L)+1] = p;
    767    }
    768    if(p == 2){ERROR("no more primes");}
    769    for(i = 2; i <= n; i++)
    770    {
    771       p = prime(p-1);
    772       L[size(L)+1] = p;
    773    }
    774    return(L);
    775 }
    776 example
    777 { "EXAMPLE:"; echo = 2;
    778    intvec L=primeList(10);
    779    size(L);
    780    L[size(L)];
    781    L=primeList(5,L);
    782    size(L);
    783    L[size(L)];
    784 }
    785 
    786 ////////////////////////////////////////////////////////////////////////////////
    787 
    788742static proc zeroRadP(ideal I, int p)
    789743{
    … …  
    874828
    875829//Test for zeroR
    876 ring R = 0, (x,y), dp;
     830ring R = 0,(x,y),dp;
    877831ideal I = xy4-2xy2+x, x2-x, y4-2y2+1;
    878832
    879833//Cyclic 6
    880 ring R=0,x(1..6),dp;
    881 ideal I=cyclic(6);
     834ring R = 0,x(1..6),dp;
     835ideal I = cyclic(6);
    882836
    883837//Amrhei
    884 ring R=0,(a,b,c,d,e,f),dp;
    885 ideal I=
    886 2fb+2ec+d2+a2+a,
    887 2fc+2ed+2ba+b,
    888 2fd+e2+2ca+c+b2,
    889 2fe+2da+d+2cb,
    890 f2+2ea+e+2db+c2,
    891 2fa+f+2eb+2dc;
     838ring R = 0,(a,b,c,d,e,f),dp;
     839ideal I = 2fb+2ec+d2+a2+a,
     840          2fc+2ed+2ba+b,
     841          2fd+e2+2ca+c+b2,
     842          2fe+2da+d+2cb,
     843          f2+2ea+e+2db+c2,
     844          2fa+f+2eb+2dc;
    892845
    893846
    894847//Becker-Niermann
    895848ring R = 0,(x,y,z),dp;
    896 ideal I=
    897 x2+xy2z-2xy+y4+y2+z2,
    898 -x3y2+xy2z+xyz3-2xy+y4,
    899 -2x2y+xy4+yz4-3;
     849ideal I = x2+xy2z-2xy+y4+y2+z2,
     850          -x3y2+xy2z+xyz3-2xy+y4,
     851          -2x2y+xy4+yz4-3;
    900852
    901853//Roczen
    902 ring R=0,(a,b,c,d,e,f,g,h,k,o),dp;
    903 ideal I=
    904 o+1,
    905 k4+k,
    906 hk,
    907 h4+h,
    908 gk,
    909 gh,
    910 g3+h3+k3+1,
    911 fk,
    912 f4+f,
    913 eh,
    914 ef,
    915 f3h3+e3k3+e3+f3+h3+k3+1,
    916 e3g+f3g+g,
    917 e4+e,
    918 dh3+dk3+d,
    919 dg,
    920 df,
    921 de,
    922 d3+e3+f3+1,
    923 e2g2+d2h2+c,
    924 f2g2+d2k2+b,
    925 f2h2+e2k2+a;
     854ring R = 0,(a,b,c,d,e,f,g,h,k,o),dp;
     855ideal I = o+1, k4+k, hk, h4+h, gk, gh, g3+h3+k3+1,
     856          fk, f4+f, eh, ef, f3h3+e3k3+e3+f3+h3+k3+1,
     857          e3g+f3g+g, e4+e, dh3+dk3+d, dg, df, de,
     858          d3+e3+f3+1, e2g2+d2h2+c, f2g2+d2k2+b,
     859          f2h2+e2k2+a;
    926860
    927861//4 bodies with equal masses, before symmetrisation.
    928862//We are looking for the real positive solutions
    929 ring R=0,(B,D,F,b,d,f),dp;
    930 
    931 ring R=0,(B,b,D,d,F,f),dp;
    932 ideal I=(b-d)*(B-D)-2*F+2,
    933 (b-d)*(B+D-2*F)+2*(B-D),
    934 (b-d)^2-2*(b+d)+f+1,
    935 B^2*b^3-1,D^2*d^3-1,F^2*f^3-1;
    936 
    937 ring R=0,(B,b,D,d,F,f),dp;
    938 ideal I=(b-d)*(B-D)-F+1,
    939 (b-d)*(B+D-F)+(B-D),
    940 (b-d)^2-(b+d)+f+1,
    941 B^2*b^3-1,D^2*d^3-1,F^2*f^3-1;
     863ring R = 0,(B,b,D,d,F,f),dp;
     864ideal I = (b-d)*(B-D)-2*F+2,
     865          (b-d)*(B+D-2*F)+2*(B-D),
     866          (b-d)^2-2*(b+d)+f+1,
     867          B^2*b^3-1,D^2*d^3-1,F^2*f^3-1;
     868
     869ring R = 0,(B,b,D,d,F,f),dp;
     870ideal I = (b-d)*(B-D)-F+1,
     871          (b-d)*(B+D-F)+(B-D),
     872          (b-d)^2-(b+d)+f+1,
     873          B^2*b^3-1,D^2*d^3-1,F^2*f^3-1;
    942874
    943875
    944876//prime
    945 ring R=0,(x,y,z,t,u),dp;
    946 
    947 ideal I=2x2-2y2+2z2-2t2+2u2-1,
     877ring R = 0,(x,y,z,t,u),dp;
     878ideal I = 2x2-2y2+2z2-2t2+2u2-1,
    948879          2x3-2y3+2z3-2t3+2u3-1,
    949          2x4-2y4+2z4-2t4+2u4-1,
    950          2x5-2y5+2z5-2t5+2u5-1,
    951          2x6-2y6+2z6-2t6+2u6-1;
     880          2x4-2y4+2z4-2t4+2u4-1,
     881          2x5-2y5+2z5-2t5+2u5-1,
     882          2x6-2y6+2z6-2t6+2u6-1;
    952883
    953884*/

  • Singular/LIB/surfacesignature.lib

    r1492bb r2e7410  
    22category="Singularities";
    33info="
    4 LIBRARY:  surfacesignature.lib        signature of surface singularity
    5 AUTHORS:  Gerhard Pfister,            pfister@mathematik.uni-kl.de
    6           Muhammad Ahsan Banyamin,    ahsanbanyamin@gmail.com
    7 
    8 OVERVIEW:
    9  A library for computing the signature of irreducible surface singularity.
    10  It contains also a procedure for computing the newton pairs of a surface
    11  singularity.
    12 
    13 THEORY: The signature of a surface singularity is defined in
    14         [Durfee,A.:The Signature of Smoothings of Complex Surface
    15         Singularities, Math. Ann.,232,85-98 (1978)].
    16         The algorithm we use has been proposed in
    17         [Nemethi,A.:The signature of f(x,y)+z^n, Proceedings of Singularity
    18         Conference (C.T.C Wall's 60th birthday meeting),Liverpool 1996,
    19         London Math.Soc. LN 263(1999),131-149].
     4LIBRARY:  signaturesurfsing.lib       signature of surface singularity
     5
     6AUTHORS:  Gerhard Pfister             pfister@mathematik.uni-kl.de
     7@*        Muhammad Ahsan Banyamin     ahsanbanyamin@gmail.com
     8@*        Stefan Steidel              steidel@mathematik.uni-kl.de
     9
     10OVERVIEW:
     11
     12  A library for computing the signature of irreducible surface singularity.
     13  The signature of a surface singularity is defined in [Durfee, A.: The
     14  Signature of Smoothings of Complex Surface Singularities, Math. Ann., 232,
     15  85-98 (1978)]. The algorithm we use has been proposed in [Nemethi, A.: The
     16  signature of f(x,y)+z^N, Proceedings of Singularity Conference (C.T.C. Wall's
     17  60th birthday meeting), Liverpool 1996, London Math.Soc. LN 263(1999),
     18  131-149].
    2019
    2120PROCEDURES:
    22  brieskornSign(a,b,c);   signature of Brieskorn singularity x^a+y^b+z^c=0
    23  newtonpairs(f);         newton pairs of surface singularity
    24  signature(N,f);         signature of singularity z^+f(x,y)=0,f irreducible
     21 brieskornSign(a1,a2,a3);  signature of Brieskorn singularity x^a1+y^a2+z^a3
     22 signature(N,f);           signature of singularity z^N+f(x,y)=0, f irreducible
    2523";
    2624
    2725LIB "hnoether.lib";
    28 ///////////////////////////////////////////////////////////////////////////////
    29 proc signature(int N,poly f)
     26LIB "alexpoly.lib";
     27LIB "gmssing.lib";
     28
     29///////////////////////////////////////////////////////////////////////////////
     30//------- sigma(z^N + f) in terms of Puiseux pairs of f for f irreducible ----
     31
     32static proc exponentSequence(poly f)
     33//=== computes the sequence a_1,...,a_s of exponents as described in [Nemethi]
     34//=== using the Puiseux pairs (m_1, n_1),...,(m_s, n_s) of f:
     35//===  - a_1 = m_1,
     36//===  - a_i = m_i - n_i * (m_[i-1] - n_[i-1] * a_[i-1]).
     37//===
     38//=== Return: list of two intvecs:
     39//===         1st entry: A = (a_1,...,a_s)
     40//===         2nd entry: N = (n_1,...,n_s)
     41{
     42   def R = basering;
     43   ring S = 0,(x,y),dp;
     44   poly f = fetch(R,f);
     45   list puiseuxPairs = invariants(f);
     46   setring R;
     47   
     48   intvec M = puiseuxPairs[1][3];
     49   intvec N = puiseuxPairs[1][4];
     50   
     51   int i;
     52   int a = M[1];
     53   intvec A = a;
     54   for(i = 2; i <= size(M); i++)
     55   {
     56      a = M[i] - N[i] * (M[i-1] - N[i-1] * a);
     57      A[size(A)+1] = a;
     58   }
     59 
     60   return(list(A,N));
     61}
     62example
     63{ "EXAMPLE:"; echo = 2;
     64   ring r = 0,(x,y),dp;
     65   exponentSequence(y4+2x3y2+x6+x5y);
     66}
     67
     68///////////////////////////////////////////////////////////////////////////////
     69
     70proc brieskornSign(a1,a2,a3)
     71"USAGE:  brieskornSign(a1,a2,a3); a1,a2,a3 = integers
     72RETURN:  signature of Brieskorn singularity x^a1+y^a2+z^a3
     73EXAMPLE: example brieskornSign; shows an example
    3074"
    31 USAGE:   signature(N,f); N=integer, f=irreducible poly in 2 variables
    32 RETURN:  signature of surface singularity defined by z^N+f=0
    33 EXAMPLE: example signature; shows an example
     75{
     76   int a_temp, t, k1, k2, k3, s_t, sigma;
     77   number s;
     78   
     79   if(a1 > a2) { a_temp = a1; a1 = a2; a2 = a_temp; }
     80   if(a2 > a3) { a_temp = a2; a2 = a3; a3 = a_temp; }
     81   if(a1 > a2) { a_temp = a1; a1 = a2; a2 = a_temp; }
     82   
     83   for(t = 0; t <= 2; t++)
     84   {
     85      s_t = 0;
     86      for(k1 = 1; k1 <= a1-1; k1++)
     87      {
     88         for(k2 = 1; k2 <= a2-1; k2++)
     89         {
     90            for(k3 = 1; k3 <= a3-1; k3++)
     91            {
     92               s = number(k1)/a1 + number(k2)/a2 + number(k3)/a3;
     93               if(t < s)
     94               {
     95                  if(s < t+1)
     96                  {
     97                     s_t = s_t + 1;
     98                  }
     99                  else
     100                  {
     101                     break;
     102                  }
     103               }
     104            }
     105            if(k3 == 1) { break; }
     106         }
     107         if(k2 == 1) { break; }
     108      }
     109      sigma = sigma + (-1)^t * s_t;
     110   }
     111   return(sigma);
     112}
     113example
     114{ "EXAMPLE:"; echo = 2;
     115   ring R = 0,x,dp;
     116   brieskornSign(11,3,5);
     117}
     118
     119///////////////////////////////////////////////////////////////////////////////
     120
     121static proc signatureP(int N,poly f)
     122"USAGE:  signatureP(N,f); N = integer, f = irreducible poly in 2 variables
     123RETURN:  signature of surface singularity defined by z^N + f(x,y) = 0
     124EXAMPLE: example signatureP; shows an example
    34125"
    35126{
    36    def R=basering;
    37    ring S=0,(x,y),dp;
    38    poly f = fetch(R,f);
    39    list L=newtonpairs(f);
    40    setring R;
    41    if(size(L)>2)
    42    {
    43      return(Generalcase(N,L));
    44    }
    45    if(size(L)==2)
    46    {
    47      return(signtwopairs(N,L));
    48    }
    49    return(brieskornSign(L[1][2],L[1][1],N));
    50  }
    51  example
     127   int i, d, prod, sigma;
     128   list L = exponentSequence(f);
     129   int s = size(L[2]);
     130   
     131   if(s == 1)
     132   {
     133      return(brieskornSign(L[1][1], L[2][1], N));
     134   }
     135   
     136   prod = 1;
     137   sigma = brieskornSign(L[1][s], L[2][s], N);
     138   for(i = s - 1; i >= 1; i--)
     139   {
     140      prod = prod * L[2][i+1];
     141      d = gcd(N, prod);
     142      sigma = sigma + d * brieskornSign(L[1][i], L[2][i], N/d);
     143   }
     144   
     145   return(sigma);
     146}
     147example
    52148{ "EXAMPLE:"; echo = 2;
    53149   ring r = 0,(x,y),dp;
    54150   int N  = 3;
    55    poly f= x15-21x14+8x13y-6x13-16x12y+20x11y2-x12+8x11y-36x10y2
    56       +24x9y3+4x9y2-16x8y3+26x7y4-6x6y4+8x5y5+4x3y6-y8;
    57    signature(N,f);
    58 }
    59 
    60 ///////////////////////////////////////////////////////////////////////////
    61 proc newtonpairs(poly f)
    62 "USAGE:   newtonpairs(f); f= irreducible bivariate poly
    63  RETURN:  newton pairs of curve singularity f(x,y)=0
    64  EXAMPLE: example newtonpairs; shows an example
    65 "
    66 {
    67  def R=basering;
    68  ring S=0,(x,y),dp;
    69  poly f = fetch(R,f);
    70  list M=invariants(f);
    71  setring R;
    72  list L=M[1][3];
    73  list K=M[1][4];
    74  int i,j;
    75  intvec v;
    76  list N1,N2,P,T;
    77  N1[1]=M[1][4][1];
    78  N2[1]=M[1][3][1];
    79  for(i=2;i<=size(M[1][3]);i++)
    80  {
    81   N1[i]=M[1][4][i];
    82   N2[i]=M[1][3][i]-N1[i]*M[1][3][i-1];
    83  }
    84   P[1]=N1;
    85   P[2]=N2;
    86   for(j=1;j<=size(P[1]);j++)
    87   {
    88      v=P[1][j],P[2][j];
    89      T[j]=v;
    90    }
    91   return(T);
    92 }
    93 example
    94 { "EXAMPLE:"; echo = 2;
    95    ring r=0,(x,y),dp;
    96    newtonpairs(y4+2x3y2+x6+x5y);
    97 }
    98 
    99 ///////////////////////////////////////////////////////////////////////////
    100 proc brieskornSign(a,b,c)
    101 "USAGE:   brieskornSign(a,b,c);a,b,c=integers
    102  RETURN:  signature of Brieskorn singularity x^a+y^b+z^c
    103  EXAMPLE: example brieskornSign; shows an example
    104 "
    105 {
    106    int i,j,k,d;
    107    for(i=1;i<=a-1;i++)
    108    {
    109        for(j=1;j<=b-1;j++)
    110        {
    111           for(k=1;k<=c-1;k++)
    112          {
    113             d=d+signa(number(i)/a+number(j)/b+number(k)/c);
    114          }
    115        }
    116     }
    117     return(d);
    118  }
    119 example
    120 { "EXAMPLE:"; echo = 2;
    121    ring R=0,x,dp;
    122    brieskornSign(11,3,5);
    123 }
    124 ///////////////////////////////////////////////////////////////////////////
    125 static proc signsin(number n)
    126 "USAGE:   signsin(n); n=number
    127  RETURN:  the sign of sin(n) (sin(n) = the value of the sine of n)
    128  EXAMPLE: example signSin; shows an example
    129 "
    130 {
    131    def r=basering;
    132    int a;
    133    int b=10;
    134    ring s=(real,10),x,dp;
    135    number m=imap(r,n);
    136    number pi = number_pi(11);
    137    number t=m*pi;
    138    poly f=sin(b);
    139    poly h=subst(f,x,t);
    140    if(h>0){a=1;}
    141    if(h<0){a=-1;}
    142    setring r;
    143    return(a);
    144 }
    145 example
    146 { "EXAMPLE:"; echo = 2;
    147    ring r=0,x,dp;
    148    signsin(11/3);
    149 }
    150 ///////////////////////////////////////////////////////////////////////////
    151 static proc split1(number n)
    152 "USAGE:   split1(n); n=number
    153  RETURN:  integral and fractional parts of number n
    154  EXAMPLE: example split1; shows an example
    155  "
    156 {
    157    number a,b;
    158    int r;
    159    a=numerator(n);
    160    b=denominator(n);
    161    int z=int(number(a));
    162    int y=int(number(b));
    163    r=z mod y;
    164    int q=(z-r) div y;
    165    number n1=q;
    166    number n2=n-n1;
    167    list l=n1,n2;
    168    return(l);
    169 }
    170 example
    171 { "EXAMPLE:"; echo = 2;
    172      ring r=0,x,dp;
    173      split1(11/3);
    174 }
    175 ///////////////////////////////////////////////////////////////////////////
    176 static proc sin( int n)
    177 "USAGE:   sin(n); n=integer
    178  RETURN:  approximate value of the sine of n by using maclaurin series
    179  EXAMPLE:  example sin; shows an example
    180  "
    181 {
    182    def R=basering;
    183    ring S=0,x,dp;
    184    int i;
    185    poly f;
    186    int z;
    187    for(i=1;i<=n;i=i+2)
    188    {
    189      f=f+(-1)^z*x^i/factorial(i) ;
    190      z++;
    191    }
    192    setring R;
    193    map phi=S,var(1);
    194    poly f=phi(f);
    195    return(f);
    196 }
    197 example
    198 { "EXAMPLE:"; echo = 2;
    199    ring r=0,x,dp;
    200    sin(10);
    201 }
    202 ///////////////////////////////////////////////////////////////////////////
    203 static proc cos(int n)
    204 "USAGE:   cos(n); n=integer
    205  RETURN:  approximate value of the cosine of n by using maclaurin series
    206  EXAMPLE: example cos; shows an example
    207  "
    208 {
    209    def R=basering;
    210    ring S=0,x,dp;
    211    poly f;
    212    int i;
    213    int z;
    214    for(i=0;i<=n;i=i+2)
    215    {
    216      // f=f+(-1)^z*x^i/factorial(i,0);
    217      f=f+(-1)^z*x^i/factorial(i) ;
    218      z++;
    219    }
    220    setring R;
    221    map phi=S,var(1);
    222    poly f=phi(f);
    223    return(f);
    224  }
    225 example
    226 { "EXAMPLE:"; echo = 2;
    227    ring r=0,x,dp;
    228    cos(10);
    229 }
    230 ///////////////////////////////////////////////////////////////////////////
    231 static proc signcos(number n)
    232 "USAGE:   signcos(n); n=number
    233  RETURN:  the sign of cosin(n) (cosin(n) = the value of the cosine of n)
    234  EXAMPLE: example signcos; shows an example
    235 "
    236 {
    237    def r=basering;
    238    int a;
    239    int b=10;
    240    ring s=(real,10),x,dp;
    241    number m=imap(r,n);
    242    number pi = number_pi(11);
    243    number t=m*pi;
    244    poly f=cos(b);
    245    poly h=subst(f,x,t);
    246    if(h>0){a=1;}
    247    if(h<0){a=-1;}
    248    setring r;
    249    return(a);
    250 }
    251 example
    252 { "EXAMPLE:"; echo = 2;
    253    ring r=0,x,dp;
    254    signcos(11/3);
    255 }
    256 ///////////////////////////////////////////////////////////////////////////
    257 static proc  signa( number u)
    258 "USAGE:   signa(u); u=number
    259  RETURN:  the signa of a number
    260  EXAMPLE: example signa; shows an example
    261 "
    262 {
    263    list l=split1(u);
    264    int z;
    265    if( l[1] mod 2==0 )
    266    { z=signsin(l[2]); }
    267    else
    268    { z=signcos(l[1]);}
    269    return(z);
    270 }
    271 example
    272 { "EXAMPLE:"; echo = 2;
    273    ring r=0,x,dp;
    274    signa(11/3);
    275 }
    276 ///////////////////////////////////////////////////////////////////////////
    277 static proc prods(list L)
    278 "USAGE:   product(L); L=list of intvec
    279  RETURN:  product of first components of Newton pairs in L
    280  EXAMPLE: example product; shows an example
    281 "
    282 {
    283    int a=L[2][1];
    284    int i;
    285    for(i=1;i<=size(L)-2;i++)
    286    {
    287       a=a*L[i+2][1];
    288    }
    289    return(a);
    290 }
    291 example
    292 { "EXAMPLE:"; echo = 2;
    293    list L=intvec(2,3),intvec(2,1);
    294    prods(L);
    295 }
    296 ///////////////////////////////////////////////////////////////////////////
    297 static proc Generalcase(int N, list L)
    298 "USAGE:   Generalcase(N,f);N=integer,list L of intvec
    299  RETURN:  signature of surface singularity with Newton pairs in L
    300  ASSUME:  number of newton pairs greater than 2
    301  EXAMPLE: example Generalcase; shows an example
    302 "
    303 {
    304    int i,j,k,n,m,t,p;
    305    int a=L[1][2];
    306    int b=prods(L);
    307    int d=gcd(N,b);
    308    int q=d*brieskornSign(a,L[1][1],N/d);
    309    list A;
    310    list B;
    311    list S;
    312    for(i=1;i<=size(L)-1;i++)
    313    {
    314      a=L[i+1][2]+L[i+1][1]*L[i][1]*a;
    315      A[i]=a;
    316      S[i]=L[i+1][1];
    317    }
    318    for(m=1;m<=size(L)-2;m++)
    319    {
    320      B[m]=L[m+2][1];
    321    }
    322    list C;
    323    int c=B[size(B)];
    324    C[size(B)]=c;
    325    for(j=1;j<=size(B)-1;j++)
    326    {
    327       c=c*B[size(B)-j];
    328       C[size(B)-j]=c;
    329    }
    330    list D;
    331    D[size(L)-1]=1;
    332    for(k=1;k<=size(L)-2;k++)
    333    {
    334       D[k]=gcd(N,C[k]);
    335    }
    336    for(n=1;n<=size(L)-1;n++)
    337    {
    338      t=t+D[n]*brieskornSign(A[n],S[n],N/D[n]);
    339    }
    340    return(q+t);
    341 }
    342 example
    343 { "EXAMPLE:"; echo = 2;
    344    int N=2;
    345    list L=intvec(2,3),intvec(2,1),intvec(2,1);
    346    Generalcase(N,L);
    347 }
    348 ///////////////////////////////////////////////////////////////////////////
    349 static proc signtwopairs(int N,list L)
    350 "USAGE:   signtwopairs(N,f);N=integer,L=list of intvec
    351  RETURN:  signature of surface singularity with Newton pairs in L
    352  ASSUME:  number of newton pairs equal to 2
    353  EXAMPLE: example signtwopairs; shows an example
    354 "
    355 {
    356   int a1,a2,b,d1,q,t;
    357   a1=L[1][2];
    358   b=prods(L);
    359   d1=gcd(N,b);
    360   q=d1*brieskornSign(a1,L[1][1],N/d1);
    361   a2=L[2][2]+L[2][1]*L[1][1]*a1;
    362   t=brieskornSign(a2,L[2][1],N);
    363   return(q+t);
    364 }
    365 example
    366 { "EXAMPLE:"; echo = 2;
    367    int N=2;
    368    list L=intvec(2,3),intvec(2,1);
    369    signtwopairs(N,L);
    370 }
    371 ///////////////////////////////////////////////////////////////////////////
    372 static proc DedekindSum(number b, number c, int a)
     151   poly f = x15-21x14+8x13y-6x13-16x12y+20x11y2-x12+8x11y-36x10y2
     152            +24x9y3+4x9y2-16x8y3+26x7y4-6x6y4+8x5y5+4x3y6-y8;
     153   signatureP(N,f);
     154}
     155
     156///////////////////////////////////////////////////////////////////////////////
     157//------- sigma(z^N + f) in terms of the imbedded resolution graph of f -------
     158
     159static proc dedekindSum(number b, number c, int a)
    373160{
    374161   number s,d,e;
    … …  
    385172   return(s);
    386173}
    387 ///////////////////////////////////////////////////////////////////////////
     174
     175///////////////////////////////////////////////////////////////////////////////
     176
     177static proc isRupture(intvec v)
     178//=== decides whether the exceptional divisor given by the row v in the
     179//=== incidence matrix of the resolution graph intersects at least 3 other divisors
     180{
     181   int i,j;
     182   for(i=1;i<=size(v);i++)
     183   {
     184       if(v[i]<0){return(0);}
     185       if(v[i]!=0){j++;}
     186   }
     187   return(j>=4);
     188}
     189
     190///////////////////////////////////////////////////////////////////////////////
     191
     192static proc sumExcepDiv(intmat N, list M, int K, int n)
     193//=== computes part of the formulae for eta(g,K), g defining an
     194//=== isolated curve singularity
     195//=== N the incidence matrix of the resolution graph of g
     196//=== M list of total multiplicities
     197//=== n = nrows(N)
     198{
     199   int i,j,m,d;
     200   for(i=1;i<=n;i++)
     201   {
     202      if(N[i,i]>0)
     203      {
     204         m=gcd(K,M[i]);
     205         for(j=1;j<=n;j++)
     206         {
     207            if((i!=j)&&(N[i,j]!=0))
     208            {
     209               if(m==1){break;}
     210               m=gcd(m,M[j]);
     211            }
     212         }
     213         d=d+m-1;
     214      }
     215   }
     216   return(d);
     217}
     218
     219///////////////////////////////////////////////////////////////////////////////
     220
     221static proc sumEdges(intmat N, list M, int K, int n)
     222//=== computes part of the formulae for eta(g,K), g defining an
     223//=== isolated curve singularity
     224//=== N the incidence matrix of the resolution graph of g
     225//=== M list of total multiplicities
     226//=== n = nrows(N)
     227{
     228   int i,j,d;
     229   for(i=1;i<=n-1;i++)
     230   {
     231      for(j=i+1;j<=n;j++)
     232      {
     233         if(N[i,j]==1)
     234         {
     235            d=d+gcd(K,gcd(M[i],M[j]))-1;
     236         }
     237      }
     238   }
     239   return(d);
     240}
     241
     242///////////////////////////////////////////////////////////////////////////////
     243
     244static proc etaRes(list L, int K)
     245//=== L total multiplicities
     246//=== eta-invariant in terms of the imbedded resolution graph of f
     247{
     248   int i,j,d;
     249   intvec v;
     250   number e;
     251   intmat N = L[1];         // incidence matrix of the resolution graph
     252   int n = ncols(L[1]);     // number of vertices in the resolution graph
     253   int a = ncols(L[2]);     // number of branches
     254   list M;                  // total multiplicities
     255   for(i=1;i<=n;i++)
     256   {
     257      d=L[2][i,1];
     258      for(j=2;j<=a;j++)
     259      {
     260         d=d+L[2][i,j];
     261      }
     262      if(d==0){d=1;}
     263      M[i]=d;
     264   }
     265   for(i=1;i<=n;i++)
     266   {
     267      v=N[i,1..n];
     268      if(isRupture(v))    // the divisor intersects more then two others
     269      {
     270         for(j=1;j<=n;j++)
     271         {
     272            if((i!=j)&&(v[j]!=0))
     273            {
     274               e=e+dedekindSum(M[j],K,M[i]);
     275            }
     276         }
     277      }   
     278   }
     279   if(a==1)
     280   {
     281      //the irreducible case
     282      return(4*e);
     283   }
     284   return(a-1+4*e+sumEdges(N,M,K,n)-sumExcepDiv(N,M,K,n));
     285}
     286
     287///////////////////////////////////////////////////////////////////////////////
     288//------------ sigma(z^N + f) in terms of the spectral pairs of f -------------
     289
     290static proc fracPart(number n)
     291//=== computes the fractional part n2 of n
     292//=== i.e. n2 is not in Z but n-n2 is in Z
     293{
     294   number a,b;
     295   int r;
     296   a = numerator(n);
     297   b = denominator(n);
     298   int z = int(number(a));
     299   int y = int(number(b));
     300   r = z mod y;
     301   int q = (z-r) div y;
     302   number n1 = q;
     303   number n2 = n-n1;
     304   return(n2);
     305}
     306
     307///////////////////////////////////////////////////////////////////////////////
     308
     309static proc etaSpec(list L, int N)
     310//=== L spectral numbers
     311//=== eta-invariant in terms of the spectral pairs of f
     312{
     313   int i;
     314   number e, h;
     315   
     316   int n = ncols(L[1]);
     317   
     318   if((n mod 2) == 0)               
     319   // 0 is not a spectral number, thus f is irreducible
     320   {
     321      for(i = n/2+1; i <= n; i++)
     322      {
     323         e = e + (1 - 2 * fracPart(N * number(L[1][i]))) * L[3][i];
     324      }
     325      return(2*e);
     326   }
     327   else
     328   // 0 is a spectral number, thus f is reducible
     329   {
     330      // sum of Hodge numbers in eta function
     331      for(i = 1; i <= n; i++)                   
     332      {
     333         if((L[2][i] == 2)&&((denominator(leadcoef(N*L[1][i]))==1)||(denominator(leadcoef(N*L[1][i]))==-1)))
     334         {
     335            h = h + L[3][i];
     336         }
     337      }
     338     
     339      // summand coming from spectral number 0 in eta function
     340      h = h + L[3][(n+1)/2];
     341     
     342      // sum coming from non-zero spectral numbers in eta function
     343      for(i = (n+3)/2; i <= n; i++)
     344      {
     345         if(!((denominator(leadcoef(N*L[1][i]))==1)||(denominator(leadcoef(N*L[1][i]))==-1)))
     346         {
     347            e = e + (1 - 2 * fracPart(N * number(L[1][i]))) * L[3][i];
     348         }
     349      }   
     350      return(h + 2*e);
     351   }
     352}
     353
     354///////////////////////////////////////////////////////////////////////////////
     355//---------------- Consolidation of the three former variants -----------------
     356
     357proc signature(int N, poly f, list #)
     358"USAGE:  signature(N,f); N = integer, f = reduced poly in 2 variables, # empty or 1,2,3
     359@*       - if # is empty or #[1] = 2 then resolution of singularities is used
     360@*       - if #[1] = 1 then f has to be analytically irreducible and Puiseux expansions are used
     361@*       - if #[1] = 3 then spectral pairs are used
     362RETURN:  signature of surface singularity defined by z^N + f(x,y) = 0
     363EXAMPLE: example signature; shows an example
     364"
     365{
     366   if(size(#) == 0)
     367   {
     368      list L = totalmultiplicities(f);
     369      return(etaRes(L,N) - N*etaRes(L,1));
     370   }
     371   
     372   if(#[1] == 1)
     373   {
     374      return(signatureP(N,f));
     375   }
     376   
     377   if(#[1] == 2)
     378   {
     379      list L = totalmultiplicities(f);
     380      return(etaRes(L,N) - N*etaRes(L,1));
     381   }
     382   
     383   if(#[1] == 3)
     384   {
     385      def R = basering;
     386      def Rds = changeord("ds");
     387      setring Rds;
     388      poly f = imap(R,f);
     389      list L = sppairs(f);
     390      setring R;
     391      list L = imap(Rds,L);
     392      return(etaSpec(L,N) - N*etaSpec(L,1));
     393   }
     394}
     395example
     396{ "EXAMPLE:"; echo = 2;
     397   ring r = 0,(x,y),dp;
     398   int N  = 3;
     399   poly f = x15-21x14+8x13y-6x13-16x12y+20x11y2-x12+8x11y-36x10y2
     400            +24x9y3+4x9y2-16x8y3+26x7y4-6x6y4+8x5y5+4x3y6-y8;
     401   signature(N,f,1);
     402   signature(N,f,2);
     403}
     404
     405///////////////////////////////////////////////////////////////////////////////
    388406
    389407/*
    390408Further examples
    391409
    392    ring r = 0,(x,y),dp;
    393    int N;
    394    poly f;
    395 
    396    N  = 5;
    397    poly f= x15-21x14+8x13y-6x13-16x12y+20x11y2-x12+8x11y-36x10y2    //3 characteristic pairs
    398       +24x9y3+4x9y2-16x8y3+26x7y4-6x6y4+8x5y5+4x3y6-y8;
    399 
    400    N=6;
    401    f= y4+2x3y2+x6+x5y;                                             //2 characteristic pairs
    402 
    403    N=7;
    404    f=x5+y11;                                                       //1 characteristc pair
     410ring r = 0,(x,y),dp;
     411int N;
     412poly f,g,g1,g2,g3;
     413
     414
     415// irreducible polynomials
     416
     417N = 5;
     418f = x15-21x14+8x13y-6x13-16x12y+20x11y2-x12+8x11y-36x10y2   
     419    +24x9y3+4x9y2-16x8y3+26x7y4-6x6y4+8x5y5+4x3y6-y8;     
     420g = f^3 + x17y17;                                               
     421
     422N = 6;
     423f = y4+2x3y2+x6+x5y;                                             
     424g1 = f^2 + x5y5;
     425g2 = f^3 + x11y11;
     426g3 = f^3 + x17y17;
     427
     428N = 7;
     429f = x5+y11;
     430g1 = f^3 + x11y11;
     431g2 = f^3 + x17y17;
     432
     433N = 6;     
     434// k0 = 30, k1 = 35, k2 = 71
     435f = x71+6x65+15x59-630x52y6+20x53+6230x46y6+910x39y12+15x47
     436    -7530x40y6+14955x33y12-285x26y18+6x41+1230x34y6+4680x27y12
     437    +1830x20y18+30x13y24+x35-5x28y6+10x21y12-10x14y18+5x7y24-y30;
     438
     439// k0 = 16, k1 = 24, k2 = 28, k3 = 30, k4 = 31
     440f = x31-781x30+16x29y-3010x29-2464x28y+104x27y2-2805x28-7024x27y
     441    -5352x26y2+368x25y3+366x27-7136x26y-984x25y2-8000x24y3
     442    +836x23y4+34x26-320x25y-6464x24y2+6560x23y3-8812x22y4+1392x21y5
     443    -12x25+256x24y-1296x23y2-1536x22y3+4416x21y4-8864x20y5+1752x19y6
     444    -x24+16x23y-88x22y2-16x21y3-404x20y4+3056x19y5-6872x18y6+1648x17y7
     445    +8x21y2-96x20y3+524x19y4-1472x18y5+3464x17y6-3808x16y7+1290x15y8
     446    -28x18y4+240x17y5-976x16y6+2208x15y7-2494x14y8+816x13y9+56x15y6
     447    -320x14y7+844x13y8-1216x12y9+440x11y10-70x12y8+240x11y9-344x10y10
     448    +240x9y11+56x9y10-96x8y11+52x7y12-28x6y12+16x5y13+8x3y14-y16;
     449
     450
     451// reducible polynomials
     452
     453N = 12;
     454f = ((y2-x3)^2 - 4x5y - x7)*(x2-y3);
     455
     456f = 2x3y3-2y5+x4-xy2;
     457
     458f = -x3y3+x6y+xy6-x4y4;
     459*/
Note: See TracChangeset for help on using the changeset viewer.