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