polynomial.H

Go to the documentation of this file.

// polynomial arithmetic
// last change: 2004-09-13

#ifndef POLYNOMIAL_HEADER_
#define POLYNOMIAL_HEADER_

#include <gmp.h>

namespace polynomial
{

//#define DEBUG    /* for explicit DEBUGGING-sessions */
//#define VERBOSE  /* be more verbose */

// this define should be triggered by Makefile
// #define USE_DFT /* enable use of discrete fourier transform */
#if !defined(USE_DFT) && (CONTINUATION_METHOD==2)
 #define USE_DFT /* enable use of discrete fourier transform */
#endif

typedef mpz_t* TPolynom;
typedef const mpz_t* TconstPolynom;
/* A polynomial is defined by an array of mpz_t-coefficients:
   P(x)=a0*x^0+a1*x^1+...+a[k-1]*x^[k-1],
   stored as (a0,a1,a2,...,a[k-1]).
*/

/* Hints for function calls:

 TPolynom denotes a pointer to an array of coefficients
 (which thereby defines a polynomial).

 As usual in C++, take respect to the following calling conventions:

 Parameter              Restrictions/Properties
 -------------------------------------------------------------
 mpz_t* p                mutable pointer, mutable data
 const mpz_t* p          mutable pointer, constant data
 mpz_t* const p          constant pointer, mutable data
 const mpz_t* const p    constant pointer, constant data

 "const TPolynom P" and "TPolynom const P" are (unfortunately) equivalent.
 Both define a constant pointer to a mutable polynomial!

 For this reason, I have defined an auxiliary datatype
 for constant polynomials: TconstPolynom
 So keep this in mind:

 Parameter              Restrictions/Properties
 ----------------------------------------------------------------
 TPolynom P              P mutable,     coefficients mutable
 const TPolynom P        P constant,    coefficients mutable
 TPolynom const P        P constant,    coefficients mutable
 TconstPolynom P         P mutable,     coefficients constant
 const TconstPolynom P   P constant,    coefficients constant
 TconstPolynom const P   P constant,    coefficients constant

*/


class TTempPolynom : private ForbidAssignment
{
private:
  mpz_t* const data;
  const int _capacity;
public:
  explicit TTempPolynom(const int n) : data(new mpz_t[n]), _capacity(n)
   {
     //cout << "allocate polynomial of size " << n << endl;
     for (int i=0; i<n; ++i) mpz_init(data[i]);
   }
  TTempPolynom(const int n, const int estimated_operand_size) : data(new mpz_t[n]), _capacity(n)
   {
     //cout << "allocate polynomial of size " << n << endl;
     for (int i=0; i<n; ++i) mpz_init2(data[i],estimated_operand_size);
   }
  ~TTempPolynom()
   {
     //cout << "dispose polynomial of size " << _capacity << endl;
     for (int i=_capacity-1; i>=0; --i) mpz_clear(data[i]);
     delete [] data;
   }
  int capacity() const { return _capacity; }

  //const mpz_t& operator[] (const unsigned int i) const { return data[i]; }
  //mpz_t& operator[] (const unsigned int i) { return data[i]; }

  // dangerous, but helpful:
  // remember: TPolynom will not live longer than TempPolynom!
  operator TPolynom() const { return data; }
};


void print (const TconstPolynom P, int k);


void eval(mpz_t res, const TconstPolynom P, const int k, const mpz_t x, const mpz_t m);
// berechnet mit dem Hornerschema res = P(x) mod m


#ifdef USE_DFT
 // fast discrete fourier transformation
 // (will allocate internal temporary memory, which you may want to release)

void clear_dft_tempmemory() ;

#endif


int classic_mul(const TPolynom Pr, const int kr,
                const TconstPolynom P1, const int k1,
                const TconstPolynom P2, const int k2);
// Pres = P1*P2, &Pres must be different from &P1,&P2
// returns degree of the resulting polynomial
// complexity: O(k1*k2)
// (classical method, implemented alternatively to Karatsuba for small polynomials)


int classic_mul(const TPolynom Pr, const int kr,
                const TconstPolynom P1, const int k1,
                const TconstPolynom P2, const int k2, const mpz_t m);
// Pres = P1*P2, &Pres must be different from &P1,&P2
// returns degree of the resulting polynomial
// complexity: O(k1*k2)
// (classical method, implemented alternatively to Karatsuba for small polynomials)



int square(const TPolynom R, const int kR, const TconstPolynom P, const int k, const mpz_t m);
// R = P^2, &R must be different from &P
// returns degree of resulting polynomial
// complexity: O(k^1.59) -> Karatsuba


int square(const TPolynom R, const int kR, const TconstPolynom P, const int k);
// R = P^2, &R must be different from &P
// returns degree of resulting polynomial
// complexity: O(k^1.59) -> Karatsuba


int mul(const TPolynom R, const int kR,
        TconstPolynom P1, int k1,
        TconstPolynom P2, int k2, const mpz_t m);
// Pres = P1*P2, &R must be different from &P1,&P2
// returns degree of the resulting polynomial
// complexity: O(max(k1,k2)^1.59) -> Karatsuba
// resp. complexity: O((k1+k2)*ld(k1+k2)) -> using (optimal) dft


int mul(const TPolynom R, const int kR,
        TconstPolynom P1, int k1,
        TconstPolynom P2, int k2);
// Pres = P1*P2, &R must be different from &P1,&P2
// returns degree of resulting polynomial
// complexity: O(max(k1,k2)^1.59) -> Karatsuba


void reciprocal(TPolynom R, int &kR, const TconstPolynom P, const int k, const mpz_t m, const unsigned int scale=0);
// computes the reciprocal polynomial of P,
// R must provide enough memory!
// the given polynomial will be scaled using x^scale as multiplier,
// thereby shifting the coefficients by scale places


void classic_div(TPolynom Q, int &kQ, TPolynom R, int &kR,
         const TconstPolynom P1, int k1,
         const TconstPolynom P2, int k2, const mpz_t m) __attribute__ ((deprecated));


void classic_mod(TPolynom R, int &kR,
            const TconstPolynom P1, int k1,
            const TconstPolynom P2, int k2, const mpz_t m);
// Remainder of polynomial division, O(n^2)
// returns P1 mod P2 in R


void div(TPolynom Q, int &kQ,
         const TconstPolynom P1, const int k1,
         const TconstPolynom P2, const int k2, const mpz_t m);
// (fast) polynomial division by multiplying with the reciprocal polynomial


void mod(TPolynom R, int &kR,
         const TconstPolynom P1, int k1,
         const TconstPolynom P2, int k2, const mpz_t m);
// polynomial remainder, fast(er) method
// returns P1 modulo P2 in R
/* asymptotically faster than the "naive modulo", because we can make use of
   fast multiplication algorithms. */


void multipoint_eval(mpz_t* res,
                     const TconstPolynom P, const int k, 
                     const mpz_t* const array_of_arguments, const int size,
                     const mpz_t m);


int construct_polynomial_from_roots
      (TPolynom &res,
       const mpz_t* const roots_array, const int size,
       const mpz_t m);
  // task:
  // Creates a new polynomial using the zeros given in "roots_array";
  // this polynomial will be placed in "res" and its size will be returned.
  // To avoid memory leaks (due to erroneous calls), we expect that
  // "res" is initially an NULL-pointer.
  // additional remark:
  // Sure, it would be possible to return this pointer instead of using a reference
  // parameter, but what if you forget to delete it later? -- Using this technique,
  // the programmer is at least urged to provide a valid container for our data!

} // namespace polynomial

#endif /* POLYNOMIAL_HEADER_ */

---