source: git/kernel/weight.cc @ 98938c

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