Changeset bb3e63b in git


Ignore:
Timestamp:
Jan 30, 2001, 9:39:29 AM (23 years ago)
Author:
Wilfred Pohl <pohl@…>
Branches:
(u'spielwiese', '4a9821a93ffdc22a6696668bd4f6b8c9de3e6c5f')
Children:
f1c4b1ef727ee82da892d0ed67ec94c7f1d1ec61
Parents:
995ba8e72c52d5e7ce134eb888bb7682e01af581
Message:
expbound and control for det


git-svn-id: file:///usr/local/Singular/svn/trunk@5130 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
Singular
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • Singular/sparsmat.cc

    r995ba8 rbb3e63b  
    22*  Computer Algebra System SINGULAR     *
    33****************************************/
    4 /* $Id: sparsmat.cc,v 1.48 2000-12-31 15:54:47 obachman Exp $ */
     4/* $Id: sparsmat.cc,v 1.49 2001-01-30 08:39:29 pohl Exp $ */
    55
    66/*
    … …  
    160160};
    161161
    162 Exponent_t smExpBound(ideal Id)
    163 {
    164   Exponent_t max = 0;
    165   long ldeg;
    166   int dummy, i, n = IDELEMS(Id);
    167  
    168   for (i=0; i<n; i++)
    169   {
    170     if (Id->m[i] != NULL)
    171     {
    172       ldeg = pLDeg(Id->m[i], &dummy);
    173       if (ldeg > max) max = ldeg;
    174     }
    175   }
    176   return max*2*n;
    177 }
    178 
    179 // #define HOMOG_LP
     162/*
     163* estimate maximal exponent for det or minors,
     164* the module m has di vectors and maximal rank ra,
     165* estimate yields for the tXt minors,
     166* we have di,ra >= t
     167*/
     168static void smMinSelect(Exponent_t *, int, int);
     169Exponent_t smExpBound( ideal m, int di, int ra, int t)
     170{
     171  poly p;
     172  Exponent_t kr, kc;
     173  Exponent_t *r, *c;
     174  int al, bl, i, j, k;
     175
     176  al = di*sizeof(Exponent_t);
     177  c = omAlloc(al);
     178  bl = ra*sizeof(Exponent_t);
     179  r = omAlloc0(bl);
     180  for (i=di-1;i>=0;i--)
     181  {
     182    kc = 0;
     183    p = m->m[i];
     184    while(p!=NULL)
     185    {
     186      k = pGetComp(p)-1;
     187      kr = r[k];
     188      for (j=pVariables;j>0;j--)
     189      {
     190        if(pGetExp(p,j)>kc)
     191          kc=pGetExp(p,j);
     192        if(pGetExp(p,j)>kr)
     193          kr=pGetExp(p,j);
     194      }
     195      r[k] = kr;
     196      pIter(p);
     197    }
     198    c[i] = kc;
     199  }
     200  if (t<di) smMinSelect(c, t, di);
     201  if (t<ra) smMinSelect(r, t, ra);
     202  kr = kc = 0;
     203  for (j=t-1;j>=0;j--)
     204  {
     205    kr += r[j];
     206    kc += c[j];
     207  }
     208  omFreeSize((ADDRESS)c, al);
     209  omFreeSize((ADDRESS)r, bl);
     210  if (kr<kc) kc = kr;
     211  if (kr<1) kr = 1;
     212  return kr;
     213}
     214
     215static void smMinSelect(Exponent_t *c, int t, int d)
     216{
     217  Exponent_t m;
     218  int pos, i;
     219  do
     220  {
     221    d--;
     222    pos = d;
     223    m = c[pos];
     224    for (i=d-1;i>=0;i--)
     225    {
     226      if(c[i]<m)
     227      {
     228        pos = i;
     229        m = c[i];
     230      }
     231    }
     232    for (i=pos;i<d;i++) c[i] = c[i+1];
     233  } while (d>t);
     234}
     235
    180236/* ----------------- ops with rings ------------------ */
    181 ideal smRingCopy(ideal I, ring *ri, sip_sring &tmpR)
    182 {
    183   ring origR =NULL;
    184   ideal II;
    185   origR =currRing;
    186   tmpR=*origR;
     237void smRingChange(ring *origR, sip_sring &tmpR, Exponent_t bound)
     238{
     239  *origR =currRing;
     240  tmpR=*currRing;
    187241  int *ord=(int*)omAlloc0(3*sizeof(int));
    188242  int *block0=(int*)omAlloc(3*sizeof(int));
    189243  int *block1=(int*)omAlloc(3*sizeof(int));
    190244  ord[0]=ringorder_c;
    191 #ifdef HOMOG_LP
    192   if (idHomIdeal(I))
    193     ord[1] = ringorder_lp;
    194   else
    195 #endif
    196     ord[1]=ringorder_dp;
     245  ord[1]=ringorder_dp;
    197246  tmpR.order=ord;
    198247  tmpR.OrdSgn=1;
    … …  
    201250  block1[1]=tmpR.N;
    202251  tmpR.block1=block1;
    203 
    204   tmpR.bitmask = smExpBound(I);
    205 
    206   // unfortunately, we can not work (yet) with r->N == 0
    207   if (tmpR.bitmask < 1) tmpR.bitmask = 1;
    208   if (tmpR.bitmask > currRing->bitmask) tmpR.bitmask = currRing->bitmask;
     252  tmpR.bitmask = bound;
     253
     254// ???
     255//  if (tmpR.bitmask > currRing->bitmask) tmpR.bitmask = currRing->bitmask;
    209256
    210257  rComplete(&tmpR,1);
    … …  
    212259  if (TEST_OPT_PROT)
    213260    Print("[%d:%d]", (long) tmpR.bitmask, tmpR.ExpL_Size);
    214   // fetch data from the old ring
    215   II = idrCopyR(I, origR);
    216   idTest(II);
    217   *ri = origR;
    218   return II;
    219261}
    220262
    221263void smRingClean(ring origR, ip_sring &tmpR)
    222264{
    223   rChangeCurrRing(origR);
    224265  rUnComplete(&tmpR);
    225266  omFreeSize((ADDRESS)tmpR.order,3*sizeof(int));
    … …  
    228269}
    229270
     271/*2
     272* Bareiss or Chinese remainder ?
     273* I is dXd
     274* sw = TRUE  -> I Matrix
     275*      FALSE -> I Module
     276* return True  -> change Type
     277*        FALSE -> same Type
     278*/
     279BOOLEAN smCheckDet(ideal I, int d, BOOLEAN sw)
     280{
     281  int s,t,i;
     282  poly p;
     283
     284  if ((d>100) || (currRing->parameter!=NULL))
     285    return sw;
     286  if (rField_is_Q(currRing)==FALSE)
     287    return sw;
     288  s = t = 0;
     289  if (sw)
     290  {
     291    for(i=IDELEMS(I)-1;i>=0;i--)
     292    {
     293      p=I->m[i];
     294      if (p!=NULL)
     295      {
     296        if(!pIsConstant(p))
     297          return sw;
     298        s++;
     299        t+=nSize(pGetCoeff(p));
     300      }
     301    }
     302  }
     303  else
     304  {
     305    for(i=IDELEMS(I)-1;i>=0;i--)
     306    {
     307      p=I->m[i];
     308      if (!pIsConstantPoly(p))
     309        return sw;
     310      while (p!=NULL)
     311      {
     312        s++;
     313        t+=nSize(pGetCoeff(p));
     314        pIter(p);
     315      }
     316    }
     317  }
     318  s*=15;
     319  if (t>s)
     320    return !sw;
     321  else
     322    return sw;
     323}
     324
    230325/* ----------------- basics (used from 'C') ------------------ */
    231326/*2
    … …  
    245340    return NULL;
    246341  }
     342  Exponent_t bound=smExpBound(I,r,r,r);
    247343  number diag,h=nInit(1);
    248   poly res,save;
     344  poly res;
    249345  ring origR;
    250346  sip_sring tmpR;
    251347  sparse_mat *det;
    252   ideal II=smRingCopy(I,&origR,tmpR);
    253 
     348  ideal II;
     349
     350  smRingChange(&origR,tmpR,bound);
     351  II = idrCopyR(I, origR);
    254352  diag = smCleardenom(II);
    255353  det = new sparse_mat(II);
    … …  
    258356  {
    259357    delete det;
    260     if (origR!=NULL) smRingClean(origR,tmpR);
     358    rChangeCurrRing(origR);
     359    smRingClean(origR,tmpR);
    261360    return NULL;
    262361  }
    … …  
    264363  if(det->smGetSign()<0) res=pNeg(res);
    265364  delete det;
    266   if (origR!=NULL)
    267   {
    268     rChangeCurrRing(origR);
    269     save = res;
    270     res = prCopyR( save, &tmpR);
    271     rChangeCurrRing(&tmpR);
    272     pDelete(&save);
    273     smRingClean(origR,tmpR);
    274   }
     365  rChangeCurrRing(origR);
     366  res = prMoveR(res, &tmpR);
     367  smRingClean(origR,tmpR);
    275368  if (nEqual(diag,h) == FALSE)
    276369  {
    … …  
    285378lists smCallBareiss(ideal I, int x, int y)
    286379{
     380  lists res=(lists)omAllocBin(slists_bin);
     381  int r=idRankFreeModule(I),t=r;
     382  int c=IDELEMS(I),s=c;
     383  Exponent_t bound;
    287384  ring origR;
    288385  sip_sring tmpR;
    289   lists res=(lists)omAllocBin(slists_bin);
    290   ideal II = smRingCopy(I,&origR,tmpR);
    291   sparse_mat *bareiss = new sparse_mat(II);
     386  ideal II;
     387  sparse_mat *bareiss;
    292388  ideal mm=II;
    293389  intvec *v;
    294   ideal m;
    295 
     390
     391  res->Init(2);
     392  res->m[0].rtyp=MODUL_CMD;
     393  res->m[1].rtyp=INTVEC_CMD;
     394  if ((x>0) && (x<t))
     395    t-=x;
     396  if ((y>1) && (y<s))
     397    s-=y;
     398  if (t>s) t=s;
     399  bound=2*smExpBound(I,c,r,t);
     400  smRingChange(&origR,tmpR,bound);
     401  II = idrCopyR(I, origR);
     402  bareiss = new sparse_mat(II);
    296403  if (bareiss->smGetAct() == NULL)
    297404  {
    298405    delete bareiss;
    299     if (origR!=NULL) smRingClean(origR,tmpR);
    300406    v=new intvec(1,pVariables);
     407    rChangeCurrRing(origR);
    301408  }
    302409  else
    … …  
    304411    idDelete(&II);
    305412    bareiss->smBareiss(x, y);
    306     m = bareiss->smRes2Mod();
     413    II = bareiss->smRes2Mod();
    307414    v = new intvec(bareiss->smGetRed());
    308415    bareiss->smToIntvec(v);
    309416    delete bareiss;
    310     if (origR!=NULL)
    311     {
    312       rChangeCurrRing(origR);
    313       mm=idInit(IDELEMS(m),m->rank);
    314       int k;
    315       for (k=0;k<IDELEMS(m);k++) mm->m[k] = prCopyR( m->m[k], &tmpR);
    316       rChangeCurrRing(&tmpR);
    317       idDelete(&m);
    318       smRingClean(origR,tmpR);
    319     }
    320     else
    321     {
    322       mm=m;
    323     }
    324   }
     417    rChangeCurrRing(origR);
     418    II = idrMoveR(II,&tmpR);
     419  }
     420  smRingClean(origR,tmpR);
     421  res->m[0].data=(void *)II;
     422  res->m[1].data=(void *)v;
     423  return res;
     424}
     425
     426lists smCallNewBareiss(ideal I, int x, int y)
     427{
     428  lists res=(lists)omAllocBin(slists_bin);
     429  int r=idRankFreeModule(I),t=r;
     430  int c=IDELEMS(I),s=c;
     431  Exponent_t bound;
     432  ring origR;
     433  sip_sring tmpR;
     434  ideal II;
     435  sparse_mat *bareiss;
     436  ideal mm=II;
     437  intvec *v;
     438
    325439  res->Init(2);
    326440  res->m[0].rtyp=MODUL_CMD;
    327   res->m[0].data=(void *)mm;
    328441  res->m[1].rtyp=INTVEC_CMD;
    329   res->m[1].data=(void *)v;
    330   return res;
    331 }
    332 
    333 lists smCallNewBareiss(ideal I, int x, int y)
    334 {
    335   ring origR;
    336   sip_sring tmpR;
    337   lists res=(lists)omAllocBin(slists_bin);
    338   ideal II=smRingCopy(I,&origR,tmpR);
    339   sparse_mat *bareiss = new sparse_mat(II);
    340   ideal mm=II;
    341   intvec *v;
    342   ideal m;
    343 
     442  if ((x>0) && (x<t))
     443    t-=x;
     444  if ((y>1) && (y<s))
     445    s-=y;
     446  if (t>s) t=s;
     447  bound=smExpBound(I,c,r,t);
     448  smRingChange(&origR,tmpR,bound);
     449  II = idrCopyR(I, origR);
     450  bareiss = new sparse_mat(II);
    344451  if (bareiss->smGetAct() == NULL)
    345452  {
    346453    delete bareiss;
    347     if (origR!=NULL) smRingClean(origR,tmpR);
    348454    v=new intvec(1,pVariables);
     455    rChangeCurrRing(origR);
    349456  }
    350457  else
    … …  
    352459    idDelete(&II);
    353460    bareiss->smNewBareiss(x, y);
    354     m = bareiss->smRes2Mod();
     461    II = bareiss->smRes2Mod();
    355462    v = new intvec(bareiss->smGetRed());
    356463    bareiss->smToIntvec(v);
    357464    delete bareiss;
    358     if (origR!=NULL)
    359     {
    360       rChangeCurrRing(origR);
    361       mm=idInit(IDELEMS(m),m->rank);
    362       int k;
    363       for (k=0;k<IDELEMS(m);k++) mm->m[k] = prCopyR( m->m[k], &tmpR);
    364       rChangeCurrRing(&tmpR);
    365       idDelete(&m);
    366       smRingClean(origR,tmpR);
    367     }
    368     else
    369     {
    370       mm=m;
    371     }
    372   }
    373   res->Init(2);
    374   res->m[0].rtyp=MODUL_CMD;
    375   res->m[0].data=(void *)mm;
    376   res->m[1].rtyp=INTVEC_CMD;
     465    rChangeCurrRing(origR);
     466    II = idrMoveR(II,&tmpR);
     467  }
     468  smRingClean(origR,tmpR);
     469  res->m[0].data=(void *)II;
    377470  res->m[1].data=(void *)v;
    378471  return res;
    … …  
    24912584{
    24922585  sparse_number_mat *linsolv;
    2493   int k;
    24942586  ring origR;
    24952587  sip_sring tmpR;
    2496   ideal rr, ss;
     2588  ideal rr;
    24972589  lists res;
    24982590
    … …  
    25052597  if (smCheckSolv(I)) return NULL;
    25062598  res=(lists)omAllocBin(slists_bin);
    2507   rr=smRingCopy(I,&origR,tmpR);
    2508   ss = NULL;
     2599  smRingChange(&origR,tmpR,1);
     2600  rr=idrCopyR(I,origR);
    25092601  linsolv = new sparse_number_mat(rr);
     2602  rr=NULL;
    25102603  linsolv->smTriangular();
    25112604  if (linsolv->smIsSing() == 0)
    25122605  {
    25132606    linsolv->smSolv();
    2514     ss = linsolv->smRes2Ideal();
     2607    rr = linsolv->smRes2Ideal();
    25152608  }
    25162609  else
    25172610    WerrorS("singular problem for linsolv");
    25182611  delete linsolv;
    2519   if ((origR!=NULL) && (ss!=NULL))
    2520   {
    2521     rChangeCurrRing(origR);
    2522     rr = idInit(IDELEMS(ss), 1);
    2523     for (k=0;k<IDELEMS(ss);k++)
    2524       rr->m[k] = prCopyR(ss->m[k], &tmpR);
    2525     rChangeCurrRing(&tmpR);
    2526     idDelete(&ss);
    2527     ss = rr;
    2528   }
    2529   if(origR!=NULL)
    2530     smRingClean(origR,tmpR);
     2612  rChangeCurrRing(origR);
     2613  if (rr!=NULL)
     2614    rr = idrMoveR(rr,&tmpR);
     2615  smRingClean(origR,tmpR);
    25312616  res->Init(1);
    25322617  res->m[0].rtyp=IDEAL_CMD;
    2533   res->m[0].data=(void *)ss;
     2618  res->m[0].data=(void *)rr;
    25342619  return res;
    25352620}

  • Singular/sparsmat.h

    r995ba8 rbb3e63b  
    88 *
    99 *******************************************************************/
    10 /* $Id: sparsmat.h,v 1.8 2000-07-06 13:24:21 pohl Exp $ */
     10/* $Id: sparsmat.h,v 1.9 2001-01-30 08:39:29 pohl Exp $ */
    1111
    1212
    … …  
    3131lists smCallSolv(ideal I);
    3232
     33void smRingChange(ring *, sip_sring &, Exponent_t);
     34void smRingClean(ring, ip_sring &);
     35Exponent_t smExpBound(ideal, int, int, int);
     36BOOLEAN smCheckDet(ideal, int, BOOLEAN);
    3337#endif
Note: See TracChangeset for help on using the changeset viewer.