source: git/factory/ffreval.cc @ b1dfaf

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