source: git/Singular/weight.cc @ 666c90

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