Changeset a24e0a in git
 Timestamp:
 Oct 18, 2011, 2:06:44 PM (12 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '0604212ebb110535022efecad887940825b97c3f')
 Children:
 48c2dff0cd110870e54b65568f912ab293a5948b
 Parents:
 de71f84a62bde9a34b88ab2bc6de78a67ea87957
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

kernel/sparsmat.cc
rde71f84 ra24e0a 46 46 47 47 48 /* declare internal 'C' stuff */49 static void smExactPolyDiv(poly, poly);50 static BOOLEAN smIsNegQuot(poly, const poly, const poly);51 static void smExpMultDiv(poly, const poly, const poly);52 static void smPolyDivN(poly, const number);53 static BOOLEAN smSmaller(poly, poly);54 static void smCombineChain(poly *, poly);55 static void smFindRef(poly *, poly *, poly);56 57 static void smElemDelete(smpoly *);58 static smpoly smElemCopy(smpoly);59 static float smPolyWeight(smpoly);60 static smpoly smPoly2Smpoly(poly);61 static poly smSmpoly2Poly(smpoly);62 static BOOLEAN smHaveDenom(poly);63 static number smCleardenom(ideal);64 65 48 omBin smprec_bin = omGetSpecBin(sizeof(smprec)); 49 50 static void smExpMultDiv(poly t, const poly b, const poly c) 51 { 52 int i; 53 pTest(t); 54 pLmTest(b); 55 pLmTest(c); 56 poly bc = p_New(currRing); 57 58 p_ExpVectorDiff(bc, b, c, currRing); 59 60 while(t!=NULL) 61 { 62 pExpVectorAdd(t, bc); 63 pIter(t); 64 } 65 p_LmFree(bc, currRing); 66 } 66 67 67 68 static poly pp_Mult_Coeff_mm_DivSelect_MultDiv(poly p, int &lp, poly m, … … 336 337 337 338 /*  basics (used from 'C')  */ 339 static BOOLEAN smHaveDenom(poly a) 340 { 341 BOOLEAN sw; 342 number x; 343 344 while (a != NULL) 345 { 346 x = nGetDenom(pGetCoeff(a)); 347 sw = nIsOne(x); 348 nDelete(&x); 349 if (!sw) 350 { 351 return TRUE; 352 } 353 pIter(a); 354 } 355 return FALSE; 356 } 357 358 static number smCleardenom(ideal id) 359 { 360 poly a; 361 number x,y,res=nInit(1); 362 BOOLEAN sw=FALSE; 363 364 for (int i=0; i<IDELEMS(id); i++) 365 { 366 a = id>m[i]; 367 sw = smHaveDenom(a); 368 if (sw) break; 369 } 370 if (!sw) return res; 371 for (int i=0; i<IDELEMS(id); i++) 372 { 373 a = id>m[i]; 374 if (a!=NULL) 375 { 376 x = nCopy(pGetCoeff(a)); 377 p_Cleardenom(a, currRing); 378 y = nDiv(x,pGetCoeff(a)); 379 nDelete(&x); 380 x = nMult(res,y); 381 nNormalize(x); 382 nDelete(&res); 383 res = x; 384 } 385 } 386 return res; 387 } 388 338 389 /*2 339 390 *returns the determinant of the module I; … … 425 476 smKillModifiedRing(tmpR); 426 477 M=II; 478 } 479 480 /* 481 * from poly to smpoly 482 * do not destroy p 483 */ 484 static smpoly smPoly2Smpoly(poly q) 485 { 486 poly pp; 487 smpoly res, a; 488 long x; 489 490 if (q == NULL) 491 return NULL; 492 a = res = (smpoly)omAllocBin(smprec_bin); 493 a>pos = x = pGetComp(q); 494 a>m = q; 495 a>e = 0; 496 loop 497 { 498 pSetComp(q,0); 499 pp = q; 500 pIter(q); 501 if (q == NULL) 502 { 503 a>n = NULL; 504 return res; 505 } 506 if (pGetComp(q) != x) 507 { 508 a = a>n = (smpoly)omAllocBin(smprec_bin); 509 pNext(pp) = NULL; 510 a>pos = x = pGetComp(q); 511 a>m = q; 512 a>e = 0; 513 } 514 } 427 515 } 428 516 … … 488 576 489 577 /* 578 * from smpoly to poly 579 * destroy a 580 */ 581 static poly smSmpoly2Poly(smpoly a) 582 { 583 smpoly b; 584 poly res, pp, q; 585 long x; 586 587 if (a == NULL) 588 return NULL; 589 x = a>pos; 590 q = res = a>m; 591 loop 592 { 593 pSetComp(q,x); 594 pp = q; 595 pIter(q); 596 if (q == NULL) 597 break; 598 } 599 loop 600 { 601 b = a; 602 a = a>n; 603 omFreeBin((void *)b, smprec_bin); 604 if (a == NULL) 605 return res; 606 x = a>pos; 607 q = pNext(pp) = a>m; 608 loop 609 { 610 pSetComp(q,x); 611 pp = q; 612 pIter(q); 613 if (q == NULL) 614 break; 615 } 616 } 617 } 618 619 /* 490 620 * transform the result to a module 491 621 */ … … 647 777 648 778 /* 779 * weigth of a polynomial, for pivot strategy 780 */ 781 static float smPolyWeight(smpoly a) 782 { 783 poly p = a>m; 784 int i; 785 float res = (float)nSize(pGetCoeff(p)); 786 787 if (pNext(p) == NULL) 788 { 789 for(i=pVariables; i>0; i) 790 { 791 if (pGetExp(p,i) != 0) return res+1.0; 792 } 793 return res; 794 } 795 else 796 { 797 i = 0; 798 res = 0.0; 799 do 800 { 801 i++; 802 res += (float)nSize(pGetCoeff(p)); 803 pIter(p); 804 } 805 while (p); 806 return res+(float)i; 807 } 808 } 809 810 /* 649 811 * prepare smPivot, compute weights for rows and columns 650 812 * and the weight for all points … … 837 999 838 1000 /*  elimination  */ 1001 static void smElemDelete(smpoly *r) 1002 { 1003 smpoly a = *r, b = a>n; 1004 1005 pDelete(&a>m); 1006 omFreeBin((void *)a, smprec_bin); 1007 *r = b; 1008 } 1009 1010 static smpoly smElemCopy(smpoly a) 1011 { 1012 smpoly r = (smpoly)omAllocBin(smprec_bin); 1013 memcpy(r, a, sizeof(smprec)); 1014 /* r>m = pCopy(r>m); */ 1015 return r; 1016 } 839 1017 840 1018 /* first step of elimination */ … … 1653 1831 1654 1832 1833 static BOOLEAN smIsNegQuot(poly a, const poly b, const poly c) 1834 { 1835 if (pLmDivisibleByNoComp(c, b)) 1836 { 1837 pExpVectorDiff(a, b, c); 1838 // Hmm: here used to be a pSetm(a): but it is unnecessary, 1839 // if b and c are correct 1840 return FALSE; 1841 } 1842 else 1843 { 1844 int i; 1845 for (i=pVariables; i>0; i) 1846 { 1847 if(pGetExp(c,i) > pGetExp(b,i)) 1848 pSetExp(a,i,pGetExp(c,i)pGetExp(b,i)); 1849 else 1850 pSetExp(a,i,0); 1851 } 1852 // here we actually might need a pSetm, if a is to be used in 1853 // comparisons 1854 return TRUE; 1855 } 1856 } 1857 1858 static BOOLEAN smSmaller(poly a, poly b) 1859 { 1860 loop 1861 { 1862 pIter(b); 1863 if (b == NULL) return TRUE; 1864 pIter(a); 1865 if (a == NULL) return FALSE; 1866 } 1867 } 1868 1869 static void smCombineChain(poly *px, poly r) 1870 { 1871 poly pa = *px, pb; 1872 number x; 1873 int i; 1874 1875 loop 1876 { 1877 pb = pNext(pa); 1878 if (pb == NULL) 1879 { 1880 pa = pNext(pa) = r; 1881 break; 1882 } 1883 i = pLmCmp(pb, r); 1884 if (i > 0) 1885 pa = pb; 1886 else 1887 { 1888 if (i == 0) 1889 { 1890 x = nAdd(pGetCoeff(pb), pGetCoeff(r)); 1891 pLmDelete(&r); 1892 if (nIsZero(x)) 1893 { 1894 pLmDelete(&pb); 1895 pNext(pa) = pAdd(pb,r); 1896 } 1897 else 1898 { 1899 pa = pb; 1900 pSetCoeff(pa,x); 1901 pNext(pa) = pAdd(pNext(pa), r); 1902 } 1903 } 1904 else 1905 { 1906 pa = pNext(pa) = r; 1907 pNext(pa) = pAdd(pb, pNext(pa)); 1908 } 1909 break; 1910 } 1911 } 1912 *px = pa; 1913 } 1914 1915 static void smFindRef(poly *ref, poly *px, poly r) 1916 { 1917 number x; 1918 int i; 1919 poly pa = *px, pp = NULL; 1920 1921 loop 1922 { 1923 i = pLmCmp(pa, r); 1924 if (i > 0) 1925 { 1926 pp = pa; 1927 pIter(pa); 1928 if (pa==NULL) 1929 { 1930 pNext(pp) = r; 1931 break; 1932 } 1933 } 1934 else 1935 { 1936 if (i == 0) 1937 { 1938 x = nAdd(pGetCoeff(pa), pGetCoeff(r)); 1939 pLmDelete(&r); 1940 if (nIsZero(x)) 1941 { 1942 pLmDelete(&pa); 1943 if (pp!=NULL) 1944 pNext(pp) = pAdd(pa,r); 1945 else 1946 *px = pAdd(pa,r); 1947 } 1948 else 1949 { 1950 pp = pa; 1951 pSetCoeff(pp,x); 1952 pNext(pp) = pAdd(pNext(pp), r); 1953 } 1954 } 1955 else 1956 { 1957 if (pp!=NULL) 1958 pp = pNext(pp) = r; 1959 else 1960 *px = pp = r; 1961 pNext(pp) = pAdd(pa, pNext(r)); 1962 } 1963 break; 1964 } 1965 } 1966 *ref = pp; 1967 } 1968 1655 1969 //disable that, as it fails with coef buckets 1656 1970 //#define X_MAS … … 1874 2188 } 1875 2189 #endif 1876 /*n1877 * exact division a/b1878 * a is a result of smMultDiv1879 * a destroyed, b NOT destroyed1880 */1881 void smSpecialPolyDiv(poly a, poly b)1882 {1883 if (pNext(b) == NULL)1884 {1885 smPolyDivN(a, pGetCoeff(b));1886 return;1887 }1888 smExactPolyDiv(a, b);1889 }1890 1891 1892 2190 /*  internals arithmetic  */ 1893 2191 static void smExactPolyDiv(poly a, poly b) … … 1945 2243 } 1946 2244 1947 // obachman > Wilfried: check the following1948 static BOOLEAN smIsNegQuot(poly a, const poly b, const poly c)1949 {1950 if (pLmDivisibleByNoComp(c, b))1951 {1952 pExpVectorDiff(a, b, c);1953 // Hmm: here used to be a pSetm(a): but it is unnecessary,1954 // if b and c are correct1955 return FALSE;1956 }1957 else1958 {1959 int i;1960 for (i=pVariables; i>0; i)1961 {1962 if(pGetExp(c,i) > pGetExp(b,i))1963 pSetExp(a,i,pGetExp(c,i)pGetExp(b,i));1964 else1965 pSetExp(a,i,0);1966 }1967 // here we actually might need a pSetm, if a is to be used in1968 // comparisons1969 return TRUE;1970 }1971 }1972 1973 static void smExpMultDiv(poly t, const poly b, const poly c)1974 {1975 int i;1976 pTest(t);1977 pLmTest(b);1978 pLmTest(c);1979 poly bc = p_New(currRing);1980 1981 p_ExpVectorDiff(bc, b, c, currRing);1982 1983 while(t!=NULL)1984 {1985 pExpVectorAdd(t, bc);1986 pIter(t);1987 }1988 p_LmFree(bc, currRing);1989 }1990 1991 1992 2245 static void smPolyDivN(poly a, const number x) 1993 2246 { … … 2003 2256 } 2004 2257 2005 static BOOLEAN smSmaller(poly a, poly b) 2006 { 2007 loop 2008 { 2009 pIter(b); 2010 if (b == NULL) return TRUE; 2011 pIter(a); 2012 if (a == NULL) return FALSE; 2013 } 2014 } 2015 2016 static void smCombineChain(poly *px, poly r) 2017 { 2018 poly pa = *px, pb; 2019 number x; 2020 int i; 2021 2022 loop 2023 { 2024 pb = pNext(pa); 2025 if (pb == NULL) 2026 { 2027 pa = pNext(pa) = r; 2028 break; 2029 } 2030 i = pLmCmp(pb, r); 2031 if (i > 0) 2032 pa = pb; 2033 else 2034 { 2035 if (i == 0) 2036 { 2037 x = nAdd(pGetCoeff(pb), pGetCoeff(r)); 2038 pLmDelete(&r); 2039 if (nIsZero(x)) 2040 { 2041 pLmDelete(&pb); 2042 pNext(pa) = pAdd(pb,r); 2043 } 2044 else 2045 { 2046 pa = pb; 2047 pSetCoeff(pa,x); 2048 pNext(pa) = pAdd(pNext(pa), r); 2049 } 2050 } 2051 else 2052 { 2053 pa = pNext(pa) = r; 2054 pNext(pa) = pAdd(pb, pNext(pa)); 2055 } 2056 break; 2057 } 2058 } 2059 *px = pa; 2060 } 2061 2062 2063 static void smFindRef(poly *ref, poly *px, poly r) 2064 { 2065 number x; 2066 int i; 2067 poly pa = *px, pp = NULL; 2068 2069 loop 2070 { 2071 i = pLmCmp(pa, r); 2072 if (i > 0) 2073 { 2074 pp = pa; 2075 pIter(pa); 2076 if (pa==NULL) 2077 { 2078 pNext(pp) = r; 2079 break; 2080 } 2081 } 2082 else 2083 { 2084 if (i == 0) 2085 { 2086 x = nAdd(pGetCoeff(pa), pGetCoeff(r)); 2087 pLmDelete(&r); 2088 if (nIsZero(x)) 2089 { 2090 pLmDelete(&pa); 2091 if (pp!=NULL) 2092 pNext(pp) = pAdd(pa,r); 2093 else 2094 *px = pAdd(pa,r); 2095 } 2096 else 2097 { 2098 pp = pa; 2099 pSetCoeff(pp,x); 2100 pNext(pp) = pAdd(pNext(pp), r); 2101 } 2102 } 2103 else 2104 { 2105 if (pp!=NULL) 2106 pp = pNext(pp) = r; 2107 else 2108 *px = pp = r; 2109 pNext(pp) = pAdd(pa, pNext(r)); 2110 } 2111 break; 2112 } 2113 } 2114 *ref = pp; 2115 } 2116 2117 /*  internal 'C' stuff  */ 2118 2119 static void smElemDelete(smpoly *r) 2120 { 2121 smpoly a = *r, b = a>n; 2122 2123 pDelete(&a>m); 2124 omFreeBin((void *)a, smprec_bin); 2125 *r = b; 2126 } 2127 2128 static smpoly smElemCopy(smpoly a) 2129 { 2130 smpoly r = (smpoly)omAllocBin(smprec_bin); 2131 memcpy(r, a, sizeof(smprec)); 2132 /* r>m = pCopy(r>m); */ 2133 return r; 2134 } 2135 2136 /* 2137 * from poly to smpoly 2138 * do not destroy p 2139 */ 2140 static smpoly smPoly2Smpoly(poly q) 2141 { 2142 poly pp; 2143 smpoly res, a; 2144 long x; 2145 2146 if (q == NULL) 2147 return NULL; 2148 a = res = (smpoly)omAllocBin(smprec_bin); 2149 a>pos = x = pGetComp(q); 2150 a>m = q; 2151 a>e = 0; 2152 loop 2153 { 2154 pSetComp(q,0); 2155 pp = q; 2156 pIter(q); 2157 if (q == NULL) 2158 { 2159 a>n = NULL; 2160 return res; 2161 } 2162 if (pGetComp(q) != x) 2163 { 2164 a = a>n = (smpoly)omAllocBin(smprec_bin); 2165 pNext(pp) = NULL; 2166 a>pos = x = pGetComp(q); 2167 a>m = q; 2168 a>e = 0; 2169 } 2170 } 2171 } 2172 2173 /* 2174 * from smpoly to poly 2175 * destroy a 2176 */ 2177 static poly smSmpoly2Poly(smpoly a) 2178 { 2179 smpoly b; 2180 poly res, pp, q; 2181 long x; 2182 2183 if (a == NULL) 2184 return NULL; 2185 x = a>pos; 2186 q = res = a>m; 2187 loop 2188 { 2189 pSetComp(q,x); 2190 pp = q; 2191 pIter(q); 2192 if (q == NULL) 2193 break; 2194 } 2195 loop 2196 { 2197 b = a; 2198 a = a>n; 2199 omFreeBin((void *)b, smprec_bin); 2200 if (a == NULL) 2201 return res; 2202 x = a>pos; 2203 q = pNext(pp) = a>m; 2204 loop 2205 { 2206 pSetComp(q,x); 2207 pp = q; 2208 pIter(q); 2209 if (q == NULL) 2210 break; 2211 } 2212 } 2213 } 2214 2215 /* 2216 * weigth of a polynomial, for pivot strategy 2217 */ 2218 static float smPolyWeight(smpoly a) 2219 { 2220 poly p = a>m; 2221 int i; 2222 float res = (float)nSize(pGetCoeff(p)); 2223 2224 if (pNext(p) == NULL) 2225 { 2226 for(i=pVariables; i>0; i) 2227 { 2228 if (pGetExp(p,i) != 0) return res+1.0; 2229 } 2230 return res; 2231 } 2232 else 2233 { 2234 i = 0; 2235 res = 0.0; 2236 do 2237 { 2238 i++; 2239 res += (float)nSize(pGetCoeff(p)); 2240 pIter(p); 2241 } 2242 while (p); 2243 return res+(float)i; 2244 } 2245 } 2246 2247 static BOOLEAN smHaveDenom(poly a) 2248 { 2249 BOOLEAN sw; 2250 number x; 2251 2252 while (a != NULL) 2253 { 2254 x = nGetDenom(pGetCoeff(a)); 2255 sw = nIsOne(x); 2256 nDelete(&x); 2257 if (!sw) 2258 { 2259 return TRUE; 2260 } 2261 pIter(a); 2262 } 2263 return FALSE; 2264 } 2265 2266 static number smCleardenom(ideal id) 2267 { 2268 poly a; 2269 number x,y,res=nInit(1); 2270 BOOLEAN sw=FALSE; 2271 2272 for (int i=0; i<IDELEMS(id); i++) 2273 { 2274 a = id>m[i]; 2275 sw = smHaveDenom(a); 2276 if (sw) break; 2277 } 2278 if (!sw) return res; 2279 for (int i=0; i<IDELEMS(id); i++) 2280 { 2281 a = id>m[i]; 2282 if (a!=NULL) 2283 { 2284 x = nCopy(pGetCoeff(a)); 2285 p_Cleardenom(a, currRing); 2286 y = nDiv(x,pGetCoeff(a)); 2287 nDelete(&x); 2288 x = nMult(res,y); 2289 nNormalize(x); 2290 nDelete(&res); 2291 res = x; 2292 } 2293 } 2294 return res; 2258 /*n 2259 * exact division a/b 2260 * a is a result of smMultDiv 2261 * a destroyed, b NOT destroyed 2262 */ 2263 void smSpecialPolyDiv(poly a, poly b) 2264 { 2265 if (pNext(b) == NULL) 2266 { 2267 smPolyDivN(a, pGetCoeff(b)); 2268 return; 2269 } 2270 smExactPolyDiv(a, b); 2295 2271 } 2296 2272
Note: See TracChangeset
for help on using the changeset viewer.