source: git/libpolys/polys/pDebug.cc @ 9d68fd

spielwiese
Last change on this file since 9d68fd was 9d68fd, checked in by Hans Schoenemann <hannes@…>, 11 years ago
fix: better tests and normalization for subst
  • Property mode set to 100644
File size: 10.8 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/***************************************************************
5 *  File:    pDebug.h
6 *  Purpose: implementation of debug related poly routines
7 *  Author:  obachman (Olaf Bachmann)
8 *  Created: 8/00
9 *******************************************************************/
10
11#ifndef PDEBUG_CC
12#define PDEBUG_CC
13
14#include <stdarg.h>
15#include <stdio.h>
16
17
18#include "config.h"
19#include <misc/auxiliary.h>
20
21
22#ifdef PDEBUG
23
24// do the following to always enforce checking of pSetm
25// #undef PDEBUG
26// #define PDEBUG 2
27
28#include <omalloc/omalloc.h>
29
30#include <polys/monomials/ring.h>
31#include <polys/monomials/p_polys.h>
32
33#include <coeffs/coeffs.h>
34
35/***************************************************************
36 *
37 * Error reporting
38 *
39 ***************************************************************/
40// avoid recursive calls
41static BOOLEAN d_poly_error_reporting = FALSE;
42BOOLEAN dPolyReportError(poly p, ring r, const char* fmt, ...)
43{
44  if (d_poly_error_reporting) return FALSE;
45  d_poly_error_reporting = TRUE;
46  va_list ap;
47  va_start(ap, fmt);
48
49  fprintf(stderr, "\n// ***dPolyError: ");
50  vfprintf(stderr, fmt, ap);
51  fprintf(stderr, "\n occured at\n");
52  omPrintCurrentBackTraceMax(stderr, 8);
53  if (p != NULL)
54  {
55    fprintf(stderr, " occured for poly: ");
56    p_wrp(p, r);
57    omPrintAddrInfo(stderr, p, " ");
58  }
59  dErrorBreak();
60  d_poly_error_reporting = FALSE;
61  return FALSE;
62}
63
64/***************************************************************
65 *
66 * checking for ring stuff
67 *
68 ***************************************************************/
69BOOLEAN p_LmCheckIsFromRing(poly p, ring r)
70{
71  if (p != NULL)
72  {
73    #if (OM_TRACK > 0) && defined(OM_TRACK_CUSTOM)
74    void* custom = omGetCustomOfAddr(p);
75    if (custom != NULL)
76    {
77      pPolyAssumeReturnMsg(custom == r ||
78                           // be more sloppy for qrings
79                           (r->qideal != NULL &&
80                            omIsBinPageAddr(p) &&
81                            omSizeWOfAddr(p)==omSizeWOfBin(r->PolyBin)) ||
82                           rSamePolyRep((ring) custom, r),
83                           "monomial not from specified ring",p,r);
84      return TRUE;
85    }
86    else
87    #endif
88    #ifndef X_OMALLOC
89    {
90      _pPolyAssumeReturn(omIsBinPageAddr(p) && omSizeWOfAddr(p)==omSizeWOfBin(r->PolyBin),p,r);
91      return TRUE;
92    }
93    return FALSE;
94    #endif
95  }
96  return TRUE;
97}
98
99BOOLEAN p_CheckIsFromRing(poly p, ring r)
100{
101  while (p!=NULL)
102  {
103    pFalseReturn(p_LmCheckIsFromRing(p, r));
104    pIter(p);
105  }
106  return TRUE;
107}
108
109BOOLEAN p_CheckPolyRing(poly p, ring r)
110{
111  #ifndef X_OMALLOC
112  pAssumeReturn(r != NULL && r->PolyBin != NULL);
113  #endif
114  return p_CheckIsFromRing(p, r);
115}
116
117BOOLEAN p_LmCheckPolyRing(poly p, ring r)
118{
119  #ifndef X_OMALLOC
120  pAssumeReturn(r != NULL && r->PolyBin != NULL);
121  #endif
122  pAssumeReturn(p != NULL);
123  return p_LmCheckIsFromRing(p, r);
124}
125BOOLEAN p_CheckRing(ring r)
126{
127  #ifndef X_OMALLOC
128  pAssumeReturn(r != NULL && r->PolyBin != NULL);
129  #endif
130  return TRUE;
131}
132
133/***************************************************************
134 *
135 * Debugging/statistics of pDivisibleBy
136 *
137 ***************************************************************/
138BOOLEAN p_DebugLmDivisibleByNoComp(poly a, poly b, ring r)
139{
140  int i=r->N;
141
142  do
143  {
144    if (p_GetExp(a,i,r) > p_GetExp(b,i,r))
145      return FALSE;
146    i--;
147  }
148  while (i);
149#ifdef HAVE_RINGS
150  return n_DivBy(pGetCoeff(b), pGetCoeff(a), r->cf);
151#else
152  return TRUE;
153#endif
154     }
155
156
157/***************************************************************
158 *
159 * Misc things helpful for debugging
160 *
161 ***************************************************************/
162BOOLEAN pIsMonomOf(poly p, poly m)
163{
164  if (m == NULL) return TRUE;
165  while (p != NULL)
166  {
167    if (p == m) return TRUE;
168    pIter(p);
169  }
170  return FALSE;
171}
172BOOLEAN pHaveCommonMonoms(poly p, poly q)
173{
174  while (p != NULL)
175  {
176    if (pIsMonomOf(q, p))
177    {
178      return TRUE;
179    }
180    pIter(p);
181  }
182  return FALSE;
183}
184
185/***************************************************************
186 *
187 * Testing of polys
188 *
189 ***************************************************************/
190extern void p_Setm_General(poly p, ring r);
191
192static poly p_DebugInit(poly p, ring src_ring, ring dest_ring)
193{
194  poly d_p = p_Init(dest_ring);
195  int i;
196  assume(dest_ring->N == src_ring->N);
197
198  for (i=1; i<= src_ring->N; i++)
199  {
200    p_SetExp(d_p, i, p_GetExp(p, i, src_ring), dest_ring);
201  }
202  if (rRing_has_Comp(dest_ring))
203    p_SetComp(d_p, p_GetComp(p, src_ring), dest_ring);
204
205  p_Setm_General(d_p, dest_ring);
206  return d_p;
207}
208
209BOOLEAN _p_Test(poly p, ring r, int level)
210{
211  assume(r->cf !=NULL);
212
213  if (PDEBUG > level) level = PDEBUG;
214  if (level < 0 || p == NULL) return TRUE;
215
216  poly p_prev = NULL;
217
218  #ifndef OM_NDEBUG
219  #ifndef X_OMALLOC
220  // check addr with level+1 so as to check bin/page of addr
221  _pPolyAssumeReturnMsg(omTestBinAddrSize(p, (omSizeWOfBin(r->PolyBin))*SIZEOF_LONG, level+1)
222                        == omError_NoError, "memory error",p,r);
223  #endif
224  #endif
225
226  pFalseReturn(p_CheckRing(r));
227
228  // this checks that p does not contain a loop: rather expensive O(length^2)
229  #ifndef OM_NDEBUG
230  if (level > 1)
231    pFalseReturn(omTestList(p, level) == omError_NoError);
232  #endif
233
234  int ismod = p_GetComp(p, r) != 0;
235
236  while (p != NULL)
237  {
238    // ring check
239    pFalseReturn(p_LmCheckIsFromRing(p, r));
240    #ifndef OM_NDEBUG
241    #ifndef X_OMALLOC
242    // omAddr check
243    _pPolyAssumeReturnMsg(omTestBinAddrSize(p, (omSizeWOfBin(r->PolyBin))*SIZEOF_LONG, 1)
244                     == omError_NoError, "memory error",p,r);
245    #endif
246    #endif
247    // number/coef check
248    _pPolyAssumeReturnMsg(p->coef != NULL || (n_GetChar(r->cf) >= 2), "NULL coef",p,r);
249    #ifdef LDEBUG
250    _pPolyAssumeReturnMsg(r->cf->cfDBTest(p->coef,__FILE__,__LINE__,r->cf),"coeff err",p,r);
251    #endif
252    _pPolyAssumeReturnMsg(!n_IsZero(p->coef, r->cf), "Zero coef",p,r);
253
254    // check for valid comp
255    _pPolyAssumeReturnMsg(p_GetComp(p, r) >= 0 && (p_GetComp(p, r)<65000), "component out of range ?",p,r);
256    // check for mix poly/vec representation
257    _pPolyAssumeReturnMsg(ismod == (p_GetComp(p, r) != 0), "mixed poly/vector",p,r);
258
259    // special check for ringorder_s/S
260    if ((r->typ!=NULL) && (r->typ[0].ord_typ == ro_syzcomp))
261    {
262      long c1, cc1, ccc1, ec1;
263      sro_ord* o = &(r->typ[0]);
264
265      c1 = p_GetComp(p, r);
266      if (o->data.syzcomp.Components!=NULL)
267      {
268        cc1 = o->data.syzcomp.Components[c1];
269        ccc1 = o->data.syzcomp.ShiftedComponents[cc1];
270      }
271      else { cc1=0; ccc1=0; }
272      _pPolyAssumeReturnMsg(c1 == 0 || cc1 != 0, "Component <-> TrueComponent zero mismatch",p,r);
273      _pPolyAssumeReturnMsg(c1 == 0 || ccc1 != 0,"Component <-> ShiftedComponent zero mismatch",p,r);
274      ec1 = p->exp[o->data.syzcomp.place];
275      //pPolyAssumeReturnMsg(ec1 == ccc1, "Shifted comp out of sync. should %d, is %d");
276      if (ec1 != ccc1)
277      {
278        dPolyReportError(p,r,"Shifted comp out of sync. should %d, is %d",ccc1,ec1);
279        return FALSE;
280      }
281    }
282
283    // check that p_Setm works ok
284    if (level > 0)
285    {
286      poly p_should_equal = p_DebugInit(p, r, r);
287      _pPolyAssumeReturnMsg(p_ExpVectorEqual(p, p_should_equal, r), "p_Setm field(s) out of sync",p,r);
288      p_LmFree(p_should_equal, r);
289    }
290
291    // check order
292    if (p_prev != NULL)
293    {
294      int cmp = p_LmCmp(p_prev, p, r);
295      if (cmp == 0)
296      {
297        _pPolyAssumeReturnMsg(0, "monoms p and p->next are equal", p_prev, r);
298      }
299      else
300        _pPolyAssumeReturnMsg(p_LmCmp(p_prev, p, r) == 1, "wrong order", p_prev, r);
301
302      // check that compare worked sensibly
303      if (level > 1 && p_GetComp(p_prev, r) == p_GetComp(p, r))
304      {
305        int i;
306        for (i=r->N; i>0; i--)
307        {
308          if (p_GetExp(p_prev, i, r) != p_GetExp(p, i, r)) break;
309        }
310        _pPolyAssumeReturnMsg(i > 0, "Exponents equal but compare different", p_prev, r);
311      }
312    }
313    p_prev = p;
314    pIter(p);
315  }
316  return TRUE;
317}
318
319BOOLEAN _p_LmTest(poly p, ring r, int level)
320{
321  if (level < 0 || p == NULL) return TRUE;
322  poly pnext = pNext(p);
323  pNext(p) = NULL;
324  BOOLEAN test_res = _p_Test(p, r, level);
325  pNext(p) = pnext;
326  return test_res;
327}
328
329BOOLEAN _pp_Test(poly p, ring lmRing, ring tailRing, int level)
330{
331  if (PDEBUG > level) level = PDEBUG;
332  if (level < 0 || p == NULL) return TRUE;
333  if (pNext(p) == NULL || lmRing == tailRing) return _p_Test(p, lmRing, level);
334
335  pFalseReturn(_p_LmTest(p, lmRing, level));
336  pFalseReturn(_p_Test(pNext(p), tailRing, level));
337
338  // check that lm > Lm(tail)
339  if (level > 1)
340  {
341    poly lm = p;
342    poly tail = p_DebugInit(pNext(p), tailRing, lmRing);
343    poly pnext = pNext(lm);
344    pNext(lm) = tail;
345    BOOLEAN cmp = p_LmCmp(lm, tail, lmRing);
346    if (cmp != 1)
347      dPolyReportError(lm, lmRing, "wrong order: lm <= Lm(tail)");
348    p_LmFree(tail, lmRing);
349    pNext(lm) = pnext;
350    return (cmp == 1);
351  }
352  return TRUE;
353}
354
355#endif // PDEBUG
356
357#if defined(PDEBUG) || defined(PDIV_DEBUG)
358static unsigned long pDivisibleBy_number = 1;
359static unsigned long pDivisibleBy_FALSE = 1;
360static unsigned long pDivisibleBy_ShortFalse = 1;
361
362BOOLEAN pDebugLmShortDivisibleBy(poly p1, unsigned long sev_1, ring r_1,
363                               poly p2, unsigned long not_sev_2, ring r_2)
364{
365  _pPolyAssume(p_GetShortExpVector(p1, r_1) == sev_1, p1, r_1);
366  _pPolyAssume(p_GetShortExpVector(p2, r_2) == ~ not_sev_2, p2, r_2);
367
368  pDivisibleBy_number++;
369  BOOLEAN ret;
370  if (r_1 == r_2)
371    ret = p_LmDivisibleBy(p1, p2, r_1);
372  else
373    ret = p_LmDivisibleBy(p1, r_1, p2, r_2);
374
375  if (! ret) pDivisibleBy_FALSE++;
376  if (sev_1 & not_sev_2)
377  {
378    pDivisibleBy_ShortFalse++;
379    if (ret)
380      dReportError("p1 divides p2, but sev's are wrong");
381  }
382  return ret;
383}
384
385BOOLEAN pDebugLmShortDivisibleByNoComp(poly p1, unsigned long sev_1, ring r_1,
386                                       poly p2, unsigned long not_sev_2, ring r_2)
387{
388//  _pPolyAssume((p_GetComp(p1, r_1) == p_GetComp(p2, r_2)) || (p_GetComp(p1, r_1) == 0));
389  _pPolyAssume(p_GetShortExpVector(p1, r_1) == sev_1, p1, r_1);
390  _pPolyAssume(p_GetShortExpVector(p2, r_2) == ~ not_sev_2, p2, r_2);
391
392  pDivisibleBy_number++;
393  BOOLEAN ret;
394  if (r_1 == r_2)
395    ret = p_LmDivisibleByNoComp(p1, p2, r_1);
396  else
397    ret = p_LmDivisibleByNoComp(p1, r_1, p2, r_2);
398
399  if (! ret) pDivisibleBy_FALSE++;
400  if (sev_1 & not_sev_2)
401  {
402    pDivisibleBy_ShortFalse++;
403    if (ret)
404      dReportError("p1 divides p2, but sev's are wrong");
405  }
406  return ret;
407}
408
409void pPrintDivisbleByStat()
410{
411  Print("#Tests: %ld; #FALSE %ld(%ld); #SHORT %ld(%ld)\n",
412        pDivisibleBy_number,
413        pDivisibleBy_FALSE, (unsigned long) ((double)pDivisibleBy_FALSE*((double) 100)/(double)pDivisibleBy_number),
414        pDivisibleBy_ShortFalse, (unsigned long) ((double)pDivisibleBy_ShortFalse*((double)100)/(double)pDivisibleBy_FALSE));
415}
416
417#endif // PDEBUG
418
419#endif // PDEBUG_CC
Note: See TracBrowser for help on using the repository browser.