[35aab3] | 1 | /**************************************** |
---|
| 2 | * Computer Algebra System SINGULAR * |
---|
| 3 | ****************************************/ |
---|
[341696] | 4 | /* $Id$ */ |
---|
[35aab3] | 5 | |
---|
| 6 | /* |
---|
| 7 | * ABSTRACT: |
---|
| 8 | */ |
---|
| 9 | |
---|
| 10 | #include <math.h> |
---|
[599326] | 11 | #include <kernel/mod2.h> |
---|
| 12 | #include <kernel/options.h> |
---|
[b1dfaf] | 13 | #include <omalloc/omalloc.h> |
---|
[599326] | 14 | #include <kernel/polys.h> |
---|
| 15 | #include <kernel/intvec.h> |
---|
| 16 | #include <kernel/febase.h> |
---|
| 17 | #include <kernel/ideals.h> |
---|
| 18 | #include <kernel/ring.h> |
---|
| 19 | #include <kernel/weight.h> |
---|
[35aab3] | 20 | |
---|
| 21 | /*0 implementation*/ |
---|
| 22 | extern "C" double (*wFunctional)(int *degw, int *lpol, int npol, |
---|
[6aacb6] | 23 | double *rel, double wx, double wNsqr); |
---|
[35aab3] | 24 | extern "C" double wFunctionalMora(int *degw, int *lpol, int npol, |
---|
[6aacb6] | 25 | double *rel, double wx, double wNsqr); |
---|
[35aab3] | 26 | extern "C" double wFunctionalBuch(int *degw, int *lpol, int npol, |
---|
[6aacb6] | 27 | double *rel, double wx, double wNsqr); |
---|
[35aab3] | 28 | extern "C" void wAdd(int *A, int mons, int kn, int xx); |
---|
| 29 | extern "C" void wNorm(int *degw, int *lpol, int npol, double *rel); |
---|
| 30 | extern "C" void wFirstSearch(int *A, int *x, int mons, |
---|
[6aacb6] | 31 | int *lpol, int npol, double *rel, double *fopt, double wNsqr); |
---|
[35aab3] | 32 | extern "C" void wSecondSearch(int *A, int *x, int *lpol, |
---|
[6aacb6] | 33 | int npol, int mons, double *rel, double *fk, double wNsqr); |
---|
[35aab3] | 34 | extern "C" void wGcd(int *x, int n); |
---|
| 35 | |
---|
| 36 | static void wDimensions(polyset s, int sl, int *lpol, int *npol, int *mons) |
---|
| 37 | { |
---|
| 38 | int i, i1, j, k; |
---|
| 39 | poly p, q; |
---|
| 40 | |
---|
| 41 | i1 = j = 0; |
---|
| 42 | for (i = 0; i <= sl; i++) |
---|
| 43 | { |
---|
| 44 | p = s[i]; |
---|
| 45 | if (p!=NULL) |
---|
| 46 | { |
---|
| 47 | k = 1; |
---|
| 48 | q = pNext(p); |
---|
| 49 | while (q!=NULL) |
---|
| 50 | { |
---|
| 51 | k++; |
---|
| 52 | q = pNext(q); |
---|
| 53 | } |
---|
| 54 | if (k > 1) |
---|
| 55 | { |
---|
| 56 | lpol[i1++] = k; |
---|
| 57 | j += k; |
---|
| 58 | } |
---|
| 59 | } |
---|
| 60 | } |
---|
| 61 | *npol = i1; |
---|
| 62 | *mons = j; |
---|
| 63 | } |
---|
| 64 | |
---|
| 65 | |
---|
| 66 | static void wInit(polyset s, int sl, int mons, int *A) |
---|
| 67 | { |
---|
| 68 | int n, a, i, j, *B, *C; |
---|
| 69 | poly p, q; |
---|
| 70 | int *pl; |
---|
| 71 | |
---|
| 72 | B = A; |
---|
| 73 | n = pVariables; |
---|
| 74 | a = (n + 1) * sizeof(int); |
---|
| 75 | pl = (int *)omAlloc(a); |
---|
| 76 | for (i = 0; i <= sl; i++) |
---|
| 77 | { |
---|
| 78 | p = s[i]; |
---|
| 79 | if (p!=NULL) |
---|
| 80 | { |
---|
| 81 | q = pNext(p); |
---|
| 82 | if (q!=NULL) |
---|
| 83 | { |
---|
| 84 | C = B; |
---|
| 85 | B++; |
---|
| 86 | pGetExpV(p, pl); |
---|
| 87 | for (j = 0; j < n; j++) |
---|
| 88 | { |
---|
| 89 | *C = pl[j+1]; |
---|
| 90 | C += mons; |
---|
| 91 | } |
---|
| 92 | } |
---|
| 93 | while (q!=NULL) |
---|
| 94 | { |
---|
| 95 | C = B; |
---|
| 96 | B++; |
---|
| 97 | pGetExpV(q, pl); |
---|
| 98 | for (j = 0; j < n; j++) |
---|
| 99 | { |
---|
| 100 | *C = pl[j+1]; |
---|
| 101 | C += mons; |
---|
| 102 | } |
---|
| 103 | pIter(q); |
---|
| 104 | } |
---|
| 105 | } |
---|
| 106 | } |
---|
| 107 | omFreeSize((ADDRESS)pl, a); |
---|
| 108 | } |
---|
| 109 | |
---|
[6aacb6] | 110 | void wCall(polyset s, int sl, int *x, double wNsqr) |
---|
[35aab3] | 111 | { |
---|
| 112 | int n, q, npol, mons, i; |
---|
| 113 | int *A, *xopt, *lpol, *degw; |
---|
| 114 | double f1, fx, eps, *rel; |
---|
| 115 | void *adr; |
---|
| 116 | |
---|
| 117 | n = pVariables; |
---|
| 118 | lpol = (int * )omAlloc((sl + 1) * sizeof(int)); |
---|
| 119 | wDimensions(s, sl, lpol, &npol, &mons); |
---|
| 120 | xopt = x + (n + 1); |
---|
| 121 | for (i = n; i!=0; i--) |
---|
| 122 | xopt[i] = 1; |
---|
| 123 | if (mons==0) |
---|
| 124 | { |
---|
| 125 | omFreeSize((ADDRESS)lpol, (sl + 1) * sizeof(int)); |
---|
| 126 | return; |
---|
| 127 | } |
---|
| 128 | adr = (void * )omAllocAligned(npol * sizeof(double)); |
---|
| 129 | rel = (double*)adr; |
---|
| 130 | q = (n + 1) * mons * sizeof(int); |
---|
| 131 | A = (int * )omAlloc(q); |
---|
| 132 | wInit(s, sl, mons, A); |
---|
| 133 | degw = A + (n * mons); |
---|
| 134 | memset(degw, 0, mons * sizeof(int)); |
---|
| 135 | for (i = n; i!=0; i--) |
---|
| 136 | wAdd(A, mons, i, 1); |
---|
| 137 | wNorm(degw, lpol, npol, rel); |
---|
[6aacb6] | 138 | f1 = (*wFunctional)(degw, lpol, npol, rel, (double)1.0, wNsqr); |
---|
[35aab3] | 139 | if (TEST_OPT_PROT) Print("// %e\n",f1); |
---|
| 140 | eps = f1; |
---|
| 141 | fx = (double)2.0 * eps; |
---|
| 142 | memset(x, 0, (n + 1) * sizeof(int)); |
---|
[6aacb6] | 143 | wFirstSearch(A, x, mons, lpol, npol, rel, &fx, wNsqr); |
---|
[35aab3] | 144 | if (TEST_OPT_PROT) Print("// %e\n",fx); |
---|
| 145 | memcpy(x + 1, xopt + 1, n * sizeof(int)); |
---|
| 146 | memset(degw, 0, mons * sizeof(int)); |
---|
| 147 | for (i = n; i!=0; i--) |
---|
| 148 | { |
---|
| 149 | x[i] *= 16; |
---|
| 150 | wAdd(A, mons, i, x[i]); |
---|
| 151 | } |
---|
[6aacb6] | 152 | wSecondSearch(A, x, lpol, npol, mons, rel, &fx, wNsqr); |
---|
[35aab3] | 153 | if (TEST_OPT_PROT) Print("// %e\n",fx); |
---|
| 154 | if (fx >= eps) |
---|
| 155 | { |
---|
| 156 | for (i = n; i!=0; i--) |
---|
| 157 | xopt[i] = 1; |
---|
| 158 | } |
---|
| 159 | else |
---|
| 160 | { |
---|
| 161 | wGcd(xopt, n); |
---|
| 162 | // if (BTEST1(22)) |
---|
| 163 | // { |
---|
| 164 | // f1 = fx + (double)0.1 * (f1 - fx); |
---|
| 165 | // wSimple(x, n); |
---|
| 166 | // memset(degw, 0, mons * sizeof(int)); |
---|
| 167 | // for (i = n; i!=0; i--) |
---|
| 168 | // wAdd(A, mons, i, x[i]); |
---|
| 169 | // eps = wPrWeight(x, n); |
---|
| 170 | // fx = (*wFunctional)(degw, lpol, npol, rel, eps); |
---|
| 171 | // if (fx < f1) |
---|
| 172 | // { |
---|
| 173 | // if (TEST_OPT_PROT) Print("// %e\n",fx); |
---|
| 174 | // memcpy(xopt + 1, x + 1, n * sizeof(int)); |
---|
| 175 | // } |
---|
| 176 | // } |
---|
| 177 | } |
---|
| 178 | omFreeSize((ADDRESS)A, q); |
---|
| 179 | omFreeSize((ADDRESS)lpol, (sl + 1) * sizeof(int)); |
---|
| 180 | omFreeSize((ADDRESS)adr, npol * sizeof(double)); |
---|
| 181 | } |
---|
| 182 | |
---|
| 183 | |
---|
| 184 | void kEcartWeights(polyset s, int sl, short *eweight) |
---|
| 185 | { |
---|
| 186 | int n, i; |
---|
| 187 | int *x; |
---|
| 188 | |
---|
| 189 | *eweight = 0; |
---|
| 190 | n = pVariables; |
---|
[753f9d2] | 191 | if (rHasLocalOrMixedOrdering_currRing()) |
---|
[35aab3] | 192 | wFunctional = wFunctionalMora; |
---|
| 193 | else |
---|
| 194 | wFunctional = wFunctionalBuch; |
---|
| 195 | x = (int * )omAlloc(2 * (n + 1) * sizeof(int)); |
---|
[6aacb6] | 196 | wCall(s, sl, x, (double)2.0 / (double)n); |
---|
[35aab3] | 197 | for (i = n; i!=0; i--) |
---|
| 198 | eweight[i] = x[i + n + 1]; |
---|
| 199 | omFreeSize((ADDRESS)x, 2 * (n + 1) * sizeof(int)); |
---|
| 200 | } |
---|
| 201 | |
---|
| 202 | short * iv2array(intvec * iv) |
---|
| 203 | { |
---|
| 204 | short *s=(short *)omAlloc0((pVariables+1)*sizeof(short)); |
---|
| 205 | int len=0; |
---|
| 206 | if(iv!=NULL) |
---|
[48c2df] | 207 | len=si_min(iv->length(),pVariables); // usually: pVariables==length() |
---|
[35aab3] | 208 | int i; |
---|
| 209 | //for(i=pVariables;i>len;i--) s[i]=1; |
---|
| 210 | for(i=len;i>0;i--) s[i]=(*iv)[i-1]; |
---|
| 211 | return s; |
---|
| 212 | } |
---|
| 213 | |
---|
| 214 | /*2 |
---|
| 215 | *computes the degree of the leading term of the polynomial |
---|
| 216 | *with respect to given ecartWeights |
---|
| 217 | *used for Graebes method if BTEST1(31) is set |
---|
| 218 | */ |
---|
| 219 | long totaldegreeWecart(poly p, ring r) |
---|
| 220 | { |
---|
| 221 | int i; |
---|
| 222 | long j =0; |
---|
| 223 | |
---|
[5f4c317] | 224 | for (i=rVar(r); i>0; i--) |
---|
[35aab3] | 225 | j += (int)(p_GetExp(p,i,r) * ecartWeights[i]); |
---|
| 226 | return j; |
---|
| 227 | } |
---|
| 228 | |
---|
| 229 | /*2 |
---|
| 230 | *computes the degree of the leading term of the polynomial |
---|
| 231 | *with respect to given weights |
---|
| 232 | */ |
---|
| 233 | long totaldegreeWecart_IV(poly p, ring r, const short *w) |
---|
| 234 | { |
---|
| 235 | int i; |
---|
| 236 | long j =0; |
---|
| 237 | |
---|
[5f4c317] | 238 | for (i=rVar(r); i>0; i--) |
---|
[db1fa7] | 239 | j += (long)((int)(p_GetExp(p,i,r) * w[i])); |
---|
[35aab3] | 240 | return j; |
---|
| 241 | } |
---|
| 242 | |
---|
| 243 | /*2 |
---|
| 244 | *computes the maximal degree of all terms of the polynomial |
---|
| 245 | *with respect to given ecartWeights and |
---|
| 246 | *computes the length of the polynomial |
---|
| 247 | *used for Graebes method if BTEST1(31) is set |
---|
| 248 | */ |
---|
| 249 | long maxdegreeWecart(poly p,int *l, ring r) |
---|
| 250 | { |
---|
| 251 | short k=p_GetComp(p, r); |
---|
| 252 | int ll=1; |
---|
| 253 | long t,max; |
---|
| 254 | |
---|
| 255 | max=totaldegreeWecart(p, r); |
---|
| 256 | pIter(p); |
---|
| 257 | while ((p!=NULL) && (p_GetComp(p, r)==k)) |
---|
| 258 | { |
---|
| 259 | t=totaldegreeWecart(p, r); |
---|
| 260 | if (t>max) max=t; |
---|
| 261 | ll++; |
---|
| 262 | pIter(p); |
---|
| 263 | } |
---|
| 264 | *l=ll; |
---|
| 265 | return max; |
---|
| 266 | } |
---|