1 //==============================================================================================
2 //
3 // This file is part of LiDIA --- a library for computational number theory
4 //
5 // Copyright (c) 1994--2001 the LiDIA Group. All rights reserved.
6 //
7 // See http://www.informatik.tu-darmstadt.de/TI/LiDIA/
8 //
9 //----------------------------------------------------------------------------------------------
10 //
11 // $Id$
12 //
13 // Author : Stefan Neis (SN)
14 // Changes : See CVS log
15 //
16 //==============================================================================================
17
18
19 #ifndef LIDIA_POLY_INTERN_CC_GUARD_
20 #define LIDIA_POLY_INTERN_CC_GUARD_
21
22
23
24 #ifndef LIDIA_POLY_INTERN_H_GUARD_
25 # include "LiDIA/base/poly_intern.h"
26 #endif
27 #include <cctype>
28
29
30
31 #ifdef LIDIA_NAMESPACE
32 # ifndef IN_NAMESPACE_LIDIA
33 namespace LiDIA {
34 # endif
35 #endif
36
37
38
39 #define DV_BP LDBL_UNIPOL
40 #define DM_BP "base_polynomial"
41 #define DV_P LDBL_UNIPOL+12
42 #define DM_P "polynomial< T >"
43 #define LP_ERROR poly_error_msg
44
45 // BASE_POLYNOMIAL FUNCTIONS
46
47 //
48 // constructors and destructor
49 //
50
51 template< class T >
base_polynomial()52 base_polynomial< T >::base_polynomial()
53 {
54 debug_handler_l(DM_BP, "in constructor "
55 "base_polynomial()", DV_BP);
56 this->deg = -1;
57 this->coeff = NULL;
58 }
59
60
61
62 template< class T >
base_polynomial(const T & x)63 base_polynomial< T >::base_polynomial(const T & x)
64 {
65 debug_handler_l(DM_BP, "in constructor "
66 "base_polynomial(const T &)", DV_BP);
67 if (x == 0) {
68 this->deg = -1;
69 this->coeff = NULL;
70 }
71 else {
72 this->deg = 0;
73 this->coeff = new T[1];
74 memory_handler(this->coeff, DM_BP, "base_polynomial(const T &) :: "
75 "Error in memory allocation (coeff)");
76 this->coeff[0] = x;
77 }
78 }
79
80
81
82 template< class T >
base_polynomial(const T * v,lidia_size_t d)83 base_polynomial< T >::base_polynomial(const T * v, lidia_size_t d)
84 {
85 debug_handler_l(DM_BP, "in constructor "
86 "base_polynomial(const T *, lidia_size_t)", DV_BP);
87 if (d < 0 || v == NULL)
88 lidia_error_handler_para(d, "d", "d >= 0",
89 PRT, "v", "v != NULL",
90 "base_polynomial< T >::"
91 "base_polynomial(const T * v, lidia_size_t d)",
92 DM_BP, LP_ERROR[0]);
93 this->deg = d;
94 this->coeff = new T[d + 1];
95 memory_handler(this->coeff, DM_BP, "base_polynomial(const T *, lidia_size_t)"
96 " :: Error in memory allocation (coeff)");
97 this->copy_data(this->coeff, v, this->deg);
98 this->remove_leading_zeros();
99 }
100
101
102
103 template< class T >
base_polynomial(const base_vector<T> v)104 base_polynomial< T >::base_polynomial(const base_vector< T > v)
105 {
106 debug_handler_l(DM_BP, "in constructor "
107 "base_polynomial(const base_vector< T > )", DV_BP);
108 this->deg = v.size()-1;
109 this->coeff = v.get_data(); // Note: the ownership of the array
110 // returned by v.get_data() is taken
111 // by this->coeff!
112 this->remove_leading_zeros();
113 }
114
115
116
117 template< class T >
base_polynomial(const base_polynomial<T> & p)118 base_polynomial< T >::base_polynomial(const base_polynomial< T > &p)
119 {
120 debug_handler_l(DM_BP, "in constructor "
121 "base_polynomial(const base_polynomial< T > )", DV_BP);
122 this->deg = p.deg;
123 if (this->deg < 0)
124 this->coeff = NULL;
125 else {
126 this->coeff = new T[this->deg + 1];
127 memory_handler(this->coeff, DM_BP, "base_polynomial(const T *, lidia_size_t)"
128 " :: Error in memory allocation (coeff)");
129 this->copy_data(coeff, p.coeff, this->deg);
130 }
131 }
132
133
134
135 template< class T >
~base_polynomial()136 base_polynomial< T >::~base_polynomial()
137 {
138 debug_handler_l(DM_BP, "in destructor "
139 "~base_polynomial()", DV_BP);
140 if (this->deg >= 0) {
141 delete[] this->coeff;
142 }
143 }
144
145
146
147 //
148 // comparisons
149 //
150
151 template< class T >
equal(const base_polynomial<T> & b) const152 bool base_polynomial< T >::equal(const base_polynomial< T > &b) const
153 {
154 debug_handler_l(DM_BP, "in member - function "
155 "bool equal(const base_polynomial< T > &)", DV_BP + 4);
156 const T *ap, *bp;
157 lidia_size_t i;
158 if (this->deg != b.deg)
159 return false;
160 for (i = this->deg + 1, ap = this->coeff, bp = b.coeff; i; i--, ap++, bp++)
161 if (*ap != *bp)
162 return false;
163 return true;
164 }
165
166
167
168 //
169 // member functions
170 //
171
172 template< class T >
set_data(const T * d,lidia_size_t l)173 int base_polynomial< T >::set_data (const T * d, lidia_size_t l)
174 {
175 debug_handler_l(DM_BP, "in member - function "
176 "set_data (const T *, lidia_size_t)" , DV_BP + 2);
177 if (l <= 0 || d == NULL)
178 lidia_error_handler_para(l, "l", "l > 0",
179 PRT, "d", "d != NULL",
180 "base_polynomial< T >::"
181 "set_data(const T * d, lidia_size_t l)",
182 DM_BP, LP_ERROR[0]);
183 if (this->deg >= 0)
184 delete[] this->coeff;
185 this->coeff = new T[l];
186 memory_handler(this->coeff, DM_BP, "set_data(const T *, lidia_size_t) :: "
187 "Error in memory allocation (this->coeff)");
188
189 this->deg = l-1;
190
191 for (register lidia_size_t i = 0; i < l; i++)
192 this->coeff[i] = d[i];
193 this->remove_leading_zeros();
194 return 0;
195 }
196
197
198
199 template< class T >
get_data() const200 T* base_polynomial< T >::get_data () const
201 {
202 debug_handler_l(DM_BP, "in member - function "
203 "get_data ()" , DV_BP + 2);
204
205 T * d;
206
207 if (this->deg < 0) {
208 d = new T[1];
209 memory_handler(d, DM_BP, "get_data () :: "
210 "Error in memory allocation (d)");
211 d[0] = 0;
212 return d;
213 }
214
215 d = new T [this->deg+1];
216 memory_handler(d, DM_BP, "get_data () :: "
217 "Error in memory allocation (d)");
218
219 for (register lidia_size_t i = 0; i <= this->deg; i++)
220 d[i] = this->coeff[i];
221 return d;
222 }
223
224
225
226 template< class T >
copy_data(T * d,const T * vd,lidia_size_t al)227 void base_polynomial< T >::copy_data(T * d, const T * vd, lidia_size_t al)
228 {
229 debug_handler_l(DM_BP, "in member - function "
230 "copy_data(T *, const T *, lidia_size_t)", DV_BP + 1);
231 for (lidia_size_t i = al +1; i; i--, d++, vd++)
232 (*d) = (*vd);
233 }
234
235
236
237 template< class T >
remove_leading_zeros()238 void base_polynomial< T >::remove_leading_zeros()
239 {
240 debug_handler_c(DM_BP, "in member - function remove_leading_zeros ()",
241 DV_BP + 1, std::cout << "original degree is " << deg << std::endl << std::flush);
242 T c, *np;
243 lidia_size_t d, i;
244
245 d = this->deg;
246 np = this->coeff + d;
247 #ifdef LIDIA_T_IS_BUILTIN
248 while (d >= 0 && (*np) == 0)
249 #else
250 while (d >= 0 && np->is_zero())
251 #endif
252 d--, np--;
253
254 if (d < 0) {
255 this->deg = d;
256 delete[] this->coeff;
257 this->coeff = NULL;
258 }
259 else if (d != this->deg) {
260 debug_handler_c(DM_BP, "in member - function remove_leading_zeros()",
261 DV_BP + 1, std::cout << "new degree is " << d << std::endl << std::flush);
262 this->deg = d;
263 np = new T[d + 1];
264 memory_handler(np, DM_BP, "remove_leading_zeros() :: "
265 "Error in memory allocation (np)");
266 for (i = 0; i <= d; i++)
267 np[i] = this->coeff[i];
268 delete[] this->coeff;
269 this->coeff = np;
270 }
271 }
272
273
274
275 template< class T >
set_degree(lidia_size_t d)276 void base_polynomial< T >::set_degree(lidia_size_t d)
277 {
278 debug_handler_l(DM_BP, "in member - function "
279 "set_degree(lidia_size_t)", DV_BP + 2);
280
281 if (d < -1)
282 lidia_error_handler_para(d, "d", "d >= -1",
283 "void base_polynomial< T >::"
284 "set_degree(lidia_size_t d)",
285 DM_BP, LP_ERROR[2]);
286
287 if (d == this->deg)
288 return;
289
290 if (d < 0) {
291 delete[] this->coeff; // Note: coeff != NULL, since otherwise d == deg !
292 this->coeff = NULL;
293 this->deg = d;
294 return;
295 }
296
297 T *tmp = this->coeff;
298 this->coeff = new T [d + 1];
299 memory_handler(this->coeff, DM_BP, "set_degree(lidia_size_t d) :: "
300 "Error in memory allocation (this->coeff)");
301
302 register lidia_size_t minimum = (d < this->deg)? d : this->deg;
303
304 for (register lidia_size_t i = 0; i <= minimum; i++)
305 this->coeff[i] = tmp[i];
306
307 if (tmp != NULL)
308 delete[] tmp;
309 this->deg = d;
310 }
311
312
313
314 //
315 // assignment
316 //
317
318 template< class T >
assign(const T & a)319 void base_polynomial< T >::assign(const T & a)
320 {
321 debug_handler_l(DM_BP, "in member - function "
322 "assign (const T &)", DV_BP + 3);
323 if (a == 0) {
324 if (this->deg >= 0) {
325 delete[] this->coeff;
326 this->coeff = NULL;
327 this->deg = -1;
328 }
329 return;
330 }
331 if (this->deg > 0)
332 delete[] this->coeff;
333 if (this->deg != 0) {
334 this->deg = 0;
335 this->coeff = new T[1];
336 memory_handler(this->coeff, DM_BP, "assign(const T &) :: "
337 "Error in memory allocation (coeff)");
338 }
339 this->coeff[0] = a;
340 }
341
342
343
344 template< class T >
assign(const base_polynomial<T> & a)345 void base_polynomial< T >::assign(const base_polynomial< T > &a)
346 {
347 debug_handler_l(DM_BP, "in member - function "
348 "assign (const base_polynomial< T > &)", DV_BP + 3);
349 if (this->deg != a.deg) {
350 if (this->deg >= 0)
351 delete[] this->coeff;
352 this->deg = a.deg;
353 if (this->deg >= 0) {
354 this->coeff = new T[this->deg + 1];
355 memory_handler(this->coeff, DM_BP, "assign(const base_polynomial< T > &) :: "
356 "Error in memory allocation (coeff)");
357 }
358 else this->coeff = NULL;
359 }
360 for (register lidia_size_t i = 0; i <= this->deg; i++)
361 this->coeff[i] = a.coeff[i];
362 }
363
364
365
366 template< class T >
assign_zero()367 void base_polynomial< T >::assign_zero()
368 {
369 debug_handler_l(DM_BP, "in member - function "
370 "assign_zero ()" , DV_BP + 3);
371 if (this->deg >= 0) {
372 this->deg = -1;
373 delete[] this->coeff;
374 this->coeff = NULL;
375 }
376 }
377
378
379
380 template< class T >
assign_one()381 void base_polynomial< T >::assign_one()
382 {
383 debug_handler_l(DM_BP, "in member - function "
384 "assign_one ()" , DV_BP + 3);
385 if (this->deg > 0)
386 delete[] this->coeff;
387 if (this->deg != 0) {
388 this->deg = 0;
389 this->coeff = new T[1];
390 memory_handler(this->coeff, DM_BP, "base_polynomial(const T &) :: "
391 "Error in memory allocation (coeff)");
392 }
393 this->coeff[0] = 1;
394 }
395
396
397
398 template< class T >
assign_x()399 void base_polynomial< T >::assign_x()
400 {
401 debug_handler_l(DM_BP, "in member - function "
402 "assign_x ()" , DV_BP + 3);
403 if (this->deg > 1 || this->deg == 0)
404 delete[] this->coeff;
405 if (this->deg != 1) {
406 this->deg = 1;
407 this->coeff = new T[2];
408 memory_handler(this->coeff, DM_BP, "base_polynomial(const T &) :: "
409 "Error in memory allocation (coeff)");
410 }
411 this->coeff[0] = 0;
412 this->coeff[1] = 1;
413 }
414
415
416
417 template< class T >
swap(base_polynomial<T> & b)418 void base_polynomial< T >::swap(base_polynomial< T > &b)
419 {
420 debug_handler_l(DM_BP, "in member - function "
421 "swap(const base_polynomial< T > &)" , DV_BP + 3);
422 T* tmp = this->coeff;
423 this->coeff = b.coeff;
424 b.coeff = tmp;
425 lidia_size_t deg_tmp = this->deg;
426 this->deg = b.deg;
427 b.deg = deg_tmp;
428 }
429
430
431
432 //
433 // operator overloading
434 //
435
436 template< class T >
operator ()(const T & value) const437 T base_polynomial< T >::operator() (const T & value) const
438 {
439 debug_handler_l(DM_BP, "in operator "
440 "() (const T &)", DV_BP + 5);
441
442 if (this->deg == 0)
443 return this->coeff[0];
444
445 T result;
446 result = 0;
447
448 if (this->deg < 0)
449 return result;
450
451 result = this->coeff[deg];
452 for (lidia_size_t i = this->deg - 1; i >= 0; i--) {
453 LiDIA::multiply(result, result, value);
454 LiDIA::add(result, result, this->coeff[i]);
455 }
456 return result;
457 }
458
459
460
461 //
462 // arithmetic procedures
463 //
464
465 template< class T >
negate(const base_polynomial<T> & a)466 void base_polynomial< T >::negate (const base_polynomial< T > &a)
467 {
468 debug_handler_l(DM_BP, "in member - function "
469 "negate (const base_polynomial< T > &)", DV_BP + 5);
470 register lidia_size_t d = a.deg;
471 this->set_degree(d);
472
473 const T * ap = a.coeff;
474 T * cp = this->coeff;
475
476 for (d++; d; d--, ap++, cp++)
477 LiDIA::negate(*cp, *ap);
478 }
479
480
481
482 template< class T >
add(const base_polynomial<T> & a,const base_polynomial<T> & b)483 void base_polynomial< T >::add(const base_polynomial< T > & a,
484 const base_polynomial< T > & b)
485 {
486 debug_handler_l(DM_BP, "in member - function "
487 "add (const base_polynomial< T > &, "
488 "const base_polynomial< T > &)", DV_BP + 5);
489 const T *ap, *bp;
490 T *cp;
491 lidia_size_t deg_a = a.deg, deg_b = b.deg;
492 lidia_size_t i, min_deg_ab, max_deg_ab;
493
494 if (deg_a > deg_b) {
495 max_deg_ab = deg_a;
496 min_deg_ab = deg_b;
497 }
498 else {
499 max_deg_ab = deg_b;
500 min_deg_ab = deg_a;
501 }
502
503 this->set_degree(max_deg_ab);
504 if (max_deg_ab < 0) return;
505
506 ap = a.coeff;
507 bp = b.coeff;
508 cp = this->coeff;
509
510 for (i = min_deg_ab + 1; i; i--, ap++, bp++, cp++)
511 LiDIA::add(*cp, *ap, *bp);
512
513 if (deg_a > min_deg_ab)
514 for (i = deg_a - min_deg_ab; i; i--, cp++, ap++)
515 *cp = *ap;
516 else if (deg_b > min_deg_ab)
517 for (i = deg_b - min_deg_ab; i; i--, cp++, bp++)
518 *cp = *bp;
519 else
520 this->remove_leading_zeros();
521 }
522
523
524
525 template< class T >
add(const base_polynomial<T> & a,const T & b)526 void base_polynomial< T >::add(const base_polynomial< T > & a, const T & b)
527 {
528 debug_handler_l(DM_BP, "in member - function "
529 "add (base_polynomial< T > &, "
530 "const T &)", DV_BP + 5);
531 if (a.deg < 0) {
532 if (b != 0) {
533 this->set_degree(0);
534 this->coeff[0] = b;
535 return;
536 }
537 this->set_degree(-1);
538 return;
539 }
540 this->set_degree(a.deg);
541
542 const T *ap = a.coeff;
543 T *cp = this->coeff;
544
545 LiDIA::add(*cp , *ap, b);
546 if (a.deg > 0 && this != &a)
547 for (register lidia_size_t i = a.deg; i; i--)
548 *(++cp) = *(++ap);
549 else
550 this->remove_leading_zeros();
551 }
552
553
554
555 template< class T >
add(const T & b,const base_polynomial<T> & a)556 void base_polynomial< T >::add(const T & b,
557 const base_polynomial< T > & a)
558 {
559 debug_handler_l(DM_BP, "in member - function "
560 "add (const T &, "
561 "const base_polynomial< T > &)", DV_BP + 5);
562 if (a.deg < 0) {
563 if (b != 0) {
564 this->set_degree(0);
565 this->coeff[0] = b;
566 return;
567 }
568 this->set_degree(-1);
569 return;
570 }
571 this->set_degree(a.deg);
572
573 const T *ap = a.coeff;
574 T *cp = this->coeff;
575
576 LiDIA::add(*cp, *ap, b);
577 if (a.deg > 0 && &a != this)
578 for (register lidia_size_t i = a.deg; i; i--)
579 *(++cp) = *(++ap);
580 else
581 this->remove_leading_zeros();
582 }
583
584
585
586 template< class T >
subtract(const base_polynomial<T> & a,const base_polynomial<T> & b)587 void base_polynomial< T >::subtract(const base_polynomial< T > & a,
588 const base_polynomial< T > & b)
589 {
590 debug_handler_l(DM_BP, "in member - function "
591 "subtract (const base_polynomial< T > &, "
592 "const base_polynomial< T > &)", DV_BP + 5);
593 const T *ap, *bp;
594 T *cp;
595 lidia_size_t deg_a = a.deg, deg_b = b.deg;
596 lidia_size_t i, min_deg_ab, max_deg_ab;
597
598 if (deg_a > deg_b) {
599 max_deg_ab = deg_a;
600 min_deg_ab = deg_b;
601 }
602 else {
603 max_deg_ab = deg_b;
604 min_deg_ab = deg_a;
605 }
606
607 this->set_degree(max_deg_ab);
608 if (max_deg_ab < 0) return;
609
610 ap = a.coeff;
611 bp = b.coeff;
612 cp = this->coeff;
613 for (i = min_deg_ab + 1; i; i--, ap++, bp++, cp++)
614 LiDIA::subtract(*cp, *ap, *bp);
615
616 if (deg_a > min_deg_ab && this != &a)
617 for (i = deg_a - min_deg_ab; i; i--, cp++, ap++)
618 *cp = *ap;
619 else if (deg_b > min_deg_ab)
620 for (i = deg_b - min_deg_ab; i; i--, cp++, bp++)
621 LiDIA::negate(*cp, *bp);
622 else
623 this->remove_leading_zeros();
624 }
625
626
627
628 template< class T >
subtract(const base_polynomial<T> & a,const T & b)629 void base_polynomial< T >::subtract(const base_polynomial< T > & a,
630 const T & b)
631 {
632 debug_handler_l(DM_BP, "in member - function "
633 "subtract (const base_polynomial< T > &, "
634 "const T &)", DV_BP + 5);
635 if (a.deg < 0) {
636 if (b != 0) {
637 this->set_degree(0);
638 this->coeff[0] = - b;
639 return;
640 }
641 this->set_degree(-1);
642 return;
643 }
644 this->set_degree(a.deg);
645
646 const T *ap = a.coeff;
647 T *cp = this->coeff;
648
649 LiDIA::subtract(*cp, *ap, b);
650 if (a.deg > 0 && this != &a)
651 for (register lidia_size_t i = a.deg; i; i--)
652 *(++cp) = *(++ap);
653 else
654 this->remove_leading_zeros();
655 }
656
657
658
659 template< class T >
subtract(const T & b,const base_polynomial<T> & a)660 void base_polynomial< T >::subtract(const T & b,
661 const base_polynomial< T > & a)
662 {
663 debug_handler_l(DM_BP, "in member - function "
664 "subtract (const T &, "
665 "const base_polynomial< T > &)", DV_BP + 5);
666 if (a.deg < 0) {
667 if (b != 0) {
668 this->set_degree(0);
669 this->coeff[0] = b;
670 return;
671 }
672 this->set_degree(-1);
673 return;
674 }
675 this->set_degree(a.deg);
676
677 const T *ap = a.coeff;
678 T *cp = this->coeff;
679
680 LiDIA::subtract(*cp, b, *ap);
681 if (a.deg > 0)
682 for (register lidia_size_t i = a.deg; i; i--)
683 LiDIA::negate(*(++cp), *(++ap));
684 else
685 this->remove_leading_zeros();
686 }
687
688
689
690 template< class T >
multiply(const base_polynomial<T> & a,const base_polynomial<T> & b)691 void base_polynomial< T >::multiply(const base_polynomial< T > & a,
692 const base_polynomial< T > & b)
693 {
694 debug_handler_l(DM_BP, "in member - function "
695 "multiply (const base_polynomial< T > &, "
696 "const base_polynomial< T > &)", DV_BP + 5);
697 const T *ap, *bp;
698 T * cp, temp;
699 lidia_size_t deg_a = a.deg, deg_b = b.deg;
700
701 if (deg_a < 0 || deg_b < 0) {
702 this->set_degree(-1);
703 return;
704 }
705
706 lidia_size_t i, j, deg_ab = deg_a + deg_b;
707
708 if (this->coeff == a.coeff || this->coeff == b.coeff) {
709 base_polynomial< T > c_temp;
710 c_temp.set_degree(deg_ab);
711
712 for (i = deg_ab + 1, cp = c_temp.coeff; i; i--, cp++)
713 (*cp) = 0;
714
715 for (i = 0, ap = a.coeff; i <= deg_a; i++, ap++)
716 for (j = deg_b + 1, bp = b.coeff, cp = c_temp.coeff + i;
717 j; j--, bp++, cp++) {
718 LiDIA::multiply(temp, *ap, *bp);
719 LiDIA::add(*cp, *cp, temp);
720 }
721 c_temp.remove_leading_zeros();
722 assign(c_temp);
723 }
724 else {
725 this->set_degree(deg_ab);
726
727 for (i = deg_ab + 1, cp = this->coeff; i; i--, cp++)
728 (*cp) = 0;
729
730 for (i = 0, ap = a.coeff; i <= deg_a; i++, ap++)
731 for (j = deg_b + 1, bp = b.coeff, cp = this->coeff + i;
732 j; j--, bp++, cp++) {
733 LiDIA::multiply(temp, *ap, *bp);
734 LiDIA::add(*cp, *cp, temp);
735 }
736 this->remove_leading_zeros();
737 }
738 }
739
740
741
742 template< class T >
multiply(const base_polynomial<T> & a,const T & b)743 void base_polynomial< T >::multiply(const base_polynomial< T > & a,
744 const T & b)
745 {
746 debug_handler_l(DM_BP, "in member - function "
747 "multiply (const base_polynomial< T > &, "
748 "const T &)", DV_BP + 5);
749 if (b == 0) {
750 this->set_degree(-1);
751 return;
752 }
753 const T *ap;
754 T *cp;
755 lidia_size_t deg_a = a.deg;
756
757 this->set_degree(deg_a);
758
759 register lidia_size_t i = 0;
760 for (ap = a.coeff, cp = this->coeff; i <= deg_a; i++, ap++, cp++)
761 LiDIA::multiply(*cp, *ap, b);
762 this->remove_leading_zeros(); // necessary if characteristic != 0
763 }
764
765
766
767 template< class T >
multiply(const T & b,const base_polynomial<T> & a)768 void base_polynomial< T >::multiply(const T & b,
769 const base_polynomial< T > & a)
770 {
771 debug_handler_l(DM_BP, "in member - function "
772 "multiply (const T &, "
773 "const base_polynomial< T > &)", DV_BP + 5);
774 if (b == 0) {
775 this->set_degree(-1);
776 return;
777 }
778 const T *ap;
779 T *cp;
780 lidia_size_t deg_a = a.deg;
781
782 this->set_degree(deg_a);
783
784 register lidia_size_t i = 0;
785 for (ap = a.coeff, cp = this->coeff; i <= deg_a; i++, ap++, cp++)
786 LiDIA::multiply(*cp, b, *ap);
787 this->remove_leading_zeros(); // necessary if characteristic != 0
788 }
789
790
791
792 template< class T >
power(const base_polynomial<T> & a,const bigint & b)793 void base_polynomial< T >::power(const base_polynomial< T > & a,
794 const bigint & b)
795 {
796 debug_handler_l(DM_BP, "in member - function "
797 "power (const base_polynomial< T > &, "
798 "const bigint &)", DV_BP + 5);
799 bigint exponent;
800 base_polynomial< T > multiplier;
801 if (b.is_negative())
802 assign_zero();
803 else if (b.is_zero() || a.is_one())
804 assign_one();
805 else {
806 exponent.assign(b);
807 multiplier.assign(a);
808 assign_one();
809 while (exponent.is_gt_zero()) {
810 if (!exponent.is_even())
811 LiDIA::multiply(*this, *this, multiplier);
812 LiDIA::multiply(multiplier, multiplier, multiplier);
813 exponent.divide_by_2();
814 }
815 }
816 }
817
818
819
820 //
821 // functions
822 //
823
824 template< class T >
derivative(const base_polynomial<T> & a)825 void base_polynomial< T >::derivative(const base_polynomial< T > & a)
826 {
827 debug_handler_l(DM_BP, "in member - function "
828 "derivative (const base_polynomial< T > &)", DV_BP + 5);
829
830 lidia_size_t d = a.deg;
831
832 if (d <= 0) {
833 this->set_degree(-1);
834 return;
835 }
836
837 this->set_degree(d - 1);
838 const T *ap = a.coeff + 1;
839 T *cp = this->coeff;
840 T temp;
841 for (lidia_size_t i = 1; i <= d; i++, cp++, ap++) {
842 temp = i; // necessary, since bigcomplex does not
843 // support automatic cast !!
844 LiDIA::multiply(*cp, *ap, temp);
845 }
846 this->remove_leading_zeros(); // necessary if characteristic != 0
847 }
848
849
850
851 //
852 // input / output
853 //
854
855 template< class T >
read(std::istream & s)856 void base_polynomial< T >::read(std::istream &s)
857 {
858 debug_handler_l(DM_BP, "in member - function "
859 "read(std::istream &)", DV_BP + 6);
860 char c;
861 s >> std::ws >> c;
862 s.putback(c);
863
864 if (c != '[')
865 read_verbose(s);
866 else {
867 base_vector< T > temp;
868 s >> temp;
869 temp.reverse();
870
871 this->set_degree(temp.size() - 1);
872 if (this->deg >= 0) {
873 // Note: the ownership of the array
874 // returned by temp.get_data() is taken
875 // by this->coeff!
876 this->coeff = temp.get_data();
877 }
878 }
879 this->remove_leading_zeros();
880 }
881
882
883
884 template< class T >
read_verbose(std::istream & s)885 void base_polynomial< T >::read_verbose(std::istream & s)
886 {
887 // This function reads a univariate polynomial in any variable.
888 // input format : a_n*x^n+ ... + a_1*x + a_0
889 // e.g. : 3*x^2 + 4*x - 1
890 // Monomials need not be sorted, and powers of x may even appear
891 // repeated, '*' may be omitted and coefficients may follow the variable:
892 // -4 + 8x^5 + 2 - x^2 3 + x^5 + x^5*17
893 // Note however, that the routine will work faster, if the leading monomial
894 // is read first.
895
896 debug_handler_l(DM_BP, "in member-function "
897 "void read_verbose(std::istream &)", DV_BP + 6);
898
899 register lidia_size_t n = 0;
900 lidia_size_t sz;
901 char c;
902
903 this->set_degree(8);
904 for (; n <= this->deg; n++)
905 this->coeff[n] = 0;
906
907 char variable = 0;
908 T coeff_tmp;
909 coeff_tmp = 1;
910 T tmp;
911
912 // Read a monomial, i.e. "x^k" or "- x^k"
913 // or "a*x^k" or "a x^k" or "x^k*a" or "x^k a"
914
915 do {
916 c = s.get();
917 } while (isspace(c) && c != '\n');
918 while (c != '\n' && c != EOF && s.good()) {
919 sz = 0; // Assume we read coeffizient of x^0;
920 if (c == '+') {
921 coeff_tmp = 1;
922 do {
923 c = s.get();
924 } while (isspace(c) && c != '\n');
925 }
926 if (c == '-') {
927 coeff_tmp = -1;
928 do {
929 c = s.get();
930 } while (isspace(c) && c != '\n');
931 }
932 #ifdef POLYREAD_DEBUG
933 std::cout << "\n 1) Now looking at " << c;
934 #endif
935 if (c >= '0' && c <= '9' || c == '(') {
936 s.putback(c);
937 s >> tmp;
938 coeff_tmp *= tmp;
939 do {
940 c = s.get();
941 #ifdef POLYREAD_DEBUG
942 std::cout << ", looking at " << c;
943 #endif
944 } while (isspace(c) && c != '\n');
945 if (c == '*')
946 do {
947 c = s.get();
948 } while (isspace(c) && c != '\n');
949 #ifdef POLYREAD_DEBUG
950 std::cout << "\n coeff_tmp is now " << coeff_tmp;
951 std::cout << ", looking at " << c;
952 #endif
953 }
954 #ifdef POLYREAD_DEBUG
955 std::cout << "\n 2) Now looking at " << c;
956 #endif
957 if ((c >= 'A' && c <= 'Z') || (c >= 'a' && c <= 'z')) {
958 if (variable == 0)
959 variable = c;
960 else if (variable != c)
961 lidia_error_handler_c("base_polynomial", "member function "
962 "read_verbose: The given string is not "
963 "recognized to be a univariate polynomial",
964 std::cout << "Variable name seemed to be " << variable;
965 std::cout << " and now you used " << c << "." << std::endl);
966 do {
967 c = s.get();
968 } while (isspace(c) && c != '\n');
969 #ifdef POLYREAD_DEBUG
970 std::cout << "\n 3) Now looking at " << c;
971 #endif
972
973 if (c != '^') sz = 1;
974 else {
975 s >> sz;
976 #ifdef POLYREAD_DEBUG
977 std::cout << "sz ist " << sz;
978 #endif
979 do {
980 c = s.get();
981 #ifdef POLYREAD_DEBUG
982 std::cout << "\n4') Looking at " << c;
983 #endif
984 } while (isspace(c) && c != '\n');
985 #ifdef POLYREAD_DEBUG
986 std::cout << "\n 5) Now looking at " << c;
987 #endif
988 }
989 if (c == '*') {
990 s >> tmp;
991 coeff_tmp *= tmp;
992 do {
993 c = s.get();
994 } while (isspace(c) && c != '\n');
995 }
996 #ifdef POLYREAD_DEBUG
997 std::cout << "\n 6) Now looking at " << c;
998 #endif
999
1000 if (c >= '0' && c <= '9' || c == '(') {
1001 s.putback(c);
1002 s >> tmp;
1003 #ifdef POLYREAD_DEBUG
1004 std::cout << "\n Old coeff_tmp: " << coeff_tmp;
1005 #endif
1006 coeff_tmp *= tmp;
1007 #ifdef POLYREAD_DEBUG
1008 std::cout << "\n New coeff_tmp: " << coeff_tmp;
1009 #endif
1010 do {
1011 c = s.get();
1012 } while (isspace(c) && c != '\n');
1013 }
1014 }
1015
1016 if (c != '+' && c != '-') {
1017 // No next monomial, so assume end of input is reached
1018 s.putback(c);
1019 c = '\n'; // set c to end--marker
1020 }
1021 if (sz >= n) {
1022 this->set_degree(sz);
1023 for (; n <= this->deg; n++)
1024 this->coeff[n] = 0;
1025 }
1026 this->coeff[sz] += coeff_tmp;
1027
1028 #ifdef POLYREAD_DEBUG
1029 std::cout << "\nSuccessfully read next monomial; Poly now is " << (*this);
1030 std::cout << "\n Now looking at " << c;
1031 #endif
1032 }
1033 this->remove_leading_zeros();
1034 }
1035
1036
1037
1038 // print polynomial
1039 template< class T >
print_verbose(std::ostream & s,char x) const1040 void base_polynomial< T >::print_verbose(std::ostream &s, char x) const
1041 {
1042 debug_handler_l(DM_BP, "in member - function "
1043 "print_verbose(std::ostream &, char)", DV_BP + 6);
1044 lidia_size_t d = this->deg;
1045
1046 if (d < 0)
1047 s << 0;
1048 else if (d == 0)
1049 s << this->coeff[0];
1050 else if (d == 1) {
1051 if (this->coeff[1] == 1)
1052 s << x;
1053 else
1054 s << this->coeff[1] << " * " << x;
1055 if (this->coeff[0] != 0)
1056 s << "+ " << this->coeff[0];
1057 }
1058 else {
1059 if (this->coeff[d] == 1)
1060 s << x << "^" << d;
1061 else
1062 s << this->coeff[d] << " * " << x << "^" << d;
1063 for (register lidia_size_t i = d - 1; i > 1; i--)
1064 if (this->coeff[i] == 1)
1065 s << " + " << x << "^" << i;
1066 else if (this->coeff[i] != 0)
1067 s << " + " << this->coeff[i] << " * " << x << "^" << i;
1068 if (this->coeff[1] == 1)
1069 s << " + " << x;
1070 else if (this->coeff[1] != 0)
1071 s << " + " << this->coeff[1] << " * " << x;
1072 if (this->coeff[0] != 0)
1073 s << " + " << this->coeff[0];
1074 }
1075 }
1076
1077
1078
1079 #undef DV_BP
1080 #undef DM_BP
1081 #undef DV_P
1082 #undef DM_P
1083 #undef LP_ERROR
1084
1085
1086
1087 #ifdef LIDIA_NAMESPACE
1088 # ifndef IN_NAMESPACE_LIDIA
1089 } // end of namespace LiDIA
1090 # endif
1091 #endif
1092
1093
1094
1095 #endif // LIDIA_POLY_INTERN_CC_GUARD_
1096