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