[35aab3] | 1 | /**************************************** |
---|
| 2 | * Computer Algebra System SINGULAR * |
---|
| 3 | ****************************************/ |
---|
| 4 | |
---|
| 5 | /* |
---|
| 6 | * ABSTRACT: |
---|
| 7 | */ |
---|
| 8 | |
---|
| 9 | #include <math.h> |
---|
| 10 | #include <string.h> |
---|
[1135a61] | 11 | #include <misc/auxiliary.h> |
---|
[b1dfaf] | 12 | #include <omalloc/omalloc.h> |
---|
[35aab3] | 13 | |
---|
| 14 | double wFunctionalMora(int *degw, int *lpol, int npol, |
---|
[6aacb6] | 15 | double *rel, double wx, double wwNsqr); |
---|
[35aab3] | 16 | double wFunctionalBuch(int *degw, int *lpol, int npol, |
---|
[6aacb6] | 17 | double *rel, double wx, double wwNsqr); |
---|
[1135a61] | 18 | void wAdd(int *A, int mons, int kn, int xx, int rvar); |
---|
[35aab3] | 19 | void wNorm(int *degw, int *lpol, int npol, double *rel); |
---|
| 20 | void wFirstSearch(int *A, int *x, int mons, |
---|
[1135a61] | 21 | int *lpol, int npol, double *rel, double *fopt, double wNsqr, int rvar); |
---|
[35aab3] | 22 | void wSecondSearch(int *A, int *x, int *lpol, |
---|
[1135a61] | 23 | int npol, int mons, double *rel, double *fk, double wNsqr, int rvar); |
---|
[35aab3] | 24 | void wGcd(int *x, int n); |
---|
| 25 | /*0 implementation*/ |
---|
| 26 | |
---|
| 27 | short * ecartWeights=NULL; |
---|
| 28 | |
---|
| 29 | double (*wFunctional)(int *degw, int *lpol, int npol, |
---|
[6aacb6] | 30 | double *rel, double wx, double wNsqr); |
---|
[35aab3] | 31 | |
---|
| 32 | |
---|
| 33 | double wFunctionalMora(int *degw, int *lpol, int npol, |
---|
[6aacb6] | 34 | double *rel, double wx, double wNsqr) |
---|
[35aab3] | 35 | { |
---|
| 36 | int i, j, e1, ecu, ecl, ec; |
---|
| 37 | int *ex; |
---|
| 38 | double gfmax, gecart, ghom, pfmax; |
---|
| 39 | double *r; |
---|
| 40 | |
---|
| 41 | ex = degw; |
---|
| 42 | r = rel; |
---|
| 43 | gfmax = (double)0.0; |
---|
| 44 | gecart = (double)0.4 + (double)npol; |
---|
| 45 | ghom = (double)1.0; |
---|
| 46 | for (i = 0; i < npol; i++) |
---|
| 47 | { |
---|
| 48 | ecl = ecu = e1 = *ex++; |
---|
| 49 | for (j = lpol[i] - 1; j!=0; j--) |
---|
| 50 | { |
---|
| 51 | ec = *ex++; |
---|
| 52 | if (ec > ecu) |
---|
| 53 | ecu = ec; |
---|
| 54 | else if (ec < ecl) |
---|
| 55 | ecl = ec; |
---|
| 56 | } |
---|
| 57 | pfmax = (double)ecl / (double)ecu; |
---|
| 58 | if (pfmax < ghom) |
---|
| 59 | ghom = pfmax; |
---|
| 60 | pfmax = (double)e1 / (double)ecu; |
---|
| 61 | if (pfmax > (double)0.5) |
---|
| 62 | gecart -= (pfmax * pfmax); |
---|
| 63 | else |
---|
| 64 | gecart -= (double)0.25; |
---|
| 65 | ecu = 2 * ecu - ecl; |
---|
| 66 | gfmax += (double)(ecu * ecu) * (*r++); |
---|
| 67 | } |
---|
| 68 | if (ghom > (double)0.8) |
---|
| 69 | { |
---|
| 70 | ghom *= (double)5.0; |
---|
| 71 | gecart *= ((double)5.0 - ghom); |
---|
| 72 | } |
---|
| 73 | return (gfmax * gecart) / pow(wx, wNsqr); |
---|
| 74 | } |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | double wFunctionalBuch(int *degw, int *lpol, int npol, |
---|
[6aacb6] | 78 | double *rel, double wx, double wNsqr) |
---|
[35aab3] | 79 | { |
---|
| 80 | int i, j, ecl, ecu, ec; |
---|
| 81 | int *ex; |
---|
| 82 | double gfmax, ghom, pfmax; |
---|
| 83 | double *r; |
---|
| 84 | |
---|
| 85 | ex = degw; |
---|
| 86 | r = rel; |
---|
| 87 | gfmax = (double)0.0; |
---|
| 88 | ghom = (double)1.0; |
---|
| 89 | for (i = 0; i < npol; i++) |
---|
| 90 | { |
---|
| 91 | ecu = ecl = *ex++; |
---|
| 92 | for (j = lpol[i] - 1; j!=0 ; j--) |
---|
| 93 | { |
---|
| 94 | ec = *ex++; |
---|
| 95 | if (ec < ecl) |
---|
| 96 | ecl = ec; |
---|
| 97 | else if (ec > ecu) |
---|
| 98 | ecu = ec; |
---|
| 99 | } |
---|
| 100 | pfmax = (double)ecl / (double)ecu; |
---|
| 101 | if (pfmax < ghom) |
---|
| 102 | ghom = pfmax; |
---|
| 103 | gfmax += (double)(ecu * ecu) * (*r++); |
---|
| 104 | } |
---|
| 105 | if (ghom > (double)0.5) |
---|
| 106 | gfmax *= ((double)1.0 - (ghom * ghom)) / (double)0.75; |
---|
| 107 | return gfmax / pow(wx, wNsqr); |
---|
| 108 | } |
---|
| 109 | |
---|
| 110 | |
---|
[1135a61] | 111 | static void wSub(int *A, int mons, int kn, int xx,int rvar) |
---|
[35aab3] | 112 | { |
---|
| 113 | int i, *B, *ex; |
---|
| 114 | |
---|
| 115 | B = A + ((kn - 1) * mons); |
---|
[1135a61] | 116 | ex = A + (rvar * mons); |
---|
[35aab3] | 117 | i = mons; |
---|
| 118 | if (xx == 1) |
---|
| 119 | { |
---|
| 120 | for (/* i=mons */; i!=0 ; i--) |
---|
| 121 | *ex++ -= *B++; |
---|
| 122 | } |
---|
| 123 | else |
---|
| 124 | { |
---|
| 125 | for (/* i=mons */; i!=0 ; i--) |
---|
| 126 | *ex++ -= (*B++) * xx; |
---|
| 127 | } |
---|
| 128 | } |
---|
| 129 | |
---|
| 130 | |
---|
[1135a61] | 131 | void wAdd(int *A, int mons, int kn, int xx, int rvar) |
---|
[35aab3] | 132 | { |
---|
| 133 | int i, *B, *ex; |
---|
| 134 | |
---|
| 135 | B = A + ((kn - 1) * mons); |
---|
[1135a61] | 136 | ex = A + (rvar * mons); |
---|
[35aab3] | 137 | i = mons; |
---|
| 138 | if (xx == 1) |
---|
| 139 | { |
---|
| 140 | for (/* i=mons */; i!=0 ; i--) |
---|
| 141 | *ex++ += *B++; |
---|
| 142 | } |
---|
| 143 | else |
---|
| 144 | { |
---|
| 145 | for (/* i=mons */; i!=0 ; i--) |
---|
| 146 | *ex++ += (*B++) * xx; |
---|
| 147 | } |
---|
| 148 | } |
---|
| 149 | |
---|
| 150 | |
---|
| 151 | void wFirstSearch(int *A, int *x, int mons, |
---|
[1135a61] | 152 | int *lpol, int npol, double *rel, double *fopt, double wNsqr, int rvar) |
---|
[35aab3] | 153 | { |
---|
| 154 | int a0, a, n, xn, t, xx, y1; |
---|
| 155 | int *y, *degw, *xopt; |
---|
| 156 | double fy, fmax, wx; |
---|
| 157 | double *pr; |
---|
| 158 | void *adr; |
---|
| 159 | |
---|
| 160 | fy = *fopt; |
---|
[1135a61] | 161 | n = rvar; |
---|
[35aab3] | 162 | xn = n + 6 + (21 / n); |
---|
| 163 | a0 = n * sizeof(double); |
---|
| 164 | a = n * sizeof(int); |
---|
[b21eae] | 165 | y = (int * )omAlloc((long)a); |
---|
| 166 | adr = (void * )omAllocAligned((long)a0); |
---|
[35aab3] | 167 | pr = adr; |
---|
| 168 | *pr = (double)1.0; |
---|
| 169 | *y = 0; |
---|
| 170 | degw = A + (n * mons); |
---|
| 171 | xopt = x + (n + 2); |
---|
| 172 | t = 1; |
---|
| 173 | loop |
---|
| 174 | { |
---|
| 175 | while (t < n) |
---|
| 176 | { |
---|
| 177 | xx = x[t] + 1; |
---|
| 178 | wx = pr[t-1] * (double)xx; |
---|
| 179 | y1 = y[t-1] + xx; |
---|
| 180 | if ((y1 + n - t) <= xn) |
---|
| 181 | { |
---|
| 182 | pr[t] = wx; |
---|
| 183 | y[t] = y1; |
---|
| 184 | x[t] = xx; |
---|
| 185 | if (xx > 1) |
---|
[1135a61] | 186 | wAdd(A, mons, t, 1, rvar); |
---|
[35aab3] | 187 | t++; |
---|
| 188 | } |
---|
| 189 | else |
---|
| 190 | { |
---|
| 191 | xx = x[t] - 1; |
---|
| 192 | x[t] = 0; |
---|
| 193 | if (xx!=0) |
---|
[1135a61] | 194 | wSub(A, mons, t, xx, rvar); |
---|
[35aab3] | 195 | t--; |
---|
| 196 | if (t==0) |
---|
| 197 | { |
---|
| 198 | *fopt = fy; |
---|
[b21eae] | 199 | omFreeSize((ADDRESS)y, (long)a); |
---|
| 200 | omFreeSize((ADDRESS)adr, (long)a0); |
---|
[35aab3] | 201 | return; |
---|
| 202 | } |
---|
| 203 | } |
---|
| 204 | } |
---|
| 205 | xx = xn - y[t-1]; |
---|
| 206 | wx = pr[t-1] * (double)xx; |
---|
| 207 | x[t] = xx; |
---|
| 208 | xx--; |
---|
| 209 | if (xx!=0) |
---|
[1135a61] | 210 | wAdd(A, mons, t, xx, rvar); |
---|
[6aacb6] | 211 | fmax = (*wFunctional)(degw, lpol, npol, rel, wx,wNsqr); |
---|
[35aab3] | 212 | if (xx!=0) |
---|
[1135a61] | 213 | wSub(A, mons, t, xx, rvar); |
---|
[35aab3] | 214 | if (fmax < fy) |
---|
| 215 | { |
---|
| 216 | fy = fmax; |
---|
[611871] | 217 | memcpy(xopt, x + 1, a); |
---|
[35aab3] | 218 | } |
---|
| 219 | t--; |
---|
| 220 | } /* end loop */ |
---|
| 221 | } |
---|
| 222 | |
---|
| 223 | |
---|
| 224 | static double wPrWeight(int *x, int n) |
---|
| 225 | { |
---|
| 226 | int i; |
---|
| 227 | double y; |
---|
| 228 | |
---|
| 229 | y = (double)x[n]; |
---|
| 230 | for (i = n - 1; i!=0 ; i--) |
---|
| 231 | y *= (double)x[i]; |
---|
| 232 | return y; |
---|
| 233 | } |
---|
| 234 | |
---|
| 235 | |
---|
| 236 | static void wEstimate(int *A, int *x, int *lpol, int npol, int mons, |
---|
[1135a61] | 237 | double wx, double *rel, double *fopt, int *s0, int *s1, int *s2, double wNsqr, int rvar) |
---|
[35aab3] | 238 | { |
---|
| 239 | int n, i1, i2, k0 = 0, k1 = 0, k2 = 0; |
---|
| 240 | int *degw; |
---|
| 241 | double fo1, fo2, fmax, wx1, wx2; |
---|
| 242 | |
---|
[1135a61] | 243 | n = rvar; |
---|
[35aab3] | 244 | degw = A + (n * mons); |
---|
| 245 | fo2 = fo1 = (double)1.0e10; |
---|
| 246 | for (i1 = n; i1!=0 ; i1--) |
---|
| 247 | { |
---|
| 248 | if (x[i1] > 1) |
---|
| 249 | { |
---|
[1135a61] | 250 | wSub(A, mons, i1, 1, rvar); |
---|
[35aab3] | 251 | wx1 = wx - wx / (double)x[i1]; |
---|
| 252 | x[i1]--; |
---|
[6aacb6] | 253 | fmax = (*wFunctional)(degw, lpol, npol, rel, wx1,wNsqr); |
---|
[35aab3] | 254 | if (fmax < fo1) |
---|
| 255 | { |
---|
| 256 | fo1 = fmax; |
---|
| 257 | k0 = i1; |
---|
| 258 | } |
---|
| 259 | for (i2 = i1; i2!=0 ; i2--) |
---|
| 260 | { |
---|
| 261 | if (x[i2] > 1) |
---|
| 262 | { |
---|
[1135a61] | 263 | wSub(A, mons, i2, 1, rvar); |
---|
[35aab3] | 264 | wx2 = wx1 - wx1 / (double)x[i2]; |
---|
[6aacb6] | 265 | fmax = (*wFunctional)(degw, lpol, npol, rel, wx2, wNsqr); |
---|
[35aab3] | 266 | if (fmax < fo2) |
---|
| 267 | { |
---|
| 268 | fo2 = fmax; |
---|
| 269 | k1 = i1; |
---|
| 270 | k2 = i2; |
---|
| 271 | } |
---|
[1135a61] | 272 | wAdd(A, mons, i2, 1, rvar); |
---|
[35aab3] | 273 | } |
---|
| 274 | } |
---|
[1135a61] | 275 | wAdd(A, mons, i1, 1, rvar); |
---|
[35aab3] | 276 | x[i1]++; |
---|
| 277 | } |
---|
| 278 | } |
---|
| 279 | if (fo1 < fo2) |
---|
| 280 | { |
---|
| 281 | *fopt = fo1; |
---|
| 282 | *s0 = k0; |
---|
| 283 | } |
---|
| 284 | else |
---|
| 285 | { |
---|
| 286 | *fopt = fo2; |
---|
| 287 | *s0 = 0; |
---|
| 288 | } |
---|
| 289 | *s1 = k1; |
---|
| 290 | *s2 = k2; |
---|
| 291 | } |
---|
| 292 | |
---|
| 293 | |
---|
| 294 | void wSecondSearch(int *A, int *x, int *lpol, |
---|
[1135a61] | 295 | int npol, int mons, double *rel, double *fk, double wNsqr, int rvar) |
---|
[35aab3] | 296 | { |
---|
| 297 | int n, s0, s1, s2, *xopt; |
---|
| 298 | double one, fx, fopt, wx; |
---|
| 299 | |
---|
[1135a61] | 300 | n = rvar; |
---|
[35aab3] | 301 | xopt = x + (n + 2); |
---|
| 302 | fopt = *fk * (double)0.999999999999; |
---|
| 303 | wx = wPrWeight(x, n); |
---|
| 304 | one = (double)1.0; |
---|
| 305 | loop |
---|
| 306 | { |
---|
[1135a61] | 307 | wEstimate(A, x, lpol, npol, mons, wx, rel, &fx, &s0, &s1, &s2, wNsqr, rvar); |
---|
[35aab3] | 308 | if (fx > fopt) |
---|
| 309 | { |
---|
| 310 | if (s0!=0) |
---|
| 311 | x[s0]--; |
---|
| 312 | else if (s1!=0) |
---|
| 313 | { |
---|
| 314 | x[s1]--; |
---|
| 315 | x[s2]--; |
---|
| 316 | } |
---|
| 317 | else |
---|
| 318 | break; |
---|
| 319 | } |
---|
| 320 | else |
---|
| 321 | { |
---|
| 322 | fopt = fx; |
---|
| 323 | if (s0!=0) |
---|
| 324 | { |
---|
| 325 | x[s0]--; |
---|
[611871] | 326 | memcpy(xopt, x + 1, n * sizeof(int)); |
---|
[35aab3] | 327 | if (s1==0) |
---|
| 328 | break; |
---|
| 329 | } |
---|
| 330 | else if (s1!=0) |
---|
| 331 | { |
---|
| 332 | x[s1]--; |
---|
| 333 | x[s2]--; |
---|
[611871] | 334 | memcpy(xopt, x + 1, n * sizeof(int)); |
---|
[35aab3] | 335 | } |
---|
| 336 | else |
---|
| 337 | break; |
---|
| 338 | } |
---|
| 339 | if (s0!=0) |
---|
[1135a61] | 340 | wSub(A, mons, s0, 1, rvar); |
---|
[35aab3] | 341 | else |
---|
| 342 | { |
---|
[1135a61] | 343 | wSub(A, mons, s1, 1, rvar); |
---|
| 344 | wSub(A, mons, s2, 1, rvar); |
---|
[35aab3] | 345 | } |
---|
| 346 | wx = wPrWeight(x, n); |
---|
| 347 | } |
---|
| 348 | *fk = fopt; |
---|
| 349 | } |
---|
| 350 | |
---|
| 351 | |
---|
| 352 | void wGcd(int *x, int n) |
---|
| 353 | { |
---|
| 354 | int i, b, a, h; |
---|
| 355 | |
---|
| 356 | i = n; |
---|
| 357 | b = x[i]; |
---|
| 358 | loop |
---|
| 359 | { |
---|
| 360 | i--; |
---|
| 361 | if (i==0) |
---|
| 362 | break; |
---|
| 363 | a = x[i]; |
---|
| 364 | if (a < b) |
---|
| 365 | { |
---|
| 366 | h = a; |
---|
| 367 | a = b; |
---|
| 368 | b = h; |
---|
| 369 | } |
---|
| 370 | do |
---|
| 371 | { |
---|
| 372 | h = a % b; |
---|
| 373 | a = b; |
---|
| 374 | b = h; |
---|
| 375 | } |
---|
| 376 | while (b!=0); |
---|
| 377 | b = a; |
---|
| 378 | if (b == 1) |
---|
| 379 | return; |
---|
| 380 | } |
---|
| 381 | for (i = n; i!=0 ; i--) |
---|
| 382 | x[i] /= b; |
---|
| 383 | } |
---|
| 384 | |
---|
| 385 | |
---|
| 386 | static void wSimple(int *x, int n) |
---|
| 387 | { |
---|
| 388 | int g, min, c, d, f, kopt, k, i; |
---|
| 389 | int *xopt; |
---|
| 390 | double sopt, s1, s2; |
---|
| 391 | |
---|
| 392 | xopt = x + (n + 1); |
---|
| 393 | kopt = k = g = 0; |
---|
| 394 | min = 1000000; |
---|
| 395 | for (i = n; i!=0 ; i--) |
---|
| 396 | { |
---|
| 397 | c = xopt[i]; |
---|
| 398 | if (c > 1) |
---|
| 399 | { |
---|
| 400 | if (c < min) |
---|
| 401 | min = c; |
---|
| 402 | if (c > k) |
---|
| 403 | k = c; |
---|
| 404 | } |
---|
| 405 | else |
---|
| 406 | g = 1; |
---|
| 407 | } |
---|
| 408 | k -= min; |
---|
| 409 | if ((g==0) && (k < 4)) |
---|
| 410 | return; |
---|
| 411 | if (k < min) |
---|
| 412 | min = k+1; |
---|
| 413 | sopt = (double)1.0e10; |
---|
| 414 | for (k = min; k > 1; k--) |
---|
| 415 | { |
---|
| 416 | s2 = s1 = (double)0.0; |
---|
| 417 | for(i = n; i!=0 ; i--) |
---|
| 418 | { |
---|
| 419 | c = xopt[i]; |
---|
| 420 | if (c > 1) |
---|
| 421 | { |
---|
| 422 | d = c / k; |
---|
| 423 | d *= k; |
---|
| 424 | f = d = c - d; |
---|
| 425 | if (f!=0) |
---|
| 426 | { |
---|
| 427 | f = k - f; |
---|
| 428 | if (f < d) |
---|
| 429 | s2 += (double)f / (double)c; |
---|
| 430 | else |
---|
| 431 | s1 += (double)d / (double)c; |
---|
| 432 | } |
---|
| 433 | } |
---|
| 434 | } |
---|
| 435 | s1 += s2 + sqrt(s1 * s2); |
---|
| 436 | s1 -= (double)0.01 * sqrt((double)k); |
---|
| 437 | if (s1 < sopt) |
---|
| 438 | { |
---|
| 439 | sopt = s1; |
---|
| 440 | kopt = k; |
---|
| 441 | } |
---|
| 442 | } |
---|
| 443 | for(i = n; i!=0 ; i--) |
---|
| 444 | { |
---|
| 445 | x[i] = 1; |
---|
| 446 | c = xopt[i]; |
---|
| 447 | if (c > 1) |
---|
| 448 | { |
---|
| 449 | d = c / kopt; |
---|
| 450 | d *= kopt; |
---|
| 451 | x[i] = d; |
---|
| 452 | d = c - d; |
---|
| 453 | if ((d!=0) && (kopt < 2 * d)) |
---|
| 454 | x[i] += kopt; |
---|
| 455 | } |
---|
| 456 | } |
---|
| 457 | if (g==0) |
---|
| 458 | wGcd(x, n); |
---|
| 459 | } |
---|
| 460 | |
---|
| 461 | |
---|
| 462 | void wNorm(int *degw, int *lpol, int npol, double *rel) |
---|
| 463 | { |
---|
| 464 | int i, j, ecu, ec; |
---|
| 465 | int *ex; |
---|
| 466 | double *r; |
---|
| 467 | |
---|
| 468 | ex = degw; |
---|
| 469 | r = rel; |
---|
| 470 | for (i = 0; i < npol; i++) |
---|
| 471 | { |
---|
| 472 | ecu = *ex++; |
---|
| 473 | for (j = lpol[i] - 1; j!=0 ; j--) |
---|
| 474 | { |
---|
| 475 | ec = *ex++; |
---|
| 476 | if (ec > ecu) |
---|
| 477 | ecu = ec; |
---|
| 478 | } |
---|
| 479 | *r = (double)1.0 / (double)(ecu * ecu); |
---|
| 480 | r++; |
---|
| 481 | } |
---|
| 482 | } |
---|