File: poly.h

package info (click to toggle)
form 5.0.0-1
  • links: PTS, VCS
  • area: main
  • in suites: experimental
  • size: 8,312 kB
  • sloc: ansic: 110,546; cpp: 20,395; sh: 5,874; makefile: 545
file content (263 lines) | stat: -rw-r--r-- 7,802 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
#pragma once
/* #[ License : */
/*
 *   Copyright (C) 1984-2026 J.A.M. Vermaseren
 *   When using this file you are requested to refer to the publication
 *   J.A.M.Vermaseren "New features of FORM" math-ph/0010025
 *   This is considered a matter of courtesy as the development was paid
 *   for by FOM the Dutch physics granting agency and we would like to
 *   be able to track its scientific use to convince FOM of its value
 *   for the community.
 *
 *   This file is part of FORM.
 *
 *   FORM is free software: you can redistribute it and/or modify it under the
 *   terms of the GNU General Public License as published by the Free Software
 *   Foundation, either version 3 of the License, or (at your option) any later
 *   version.
 *
 *   FORM is distributed in the hope that it will be useful, but WITHOUT ANY
 *   WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 *   FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
 *   details.
 *
 *   You should have received a copy of the GNU General Public License along
 *   with FORM.  If not, see <http://www.gnu.org/licenses/>.
 */
/* #] License : */ 

extern "C" {
#include "form3.h"
}

#include <iostream>
#include <string>
#include <vector>

// for debugging
// #define POLY_MOVE_DEBUG

// macros for tform
#ifndef WITHPTHREADS
#define POLY_GETIDENTITY(X)
#define POLY_STOREIDENTITY
#else
#define POLY_GETIDENTITY(X) ALLPRIVATES *B = (X).Bpointer
#define POLY_STOREIDENTITY Bpointer = B
#endif

// maximum size of the hash table used for multiplication and division

const int POLY_MAX_HASH_SIZE = MiN(1<<20, MAXPOSITIVE);

class poly {

public:

	// variables
#ifdef WITHPTHREADS
	ALLPRIVATES *Bpointer;
#endif

	WORD *terms;
	LONG size_of_terms;
	WORD modp,modn;

	// constructors/destructor
	poly (PHEAD int, WORD=-1, WORD=1);
	poly (PHEAD const UWORD *, WORD, WORD=-1, WORD=1);
	poly (const poly &, WORD=-1, WORD=1);
	~poly ();

	// operators
	poly& operator+= (const poly&);
	poly& operator-= (const poly&);
	poly& operator*= (const poly&);
	poly& operator/= (const poly&);
	poly& operator%= (const poly&);

	const poly operator+ (const poly&) const;
	const poly operator- (const poly&) const;
	const poly operator* (const poly&) const;
	const poly operator/ (const poly&) const;
	const poly operator% (const poly&) const;

	bool operator== (const poly&) const;
	bool operator!= (const poly&) const;
	poly& operator= (const poly &);
	WORD& operator[] (int);
	const WORD& operator[] (int) const;

#if __cplusplus >= 201103L
	// support move semantics
	poly (poly &&, WORD=-1, WORD=1) noexcept;
	poly& operator= (poly &&) noexcept;
#endif

	// memory management
	void termscopy (const WORD *, int, int);
	void check_memory(int);
	void expand_memory(int);

	// check type of polynomial
	bool is_zero () const;
	bool is_one () const;
	bool is_integer () const;
	bool is_monomial () const;
	int is_dense_univariate () const;

	// properties
	int sign () const;
	int degree (int) const;
	int total_degree () const;
	int first_variable () const;
	int number_of_terms () const;
	const std::vector<int> all_variables () const;	
	const poly integer_lcoeff () const;
	const poly lcoeff_univar (int) const;
	const poly lcoeff_multivar (int) const;
	const poly coefficient (int, int) const;
	const poly derivative (int) const;
	
	// modulo calculus
	void setmod(WORD, WORD=1);
	void coefficients_modulo (UWORD *, WORD, bool);

	// simple polynomials
	static const poly simple_poly (PHEAD int, int=0, int=1, int=0, int=1);
	static const poly simple_poly (PHEAD int, const poly&, int=1, int=0, int=1);

	// conversion from/to form notation
	static void get_variables (PHEAD const std::vector<WORD *> &, bool, bool);
	static const poly argument_to_poly (PHEAD WORD *, bool, bool, poly *den=NULL);
	static void poly_to_argument (const poly &, WORD *, bool);
	static void poly_to_argument_with_den (const poly &, WORD, const UWORD *, WORD *, bool);
	int size_of_form_notation () const;
	int size_of_form_notation_with_den (WORD) const;
	const poly & normalize ();

