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