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