source: git/factory/cf_algorithm.cc @ a8af6a8

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