source: git/factory/cfModResultant.cc @ 6caa2a6

spielwiese
Last change on this file since 6caa2a6 was e4fe2b, checked in by Oleksandr Motsak <motsak@…>, 13 years ago
FIX: Fixed huge BUG in cf_gmp.h CHG: starting to cleanup factory
  • Property mode set to 100644
File size: 11.0 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file cfModResultant.cc
5 *
6 * modular resultant algorithm
7 *
8 * @author Martin Lee
9 *
10 * @internal @version \$Id$
11 *
12 **/
13/*****************************************************************************/
14
15#include "config.h"
16
17#include "cf_assert.h"
18
19#include "cfModResultant.h"
20#include "cf_primes.h"
21#include "templates/ftmpl_functions.h"
22#include "cf_map.h"
23#include "cf_algorithm.h"
24#include "cf_iter.h"
25
26#ifdef HAVE_NTL
27#include "NTLconvert.h"
28#endif
29
30//TODO arrange by bound= deg (F,xlevel)*deg (G,i)+deg (G,xlevel)*deg (F, i)
31static inline
32void myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
33                  CFMap & N, const Variable& x)
34{
35  int n= tmax (F.level(), G.level());
36  int * degsf= new int [n + 1];
37  int * degsg= new int [n + 1];
38
39  for (int i = 0; i <= n; i++)
40    degsf[i]= degsg[i]= 0;
41
42  degsf= degrees (F, degsf);
43  degsg= degrees (G, degsg);
44
45  int both_non_zero= 0;
46  int f_zero= 0;
47  int g_zero= 0;
48  int both_zero= 0;
49  int degsfx, degsgx;
50
51  if (x.level() != 1)
52  {
53    int xlevel= x.level();
54
55    for (int i= 1; i <= n; i++)
56    {
57      if (degsf[i] != 0 && degsg[i] != 0)
58      {
59        both_non_zero++;
60        continue;
61      }
62      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
63      {
64        f_zero++;
65        continue;
66      }
67      if (degsg[i] == 0 && degsf[i] && i <= F.level())
68      {
69        g_zero++;
70        continue;
71      }
72    }
73
74    M.newpair (Variable (xlevel), Variable (1));
75    N.newpair (Variable (1), Variable (xlevel));
76    degsfx= degsf [xlevel];
77    degsgx= degsg [xlevel];
78    degsf [xlevel]= 0;
79    degsg [xlevel]= 0;
80    if (getNumVars (F) == 2 || getNumVars (G) == 2)
81    {
82      int pos= 2;
83      for (int i= 1; i <= n; i++)
84      {
85        if (i != xlevel)
86        {
87          if (i != pos && (degsf[i] != 0 || degsg [i] != 0))
88          {
89            M.newpair (Variable (i), Variable (pos));
90            N.newpair (Variable (pos), Variable (i));
91            pos++;
92          }
93        }
94      }
95      delete [] degsf;
96      delete [] degsg;
97      return;
98    }
99
100    if (both_non_zero == 0)
101    {
102      delete [] degsf;
103      delete [] degsg;
104      return;
105    }
106
107    // map Variables which do not occur in both polynomials to higher levels
108    int k= 1;
109    int l= 1;
110    for (int i= 1; i <= n; i++)
111    {
112      if (i == xlevel)
113        continue;
114      if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
115      {
116        if (k + both_non_zero != i)
117        {
118          M.newpair (Variable (i), Variable (k + both_non_zero));
119          N.newpair (Variable (k + both_non_zero), Variable (i));
120        }
121        k++;
122      }
123      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
124      {
125        if (l + g_zero + both_non_zero != i)
126        {
127          M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
128          N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
129        }
130        l++;
131      }
132    }
133
134    int m= tmax (F.level(), G.level());
135    int min_max_deg;
136    k= both_non_zero;
137    l= 0;
138    int i= 1;
139    while (k > 0)
140    {
141      if (degsf [i] != 0 && degsg [i] != 0)
142        min_max_deg= degsgx*degsf[i] + degsfx*degsg[i];
143      else
144        min_max_deg= 0;
145      while (min_max_deg == 0)
146      {
147        i++;
148        if (degsf [i] != 0 && degsg [i] != 0)
149          min_max_deg= degsgx*degsf[i] + degsfx*degsg[i];
150        else
151          min_max_deg= 0;
152      }
153      for (int j= i + 1; j <=  m; j++)
154      {
155        if (degsgx*degsf[j] + degsfx*degsg[j] <= min_max_deg &&
156            degsf[j] != 0 && degsg [j] != 0)
157        {
158          min_max_deg= degsgx*degsf[j] + degsfx*degsg[j];
159          l= j;
160        }
161      }
162      if (l != 0)
163      {
164        if (l != k && l != xlevel && k != 1)
165        {
166          M.newpair (Variable (l), Variable(k));
167          N.newpair (Variable (k), Variable(l));
168          degsf[l]= 0;
169          degsg[l]= 0;
170          l= 0;
171        }
172        else
173        {
174          degsf[l]= 0;
175          degsg[l]= 0;
176          l= 0;
177        }
178      }
179      else if (l == 0)
180      {
181        if (i != k && i != xlevel && k != 1)
182        {
183          M.newpair (Variable (i), Variable (k));
184          N.newpair (Variable (k), Variable (i));
185          degsf[i]= 0;
186          degsg[i]= 0;
187        }
188        else
189        {
190          degsf[i]= 0;
191          degsg[i]= 0;
192        }
193        i++;
194      }
195      k--;
196    }
197  }
198  else
199  {
200    //arrange Variables such that no gaps occur
201    for (int i= 1; i <= n; i++)
202    {
203      if (degsf[i] == 0 && degsg[i] == 0)
204      {
205        both_zero++;
206        continue;
207      }
208      else
209      {
210        if (both_zero != 0 && i != 1)
211        {
212          M.newpair (Variable (i), Variable (i - both_zero));
213          N.newpair (Variable (i - both_zero), Variable (i));
214        }
215      }
216    }
217  }
218
219  delete [] degsf;
220  delete [] degsg;
221}
222
223static inline
224CanonicalForm oneNorm (const CanonicalForm& F)
225{
226  if (F.inZ())
227    return abs (F);
228
229  CanonicalForm result= 0;
230  for (CFIterator i= F; i.hasTerms(); i++)
231    result += oneNorm (i.coeff());
232  return result;
233}
234
235// if F and G are both non constant, make sure their level is equal
236static inline
237CanonicalForm uniResultant (const CanonicalForm& F, const CanonicalForm& G)
238{
239#ifdef HAVE_NTL
240  ASSERT (getCharacteristic() > 0, "characteristic > 0 expected");
241  if (F.inCoeffDomain() && G.inCoeffDomain())
242    return 1;
243  if (F.isZero() || G.isZero())
244    return 0;
245
246  zz_pBak bak;
247  bak.save();
248  zz_p::init (getCharacteristic());
249  zz_pX NTLF= convertFacCF2NTLzzpX (F);
250  zz_pX NTLG= convertFacCF2NTLzzpX (G);
251
252  zz_p NTLResult= resultant (NTLF, NTLG);
253
254  bak.restore();
255  return CanonicalForm (to_long (rep (NTLResult)));
256#else
257  return resultant (F, G, F.mvar());
258#endif
259}
260
261static inline
262void evalPoint (const CanonicalForm& F, const CanonicalForm& G,
263                CanonicalForm& FEval, CanonicalForm& GEval, int& evalPoint)
264{
265  int degF, degG;
266  Variable x= Variable (1);
267  degF= degree (F, x);
268  degG= degree (G, x);
269  do
270  {
271    evalPoint++;
272    if (evalPoint >= getCharacteristic())
273      break;
274    FEval= F (evalPoint, 2);
275    GEval= G (evalPoint, 2);
276    if (degree (FEval, 1) < degF || degree (GEval, 1) < degG)
277      continue;
278    else
279      return;
280  }
281  while (evalPoint < getCharacteristic());
282}
283
284static inline CanonicalForm
285newtonInterp (const CanonicalForm alpha, const CanonicalForm u,
286              const CanonicalForm newtonPoly, const CanonicalForm oldInterPoly,
287              const Variable & x)
288{
289  CanonicalForm interPoly;
290
291  interPoly= oldInterPoly+((u - oldInterPoly (alpha, x))/newtonPoly (alpha, x))
292                            *newtonPoly;
293  return interPoly;
294}
295
296CanonicalForm
297resultantFp (const CanonicalForm& A, const CanonicalForm& B, const Variable& x,
298             bool prob)
299{
300  ASSERT (getCharacteristic() > 0, "characteristic > 0 expected");
301
302  if (A.isZero() || B.isZero())
303    return 0;
304
305  int degAx= degree (A, x);
306  int degBx= degree (B, x);
307  if (A.level() < x.level())
308    return power (A, degBx);
309  if (B.level() < x.level())
310    return power (B, degAx);
311
312  if (degAx == 0)
313    return power (A, degBx);
314  else if (degBx == 0)
315    return power (B, degAx);
316
317  CanonicalForm F= A;
318  CanonicalForm G= B;
319
320  CFMap M, N;
321  myCompress (F, G, M, N, x);
322
323  F= M (F);
324  G= M (G);
325
326  Variable y= Variable (2);
327
328  if (F.isUnivariate() && G.isUnivariate() && F.level() == G.level())
329    return N(uniResultant (F, G));
330
331  int i= -1;
332  CanonicalForm GEval, FEval, recResult, H;
333  CanonicalForm newtonPoly= 1;
334  CanonicalForm modResult= 0;
335
336  Variable z= Variable (1);
337  int bound= degAx*degree (G, 2) + degree (F, 2)*degBx;
338
339  int count= 0;
340  int equalCount= 0;
341  do
342  {
343    evalPoint (F, G, FEval, GEval, i);
344
345    ASSERT (i < getCharacteristic(), "ran out of points");
346
347    recResult= resultantFp (FEval, GEval, z);
348
349    H= newtonInterp (i, recResult, newtonPoly, modResult, y);
350
351    if (H == modResult)
352      equalCount++;
353
354    count++;
355    if (count > bound || (prob && equalCount == 2))
356      break;
357
358    modResult= H;
359    newtonPoly *= (y - i);
360  } while (1);
361
362  return N (H);
363}
364
365static inline
366CanonicalForm
367balanceUni ( const CanonicalForm & f, const CanonicalForm & q )
368{
369    Variable x = f.mvar();
370    CanonicalForm result = 0, qh = q / 2;
371    CanonicalForm c;
372    CFIterator i;
373    for ( i = f; i.hasTerms(); i++ ) {
374        c = mod( i.coeff(), q );
375        if ( c > qh )
376            result += power( x, i.exp() ) * (c - q);
377        else
378            result += power( x, i.exp() ) * c;
379    }
380    return result;
381}
382
383static inline
384CanonicalForm
385symmetricRemainder (const CanonicalForm& f, const CanonicalForm& q)
386{
387  CanonicalForm result= 0;
388  if (f.isUnivariate() || f.inCoeffDomain())
389    return balanceUni (f, q);
390  else
391  {
392    Variable x= f.mvar();
393    for (CFIterator i= f; i.hasTerms(); i++)
394      result += power (x, i.exp())*symmetricRemainder (i.coeff(), q);
395  }
396  return result;
397}
398
399CanonicalForm
400resultantZ (const CanonicalForm& A, const CanonicalForm& B, const Variable& x,
401            bool prob)
402{
403  ASSERT (getCharacteristic() == 0, "characteristic > 0 expected");
404
405  int degAx= degree (A, x);
406  int degBx= degree (B, x);
407  if (A.level() < x.level())
408    return power (A, degBx);
409  if (B.level() < x.level())
410    return power (B, degAx);
411
412  if (degAx == 0)
413    return power (A, degBx);
414  else if (degBx == 0)
415    return power (B, degAx);
416
417  CanonicalForm F= A;
418  CanonicalForm G= B;
419
420  Variable X= x;
421  if (F.level() != x.level() || G.level() != x.level())
422  {
423    if (F.level() > G.level())
424      X= F.mvar();
425    else
426      X= G.mvar();
427    F= swapvar (F, X, x);
428    G= swapvar (G, X, x);
429  }
430
431  // now X is the main variable
432
433  CanonicalForm d= 0;
434  CanonicalForm dd= 0;
435  CanonicalForm buf;
436  for (CFIterator i= F; i.hasTerms(); i++)
437  {
438    buf= oneNorm (i.coeff());
439    d= (buf > d) ? buf : d;
440  }
441  CanonicalForm e= 0, ee= 0;
442  for (CFIterator i= G; i.hasTerms(); i++)
443  {
444    buf= oneNorm (i.coeff());
445    e= (buf > e) ? buf : e;
446  }
447  d= power (d, degBx);
448  e= power (e, degAx);
449  CanonicalForm bound= 1;
450  for (int i= degBx + degAx; i > 1; i--)
451    bound *= i;
452  bound *= d*e;
453  bound *= 2;
454
455  bool onRational= isOn (SW_RATIONAL);
456  if (onRational)
457    Off (SW_RATIONAL);
458  int i = cf_getNumBigPrimes() - 1;
459  int p;
460  CanonicalForm l= lc (F)*lc(G);
461  CanonicalForm resultModP, q (0), newResult, newQ;
462  CanonicalForm result;
463  do
464  {
465    p = cf_getBigPrime( i );
466    i--;
467    while ( i >= 0 && mod( l, p ) == 0)
468    {
469      p = cf_getBigPrime( i );
470      i--;
471    }
472
473    ASSERT (i >= 0, "ran out of primes"); //sic
474
475    setCharacteristic (p);
476
477    resultModP= resultantFp (mapinto (F), mapinto (G), X, prob);
478
479    setCharacteristic (0);
480
481    if ( q.isZero() )
482    {
483      result= mapinto(resultModP);
484      q= p;
485    }
486    else
487    {
488      chineseRemainder( result, q, mapinto (resultModP), p, newResult, newQ );
489      q= newQ;
490      result= newResult;
491      if (newQ > bound)
492      {
493        result= symmetricRemainder (result, q);
494        break;
495      }
496    }
497  } while (1);
498
499  if (onRational)
500    On (SW_RATIONAL);
501  return swapvar (result, X, x);
502}
503
Note: See TracBrowser for help on using the repository browser.