source: git/factory/cf_algorithm.cc @ 1a2d66

spielwiese
Last change on this file since 1a2d66 was 362fc67, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: remove $Id$
  • Property mode set to 100644
File size: 16.1 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2
3//{{{ docu
4//
5// cf_algorithm.cc - simple mathematical algorithms.
6//
7// Hierarchy: mathematical algorithms on canonical forms
8//
9// Developers note:
10// ----------------
11// A "mathematical" algorithm is an algorithm which calculates
12// some mathematical function in contrast to a "structural"
13// algorithm which gives structural information on polynomials.
14//
15// Compare these functions to the functions in `cf_ops.cc', which
16// are structural algorithms.
17//
18//}}}
19
20#include "config.h"
21
22#include "cf_assert.h"
23
24#include "cf_factory.h"
25#include "cf_defs.h"
26#include "canonicalform.h"
27#include "cf_algorithm.h"
28#include "variable.h"
29#include "cf_iter.h"
30#include "templates/ftmpl_functions.h"
31#include "algext.h"
32
33void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
34
35//{{{ CanonicalForm psr ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
36//{{{ docu
37//
38// psr() - return pseudo remainder of `f' and `g' with respect
39//   to `x'.
40//
41// `g' must not equal zero.
42//
43// For f and g in R[x], R an arbitrary ring, g != 0, there is a
44// representation
45//
46//   LC(g)^s*f = g*q + r
47//
48// with r = 0 or deg(r) < deg(g) and s = 0 if f = 0 or
49// s = max( 0, deg(f)-deg(g)+1 ) otherwise.
50// r = psr(f, g) and q = psq(f, g) are called "pseudo remainder"
51// and "pseudo quotient", resp.  They are uniquely determined if
52// LC(g) is not a zero divisor in R.
53//
54// See H.-J. Reiffen/G. Scheja/U. Vetter - "Algebra", 2nd ed.,
55// par. 15, for a reference.
56//
57// Type info:
58// ----------
59// f, g: Current
60// x: Polynomial
61//
62// Polynomials over prime power domains are admissible if
63// lc(LC(`g',`x')) is not a zero divisor.  This is a slightly
64// stronger precondition than mathematically necessary since
65// pseudo remainder and quotient are well-defined if LC(`g',`x')
66// is not a zero divisor.
67//
68// For example, psr(y^2, (13*x+1)*y) is well-defined in
69// (Z/13^2[x])[y] since (13*x+1) is not a zero divisor.  But
70// calculating it with Factory would fail since 13 is a zero
71// divisor in Z/13^2.
72//
73// Due to this inconsistency with mathematical notion, we decided
74// not to declare type `CurrentPP' for `f' and `g'.
75//
76// Developers note:
77// ----------------
78// This is not an optimal implementation.  Better would have been
79// an implementation in `InternalPoly' avoiding the
80// exponentiation of the leading coefficient of `g'.  In contrast
81// to `psq()' and `psqr()' it definitely seems worth to implement
82// the pseudo remainder on the internal level.  Should be done
83// soon.
84//
85//}}}
86CanonicalForm
87#if 0
88psr ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
89{
90
91    ASSERT( x.level() > 0, "type error: polynomial variable expected" );
92    ASSERT( ! g.isZero(), "math error: division by zero" );
93
94    // swap variables such that x's level is larger or equal
95    // than both f's and g's levels.
96    Variable X = tmax( tmax( f.mvar(), g.mvar() ), x );
97    CanonicalForm F = swapvar( f, x, X );
98    CanonicalForm G = swapvar( g, x, X );
99
100    // now, we have to calculate the pseudo remainder of F and G
101    // w.r.t. X
102    int fDegree = degree( F, X );
103    int gDegree = degree( G, X );
104    if ( (fDegree < 0) || (fDegree < gDegree) )
105        return f;
106    else
107    {
108        CanonicalForm xresult = (power( LC( G, X ), fDegree-gDegree+1 ) * F) ;
109        CanonicalForm result = xresult -(xresult/G)*G;
110        return swapvar( result, x, X );
111    }
112}
113#else
114psr ( const CanonicalForm &rr, const CanonicalForm &vv, const Variable & x )
115{
116  CanonicalForm r=rr, v=vv, l, test, lu, lv, t, retvalue;
117  int dr, dv, d,n=0;
118
119
120  dr = degree( r, x );
121  if (dr>0)
122  {
123    dv = degree( v, x );
124    if (dv <= dr) {l=LC(v,x); v = v -l*power(x,dv);}
125    else { l = 1; }
126    d= dr-dv+1;
127    //out_cf("psr(",rr," ");
128    //out_cf("",vv," ");
129    //printf(" var=%d\n",x.level());
130    while ( ( dv <= dr  ) && ( !r.isZero()) )
131    {
132      test = power(x,dr-dv)*v*LC(r,x);
133      if ( dr == 0 ) { r= CanonicalForm(0); }
134      else { r= r - LC(r,x)*power(x,dr); }
135      r= l*r -test;
136      dr= degree(r,x);
137      n+=1;
138    }
139    r= power(l, d-n)*r;
140  }
141  return r;
142}
143#endif
144//}}}
145
146//{{{ CanonicalForm psq ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
147//{{{ docu
148//
149// psq() - return pseudo quotient of `f' and `g' with respect
150//   to `x'.
151//
152// `g' must not equal zero.
153//
154// See `psr()' for more detailed information.
155//
156// Type info:
157// ----------
158// f, g: Current
159// x: Polynomial
160//
161// Developers note:
162// ----------------
163// This is not an optimal implementation.  Better would have been
164// an implementation in `InternalPoly' avoiding the
165// exponentiation of the leading coefficient of `g'.  It seemed
166// not worth to do so.
167//
168//}}}
169CanonicalForm
170psq ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
171{
172    ASSERT( x.level() > 0, "type error: polynomial variable expected" );
173    ASSERT( ! g.isZero(), "math error: division by zero" );
174
175    // swap variables such that x's level is larger or equal
176    // than both f's and g's levels.
177    Variable X = tmax( tmax( f.mvar(), g.mvar() ), x );
178    CanonicalForm F = swapvar( f, x, X );
179    CanonicalForm G = swapvar( g, x, X );
180
181    // now, we have to calculate the pseudo remainder of F and G
182    // w.r.t. X
183    int fDegree = degree( F, X );
184    int gDegree = degree( G, X );
185    if ( fDegree < 0 || fDegree < gDegree )
186        return 0;
187    else {
188        CanonicalForm result = (power( LC( G, X ), fDegree-gDegree+1 ) * F) / G;
189        return swapvar( result, x, X );
190    }
191}
192//}}}
193
194//{{{ void psqr ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, CanonicalForm & r, const Variable & x )
195//{{{ docu
196//
197// psqr() - calculate pseudo quotient and remainder of `f' and
198//   `g' with respect to `x'.
199//
200// Returns the pseudo quotient of `f' and `g' in `q', the pseudo
201// remainder in `r'.  `g' must not equal zero.
202//
203// See `psr()' for more detailed information.
204//
205// Type info:
206// ----------
207// f, g: Current
208// q, r: Anything
209// x: Polynomial
210//
211// Developers note:
212// ----------------
213// This is not an optimal implementation.  Better would have been
214// an implementation in `InternalPoly' avoiding the
215// exponentiation of the leading coefficient of `g'.  It seemed
216// not worth to do so.
217//
218//}}}
219void
220psqr ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, CanonicalForm & r, const Variable& x )
221{
222    ASSERT( x.level() > 0, "type error: polynomial variable expected" );
223    ASSERT( ! g.isZero(), "math error: division by zero" );
224
225    // swap variables such that x's level is larger or equal
226    // than both f's and g's levels.
227    Variable X = tmax( tmax( f.mvar(), g.mvar() ), x );
228    CanonicalForm F = swapvar( f, x, X );
229    CanonicalForm G = swapvar( g, x, X );
230
231    // now, we have to calculate the pseudo remainder of F and G
232    // w.r.t. X
233    int fDegree = degree( F, X );
234    int gDegree = degree( G, X );
235    if ( fDegree < 0 || fDegree < gDegree ) {
236        q = 0; r = f;
237    } else {
238        divrem( power( LC( G, X ), fDegree-gDegree+1 ) * F, G, q, r );
239        q = swapvar( q, x, X );
240        r = swapvar( r, x, X );
241    }
242}
243//}}}
244
245//{{{ static CanonicalForm internalBCommonDen ( const CanonicalForm & f )
246//{{{ docu
247//
248// internalBCommonDen() - recursively calculate multivariate
249//   common denominator of coefficients of `f'.
250//
251// Used by: bCommonDen()
252//
253// Type info:
254// ----------
255// f: Poly( Q )
256// Switches: isOff( SW_RATIONAL )
257//
258//}}}
259static CanonicalForm
260internalBCommonDen ( const CanonicalForm & f )
261{
262    if ( f.inBaseDomain() )
263        return f.den();
264    else {
265        CanonicalForm result = 1;
266        for ( CFIterator i = f; i.hasTerms(); i++ )
267            result = blcm( result, internalBCommonDen( i.coeff() ) );
268        return result;
269    }
270}
271//}}}
272
273//{{{ CanonicalForm bCommonDen ( const CanonicalForm & f )
274//{{{ docu
275//
276// bCommonDen() - calculate multivariate common denominator of
277//   coefficients of `f'.
278//
279// The common denominator is calculated with respect to all
280// coefficients of `f' which are in a base domain.  In other
281// words, common_den( `f' ) * `f' is guaranteed to have integer
282// coefficients only.  The common denominator of zero is one.
283//
284// Returns something non-trivial iff the current domain is Q.
285//
286// Type info:
287// ----------
288// f: CurrentPP
289//
290//}}}
291CanonicalForm
292bCommonDen ( const CanonicalForm & f )
293{
294    if ( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) {
295        // otherwise `bgcd()' returns one
296        Off( SW_RATIONAL );
297        CanonicalForm result = internalBCommonDen( f );
298        On( SW_RATIONAL );
299        return result;
300    } else
301        return CanonicalForm( 1 );
302}
303//}}}
304
305//{{{ bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
306//{{{ docu
307//
308// fdivides() - check whether `f' divides `g'.
309//
310// Returns true iff `f' divides `g'.  Uses some extra heuristic
311// to avoid polynomial division.  Without the heuristic, the test
312// essentialy looks like `divremt(g, f, q, r) && r.isZero()'.
313//
314// Type info:
315// ----------
316// f, g: Current
317//
318// Elements from prime power domains (or polynomials over such
319// domains) are admissible if `f' (or lc(`f'), resp.) is not a
320// zero divisor.  This is a slightly stronger precondition than
321// mathematically necessary since divisibility is a well-defined
322// notion in arbitrary rings.  Hence, we decided not to declare
323// the weaker type `CurrentPP'.
324//
325// Developers note:
326// ----------------
327// One may consider the the test `fdivides( f.LC(), g.LC() )' in
328// the main `if'-test superfluous since `divremt()' in the
329// `if'-body repeats the test.  However, `divremt()' does not use
330// any heuristic to do so.
331//
332// It seems not reasonable to call `fdivides()' from `divremt()'
333// to check divisibility of leading coefficients.  `fdivides()' is
334// on a relatively high level compared to `divremt()'.
335//
336//}}}
337bool
338fdivides ( const CanonicalForm & f, const CanonicalForm & g )
339{
340    // trivial cases
341    if ( g.isZero() )
342        return true;
343    else if ( f.isZero() )
344        return false;
345
346    if ( (f.inCoeffDomain() || g.inCoeffDomain())
347         && ((getCharacteristic() == 0 && isOn( SW_RATIONAL ))
348             || (getCharacteristic() > 0 && CFFactory::gettype() != PrimePowerDomain)) )
349    {
350        // if we are in a field all elements not equal to zero are units
351        if ( f.inCoeffDomain() )
352            return true;
353        else
354            // g.inCoeffDomain()
355            return false;
356    }
357
358    // we may assume now that both levels either equal LEVELBASE
359    // or are greater zero
360    int fLevel = f.level();
361    int gLevel = g.level();
362    if ( (gLevel > 0) && (fLevel == gLevel) )
363        // f and g are polynomials in the same main variable
364        if ( degree( f ) <= degree( g )
365             && fdivides( f.tailcoeff(), g.tailcoeff() )
366             && fdivides( f.LC(), g.LC() ) )
367        {
368            CanonicalForm q, r;
369            return divremt( g, f, q, r ) && r.isZero();
370        }
371        else
372            return false;
373    else if ( gLevel < fLevel )
374        // g is a coefficient w.r.t. f
375        return false;
376    else
377    {
378        // either f is a coefficient w.r.t. polynomial g or both
379        // f and g are from a base domain (should be Z or Z/p^n,
380        // then)
381        CanonicalForm q, r;
382        return divremt( g, f, q, r ) && r.isZero();
383    }
384}
385//}}}
386
387/// same as fdivides if true returns quotient quot of g by f otherwise quot == 0
388bool
389fdivides ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm& quot )
390{
391    quot= 0;
392    // trivial cases
393    if ( g.isZero() )
394        return true;
395    else if ( f.isZero() )
396        return false;
397
398    if ( (f.inCoeffDomain() || g.inCoeffDomain())
399         && ((getCharacteristic() == 0 && isOn( SW_RATIONAL ))
400             || (getCharacteristic() > 0 && CFFactory::gettype() != PrimePowerDomain)) )
401    {
402        // if we are in a field all elements not equal to zero are units
403        if ( f.inCoeffDomain() )
404        {
405            quot= g/f;
406            return true;
407        }
408        else
409            // g.inCoeffDomain()
410            return false;
411    }
412
413    // we may assume now that both levels either equal LEVELBASE
414    // or are greater zero
415    int fLevel = f.level();
416    int gLevel = g.level();
417    if ( (gLevel > 0) && (fLevel == gLevel) )
418        // f and g are polynomials in the same main variable
419        if ( degree( f ) <= degree( g )
420             && fdivides( f.tailcoeff(), g.tailcoeff() )
421             && fdivides( f.LC(), g.LC() ) )
422        {
423            CanonicalForm q, r;
424            if (divremt( g, f, q, r ) && r.isZero())
425            {
426              quot= q;
427              return true;
428            }
429            else
430              return false;
431        }
432        else
433            return false;
434    else if ( gLevel < fLevel )
435        // g is a coefficient w.r.t. f
436        return false;
437    else
438    {
439        // either f is a coefficient w.r.t. polynomial g or both
440        // f and g are from a base domain (should be Z or Z/p^n,
441        // then)
442        CanonicalForm q, r;
443        if (divremt( g, f, q, r ) && r.isZero())
444        {
445          quot= q;
446          return true;
447        }
448        else
449          return false;
450    }
451}
452
453/// same as fdivides but handles zero divisors in Z_p[t]/(f)[x1,...,xn] for reducible f
454bool
455tryFdivides ( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm& M, bool& fail )
456{
457    fail= false;
458    // trivial cases
459    if ( g.isZero() )
460        return true;
461    else if ( f.isZero() )
462        return false;
463
464    if (f.inCoeffDomain() || g.inCoeffDomain())
465    {
466        // if we are in a field all elements not equal to zero are units
467        if ( f.inCoeffDomain() )
468        {
469            CanonicalForm inv;
470            tryInvert (f, M, inv, fail);
471            return !fail;
472        }
473        else
474        {
475            return false;
476        }
477    }
478
479    // we may assume now that both levels either equal LEVELBASE
480    // or are greater zero
481    int fLevel = f.level();
482    int gLevel = g.level();
483    if ( (gLevel > 0) && (fLevel == gLevel) )
484    {
485        if (degree( f ) > degree( g ))
486          return false;
487        bool dividestail= tryFdivides (f.tailcoeff(), g.tailcoeff(), M, fail);
488
489        if (fail || !dividestail)
490          return false;
491        bool dividesLC= tryFdivides (f.LC(),g.LC(), M, fail);
492        if (fail || !dividesLC)
493          return false;
494        CanonicalForm q,r;
495        bool divides= tryDivremt (g, f, q, r, M, fail);
496        if (fail || !divides)
497          return false;
498        return r.isZero();
499    }
500    else if ( gLevel < fLevel )
501    {
502        // g is a coefficient w.r.t. f
503        return false;
504    }
505    else
506    {
507        // either f is a coefficient w.r.t. polynomial g or both
508        // f and g are from a base domain (should be Z or Z/p^n,
509        // then)
510        CanonicalForm q, r;
511        bool divides= tryDivremt (g, f, q, r, M, fail);
512        if (fail || !divides)
513          return false;
514        return r.isZero();
515    }
516}
517
518//{{{ CanonicalForm maxNorm ( const CanonicalForm & f )
519//{{{ docu
520//
521// maxNorm() - return maximum norm of `f'.
522//
523// That is, the base coefficient of `f' with the largest absolute
524// value.
525//
526// Valid for arbitrary polynomials over arbitrary domains, but
527// most useful for multivariate polynomials over Z.
528//
529// Type info:
530// ----------
531// f: CurrentPP
532//
533//}}}
534CanonicalForm
535maxNorm ( const CanonicalForm & f )
536{
537    if ( f.inBaseDomain() )
538        return abs( f );
539    else {
540        CanonicalForm result = 0;
541        for ( CFIterator i = f; i.hasTerms(); i++ ) {
542            CanonicalForm coeffMaxNorm = maxNorm( i.coeff() );
543            if ( coeffMaxNorm > result )
544                result = coeffMaxNorm;
545        }
546        return result;
547    }
548}
549//}}}
550
551//{{{ CanonicalForm euclideanNorm ( const CanonicalForm & f )
552//{{{ docu
553//
554// euclideanNorm() - return Euclidean norm of `f'.
555//
556// Returns the largest integer smaller or equal norm(`f') =
557// sqrt(sum( `f'[i]^2 )).
558//
559// Type info:
560// ----------
561// f: UVPoly( Z )
562//
563//}}}
564CanonicalForm
565euclideanNorm ( const CanonicalForm & f )
566{
567    ASSERT( (f.inBaseDomain() || f.isUnivariate()) && f.LC().inZ(),
568            "type error: univariate poly over Z expected" );
569
570    CanonicalForm result = 0;
571    for ( CFIterator i = f; i.hasTerms(); i++ ) {
572        CanonicalForm coeff = i.coeff();
573        result += coeff*coeff;
574    }
575    return sqrt( result );
576}
577//}}}
Note: See TracBrowser for help on using the repository browser.