Changeset 232f51c in git


Ignore:
Timestamp:
Feb 24, 1999, 2:25:19 PM (25 years ago)
Author:
Wilfred Pohl <pohl@…>
Branches:
(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
Children:
34cb4bea84226d98334d714b43a4556d8f1c5487
Parents:
d8e380c4561b7282d49138294559898e2ca63423
Message:
bug removed


git-svn-id: file:///usr/local/Singular/svn/trunk@2866 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/sparsmat.cc

    rd8e380 r232f51c  
     1/****************************************
     2*  Computer Algebra System SINGULAR     *
     3****************************************/
     4/* $Id: sparsmat.cc,v 1.2 1999-02-24 13:25:19 pohl Exp $ */
     5
     6/*
     7* ABSTRACT: operations with sparse matrices (bareiss, ...)
     8*/
     9
    110#include "mod2.h"
    211#include "tok.h"
    … …  
    514#include "intvec.h"
    615#include "lists.h"
     16#include "ring.h"
    717#include "polys.h"
    818#include "ipid.h"
    919#include "ideals.h"
    1020#include "numbers.h"
    11 #include "ring.h"
    1221#include "sparsmat.h"
    1322
    1423/* ----------------- macros ------------------ */
    15 #define OLD_DIV
     24/* #define OLD_DIV */
    1625#ifdef OLD_DIV
    1726#define SM_MULT(A,B,C) smMult(A,B)
    … …  
    3544
    3645/* declare internal 'C' stuff */
     46static void smDebug() { int k=0; }
    3747static void smExactPolyDiv(poly, poly);
    3848static BOOLEAN smIsScalar(const poly);
    3949static BOOLEAN smIsNegQuot(poly, const poly, const poly);
    4050static poly smEMult(poly, const poly);
     51static BOOLEAN smCheckLead(const poly, const poly);
    4152static poly smDMult(poly, const poly);
    4253static void smPolyDivN(poly a, const number);
    4354static BOOLEAN smSmaller(poly, poly);
    4455static void smCombineChain(poly *, poly);
     56static void smFindRef(poly *, poly *, poly);
    4557
    4658static void smElemDelete(smpoly *);
    … …  
    89101  void smRowToCol();
    90102  void smFinalMult();
    91   void smFinalToCol();
    92   void smFinalClear();
    93103  void smSparseHomog();
    94104  void smWeights();
    … …  
    109119  ~sparse_mat();
    110120  smpoly * smGetAct() { return m_act; }
     121  int smGetRed() { return tored; }
    111122  ideal smRes2Mod();
    112123  void smBareiss(int, int);
    … …  
    120131{
    121132  lists res=(lists)Alloc(sizeof(slists));
    122   sparse_mat *bareiss = new sparse_mat(I);
     133  ideal II = idCopy(I);
     134  sparse_mat *bareiss = new sparse_mat(II);
    123135  intvec *v;
    124136  ideal m;
    125137
     138  idDelete(&II);
    126139  if (bareiss->smGetAct() == NULL)
    127140  {
    … …  
    131144  bareiss->smBareiss(x, y);
    132145  m = bareiss->smRes2Mod();
    133   v = new intvec(m->rank);
     146  v = new intvec(bareiss->smGetRed());
    134147  bareiss->smToIntvec(v);
    135148  delete bareiss;
    … …  
    177190  ideal m;
    178191
     192  idDelete(&II);
    179193  if (bareiss->smGetAct() == NULL)
    180194  {
    … …  
    192206  bareiss->smNewBareiss(x, y);
    193207  m = bareiss->smRes2Mod();
    194   v = new intvec(m->rank);
     208  v = new intvec(bareiss->smGetRed());
    195209  bareiss->smToIntvec(v);
    196210  delete bareiss;
    … …  
    228242{
    229243  int i;
     244  polyset pmat;
    230245
    231246  ncols = smat->ncols;
    … …  
    243258  i = tored+1;
    244259  perm = (unsigned short *)Alloc(sizeof(unsigned short)*i);
     260  m_row = (smpoly *)Alloc0(sizeof(smpoly)*i);
    245261  wrw = (float *)Alloc(sizeof(float)*i);
    246262  i = ncols+1;
    247263  wcl = (float *)Alloc(sizeof(float)*i);
    248264  m_act = (smpoly *)Alloc(sizeof(smpoly)*i);
    249   m_row = (smpoly *)Alloc0(sizeof(smpoly)*i);
    250265  m_res = (smpoly *)Alloc0(sizeof(smpoly)*i);
    251266  dumm = (smpoly)Alloc(sizeof(smprec));
    252267  m_res[0] = (smpoly)Alloc(sizeof(smprec));
    253268  m_res[0]->m = NULL;
    254   for(i=ncols; i; i--) m_act[i] = smPoly2Smpoly(smat->m[i-1]);
     269  pmat = smat->m;
     270  for(i=ncols; i; i--)
     271  {
     272    m_act[i] = smPoly2Smpoly(pmat[i-1]);
     273    pmat[i-1] = NULL;
     274  }
    255275  this->smZeroElim();
    256276  oldpiv = NULL;
    … …  
    268288  i = ncols+1;
    269289  Free((ADDRESS)m_res, sizeof(smpoly)*i);
    270   Free((ADDRESS)m_row, sizeof(smpoly)*i);
    271290  Free((ADDRESS)m_act, sizeof(smpoly)*i);
    272291  Free((ADDRESS)wcl, sizeof(float)*i);
    273292  i = nrows+1;
    274293  Free((ADDRESS)wrw, sizeof(float)*i);
     294  Free((ADDRESS)m_row, sizeof(smpoly)*i);
    275295  Free((ADDRESS)perm, sizeof(unsigned short)*i);
    276296}
    … …  
    356376        this->smCopToRes();
    357377      }
     378      else
     379        tored = crd;
    358380      return;
    359381    }
    … …  
    419441        this->smCopToRes();
    420442      }
     443      else
     444        tored = crd;
    421445      return;
    422446    }
    … …  
    10611085}
    10621086
    1063 /*
    1064 * store red, the pivot and the assosiated row in m_row
    1065 * to m_res (result):
    1066 *   piv + m_row[rows][pos(cols)] => m_res[cols][pos(rows)]
    1067 */
    1068 void sparse_mat::smFinalToCol()
    1069 {
    1070   smpoly r, a, ap, h;
    1071 
    1072   perm[crd] = rpiv;
    1073   piv->pos = crd;
    1074   ap = m_res[crd];
    1075   if (ap == NULL)
    1076     m_res[crd] = piv;
    1077   else
    1078   {
    1079     loop
    1080     {
    1081       a = ap->n;
    1082       if (a == NULL)
    1083       {
    1084         ap->n = piv;
    1085         break;
    1086       }
    1087       ap = a;
    1088     }
    1089   }
    1090   r = m_row[rpiv];
    1091   m_row[rpiv] = NULL;
    1092   while (r != NULL)
    1093   {
    1094     ap = m_res[r->pos];
    1095     loop
    1096     {
    1097       a = ap->n;
    1098       if (a == NULL)
    1099       {
    1100         ap->n = h = r;
    1101         r = r->n;
    1102         h->n = a;
    1103         h->pos = crd;
    1104         break;
    1105       }
    1106       ap = a;
    1107     }
    1108   }
    1109   r = red;
    1110   while (r != NULL)
    1111   {
    1112     ap = m_res[r->pos];
    1113     if (ap == NULL)
    1114     {
    1115       m_res[r->pos] = h = r;
    1116       r = r->n;
    1117       h->n = NULL;
    1118       h->pos = crd;
     1087/* ----------------- C++ stuff ------------------ */
     1088
     1089/*
     1090*  clean m_act from zeros (after elim)
     1091*/
     1092void sparse_mat::smZeroElim()
     1093{
     1094  int i = 0;
     1095  int j;
     1096
     1097  loop
     1098  {
     1099    i++;
     1100    if (i > act) return;
     1101    if (m_act[i] == NULL) break;
     1102  }
     1103  j = i;
     1104  loop
     1105  {
     1106    j++;
     1107    if (j > act) break;
     1108    if (m_act[j] != NULL)
     1109    {
     1110      m_act[i] = m_act[j];
     1111      i++;
     1112    }
     1113  }
     1114  act -= (j-i);
     1115  sign = 0;
     1116}
     1117
     1118/*
     1119*  clean m_act from cols not to reduced (after elim)
     1120*  put them to m_res
     1121*/
     1122void sparse_mat::smToredElim()
     1123{
     1124  int i = 0;
     1125  int j;
     1126
     1127  loop
     1128  {
     1129    i++;
     1130    if (i > act) return;
     1131    if (m_act[i]->pos > tored)
     1132    {
     1133      m_res[crd] = m_act[i];
     1134      crd++;
     1135      break;
     1136    }
     1137  }
     1138  j = i;
     1139  loop
     1140  {
     1141    j++;
     1142    if (j > act) break;
     1143    if (m_act[i]->pos > tored)
     1144    {
     1145      m_res[crd] = m_act[i];
     1146      crd++;
    11191147    }
    11201148    else
    11211149    {
     1150      m_act[i] = m_act[j];
     1151      i++;
     1152    }
     1153  }
     1154  act -= (j-i);
     1155  sign = 0;
     1156}
     1157
     1158/*
     1159*  copy m_act to m_res
     1160*/
     1161void sparse_mat::smCopToRes()
     1162{
     1163  smpoly p, r, ap, a, h;
     1164  int c, i;
     1165
     1166  if (act == 0) return;
     1167  if (act > 1)
     1168  {
     1169    Werror("Not implemented");
     1170  }
     1171  crd++;
     1172  c = crd;
     1173  p = m_res[c] = m_act[act];
     1174  do
     1175  {
     1176    r = m_row[p->pos];
     1177    m_row[p->pos] = NULL;
     1178    perm[c] = p->pos;
     1179    p->pos = c;
     1180    while (r != NULL)
     1181    {
     1182      ap = m_res[r->pos];
    11221183      loop
    11231184      {
    … …  
    11281189          r = r->n;
    11291190          h->n = a;
    1130           h->pos = crd;
     1191          h->pos = c;
    11311192          break;
    11321193        }
    … …  
    11341195      }
    11351196    }
    1136   }
    1137 }
    1138 
    1139 /*
    1140 * final clear of m_row
    1141 */
    1142 void sparse_mat::smFinalClear()
    1143 {
    1144   int i;
    1145   smpoly r, a, ap, h;
    1146 
    1147   for (i=nrows; i; i--)
    1148   {
    1149     if (m_row[i] != NULL)
    1150     {
    1151       crd++;
    1152       perm[crd] = i;
     1197    c++;
     1198    p = p->n;
     1199  } while (p != NULL);
     1200  for (i=1;i<=nrows;i++)
     1201  {
     1202    if(m_row[i] != NULL)
     1203    {
    11531204      r = m_row[i];
    1154       while (r != NULL)
     1205      m_row[i] = NULL;
     1206      do
    11551207      {
    11561208        ap = m_res[r->pos];
    1157         if (ap == NULL)
    1158         {
    1159           m_res[r->pos] = h = r;
    1160           r = r->n;
    1161           h->n = NULL;
    1162           h->pos = crd;
     1209        loop
     1210        {
     1211          a = ap->n;
     1212          if (a == NULL)
     1213          {
     1214            ap->n = h = r;
     1215            r = r->n;
     1216            h->n = a;
     1217            h->pos = c;
     1218            break;
     1219          }
     1220          ap = a;
    11631221        }
    1164         else
    1165         {
    1166           loop
    1167           {
    1168             a = ap->n;
    1169             if (a == NULL)
    1170             {
    1171               ap->n = h = r;
    1172               r = r->n;
    1173               h->n = a;
    1174               h->pos = crd;
    1175               break;
    1176             }
    1177             ap = a;
    1178           }
    1179         }
    1180       }
    1181     }
    1182   }
    1183 }
    1184 
    1185 /* ----------------- C++ stuff ------------------ */
    1186 
    1187 /*
    1188 *  clean m_act from zeros (after elim)
    1189 */
    1190 void sparse_mat::smZeroElim()
    1191 {
    1192   int i = 0;
    1193   int j;
    1194 
    1195   loop
    1196   {
    1197     i++;
    1198     if (i > act) return;
    1199     if (m_act[i] == NULL) break;
    1200   }
    1201   j = i;
    1202   loop
    1203   {
    1204     j++;
    1205     if (j > act) break;
    1206     if (m_act[j] != NULL)
    1207     {
    1208       m_act[i] = m_act[j];
    1209       i++;
    1210     }
    1211   }
    1212   act -= (j-i);
    1213   sign = 0;
    1214 }
    1215 
    1216 /*
    1217 *  clean m_act from cols not to reduced (after elim)
    1218 *  put them to m_res
    1219 */
    1220 void sparse_mat::smToredElim()
    1221 {
    1222   int i = 0;
    1223   int j;
    1224 
    1225   loop
    1226   {
    1227     i++;
    1228     if (i > act) return;
    1229     if (m_act[i]->pos > tored)
    1230     {
    1231       m_res[crd] = m_act[i];
    1232       crd++;
    1233       break;
    1234     }
    1235   }
    1236   j = i;
    1237   loop
    1238   {
    1239     j++;
    1240     if (j > act) break;
    1241     if (m_act[i]->pos > tored)
    1242     {
    1243       m_res[crd] = m_act[i];
    1244       crd++;
    1245     }
    1246     else
    1247     {
    1248       m_act[i] = m_act[j];
    1249       i++;
    1250     }
    1251   }
    1252   act -= (j-i);
    1253   sign = 0;
    1254 }
    1255 
    1256 /*
    1257 *  copy m_act to m_res
    1258 */
    1259 void sparse_mat::smCopToRes()
    1260 {
    1261   while (act != 0)
    1262   {
    1263     rpiv = m_act[act]->pos;
    1264     this->smSelectPR();
    1265     crd++;
    1266     this->smColToRow();
    1267     act--;
    1268     this->smFinalToCol();
    1269     if (act != 0) this->smZeroElim();
    1270   }
    1271 //  this->smFinalClear();
     1222      } while (r!=NULL);
     1223      c++;
     1224    }
     1225  }
     1226  tored = c-1;
    12721227}
    12731228
    … …  
    14001355        pSetm(a);
    14011356      }
    1402       y = nIntDiv(pGetCoeff(a),x);
     1357      y = nDiv(pGetCoeff(a),x);
     1358      nNormalize(y);
    14031359      pSetCoeff(a,y);
    14041360      pIter(a);
    … …  
    14121368      pSubExp(a,i,pGetExp(b,i));
    14131369    pSetm(a);
    1414     y = nIntDiv(pGetCoeff(a),x);
     1370    y = nDiv(pGetCoeff(a),x);
     1371    nNormalize(y);
    14151372    pSetCoeff(a,y);
    14161373    yn = nNeg(nCopy(y));
    … …  
    14401397{
    14411398  poly pa, e, res, r;
     1399  BOOLEAN lead;
    14421400
    14431401  if (smSmaller(a, b))
    … …  
    14661424    return res;
    14671425  }
     1426  res = NULL;
    14681427  e = pNew();
    1469   loop
     1428  lead = FALSE;
     1429  while (!lead)
    14701430  {
    14711431    pSetCoeff0(e,pGetCoeff(b));
    14721432    if (smIsNegQuot(e, b, c))
    1473       res = smDMult(a, e);
     1433    {
     1434      lead = smCheckLead(a, e);
     1435      r = smDMult(a, e);
     1436    }
    14741437    else
    1475       res = smEMult(a, e);
     1438    {
     1439      lead = TRUE;
     1440      r = smEMult(a, e);
     1441    }
     1442    if (lead)
     1443    {
     1444      if (res != NULL)
     1445      {
     1446        smFindRef(&pa, &res, r);
     1447        if (pa == NULL)
     1448          lead = FALSE;
     1449      }
     1450      else
     1451      {
     1452        pa = res = r;
     1453      }
     1454    }
     1455    else
     1456      res = pAdd(res, r);
    14761457    pIter(b);
    14771458    if (b == NULL)
    … …  
    14801461      return res;
    14811462    }
    1482     if (res != NULL) break;
    14831463  }
    14841464  do
    … …  
    14881468    {
    14891469      r = smDMult(a, e);
    1490       if (r != NULL)
    1491         res = pAdd(res,r);
     1470      if (smCheckLead(a, e))
     1471        smCombineChain(&pa, r);
     1472      else
     1473        pa = pAdd(pa,r);
    14921474    }
    14931475    else
    14941476    {
    14951477      r = smEMult(a, e);
    1496       res = pAdd(res,r);
     1478      smCombineChain(&pa, r);
    14971479    }
    14981480    pIter(b);
    … …  
    15271509  do
    15281510  {
    1529     y = nIntDiv(pGetCoeff(a), x);
     1511    y = nDiv(pGetCoeff(a), x);
     1512    nNormalize(y);
    15301513    pSetCoeff(a,y);
    15311514    yn = nNeg(nCopy(y));
    … …  
    15921575}
    15931576
     1577static BOOLEAN smCheckLead(const poly t, const poly e)
     1578{
     1579  int i;
     1580  for (i=pVariables; i; i--)
     1581  {
     1582    if ((pGetExp(e,i)+pGetExp(t,i)) < 0)
     1583      return FALSE;
     1584  }
     1585  return TRUE;
     1586}
     1587
    15941588static poly smDMult(poly t, const poly e)
    15951589{
    … …  
    16421636  do
    16431637  {
    1644     y = nIntDiv(pGetCoeff(a),x);
     1638    y = nDiv(pGetCoeff(a),x);
     1639    nNormalize(y);
    16451640    pSetCoeff(a,y);
    16461641    pIter(a);
    … …  
    17051700}
    17061701
     1702static void smFindRef(poly *ref, poly *px, poly r)
     1703{
     1704  number x;
     1705  int i;
     1706  poly pa = *px, pp = NULL;
     1707
     1708  loop
     1709  {
     1710    i = pComp0(pa, r);
     1711    if (i > 0)
     1712    {
     1713      pp = pa;
     1714      pIter(pa);
     1715      if (pa==NULL)
     1716      {
     1717        pNext(pp) = r;
     1718        break;
     1719      }
     1720    }
     1721    else
     1722    {
     1723      if (i == 0)
     1724      {
     1725        x = nAdd(pGetCoeff(pa), pGetCoeff(r));
     1726        pDelete1(&r);
     1727        if (nIsZero(x))
     1728        {
     1729          pDelete1(&pa);
     1730          if (pp!=NULL)
     1731            pNext(pp) = pAdd(pa,r);
     1732          else
     1733            *px = pAdd(pa,r);
     1734        }
     1735        else
     1736        {
     1737          pp = pa;
     1738          pSetCoeff(pp,x);
     1739          pNext(pp) = pAdd(pNext(pp), r);
     1740        }
     1741      }
     1742      else
     1743      {
     1744        if (pp!=NULL)
     1745          pp = pNext(pp) = r;
     1746        else
     1747          *px = pp = r;
     1748        pNext(pp) = pAdd(pa, pNext(r));
     1749      }
     1750      break;
     1751    }
     1752  }
     1753  *ref = pp;
     1754}
     1755
    17071756/* ----------------- internal 'C' stuff ------------------ */
    17081757
    … …  
    17291778* do not destroy p
    17301779*/
    1731 static smpoly smPoly2Smpoly(poly p)
    1732 {
    1733   poly pp, q;
     1780static smpoly smPoly2Smpoly(poly q)
     1781{
     1782  poly pp;
    17341783  smpoly res, a;
    17351784  EXPONENT_TYPE x;
    17361785
    1737   if (p == NULL)
     1786  if (q == NULL)
    17381787    return NULL;
    1739   q = pCopy(p);
    17401788  a = res = (smpoly)Alloc(sizeof(smprec));
    17411789  a->pos = x = pGetComp(q);
Note: See TracChangeset for help on using the changeset viewer.