source: git/factory/ffreval.cc @ 13be01

fieker-DuValspielwiese
Last change on this file since 13be01 was 3fd1ff, checked in by Hans Schönemann <hannes@…>, 16 years ago
*hannes: ezgcd_p git-svn-id: file:///usr/local/Singular/svn/trunk@10821 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 10.9 KB
Line 
1#include <config.h>
2
3#include "assert.h"
4#include "debug.h"
5
6#include "cf_defs.h"
7#include "canonicalform.h"
8#include "fac_util.h"
9#include "cf_algorithm.h"
10#include "cf_reval.h"
11#include "cf_random.h"
12#include "cf_primes.h"
13#include "fac_distrib.h"
14#include "ftmpl_functions.h"
15#include "ffreval.h"
16#include "cf_binom.h"
17#include "fac_iterfor.h"
18#include "cf_iter.h"
19
20void FFREvaluation::nextpoint()
21{
22  // enumerates the points stored in values
23  int n = values.max();
24  int p = getCharacteristic();
25  for( int i=values.min(); i<=n; i++ )
26  {
27    if( values[i] != p-1 )
28    {
29      values[i] += 1;
30      break;
31    }
32    values[i] += 1; // becomes 0
33  }
34}
35
36bool FFREvaluation::hasNext()
37{
38  int n = values.max();
39
40  for( int i=values.min(); i<=n; i++ )
41  {
42    if( values[i]!=start[i] )
43      return true;
44  }
45  return false;
46}
47
48FFREvaluation& FFREvaluation::operator= ( const FFREvaluation & e )
49{
50  if( this != &e )
51  {
52    if( gen != NULL )
53      delete gen;
54    values = e.values;
55    start = e.start;
56    if( e.gen == NULL )
57      gen = NULL;
58    else
59      gen = e.gen->clone();
60  }
61  return *this;
62}
63
64/* ------------------------------------------------------------------------ */
65/* forward declarations: fin_ezgcd stuff*/
66static bool findeval_P( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG, int bound, int & counter );
67static void solveF_P ( const CFArray & P, const CFArray & Q, const CFArray & S, const CFArray & T, const CanonicalForm & C, int r, CFArray & a );
68static CanonicalForm derivAndEval_P ( const CanonicalForm & f, int n, const Variable & x, const CanonicalForm & a );
69static CanonicalForm checkCombination_P ( const CanonicalForm & W, const Evaluation & a, const IteratedFor & e, int k );
70static CFArray findCorrCoeffs_P ( const CFArray & P, const CFArray & Q, const CFArray & P0, const CFArray & Q0, const CFArray & S, const CFArray & T, const CanonicalForm & C, const Evaluation & I, int r, int k, int h, int * n );
71static bool liftStep_P ( CFArray & P, int k, int r, int t, const Evaluation & A, const CFArray & lcG, const CanonicalForm & Uk, int * n, int h );
72static bool Hensel_P ( const CanonicalForm & U, CFArray & G, const CFArray & lcG, const Evaluation & A, const Variable & x );
73
74CanonicalForm fin_ezgcd( const CanonicalForm & FF, const CanonicalForm & GG )
75{
76  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG, lcD;
77  CFArray DD( 1, 2 ), lcDD( 1, 2 );
78  int degF, degG, delta, count, maxeval;
79  count = 0; // number of eval. used
80  maxeval = getCharacteristic(); // bound on the number of eval. to use
81  REvaluation b, bt;
82  bool gcdfound = false;
83  Variable x = Variable(1);
84  f = content( FF, x ); g = content( GG, x ); d = gcd( f, g );
85  F = FF / f; G = GG / g;
86  if( F.isUnivariate() && G.isUnivariate() )
87  {
88    if( F.mvar() == G.mvar() )
89      d *= gcd( F, G );
90    return d;
91  }
92  if( gcd_test_one( F, G, false ) )
93    return d;
94
95  lcF = LC( F, x ); lcG = LC( G, x );
96  lcD = gcd( lcF, lcG );
97  delta = 0;
98  degF = degree( F, x ); degG = degree( G, x );
99  b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
100  while( !gcdfound )
101  {
102    if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count ))
103    { // too many eval. used --> try another method
104      Off( SW_USE_EZGCD_P );
105      d *= gcd( F, G );
106      On( SW_USE_EZGCD_P );
107      return d;
108    }
109    delta = degree( Db );
110    if( delta == 0 )
111      return d;
112    while( true )
113    {
114      bt = b;
115      if( !findeval_P( F, G, Fbt, Gbt, Dbt, bt, delta+1, degF, degG, maxeval, count ))
116      { // too many eval. used --> try another method
117        Off( SW_USE_EZGCD_P );
118        d *= gcd( F, G );
119        On( SW_USE_EZGCD_P );
120        return d;
121      }
122      int dd = degree( Dbt );
123      if( dd == 0 )
124        return d;
125      if( dd == delta )
126        break;
127      if( dd < delta )
128      {
129        delta = dd;
130        b = bt;
131        Db = Dbt; Fb = Fbt; Gb = Gbt;
132      }
133    }
134    if( degF <= degG && delta == degF && fdivides( F, G ) )
135      return d*F;
136    if( degG < degF && delta == degG && fdivides( G, F ) )
137      return d*G;
138    if( delta != degF && delta != degG )
139    {
140      bool B_is_F;
141      CanonicalForm xxx;
142      DD[1] = Fb / Db;
143      xxx = gcd( DD[1], Db );
144      if( xxx.inCoeffDomain() )
145      {
146        B = F;
147        DD[2] = Db;
148        lcDD[1] = lcF;
149        lcDD[2] = lcF;
150        B *= lcF;
151        B_is_F = true;
152      }
153      else
154      {
155        DD[1] = Gb / Db;
156        xxx = gcd( DD[1], Db );
157        if( xxx.inCoeffDomain() )
158        {
159          B = G;
160          DD[2] = Db;
161          lcDD[1] = lcG;
162          lcDD[2] = lcG;
163          B *= lcG;
164          B_is_F = false;
165        }
166        else // special case handling
167        {
168          Off( SW_USE_EZGCD_P );
169          d *= gcd( F, G ); // try another method
170          On( SW_USE_EZGCD_P );
171          return d;
172        }
173      }
174      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
175      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
176
177      gcdfound = Hensel_P( B, DD, lcDD, b, x );
178
179      if( gcdfound )
180      {
181        CanonicalForm cand = DD[2] / content( DD[2], Variable(1) );
182        if( B_is_F )
183          gcdfound = fdivides( cand, G );
184        else
185          gcdfound = fdivides( cand, F );
186      }
187    }
188    delta++;
189  }
190  return d * DD[2] / content( DD[2],Variable(1) );
191}
192
193
194static bool findeval_P( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG, int maxeval, int & count )
195{
196  if( delta != 0 )
197  {
198    b.nextpoint();
199    if( count++ > maxeval )
200      return false; // too many eval. used
201  }
202  while( true )
203  {
204    Fb = b( F );
205    if( degree( Fb ) == degF )
206    {
207      Gb = b( G );
208      if( degree( Gb ) == degG )
209      {
210        Db = gcd( Fb, Gb );
211        if( delta > 0 )
212        {
213          if( degree( Db ) < delta )
214            return true;
215        }
216        else
217          return true;
218      }
219    }
220    b.nextpoint();
221    if( count++ > maxeval ) // too many eval. used
222      return false;
223  }
224}
225
226
227static void solveF_P ( const CFArray & P, const CFArray & Q, const CFArray & S, const CFArray & T, const CanonicalForm & C, int r, CFArray & a )
228{
229  CanonicalForm bb, b = C;
230  for ( int j = 1; j < r; j++ )
231  {
232    divrem( S[j] * b, Q[j], a[j], bb );
233    a[j] = T[j] * b + a[j] * P[j];
234    b = bb;
235  }
236  a[r] = b;
237}
238
239
240static CanonicalForm derivAndEval_P ( const CanonicalForm & f, int n, const Variable & x, const CanonicalForm & a )
241{
242  if ( n == 0 )
243    return f( a, x );
244  if ( f.degree( x ) < n )
245    return 0;
246  CFIterator i;
247  CanonicalForm sum = 0, fact;
248  int min, j;
249  Variable v = Variable( f.level() + 1 );
250  for ( i = swapvar( f, x, v); i.hasTerms() && i.exp() >= n; i++ )
251  {
252    fact = 1;
253    min = i.exp() - n;
254    for ( j = i.exp(); j > min; j-- )
255      fact *= j;
256    sum += fact * i.coeff() * power( v, min );
257  }
258  return sum( a, v );
259}
260
261
262static CFArray findCorrCoeffs_P ( const CFArray & P, const CFArray & Q, const CFArray & P0, const CFArray & Q0, const CFArray & S, const CFArray & T, const CanonicalForm & C, const Evaluation & I, int r, int k, int h, int * n )
263{
264  bool what;
265  int i, j, m;
266  CFArray A(1,r), a(1,r);
267  CanonicalForm C0, Dm, g, prodcomb;
268  C0 = I( C, 2, k-1 );
269  solveF_P( P0, Q0, S, T, 1, r, a );
270  for ( i = 1; i <= r; i++ )
271    A[i] = (a[i] * C0) % P0[i];
272  for ( m = 0; m <= h && ( m == 0 || Dm != 0 ); m++ )
273  {
274    Dm = -C;
275    prodcomb = 1;
276    for ( i = 1; i <= r; i++ )
277    {
278      Dm += prodcomb * A[i] * Q[i];
279      prodcomb *= P[i];
280    }
281    if ( Dm != 0 )
282    {
283      if ( k == 2 )
284      {
285        solveF_P( P0, Q0, S, T, Dm, r, a );
286        for ( i = 1; i <= r; i++ )
287          A[i] -= a[i];
288      }
289      else
290      {
291        IteratedFor e( 2, k-1, m+1 );
292        while ( e.iterations_left() )
293        {
294          j = 0;
295          what = true;
296          for ( i = 2; i < k; i++ )
297          {
298            j += e[i];
299            if ( e[i] > n[i] )
300            {
301              what = false;
302              break;
303            }
304          }
305          if ( what && j == m+1 )
306          {
307            g = Dm;
308            for ( i = k-1; i >= 2 && ! g.isZero(); i-- )
309            g = derivAndEval_P( g, e[i], Variable( i ), I[i] );
310            if ( ! g.isZero() )
311            {
312              prodcomb = 1;
313              for ( i = 2; i < k; i++ )
314                for ( j = 2; j <= e[i]; j++ )
315                  prodcomb *= j;
316              g /= prodcomb;
317              if( ! (g.mvar() > Variable(1)) )
318              {
319                prodcomb = 1;
320                for ( i = k-1; i > 1; i-- )
321                  prodcomb *= binomialpower( Variable(i), -I[i], e[i] );
322                solveF_P( P0, Q0, S, T, g, r, a );
323                for ( i = 1; i <= r; i++ )
324                  A[i] -= a[i] * prodcomb;
325              }
326            }
327          }
328          e++;
329        }
330      }
331    }
332  }
333  return A;
334}
335
336
337static bool liftStep_P ( CFArray & P, int k, int r, int t, const Evaluation & A, const CFArray & lcG, const CanonicalForm & Uk, int * n, int h )
338{
339  CFArray K( 1, r ), Q( 1, r ), Q0( 1, r ), P0( 1, r ), S( 1, r ), T( 1, r ), alpha( 1, r );
340  CanonicalForm Rm, C, D, xa = Variable(k) - A[k];
341  CanonicalForm xa1 = xa, xa2 = xa*xa;
342  int i, m;
343  for ( i = 1; i <= r; i++ )
344  {
345    Variable vm = Variable( t + 1 );
346    Variable v1 = Variable(1);
347    K[i] = swapvar( replaceLc( swapvar( P[i], v1, vm ), swapvar( A( lcG[i], k+1, t ), v1, vm ) ), v1, vm );
348    P[i] = A( K[i], k, t );
349  }
350  Q[r] = 1;
351  for ( i = r; i > 1; i-- )
352  {
353    Q[i-1] = Q[i] * P[i];
354    P0[i] = A( P[i], 2, k-1 );
355    Q0[i] = A( Q[i], 2, k-1 );
356    (void) extgcd( P0[i], Q0[i], S[i], T[i] );
357  }
358  P0[1] = A( P[1], 2, k-1 );
359  Q0[1] = A( Q[1], 2, k-1 );
360  (void) extgcd( P0[1], Q0[1], S[1], T[1] );
361  for ( m = 1; m <= n[k]+1; m++ )
362  {
363    Rm = prod( K ) - Uk;
364    for ( i = 2; i < k; i++ )
365      Rm.mod( binomialpower( Variable(i), -A[i], n[i]+1 ) );
366    if ( mod( Rm, xa2 ) != 0 )
367    {
368      C = derivAndEval_P( Rm, m, Variable( k ), A[k] );
369      D = 1;
370      for ( i = 2; i <= m; i++ ) D *= i;
371        C /= D;
372      alpha = findCorrCoeffs_P( P, Q, P0, Q0, S, T, C, A, r, k, h, n );
373      for ( i = 1; i <= r; i++ )
374        K[i] -= alpha[i] * xa1;
375    }
376    xa1 = xa2;
377    xa2 *= xa;
378  }
379  for ( i = 1; i <= r; i++ )
380    P[i] = K[i];
381  return prod( K ) - Uk == 0;
382}
383
384
385static bool Hensel_P ( const CanonicalForm & U, CFArray & G, const CFArray & lcG, const Evaluation & A, const Variable & x )
386{
387  int k, i, h, t = A.max();
388  bool goodeval = true;
389  CFArray Uk( A.min(), A.max() );
390  int * n = new int[t+1];
391  Uk[t] = U;
392  for ( k = t-1; k > 1; k-- )
393  {
394    Uk[k] = Uk[k+1]( A[k+1], Variable( k+1 ) );
395    n[k] = degree( Uk[k], Variable( k ) );
396  }
397  for ( k = A.min(); goodeval && (k <= t); k++ )
398  {
399    h = totaldegree( Uk[k], Variable(A.min()), Variable(k-1) );
400    for ( i = A.min(); i <= k; i++ )
401      n[i] = degree( Uk[k], Variable(i) );
402    goodeval = liftStep_P( G, k, G.max(), t, A, lcG, Uk[k], n, h );
403  }
404  delete[] n;
405  return goodeval;
406}
407
Note: See TracBrowser for help on using the repository browser.