Changeset ebdaa1 in git for libpolys/polys/sparsmat.cc
 Timestamp:
 Nov 9, 2011, 4:24:50 PM (12 years ago)
 Branches:
 (u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
 Children:
 b2929e71ee342c64f5189c222795ac906cb21ee9
 Parents:
 8f40312428677d88ab37a31824f2ec53f144a18a
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

libpolys/polys/sparsmat.cc
r8f4031 rebdaa1 60 60 61 61 62 /* declare internal 'C' stuff */ 63 static void sm_ExactPolyDiv(poly, poly, const ring); 64 static BOOLEAN sm_IsNegQuot(poly, const poly, const poly, const ring); 65 static void sm_ExpMultDiv(poly, const poly, const poly, const ring); 66 static void sm_PolyDivN(poly, const number, const ring); 67 static BOOLEAN smSmaller(poly, poly); 68 static void sm_CombineChain(poly *, poly, const ring); 69 static void sm_FindRef(poly *, poly *, poly, const ring); 70 71 static void sm_ElemDelete(smpoly *, const ring); 72 static smpoly smElemCopy(smpoly); 73 static float sm_PolyWeight(smpoly, const ring); 74 static smpoly sm_Poly2Smpoly(poly, const ring); 75 static poly sm_Smpoly2Poly(smpoly, const ring); 76 static BOOLEAN sm_HaveDenom(poly, const ring); 77 static number sm_Cleardenom(ideal, const ring); 78 62 79 omBin smprec_bin = omGetSpecBin(sizeof(smprec)); 63 64 static void smExpMultDiv(poly t, const poly b, const poly c)65 {66 int i;67 pTest(t);68 pLmTest(b);69 pLmTest(c);70 poly bc = p_New(currRing);71 72 p_ExpVectorDiff(bc, b, c, currRing);73 74 while(t!=NULL)75 {76 pExpVectorAdd(t, bc);77 pIter(t);78 }79 p_LmFree(bc, currRing);80 }81 80 82 81 static poly pp_Mult_Coeff_mm_DivSelect_MultDiv(poly p, int &lp, poly m, … … 353 352 354 353 /*  basics (used from 'C')  */ 355 static BOOLEAN smHaveDenom(poly a)356 {357 BOOLEAN sw;358 number x;359 360 while (a != NULL)361 {362 x = nGetDenom(pGetCoeff(a));363 sw = nIsOne(x);364 nDelete(&x);365 if (!sw)366 {367 return TRUE;368 }369 pIter(a);370 }371 return FALSE;372 }373 374 static number smCleardenom(ideal id)375 {376 poly a;377 number x,y,res=nInit(1);378 BOOLEAN sw=FALSE;379 380 for (int i=0; i<IDELEMS(id); i++)381 {382 a = id>m[i];383 sw = smHaveDenom(a);384 if (sw) break;385 }386 if (!sw) return res;387 for (int i=0; i<IDELEMS(id); i++)388 {389 a = id>m[i];390 if (a!=NULL)391 {392 x = nCopy(pGetCoeff(a));393 p_Cleardenom(a, currRing);394 y = nDiv(x,pGetCoeff(a));395 nDelete(&x);396 x = nMult(res,y);397 nNormalize(x);398 nDelete(&res);399 res = x;400 }401 }402 return res;403 }404 405 354 /*2 406 355 *returns the determinant of the module I; … … 486 435 sm_KillModifiedRing(tmpR); 487 436 M=II; 488 }489 490 /*491 * from poly to smpoly492 * do not destroy p493 */494 static smpoly smPoly2Smpoly(poly q)495 {496 poly pp;497 smpoly res, a;498 long x;499 500 if (q == NULL)501 return NULL;502 a = res = (smpoly)omAllocBin(smprec_bin);503 a>pos = x = pGetComp(q);504 a>m = q;505 a>e = 0;506 loop507 {508 pSetComp(q,0);509 pp = q;510 pIter(q);511 if (q == NULL)512 {513 a>n = NULL;514 return res;515 }516 if (pGetComp(q) != x)517 {518 a = a>n = (smpoly)omAllocBin(smprec_bin);519 pNext(pp) = NULL;520 a>pos = x = pGetComp(q);521 a>m = q;522 a>e = 0;523 }524 }525 437 } 526 438 … … 587 499 588 500 /* 589 * from smpoly to poly590 * destroy a591 */592 static poly smSmpoly2Poly(smpoly a)593 {594 smpoly b;595 poly res, pp, q;596 long x;597 598 if (a == NULL)599 return NULL;600 x = a>pos;601 q = res = a>m;602 loop603 {604 pSetComp(q,x);605 pp = q;606 pIter(q);607 if (q == NULL)608 break;609 }610 loop611 {612 b = a;613 a = a>n;614 omFreeBin((void *)b, smprec_bin);615 if (a == NULL)616 return res;617 x = a>pos;618 q = pNext(pp) = a>m;619 loop620 {621 pSetComp(q,x);622 pp = q;623 pIter(q);624 if (q == NULL)625 break;626 }627 }628 }629 630 /*631 501 * transform the result to a module 632 502 */ … … 789 659 790 660 /* 791 * weigth of a polynomial, for pivot strategy792 */793 static float smPolyWeight(smpoly a)794 {795 poly p = a>m;796 int i;797 float res = (float)nSize(pGetCoeff(p));798 799 if (pNext(p) == NULL)800 {801 for(i=pVariables; i>0; i)802 {803 if (pGetExp(p,i) != 0) return res+1.0;804 }805 return res;806 }807 else808 {809 i = 0;810 res = 0.0;811 do812 {813 i++;814 res += (float)nSize(pGetCoeff(p));815 pIter(p);816 }817 while (p);818 return res+(float)i;819 }820 }821 822 /*823 661 * prepare smPivot, compute weights for rows and columns 824 662 * and the weight for all points … … 1011 849 1012 850 /*  elimination  */ 1013 static void smElemDelete(smpoly *r)1014 {1015 smpoly a = *r, b = a>n;1016 1017 pDelete(&a>m);1018 omFreeBin((void *)a, smprec_bin);1019 *r = b;1020 }1021 1022 static smpoly smElemCopy(smpoly a)1023 {1024 smpoly r = (smpoly)omAllocBin(smprec_bin);1025 memcpy(r, a, sizeof(smprec));1026 /* r>m = pCopy(r>m); */1027 return r;1028 }1029 851 1030 852 /* first step of elimination */ … … 1843 1665 1844 1666 1845 static BOOLEAN smIsNegQuot(poly a, const poly b, const poly c)1846 {1847 if (pLmDivisibleByNoComp(c, b))1848 {1849 pExpVectorDiff(a, b, c);1850 // Hmm: here used to be a pSetm(a): but it is unnecessary,1851 // if b and c are correct1852 return FALSE;1853 }1854 else1855 {1856 int i;1857 for (i=pVariables; i>0; i)1858 {1859 if(pGetExp(c,i) > pGetExp(b,i))1860 pSetExp(a,i,pGetExp(c,i)pGetExp(b,i));1861 else1862 pSetExp(a,i,0);1863 }1864 // here we actually might need a pSetm, if a is to be used in1865 // comparisons1866 return TRUE;1867 }1868 }1869 1870 static BOOLEAN smSmaller(poly a, poly b)1871 {1872 loop1873 {1874 pIter(b);1875 if (b == NULL) return TRUE;1876 pIter(a);1877 if (a == NULL) return FALSE;1878 }1879 }1880 1881 static void smCombineChain(poly *px, poly r)1882 {1883 poly pa = *px, pb;1884 number x;1885 int i;1886 1887 loop1888 {1889 pb = pNext(pa);1890 if (pb == NULL)1891 {1892 pa = pNext(pa) = r;1893 break;1894 }1895 i = pLmCmp(pb, r);1896 if (i > 0)1897 pa = pb;1898 else1899 {1900 if (i == 0)1901 {1902 x = nAdd(pGetCoeff(pb), pGetCoeff(r));1903 pLmDelete(&r);1904 if (nIsZero(x))1905 {1906 pLmDelete(&pb);1907 pNext(pa) = pAdd(pb,r);1908 }1909 else1910 {1911 pa = pb;1912 pSetCoeff(pa,x);1913 pNext(pa) = pAdd(pNext(pa), r);1914 }1915 }1916 else1917 {1918 pa = pNext(pa) = r;1919 pNext(pa) = pAdd(pb, pNext(pa));1920 }1921 break;1922 }1923 }1924 *px = pa;1925 }1926 1927 static void smFindRef(poly *ref, poly *px, poly r)1928 {1929 number x;1930 int i;1931 poly pa = *px, pp = NULL;1932 1933 loop1934 {1935 i = pLmCmp(pa, r);1936 if (i > 0)1937 {1938 pp = pa;1939 pIter(pa);1940 if (pa==NULL)1941 {1942 pNext(pp) = r;1943 break;1944 }1945 }1946 else1947 {1948 if (i == 0)1949 {1950 x = nAdd(pGetCoeff(pa), pGetCoeff(r));1951 pLmDelete(&r);1952 if (nIsZero(x))1953 {1954 pLmDelete(&pa);1955 if (pp!=NULL)1956 pNext(pp) = pAdd(pa,r);1957 else1958 *px = pAdd(pa,r);1959 }1960 else1961 {1962 pp = pa;1963 pSetCoeff(pp,x);1964 pNext(pp) = pAdd(pNext(pp), r);1965 }1966 }1967 else1968 {1969 if (pp!=NULL)1970 pp = pNext(pp) = r;1971 else1972 *px = pp = r;1973 pNext(pp) = pAdd(pa, pNext(r));1974 }1975 break;1976 }1977 }1978 *ref = pp;1979 }1980 1981 1667 //disable that, as it fails with coef buckets 1982 1668 //#define X_MAS … … 1994 1680 if ((c == NULL)  pLmIsConstantComp(c)) 1995 1681 { 1996 return pp Mult_qq(a, b);1682 return pp_Mult_qq(a, b); 1997 1683 } 1998 1684 … … 2618 2304 } 2619 2305 } 2620 smExactPolyDiv(a, b);2306 return res; 2621 2307 } 2622 2308
Note: See TracChangeset
for help on using the changeset viewer.