	// operations for coefficient lists
	static const poly from_coefficient_list (PHEAD const std::vector<WORD> &, int, WORD);
	static const std::vector<WORD> to_coefficient_list (const poly &);
	static const std::vector<WORD> coefficient_list_divmod (const std::vector<WORD> &, const std::vector<WORD> &, WORD, int);

	// string output for debugging
	const std::string to_string() const;

	// monomials
	static int monomial_compare (PHEAD const WORD *, const WORD *);
	WORD last_monomial_index () const;
	WORD* last_monomial () const;
	int compare_degree_vector(const poly &) const;
	std::vector<int> degree_vector() const;
	int compare_degree_vector(const std::vector<int> &) const;

	// (internal) mathematical operations
	static void add (const poly &, const poly &, poly &);
	static void sub (const poly &, const poly &, poly &);
	static void mul (const poly &, const poly &, poly &);
	static void div (const poly &, const poly &, poly &);
	static void mod (const poly &, const poly &, poly &);
	static void divmod (const poly &, const poly &, poly &, poly &, bool only_divides);
	static bool divides (const poly &, const poly &);
	
	static void mul_one_term (const poly &, const poly &, poly &);
	static void mul_univar (const poly &, const poly &, poly &, int);
	static void mul_heap (const poly &, const poly &, poly &);

	static void divmod_one_term (const poly &, const poly &, poly &, poly &, bool);
	static void divmod_univar (const poly &, const poly &, poly &, poly &, int, bool);
	static void divmod_heap (const poly &, const poly &, poly &, poly &, bool, bool, bool&);
	
	static void push_heap (PHEAD WORD **, int);
	static void pop_heap (PHEAD WORD **, int);
};

// comparison class for monomials (for std::sort)
class monomial_larger {

public:
	
#ifndef WITHPTHREADS
	monomial_larger() {}
#else
	ALLPRIVATES *B;
	monomial_larger(ALLPRIVATES *b): B(b) {}
#endif
	
	bool operator()(const WORD *a, const WORD *b) {
		return poly::monomial_compare(BHEAD a, b) > 0;
	}
};

// stream operator
std::ostream& operator<< (std::ostream &, const poly &);

// inline function definitions

#if __cplusplus >= 201103L

inline poly::poly (poly &&a, WORD modp, WORD modn) noexcept {
#ifdef POLY_MOVE_DEBUG
	std::cout << "poly move ctor" << std::endl;
#endif
	POLY_GETIDENTITY(a);
	POLY_STOREIDENTITY;
	terms = a.terms;
	size_of_terms = a.size_of_terms;
	this->modp = a.modp;
	this->modn = a.modn;
	if (modp != -1) {
		setmod(modp,modn);
	}
	a.terms = nullptr;
	a.size_of_terms = -1;
}

inline poly& poly::operator= (poly &&a) noexcept {
#ifdef POLY_MOVE_DEBUG
	std::cout << "poly move assign" << std::endl;
#endif
	if (this != &a) {
		WORD *old_terms = terms;
		LONG old_size_of_terms = size_of_terms;
		terms = a.terms;
		size_of_terms = a.size_of_terms;
		modp = a.modp;
		modn = a.modn;
		a.terms = old_terms;
		a.size_of_terms = old_size_of_terms;
	}
	return *this;
}

#endif

/*   Checks whether the terms array is large enough to add another
 *   term (of size AM.MaxTal) to the polynomials. In case not, it is
 *   expanded.
 */
inline void poly::check_memory (int i) {
	POLY_GETIDENTITY(*this);
	if (i + 3 + AN.poly_num_vars + AM.MaxTal >= size_of_terms) expand_memory(i + AM.MaxTal);
//  Used to be i+2 but there should also be space for a trailing zero
}

// indexing operators
inline WORD& poly::operator[] (int i) {
	return terms[i];
}

inline const WORD& poly::operator[] (int i) const {
	return terms[i];
}

/*   Copies "num" WORD-sized terms from the pointer "source" to the
 *   current polynomial at index "dest"
 */
inline void poly::termscopy (const WORD *source, int dest, int num) {
#ifdef POLY_MOVE_DEBUG
	std::cout << "poly termscopy: num=" << num << std::endl;
#endif
	memcpy (terms+dest, source, num*sizeof(WORD));
}