Changeset f2f858 in git


Ignore:
Timestamp:
Jun 2, 2014, 3:52:56 PM (9 years ago)
Author:
Andreas Steenpass <steenpass@…>
Branches:
(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
Children:
9f1bf0bd7d089909777bd2491baed90ca7fcfeaf
Parents:
175eb21f4e380c8139d078ea81f56e1483e54df4
git-author:
Andreas Steenpass <steenpass@mathematik.uni-kl.de>2014-06-02 15:52:56+02:00
git-committer:
Andreas Steenpass <steenpass@mathematik.uni-kl.de>2014-06-03 10:42:24+02:00
Message:
chg: rewrite assprimeszerodim.lib based on modular.lib

(cherry picked from commit eeb2187019e8f379569be07360d245edfba0903a)
Signed-off-by: Andreas Steenpass <steenpass@mathematik.uni-kl.de>

Conflicts:
	Singular/LIB/assprimeszerodim.lib
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/assprimeszerodim.lib

    r175eb21 rf2f858  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="version assprimeszerodim.lib 4.0.0.0 Jun_2013 "; // $Id$
    3 category = "Commutative Algebra";
     2version="version assprimeszerodim.lib 4.0.0.0 Jun_2014 "; // $Id$
     3category="Commutative Algebra";
    44info="
    55LIBRARY:  assprimeszerodim.lib   associated primes of a zero-dimensional ideal
    66
    77AUTHORS:  N. Idrees       nazeranjawwad@gmail.com
    8 @*        G. Pfister      pfister@mathematik.uni-kl.de
    9 @*        S. Steidel      steidel@mathematik.uni-kl.de
     8          G. Pfister      pfister@mathematik.uni-kl.de
     9          A. Steenpass    steenpass@mathematik.uni-kl.de
     10          S. Steidel      steidel@mathematik.uni-kl.de
    1011
    1112OVERVIEW:
    … …  
    2829
    2930proc zeroRadical(ideal I, list #)
    30 "USAGE:  zeroRadical(I,[n]); I ideal, optional: n number of processors (for
    31          parallel computing)
     31"USAGE:  zeroRadical(I[, exactness]); I ideal, exactness int
    3232ASSUME:  I is zero-dimensional in Q[variables]
    3333RETURN:  the radical of I
     34NOTE:    A final test is applied to the result if exactness != 0 (default),
     35         otherwise no final test is done.
    3436EXAMPLE: example zeroRadical; shows an example
    3537"
    3638{
    37    return(zeroR(modStd(I,#),#));
     39   /* read optional parameter */
     40   int exactness = 1;
     41   if(size(#) > 0)
     42   {
     43      if(size(#) > 1 || typeof(#[1]) != "int")
     44      {
     45         ERROR("wrong optional parameter");
     46      }
     47      exactness = #[1];
     48   }
     49
     50   /* compute a standard basis if necessary */
     51   if (!attrib(I, "isSB"))
     52   {
     53      I = modStd(I, exactness);
     54   }
     55
     56   /* call modular() */
     57   // TODO: write deleteUnluckyPrimes_zeroRadical()
     58   if(exactness)
     59   {
     60      ideal F = modular("Assprimeszerodim::zeroRadP", list(I),
     61         Modstd::primeTest_std, Modular::deleteUnluckyPrimes_default,
     62         pTest_zeroRadical, finalTest_zeroRadical);
     63   }
     64   else
     65   {
     66      ideal F = modular("Assprimeszerodim::zeroRadP", list(I),
     67         Modstd::primeTest_std, Modular::deleteUnluckyPrimes_default,
     68         pTest_zeroRadical);
     69   }
     70
     71   /* compute the squarefree parts */
     72   poly f;
     73   int k;
     74   int i;
     75   for(i = nvars(basering); i > 0; i--)
     76   {
     77      f = gcd(F[i], diff(F[i], var(i)));
     78      k = k + deg(f);
     79      F[i] = F[i]/f;
     80   }
     81
     82   /* return the result */
     83   if(k == 0)
     84   {
     85      return(I);
     86   }
     87   else
     88   {
     89      return(modStd(I + F, exactness));
     90   }
    3891}
    3992example
    … …  
    4699////////////////////////////////////////////////////////////////////////////////
    47100
    48 static proc zeroR(ideal I, list #)
    49 // compute the radical of I provided that I is zero-dimensional in Q[variables]
    50 // and a standard basis
    51 {
    52    attrib(I,"isSB",1);
    53    int i, k;
    54    int j = 1;
    55    int index = 1;
    56    int crit;
    57 
    58    list CO1, CO2, P;
    59    ideal G, F;
    60    bigint N;
    61    poly f;
    62 
    63 //---------------------  Initialize optional parameter  ------------------------
    64    if(size(#) > 0)
    65    {
    66       int n1 = #[1];
    67       if(n1 >= 10)
    68       {
    69          int n2 = n1 + 1;
    70          int n3 = n1;
    71       }
    72       else
    73       {
    74          int n2 = 10;
    75          int n3 = 10;
    76       }
     101/* The pTest for zeroRadical(), to be used in modular(). */
     102static proc pTest_zeroRadical(string command, list args, ideal result, int p)
     103{
     104   /* change to characteristic p */
     105   def br = basering;
     106   list lbr = ringlist(br);
     107   if(typeof(lbr[1]) == "int")
     108   {
     109      lbr[1] = p;
    77110   }
    78111   else
    79112   {
    80       int n1 = 1;
    81       int n2 = 10;
    82       int n3 = 10;
    83    }
    84 
    85 //--------------------  Initialize the list of primes  -------------------------
    86    intvec L = primeList(I,n2);
    87    L[5] = prime(random(100000000,536870912));
    88 
    89    if(n1 > 1)
    90    {
    91 
    92 //-----  Create n links l(1),...,l(n1), open all of them and compute  ----------
    93 //-----  polynomial F for the primes L[2],...,L[n1 + 1].              ----------
    94 
    95       for(i = 1; i <= n1; i++)
    96       {
    97          //link l(i) = "MPtcp:fork";
    98          link l(i) = "ssi:fork";
    99          open(l(i));
    100          write(l(i), quote(zeroRadP(eval(I), eval(L[i + 1]))));
    101       }
    102 
    103       int t = timer;
    104       P = zeroRadP(I, L[1]);
    105       t = timer - t;
    106       if(t > 60) { t = 60; }
    107       int i_sleep = system("sh", "sleep "+string(t));
    108       CO1[index] = P[1];
    109       CO2[index] = bigint(P[2]);
    110       index++;
    111 
    112       j = j + n1 + 1;
    113    }
    114 
    115 //---------  Main computations in positive characteristic start here  ----------
    116 
    117    while(!crit)
    118    {
    119       if(n1 > 1)
    120       {
    121          while(j <= size(L) + 1)
    122          {
    123             for(i = 1; i <= n1; i++)
    124             {
    125                if(status(l(i), "read", "ready"))
    126                {
    127                   //--- read the result from l(i) ---
    128                   P = read(l(i));
    129                   CO1[index] = P[1];
    130                   CO2[index] = bigint(P[2]);
    131                   index++;
    132 
    133                   if(j <= size(L))
    134                   {
    135                      write(l(i), quote(zeroRadP(eval(I), eval(L[j]))));
    136                      j++;
    137                   }
    138                   else
    139                   {
    140                      k++;
    141                      close(l(i));
    142                   }
    143                }
    144             }
    145             //--- k describes the number of closed links ---
    146             if(k == n1)
    147             {
    148                j++;
    149             }
    150             //--- sleep for t seconds ---
    151             i_sleep = system("sh", "sleep "+string(t));
    152          }
    153       }
    154       else
    155       {
    156          while(j <= size(L))
    157          {
    158             P = zeroRadP(I, L[j]);
    159             CO1[index] = P[1];
    160             CO2[index] = bigint(P[2]);
    161             index++;
    162             j++;
    163          }
    164       }
    165 
    166       // insert deleteUnluckyPrimes
    167       G = chinrem(CO1,CO2);
    168       N = CO2[1];
    169       for(i = 2; i <= size(CO2); i++){ N = N*CO2[i]; }
    170       F = farey(G,N);
    171 
    172       crit = 1;
    173       for(i = 1; i <= nvars(basering); i++)
    174       {
    175          if(reduce(F[i],I) != 0) { crit = 0; break; }
    176       }
    177 
    178       if(!crit)
    179       {
    180          CO1 = G;
    181          CO2 = N;
    182          index = 2;
    183 
    184          j = size(L) + 1;
    185          L = primeList(I,n3,L);
    186 
    187          if(n1 > 1)
    188          {
    189             for(i = 1; i <= n1; i++)
    190             {
    191                open(l(i));
    192                write(l(i), quote(zeroRadP(eval(I), eval(L[j+i-1]))));
    193             }
    194             j = j + n1;
    195             k = 0;
    196          }
    197       }
    198    }
    199 
    200    k = 0;
    201    for(i = 1; i <= nvars(basering); i++)
    202    {
    203       f = gcd(F[i],diff(F[i],var(i)));
    204       k = k + deg(f);
    205       F[i] = F[i]/f;
    206    }
    207 
    208    if(k == 0) { return(I); }
    209    else { return(modStd(I + F, n1)); }
    210 }
    211 
    212 ////////////////////////////////////////////////////////////////////////////////
    213 
    214 proc assPrimes(list #)
    215 "USAGE:  assPrimes(I,[n],[a]); I ideal or module,
    216          optional: int n: number of processors (for parallel computing), int a:
    217          - a = 1: method of Eisenbud/Hunecke/Vasconcelos
    218          - a = 2: method of Gianni/Trager/Zacharias
    219          - a = 3: method of Monico
    220          assPrimes(I) chooses n = a = 1
     113      lbr[1][1] = p;
     114   }
     115   def rp = ring(lbr);
     116   setring(rp);
     117   ideal Ip = fetch(br, args)[1];
     118   ideal Fp_result = fetch(br, result);
     119
     120   /* run the command and compare with given result */
     121   execute("ideal Fp = "+command+"(Ip);");
     122   int i;
     123   for(i = nvars(br); i > 0; i--)
     124   {
     125      if(Fp[i] != Fp_result[i])
     126      {
     127         setring(br);
     128         return(0);
     129      }
     130   }
     131   setring(br);
     132   return(1);
     133}
     134
     135////////////////////////////////////////////////////////////////////////////////
     136
     137/* The finalTest for zeroRadical, to be used in modular(). */
     138static proc finalTest_zeroRadical(string command, list args, ideal F)
     139{
     140   int i;
     141   for(i = nvars(basering); i > 0; i--)
     142   {
     143      if(reduce(F[i], args[1]) != 0) { return(0); }
     144   }
     145   return(1);
     146}
     147
     148////////////////////////////////////////////////////////////////////////////////
     149
     150proc assPrimes(def I, list #)
     151"USAGE:  assPrimes(I[, alg, exactness]); I ideal or module,
     152         alg string (optional), exactness int (optional)
     153         - alg = "GTZ":    method of Gianni/Trager/Zacharias (default)
     154         - alg = "EHV":    method of Eisenbud/Hunecke/Vasconcelos
     155         - alg = "Monico": method of Monico
    221156ASSUME:  I is zero-dimensional over Q[variables]
    222 RETURN:  a list Re of associated primes of I:
     157RETURN:  a list of the associated primes of I
     158NOTE:    A final test is applied to the result if exactness != 0 (default),
     159         otherwise no final test is done.
    223160EXAMPLE: example assPrimes; shows an example
    224161"
    225162{
    226    ideal I;
    227    if(typeof(#[1]) == "ideal")
    228    {
    229       I = #[1];
    230    }
    231    else
    232    {
    233       module M = #[1];
    234       I = Ann(M);
    235    }
    236 
    237 //---------------------  Initialize optional parameter  ------------------------
    238    if(size(#) > 1)
    239    {
    240       if(size(#) == 2)
    241       {
    242          int n1 = #[2];
    243          int alg = 1;
    244          if(n1 >= 10)
    245          {
    246             int n2 = n1 + 1;
    247             int n3 = n1;
    248          }
    249          else
    250          {
    251             int n2 = 10;
    252             int n3 = 10;
    253          }
    254       }
    255       if(size(#) == 3)
    256       {
    257          int n1 = #[2];
    258          int alg = #[3];
    259          if(n1 >= 10)
    260          {
    261             int n2 = n1 + 1;
    262             int n3 = n1;
    263          }
    264          else
    265          {
    266             int n2 = 10;
    267             int n3 = 10;
    268          }
    269       }
    270    }
    271    else
    272    {
    273       int n1 = 1;
    274       int alg = 1;
    275       int n2 = 10;
    276       int n3 = 10;
    277    }
    278 
    279    if(printlevel >= 10)
    280    {
    281       "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3);
    282    }
    283 
    284    int T = timer;
    285    int RT = rtimer;
    286    int TT;
    287    int t_sleep;
    288 
    289    def SPR = basering;
    290    map phi;
    291    list H = ideal(0);
    292    ideal F;
    293    poly F1;
    294 
     163   /* read input */
     164   if(typeof(I) != "ideal")
     165   {
     166      if(typeof(I) != "module")
     167      {
     168         ERROR("The first argument must be of type 'ideal' or 'module'.");
     169      }
     170      module M = I;
     171      kill I;
     172      ideal I = Ann(M);
     173      kill M;
     174   }
     175
     176   /* read optional parameters */
     177   list defaults = list("GTZ", 1);
     178   int i;
     179   for(i = 1; i <= size(defaults); i++)
     180   {
     181      if(typeof(#[i]) != typeof(defaults[i]))
     182      {
     183         # = insert(#, defaults[i], i-1);
     184      }
     185   }
     186   if(size(#) != size(defaults))
     187   {
     188      ERROR("wrong optional parameters");
     189   }
     190   string alg = #[1];
     191   int exactness = #[2];
     192   int a;
     193   if(alg == "GTZ")
     194   {
     195      a = 1;
     196   }
     197   if(alg == "EHV")
     198   {
     199      a = 2;
     200   }
     201   if(alg == "Monico")
     202   {
     203      a = 3;
     204   }
     205   if(a == 0)   // alg != "GTZ" && alg != "EHV" && alg != "Monico"
     206   {
     207      ERROR("unknown method");
     208   }
     209
     210   /* compute a standard basis if necessary */
    295211   if(printlevel >= 10) { "========== Start modStd =========="; }
    296    I = modStd(I);
     212   if (!attrib(I, "isSB")) {
     213       I = modStd(I, exactness);
     214   }
    297215   if(printlevel >= 10) { "=========== End modStd ==========="; }
    298    if(printlevel >= 9) { "modStd takes "+string(rtimer-RT)+" seconds."; }
    299216   int d = vdim(I);
    300217   if(d == -1) { ERROR("Ideal is not zero-dimensional."); }
    301218   if(homog(I) == 1) { return(list(maxideal(1))); }
    302    poly f = findGen(I);
    303    if(printlevel >= 9) { "Coordinate change: "+string(f); }
    304 
    305    if(size(f) == nvars(SPR))
    306    {
    307       TT = timer;
    308       int spT = pTestRad(d,f,I);
    309       if(printlevel >= 9)
    310       {
    311          "pTestRad(d,f,I) = "+string(spT)+" and takes "
    312          +string(timer-TT)+" seconds.";
    313       }
    314       if(!spT)
    315       {
    316          if(typeof(attrib(#[1],"isRad")) == "int")
    317          {
    318             if(attrib(#[1],"isRad") == 0)
    319             {
    320                TT = timer;
    321                I = zeroR(I,n1);
    322                if(printlevel >= 9)
    323                {
    324                   "zeroR(I,n1) takes "+string(timer-TT)+" seconds.";
    325                }
    326                TT = timer;
    327                I = modStd(I);
    328                if(printlevel >= 9)
    329                {
    330                   "modStd(I) takes "+string(timer-TT)+" seconds.";
    331                }
    332                d = vdim(I);
    333                f = findGen(I);
    334             }
    335          }
    336          else
    337          {
    338             TT = timer;
    339             I = zeroR(I,n1);
    340             if(printlevel >= 9)
    341             {
    342                "zeroR(I,n1) takes "+string(timer-TT)+" seconds.";
    343             }
    344             TT = timer;
    345             I = modStd(I);
    346             if(printlevel >= 9)
    347             {
    348                "modStd(I) takes "+string(timer-TT)+" seconds.";
    349             }
    350             d = vdim(I);
    351             f = findGen(I);
    352          }
    353       }
    354    }
    355    if(printlevel >= 9)
    356    {
    357       "Real-time for radical-check is "+string(rtimer - RT)+" seconds.";
    358       "CPU-time for radical-check is "+string(timer - T)+" seconds.";
    359    }
    360 
    361    export(SPR);
    362    poly f_for_fork = f;
    363    export(f_for_fork);           // f available for each link
    364    ideal I_for_fork = I;
    365    export(I_for_fork);           // I available for each link
    366 
    367 //--------------------  Initialize the list of primes  -------------------------
    368    intvec L = primeList(I,n2);
    369    L[5] = prime(random(1000000000,2134567879));
    370 
    371    list Re;
    372 
    373    ring rHelp = 0,T,dp;
    374    list CO1,CO2,P,H;
    375    ideal F,G,testF;
    376    bigint N;
    377 
    378    list ringL = ringlist(SPR);
    379    int i,k,e,int_break,s;
    380    int j = 1;
    381    int index = 1;
    382 
    383 //-----  If there is more than one processor available, we parallelize the  ----
    384 //-----  main standard basis computations in positive characteristic        ----
    385 
    386    if(n1 > 1)
    387    {
    388 
    389 //-----  Create n1 links l(1),...,l(n1), open all of them and compute  ---------
    390 //-----  standard basis for the primes L[2],...,L[n1 + 1].             ---------
    391 
    392       for(i = 1; i <= n1; i++)
    393       {
    394          //link l(i) = "MPtcp:fork";
    395          link l(i) = "ssi:fork";
    396          open(l(i));
    397          write(l(i), quote(modpSpecialAlgDep(eval(ringL), eval(L[i + 1]),
    398                                                           eval(alg))));
    399       }
    400 
    401       int t = timer;
    402       P = modpSpecialAlgDep(ringL, L[1], alg);
    403       t = timer - t;
    404       if(t > 60) { t = 60; }
    405       int i_sleep = system("sh", "sleep "+string(t));
    406       CO1[index] = P[1];
    407       CO2[index] = bigint(P[2]);
    408       index++;
    409 
    410       j = j + n1 + 1;
    411    }
    412 
    413 //---------  Main computations in positive characteristic start here  ----------
    414 
    415    int tt = timer;
    416    int rt = rtimer;
    417 
    418    while(1)
    419    {
    420       tt = timer;
    421       rt = rtimer;
    422 
    423       if(printlevel >= 9) { "size(L) = "+string(size(L)); }
    424 
    425       if(n1 > 1)
    426       {
    427          while(j <= size(L) + 1)
    428          {
    429             for(i = 1; i <= n1; i++)
    430             {
    431                //--- ask if link l(i) is ready otherwise sleep for t seconds ---
    432                if(status(l(i), "read", "ready"))
    433                {
    434                   //--- read the result from l(i) ---
    435                   P = read(l(i));
    436                   CO1[index] = P[1];
    437                   CO2[index] = bigint(P[2]);
    438                   index++;
    439 
    440                   if(j <= size(L))
    441                   {
    442                      write(l(i), quote(modpSpecialAlgDep(eval(ringL),
    443                                                          eval(L[j]),
    444                                                          eval(alg))));
    445                      j++;
    446                   }
    447                   else
    448                   {
    449                      k++;
    450                      close(l(i));
    451                   }
    452                }
    453             }
    454             //--- k describes the number of closed links ---
    455             if(k == n1)
    456             {
    457                j++;
    458             }
    459             i_sleep = system("sh", "sleep "+string(t));
    460          }
    461       }
    462       else
    463       {
    464          while(j <= size(L))
    465          {
    466             P = modpSpecialAlgDep(ringL, L[j], alg);
    467             CO1[index] = P[1];
    468             CO2[index] = bigint(P[2]);
    469             index++;
    470             j++;
    471          }
    472       }
    473 
    474       if(printlevel >= 9)
    475       {
    476          "Real-time for computing list in assPrimes is "+string(rtimer - rt)+
    477          " seconds.";
    478          "CPU-time for computing list in assPrimes is "+string(timer - tt)+
    479          " seconds.";
    480       }
    481 
    482 //-------------------  Lift results to basering via farey ----------------------
    483 
    484       tt = timer;
    485       G = chinrem(CO1,CO2);
    486       N = CO2[1];
    487       for(j = 2; j <= size(CO2); j++){ N = N*CO2[j]; }
    488       F = farey(G,N);
    489       if(printlevel >= 10) { "Lifting-process takes "+string(timer - tt)
    490                              +" seconds"; }
    491 
    492       if(pTestPoly(F[1], ringL, alg, L))
    493       {
    494          F = cleardenom(F[1]);
    495 
    496          e = deg(F[1]);
    497          if(e == d)
    498          {
    499             H = factorize(F[1]);
    500 
    501             s = size(H[1]);
    502             for(i = 1; i <= s; i++)
    503             {
    504                if(H[2][i] != 1)
    505                {
    506                   int_break = 1;
    507                }
    508             }
    509 
    510             if(int_break == 0)
    511             {
    512                setring SPR;
    513                phi = rHelp,var(nvars(SPR));
    514                H = phi(H);
    515 
    516                if(printlevel >= 9)
    517                {
    518                   "Real-time without test is "+string(rtimer - RT)+" seconds.";
    519                   "CPU-time without test is "+string(timer - T)+" seconds.";
    520                }
    521 
    522                T = timer;
    523                RT = rtimer;
    524 
    525                F = phi(F);
    526 
    527                if(n1 > 1)
    528                {
    529                   open(l(1));
    530                   write(l(1), quote(quickSubst(eval(F[1]), eval(f), eval(I))));
    531                   t_sleep = timer;
    532                }
    533                else
    534                {
    535                   F1 = quickSubst(F[1],f,I);
    536                   if(F1 != 0) { int_break = 1; }
    537                }
    538 
    539                if(int_break == 0)
    540                {
    541                   for(i = 2; i <= s; i++)
    542                   {
    543                      H[1][i] = quickSubst(H[1][i],f,I);
    544                      Re[i-1] = I + ideal(H[1][i]);
    545                   }
    546 
    547                   if(n1 > 1)
    548                   {
    549                      t_sleep = timer - t_sleep;
    550                      if(t_sleep > 5) { t_sleep = 5; }
    551 
    552                      while(1)
    553                      {
    554                         if(status(l(1), "read", "ready"))
    555                         {
    556                            F1 = read(l(1));
    557                            close(l(1));
    558                            break;
    559                         }
    560                         i_sleep = system("sh", "sleep "+string(t_sleep));
    561                      }
    562                      if(F1 != 0) { int_break = 1; }
    563                   }
    564                   if(printlevel >= 9)
    565                   {
    566                      "Real-time for test is "+string(rtimer - RT)+" seconds.";
    567                      "CPU-time for test is "+string(timer - T)+" seconds.";
    568                   }
    569                   if(int_break == 0)
    570                   {
    571                      kill f_for_fork;
    572                      kill I_for_fork;
    573                      kill SPR;
    574                      return(Re);
    575                   }
    576                }
    577             }
    578          }
    579       }
    580 
    581       int_break = 0;
    582       setring rHelp;
    583       testF = F;
    584       CO1 = G;
    585       CO2 = N;
    586       index = 2;
    587 
    588       j = size(L) + 1;
    589 
    590       setring SPR;
    591       L = primeList(I,n3,L);
    592       setring rHelp;
    593 
    594       if(n1 > 1)
    595       {
    596          for(i = 1; i <= n1; i++)
    597          {
    598             open(l(i));
    599             write(l(i), quote(modpSpecialAlgDep(eval(ringL), eval(L[j+i-1]),
    600                                                              eval(alg))));
    601          }
    602          j = j + n1;
    603          k = 0;
    604       }
    605    }
     219
     220   /* compute the radical if necessary */
     221   ideal J = I;
     222   int isRad;
     223   poly f;
     224   isRad, f = pTestRad(I, d);
     225   while(!isRad)
     226   {
     227      J = zeroRadical(I, exactness);
     228      J = modStd(J, exactness);
     229      d = vdim(J);
     230      isRad, f = pTestRad(J, d);
     231   }
     232   I = J;
     233   kill J;
     234
     235   /* call modular() */
     236   if(exactness)
     237   {
     238      ideal F = modular("Assprimeszerodim::modpSpecialAlgDep",
     239         list(I, f, d, a),
     240         Modstd::primeTest_std, Modular::deleteUnluckyPrimes_default,
     241         pTest_assPrimes, finalTest_assPrimes);
     242   }
     243   else
     244   {
     245      ideal F = modular("Assprimeszerodim::modpSpecialAlgDep",
     246         list(I, f, d, a),
     247         Modstd::primeTest_std, Modular::deleteUnluckyPrimes_default,
     248         pTest_assPrimes);
     249   }
     250
     251   /* compute the components */
     252   list result;
     253   list H = factorize(F[1]);
     254   for(i = size(H[1])-1; i > 0; i--)
     255   {
     256      result[i] = I + ideal(quickSubst(H[1][i+1], f, I));
     257   }
     258
     259   /* return the result */
     260   return(result);
    606261}
    607262example
    … …  
    620275////////////////////////////////////////////////////////////////////////////////
    621276
    622 static proc specialAlgDepEHV(poly p, ideal I)
    623 {
    624 //=== computes a poly F in Q[T] such that <F>=kernel(Q[T]--->basering)
    625 //=== mapping T to p
     277/* Computes a poly F in Q[T] such that
     278 * <F> = kernel(Q[T] --> basering, T |-> f),
     279 * T := last variable in the basering.
     280 */
     281static proc specialAlgDepEHV(ideal I, poly f)
     282{
    626283   def R = basering;
    627    execute("ring Rhelp="+charstr(R)+",T,dp;");
    628    setring R;
    629    map phi = Rhelp,p;
    630    setring Rhelp;
    631    ideal F = preimage(R,phi,I); //corresponds to std(I,p-T) in dp(n),dp(1)
    632    export(F);
    633    setring R;
    634    list L = Rhelp;
    635    return(L);
    636 }
    637 
    638 ////////////////////////////////////////////////////////////////////////////////
    639 
    640 static proc specialAlgDepGTZ(poly p, ideal I)
    641 {
    642 //=== assume I is zero-dimensional
    643 //=== computes a poly F in Q[T] such that <F>=kernel(Q[T]--->basering)
    644 //=== mapping T to p
     284   execute("ring QT = ("+charstr(R)+"), "+varstr(R, nvars(R))+", dp;");
     285   setring(R);
     286   map phi = QT, f;
     287   setring QT;
     288   ideal F = preimage(R, phi, I); // corresponds to std(I, f-T) in dp(n),dp(1)
     289   setring(R);
     290   ideal F = imap(QT, F);
     291   return(F);
     292}
     293
     294////////////////////////////////////////////////////////////////////////////////
     295
     296/* Assume I is zero-dimensional.
     297 * Computes a poly F in Q[T] such that
     298 * <F> = kernel(Q[T] --> basering, T |-> f),
     299 * T := last variable in the basering.
     300 */
     301static proc specialAlgDepGTZ(ideal I, poly f)
     302{
    645303   def R = basering;
    646    execute("ring Rhelp = "+charstr(R)+",T,dp;");
    647    setring R;
    648    map phi = Rhelp,p;
    649    def Rlp = changeord(list(list("dp",1:(nvars(R)-1)),list("dp",1:1)));
    650    setring Rlp;
    651    poly p = imap(R,p);
     304   if(nvars(R) > 1)
     305   {
     306      def Rlp = changeord(list(list("dp", 1:(nvars(R)-1)), list("dp", 1:1)));
     307      setring(Rlp);
     308      poly f = imap(R, f);
     309      ideal I;
     310   }
    652311   ideal K = maxideal(1);
    653    K[nvars(R)] = 2*var(nvars(R))-p;
    654    map phi = R,K;
    655    ideal I = phi(I);
     312   K[nvars(R)] = 2*var(nvars(R))-f;
     313   map phi = R, K;
     314   I = phi(I);
    656315   I = std(I);
    657    poly q = subst(I[1],var(nvars(R)),var(1));
    658    setring Rhelp;
    659    map psi=Rlp,T;
    660    ideal F=psi(q);
    661    export(F);
    662    setring R;
    663    list L=Rhelp;
    664    return(L);
    665 }
    666 
    667 ////////////////////////////////////////////////////////////////////////////////
    668 
    669 static proc specialAlgDepMonico(poly p, ideal I)
    670 {
    671 //=== assume I is zero-dimensional
    672 //=== computes a poly F in Q[T], the characteristic polynomial of the map
    673 //=== basering/I ---> baserng/I  defined by the multiplication with p
    674 //=== in case I is radical it is the same poly as in specialAlgDepEHV
     316   ideal F = I[1];
     317   if(nvars(R) > 1)
     318   {
     319      setring(R);
     320      ideal F = imap(Rlp, F);
     321   }
     322   return(F);
     323}
     324
     325////////////////////////////////////////////////////////////////////////////////
     326
     327/* Assume I is zero-dimensional.
     328 * Computes a poly F in Q[T], the characteristic polynomial of the map
     329 * basering/I ---> basering/I  defined by the multiplication with f,
     330 * T := last variable in the basering.
     331 * In case I is radical, it is the same polynomial as in specialAlgDepEHV.
     332 */
     333static proc specialAlgDepMonico(ideal I, poly f, int d)
     334{
    675335   def R = basering;
    676    execute("ring Rhelp = "+charstr(R)+",T,dp;");
    677    setring R;
    678    map phi = Rhelp,p;
    679    poly q;
    680336   int j;
    681    matrix m ;
    682    poly va = var(1);
     337   matrix M[d][d];
    683338   ideal J = std(I);
    684    ideal ba = kbase(J);
    685    int d = vdim(J);
    686    matrix n[d][d];
    687    for(j = 2; j <= nvars(R); j++)
    688    {
    689      va = va*var(j);
     339   ideal basis = kbase(J);
     340   poly vars = var(nvars(R));
     341   for(j = nvars(R)-1; j > 0; j--)
     342   {
     343     vars = var(j)*vars;
    690344   }
    691345   for(j = 1; j <= d; j++)
    692346   {
    693      q = reduce(p*ba[j],J);
    694      m = coeffs(q,ba,va);
    695      n[j,1..d] = m[1..d,1];
    696    }
    697    setring Rhelp;
    698    matrix n = imap(R,n);
    699    ideal F = det(n-T*freemodule(d));
    700    export(F);
    701    setring R;
    702    list L = Rhelp;
    703    return(L);
     347     M[1..d, j] = coeffs(reduce(f*basis[j], J), basis, vars);
     348   }
     349   execute("ring QT = ("+charstr(R)+"), "+varstr(R, nvars(R))+", dp;");
     350   matrix M = imap(R, M);
     351   ideal F = det(M-var(1)*freemodule(d));
     352   setring(R);
     353   ideal F = imap(QT, F);
     354   return(F);
    704355}
    705356
    … …  
    711362//=== mapping T to p and test if d=deg(F)
    712363   def R = basering;
    713    execute("ring Rhelp="+charstr(R)+",T,dp;");
     364   execute("ring Rhelp = ("+charstr(R)+"), T, dp;");
    714365   setring R;
    715366   map phi = Rhelp,p;
    … …  
    718369   int e=deg(F[1]);
    719370   setring R;
    720    return((e==d));
    721 }
    722 
    723 ////////////////////////////////////////////////////////////////////////////////
    724 
    725 static proc findGen(ideal J, list #)
    726 {
    727 //=== try to find a sparse linear form r such that
    728 //=== vector space dim(basering/J)=deg(F),
    729 //=== F a poly in Q[T] such that <F>=kernel(Q[T]--->basering) mapping T to r
    730 //=== if not found returns a generic (randomly chosen) r
    731    int d = vdim(J);
     371   return((e==d), fetch(Rhelp, F)[1]);
     372}
     373
     374////////////////////////////////////////////////////////////////////////////////
     375
     376/* Assume d = vector space dim(basering/J).
     377 * Tries to find a (sparse) linear form r such that d = deg(F), where
     378 * F is a poly in Q[T] such that <F> = kernel(Q[T]-->basering) mapping T to r.
     379 * If found, returns (1, r, F). If not found, returns (0, 0, 0).
     380 */
     381static proc findGen(ideal J, int d)
     382{
    732383   def R = basering;
    733384   int n = nvars(R);
    734    list rl = ringlist(R);
    735    if(size(#) > 0) { int p = #[1]; }
    736    else { int p = prime(random(1000000000,2134567879)); }
    737    rl[1] = p;
    738    def @R = ring(rl);
    739    setring @R;
    740    ideal J = imap(R,J);
     385   int okay;
     386   poly F;
     387
     388   /* try trivial transformation */
    741389   poly r = var(n);
    742    int i,k;
    743    k = specialTest(d,r,J);
    744    if(!k)
     390   okay, F = specialTest(d, r, J);
     391   if(okay)
     392   {
     393      return(1, r, F);
     394   }
     395
     396   /* try transformations of the form var(n) + var(i) */
     397   int i;
     398   for(i = 1; i < n; i++)
     399   {
     400      okay, F = specialTest(d, r+var(i), J);
     401      if(okay)
     402      {
     403         return(1, r+var(i), F);
     404      }
     405   }
     406
     407   /* try transformations of the form var(n) + \sum var(i) */
     408   if(n > 2)
    745409   {
    746410      for(i = 1; i < n; i++)
    747411      {
    748          k = specialTest(d,r+var(i),J);
    749          if(k){ r = r + var(i); break; }
    750       }
    751    }
    752    if((!k) && (n > 2))
    753    {
    754       for(i = 1; i < n; i++)
    755       {
    756412         r = r + var(i);
    757          k = specialTest(d,r,J);
    758          if(k){ break; }
    759       }
    760    }
    761    setring R;
    762    poly r = randomLast(100)[nvars(R)];
    763    if(k){ r = imap(@R,r); }
    764    return(r);
    765 }
    766 
    767 ////////////////////////////////////////////////////////////////////////////////
    768 
    769 static proc pTestRad(int d, poly p1, ideal I)
    770 {
    771 //=== computes a poly F in Z/q1[T] such that
    772 //===                    <F> = kernel(Z/q1[T]--->Z/q1[vars(basering)])
    773 //=== mapping T to p1 and test if d=deg(squarefreepart(F)), q1 a prime randomly
    774 //=== chosen
    775 //=== If not choose randomly another prime q2 and another linear form p2 and
    776 //=== computes a poly F in Z/q2[T] such that
    777 //===                    <F> = kernel(Z/q2[T]--->Z/q2[vars(basering)])
    778 //=== mapping T to p2 and test if d=deg(squarefreepart(F))
    779 //=== if the test is positive then I is radical
     413         okay, F = specialTest(d, r, J);
     414         if(okay)
     415         {
     416            return(1, r, F);
     417         }
     418      }
     419   }
     420
     421   /* try random transformations */
     422   int N = 2;   // arbitrarily chosen
     423   for(i = N; i > 0; i--)
     424   {
     425      r = randomLast(100)[n];
     426      okay, F = specialTest(d, r, J);
     427      if(okay)
     428      {
     429         return(1, r, F);
     430      }
     431   }
     432
     433   /* not found */
     434   return(0, 0, 0);
     435}
     436
     437////////////////////////////////////////////////////////////////////////////////
     438
     439/* Assume d = vector space dim(basering/I).
     440 * Tests if I is radical over F_p, where p is some randomly chosen prime.
     441 * If yes, chooses a linear form r such that d = deg(squarefreepart(F)), where
     442 * F is a poly in Z/p[T] such that <F> = kernel(Z/p[T]-->Z/p[vars(basering)])
     443 * mapping T to r.
     444 * Returns (1, r), if I is radical over F_p, and (0, 0) otherwise.
     445 */
     446static proc pTestRad(ideal I, int d)
     447{
     448   int N = 2;   // Try N random primes. Value of N can be chosen arbitrarily.
    780449   def R = basering;
    781450   list rl = ringlist(R);
    782    int q1 = prime(random(100000000,536870912));
    783    rl[1] = q1;
    784    ring Shelp1 = q1,T,dp;
    785    setring R;
    786    def Rhelp1 = ring(rl);
    787    setring Rhelp1;
    788    poly p1 = imap(R,p1);
    789    ideal I = imap(R,I);
    790    map phi = Shelp1,p1;
    791    setring Shelp1;
    792    ideal F = preimage(Rhelp1,phi,I);
    793    poly f = gcd(F[1],diff(F[1],var(1)));
    794    int e = deg(F[1]/f);
    795    setring R;
    796    if(e != d)
    797    {
    798       poly p2 = findGen(I,q1);
    799       setring Rhelp1;
    800       poly p2 = imap(R,p2);
    801       phi = Shelp1,p2;
    802       setring Shelp1;
    803       F = preimage(Rhelp1,phi,I);
    804       f = gcd(F[1],diff(F[1],var(1)));
    805       e = deg(F[1]/f);
    806       setring R;
    807       if(e == d){ return(1); }
    808       if(e != d)
    809       {
    810          int q2 = prime(random(100000000,536870912));
    811          rl[1] = q2;
    812          ring Shelp2 = q2,T,dp;
    813          setring R;
    814          def Rhelp2 = ring(rl);
    815          setring Rhelp2;
    816          poly p1 = imap(R,p1);
    817          ideal I = imap(R,I);
    818          map phi = Shelp2,p1;
    819          setring Shelp2;
    820          ideal F = preimage(Rhelp2,phi,I);
    821          poly f = gcd(F[1],diff(F[1],var(1)));
    822          e = deg(F[1]/f);
    823          setring R;
    824          if(e == d){ return(1); }
    825       }
    826    }
    827    return((e==d));
    828 }
    829 
    830 ////////////////////////////////////////////////////////////////////////////////
    831 
    832 static proc zeroRadP(ideal I, int p)
    833 {
    834 //=== computes F=(F_1,...,F_n) such that <F_i>=IZ/p[x_1,...,x_n] intersected
    835 //=== with Z/p[x_i], F_i monic
    836    def R0 = basering;
    837    list ringL = ringlist(R0);
    838    ringL[1] = p;
    839    def @r = ring(ringL);
    840    setring @r;
    841    ideal I = fetch(R0,I);
     451   int p;
     452   int okay;
     453   int i;
     454   for(i = N; i > 0; i--)
     455   {
     456      p = prime(random(100000000,536870912));
     457
     458      // change to characteristic p
     459      if(typeof(rl[1]) == "int")
     460      {
     461         rl[1] = p;
     462      }
     463      else
     464      {
     465         rl[1][1] = p;
     466      }
     467      def Rp(i) = ring(rl);
     468      setring Rp(i);
     469      ideal I = imap(R, I);
     470
     471      // find and test transformation
     472      poly r;
     473      poly F;
     474      okay, r, F = findGen(I, d);
     475      if(okay)
     476      {
     477         poly f = gcd(F, diff(F, var(1)));
     478         if(d == deg(F/f))   // F squarefree?
     479         {
     480            setring(R);
     481            return(1, imap(Rp(i), r));
     482         }
     483      }
     484      setring(R);
     485   }
     486   return(0, 0);
     487}
     488
     489////////////////////////////////////////////////////////////////////////////////
     490
     491/* Computes an ideal F such that ncols(F) = nvars(basering),
     492 * < F[i] > = (I intersected with K[var(i)]), and F[i] is monic.
     493 */
     494static proc zeroRadP(ideal I)
     495{
     496   intvec opt = option(get);
    842497   option(redSB);
    843498   I = std(I);
    844    ideal F = finduni(I); //F[i] generates I intersected with K[var(i)]
    845    int i;
    846    for(i = 1; i <= size(F); i++){ F[i] = simplify(F[i],1); }
    847    setring R0;
    848    return(list(fetch(@r,F),p));
     499   ideal F = finduni(I);   // F[i] generates I intersected with K[var(i)]
     500   F = simplify(F, 1);
     501   option(set, opt);
     502   return(F);
    849503}
    850504
    … …  
    883537////////////////////////////////////////////////////////////////////////////////
    884538
    885 static proc modpSpecialAlgDep(list ringL, int p, list #)
    886 {
    887 //=== prepare parallel computing
    888 //===  #=1: method of Eisenbud/Hunecke/Vasconcelos
    889 //===  #=2: method of Gianni/Trager/Zacharias
    890 //===  #=3: method of Monico
    891 
    892    def R0 = basering;
    893 
    894    ringL[1] = p;
    895    def @r = ring(ringL);
    896    setring @r;
    897    poly f = fetch(SPR,f_for_fork);
    898    ideal I = fetch(SPR,I_for_fork);
    899    if(size(#) > 0)
    900    {
    901       if(#[1] == 1) { list M = specialAlgDepEHV(f,I); }
    902       if(#[1] == 2) { list M = specialAlgDepGTZ(f,I); }
    903       if(#[1] == 3) { list M = specialAlgDepMonico(f,I); }
     539/* Simple switch for specialAlgDepEHV, specialAlgDepGTZ, and
     540 * specialAlgDepMonico.
     541 */
     542static proc modpSpecialAlgDep(ideal I, poly f, int d, int alg)
     543{
     544   ideal F;
     545   if(alg == 1) { F = specialAlgDepEHV(I, f); }
     546   if(alg == 2) { F = specialAlgDepGTZ(I, f); }
     547   if(alg == 3) { F = specialAlgDepMonico(I, f, d); }
     548   F = simplify(F, 1);
     549   return(F);
     550}
     551
     552////////////////////////////////////////////////////////////////////////////////
     553
     554/* The pTest for assPrimes(), to be used in modular(). */
     555static proc pTest_assPrimes(string command, list args, ideal F, int p)
     556{
     557   def br = basering;
     558   list lbr = ringlist(br);
     559   if(typeof(lbr[1]) == "int")
     560   {
     561      lbr[1] = p;
    904562   }
    905563   else
    906564   {
    907       list M = specialAlgDepEHV(f,I);
    908    }
    909    def @S = M[1];
    910 
    911    setring R0;
    912    return(list(imap(@S,F),p));
    913 }
    914 
    915 ////////////////////////////////////////////////////////////////////////////////
    916 
    917 static proc pTestPoly(poly testF, list ringL, int alg, list L)
    918 {
    919    int i,j,p;
    920    def R0 = basering;
    921 
    922    while(!i)
    923    {
    924       i = 1;
    925       p = prime(random(1000000000,2134567879));
    926       for(j = 1; j <= size(L); j++)
    927       {
    928          if(p == L[j]) { i = 0; break; }
    929       }
    930    }
    931 
    932    ringL[1] = p;
    933    def @R = ring(ringL);
    934    setring @R;
    935    poly f = fetch(SPR,f_for_fork);
    936    ideal I = fetch(SPR,I_for_fork);
    937    if(alg == 1) { list M = specialAlgDepEHV(f,I); }
    938    if(alg == 2) { list M = specialAlgDepGTZ(f,I); }
    939    if(alg == 3) { list M = specialAlgDepMonico(f,I); }
    940    def @S = M[1];
    941    setring @S;
    942    poly testF = fetch(R0,testF);
    943    int k = (testF == F);
    944 
    945    setring R0;
     565      lbr[1][1] = p;
     566   }
     567   def rp = ring(lbr);
     568   setring(rp);
     569   list args_p = fetch(br, args);
     570   ideal F = fetch(br, F);
     571   execute("ideal Fp = "+command+"("
     572      +Tasks::argsToString("args_p", size(args_p))+");");
     573   int k = (Fp[1] == F[1]);
     574   setring br;
    946575   return(k);
     576}
     577
     578////////////////////////////////////////////////////////////////////////////////
     579
     580/* The finalTest for assPrimes(), to be used in modular(). */
     581static proc finalTest_assPrimes(string command, alias list args, ideal F)
     582{
     583   F = cleardenom(F[1]);
     584   if(deg(F[1]) != args[3]) { return(0); }
     585   if(gcd(F[1], diff(F[1], var(nvars(basering)))) != 1) { return(0); };
     586   if(quickSubst(F[1], args[2], args[1]) != 0) { return(0); }
     587   return(1);
    947588}
    948589
Note: See TracChangeset for help on using the changeset viewer.