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	: Michael Jacobson (MJ)
14 //	Changes	: See CVS log
15 //
16 //==============================================================================================
17 
18 
19 #ifdef HAVE_CONFIG_H
20 # include	"config.h"
21 #endif
22 #include	"LiDIA/quadratic_order.h"
23 #include	"LiDIA/qi_class.h"
24 #include	"LiDIA/qi_class_real.h"
25 #include	"LiDIA/number_fields/qo_list.h"
26 
27 
28 
29 #ifdef LIDIA_NAMESPACE
30 namespace LiDIA {
31 #endif
32 
33 
34 
35 //<TPf>
36 qo_node *qi_class::current_order = NULL;
37 bigint qi_class::Delta = 0;
38 bigint qi_class::rootD = 0;
39 xbigfloat qi_class::rd = 0.0;
40 bigint qi_class::PEA_L = 0;
41 bigint qi_class::ordmult = 0;
42 rational_factorization qi_class::omfact;
43 int qi_class::info = 0;
44 int qi_class::do_verify = 0;
45 long qi_class::xprec = 0;
46 //</TPf>
47 
48 
49 quadratic_order quadratic_order::zero_QO;
50 
qo_l()51 qo_list& quadratic_order::qo_l() {
52     static qo_list quadratic_order_list;
53 
54     return quadratic_order_list;
55 }
56 
57 int quadratic_order::info = 0;
58 int quadratic_order::do_verify = 0;
59 
60 // FOR TABLE GENERATION
61 bool quadratic_order::special2 = false;
62 
63 
64 const long OQvals[] = {1663, 4177, 7523, 11657, 16477, 21991, 28151, 34913, 42293,
65 		       50261, 58787, 67883, 77527, 87719, 98419, 109661, 121403, 133649, 146407, 159667};
66 
67 const long OQvals_cnum[] = {2269, 5741, 10427, 16183, 22901, 30631,
68 			    39209, 48731, 59063, 70237, 82223, 95009, 108571, 122921, 137983, 153817, 170341,
69 			    187631, 205589, 224261};
70 
71 const long OQvals_cfunc[] = {630000, 825000, 995000, 1145000,
72 			     1290000, 1425000, 1555000, 1675000, 1795000, 1910000, 2020000, 2125000, 2230000,
73 			     2335000, 2435000, 2530000, 2625000, 2720000, 2810000, 2900000};
74 
75 
76 
77 //
78 // constructor
79 //     - sets resets all the class members which are vectors
80 //
81 
quadratic_order()82 quadratic_order::quadratic_order()
83 {
84 	debug_handler("quadratic_order", "entering quadratic_order()");
85 
86 	SPL.resize(2, 70000);
87 	CL.set_mode(EXPAND);
88 	gens.set_mode(EXPAND);
89 	fact_base.set_mode(EXPAND);
90 	contributors.set_mode(EXPAND);
91 	minima_qn.set_mode(EXPAND);
92 	CL.reset();
93 	gens.reset();
94 	fact_base.reset();
95 	contributors.reset();
96 	minima_qn.reset();
97 	subused = false;
98 	U.resize(1, 1);
99 	UINV.resize(1, 1);
100 
101 	R_xbig.assign_zero(); // MM
102 
103 	if (!qi_class::get_current_order_ptr())
104 		qi_class::set_current_order(zero_QO);
105 
106 	debug_handler("quadratic_order", "leaving quadratic_order()");
107 }
108 
109 
110 
111 //
112 // constructor
113 //     - sets resets all the class members which are vectors and initializes
114 //       with discriminant D.  If D is not a quadratic discriminant, 0 will
115 //       be used.
116 //
117 
quadratic_order(const long D)118 quadratic_order::quadratic_order(const long D)
119 {
120 	debug_handler("quadratic_order", "entering quadratic_order(long)");
121 
122 	long oprec = bigfloat::get_precision();
123 	bigfloat temp;
124 
125 	// initialization
126 	SPL.resize(2, 70000);
127 	CL.set_mode(EXPAND);
128 	gens.set_mode(EXPAND);
129 	fact_base.set_mode(EXPAND);
130 	contributors.set_mode(EXPAND);
131 	minima_qn.set_mode(EXPAND);
132 	CL.reset();
133 	gens.reset();
134 	fact_base.reset();
135 	contributors.reset();
136 	minima_qn.reset();
137 	DET = 1;
138 	AMAT.kill();
139 	AMAT.set_zero_element(0);
140 	AMAT.set_representation(matrix_flags::sparse_representation);
141 	AMAT.set_orientation(matrix_flags::row_oriented);
142 	subused = false;
143 
144 	RHNF.kill();
145 	RHNF.set_representation(matrix_flags::sparse_representation);
146 	U.kill();
147 	UINV.kill();
148 	U.resize(1, 1);
149 	UINV.resize(1, 1);
150 
151 	if (is_quadratic_discriminant(D))
152 		Delta.assign(D);
153 	else
154 		Delta.assign_zero();
155 
156 	prec = static_cast<long>(decimal_length(Delta) >> 1) + 10;
157 	xprec = static_cast<long>(std::ceil(static_cast<double>(prec)*std::log(10.0L) / std::log(2.0L)));
158 
159 	bigfloat::set_precision(prec);
160 
161 	R.assign_zero();
162 	R_xbig.assign_zero(); // MM
163 	h.assign_zero();
164 	L.assign_zero();
165 	Cfunc.assign_zero();
166 	FI.assign_zero();
167 	Q = 0;
168 
169 	if (Delta.is_le_zero())
170 		newton_root(nucomp_bound, abs(Delta >> 2), 4);
171 	else {
172 		nucomp_bound.assign_zero();
173 		min_log = log((sqrt(bigfloat(Delta-4)) + sqrt(bigfloat(Delta))) /
174 			      bigfloat(2.0)) / bigfloat(2.0);
175 	}
176 
177 	bigfloat::set_precision(oprec);
178 
179 	if (!qi_class::get_current_order_ptr())
180 		qi_class::set_current_order(zero_QO);
181 
182 	qo_l().set_last(this);
183 
184 	debug_handler("quadratic_order", "leaving quadratic_order(long)");
185 }
186 
187 
188 
189 //
190 // constructor
191 //     - sets resets all the class members which are vectors and initializes
192 //       with discriminant D.  If D is not a quadratic discriminant, 0 will
193 //       be used.
194 //
195 
quadratic_order(const bigint & D)196 quadratic_order::quadratic_order(const bigint & D)
197 {
198 	debug_handler("quadratic_order", "entering quadratic_order(bigint)");
199 
200 	long oprec = bigfloat::get_precision();
201 	bigfloat temp;
202 
203 	// initialization
204 	SPL.resize(2, 70000);
205 	CL.set_mode(EXPAND);
206 	gens.set_mode(EXPAND);
207 	fact_base.set_mode(EXPAND);
208 	contributors.set_mode(EXPAND);
209 	minima_qn.set_mode(EXPAND);
210 
211 	CL.reset();
212 	gens.reset();
213 	fact_base.reset();
214 	contributors.reset();
215 	minima_qn.reset();
216 	DET = 1;
217 	AMAT.kill();
218 	AMAT.set_zero_element(0);
219 	AMAT.set_representation(matrix_flags::sparse_representation);
220 	AMAT.set_orientation(matrix_flags::row_oriented);
221 	subused = false;
222 
223 	RHNF.kill();
224 	RHNF.set_representation(matrix_flags::sparse_representation);
225 	U.kill();
226 	UINV.kill();
227 	U.resize(1, 1);
228 	UINV.resize(1, 1);
229 
230 	if (is_quadratic_discriminant(D))
231 		Delta.assign(D);
232 	else
233 		Delta.assign_zero();
234 
235 	prec = static_cast<long>(decimal_length(Delta) >> 1) + 10;
236 	xprec = static_cast<long>(std::ceil(static_cast<double>(prec)*std::log(10.0L) / std::log(2.0L)));
237 
238 	bigfloat::set_precision(prec);
239 
240 	R.assign_zero();
241 	R_xbig.assign_zero(); // MM
242 	h.assign_zero();
243 	L.assign_zero();
244 	Cfunc.assign_zero();
245 	FI.assign_zero();
246 	Q = 0;
247 
248 	if (Delta.is_le_zero())
249 		newton_root(nucomp_bound, abs(Delta >> 2), 4);
250 	else {
251 		nucomp_bound.assign_zero();
252 		min_log = log((sqrt(bigfloat(Delta-4)) + sqrt(bigfloat(Delta))) /
253 			      bigfloat(2.0)) / bigfloat(2.0);
254 	}
255 
256 	bigfloat::set_precision(oprec);
257 
258 	if (!qi_class::get_current_order_ptr())
259 		qi_class::set_current_order(zero_QO);
260 
261 	qo_l().set_last(this);
262 
263 	debug_handler("quadratic_order", "leaving quadratic_order(bigint)");
264 }
265 
266 
267 
268 //
269 // constructor
270 //     - sets resets all the class members which are vectors and initializes
271 //       with a copy of QO.
272 //
273 
quadratic_order(const quadratic_order & QO)274 quadratic_order::quadratic_order(const quadratic_order & QO)
275 {
276 	debug_handler("quadratic_order", "quadratic_order(quadratic_order)");
277 
278 	long oprec = bigfloat::get_precision();
279 
280 	SPL.resize(2, 70000);
281 	CL.set_mode(EXPAND);
282 	gens.set_mode(EXPAND);
283 	fact_base.set_mode(EXPAND);
284 	contributors.set_mode(EXPAND);
285 	minima_qn.set_mode(EXPAND);
286 
287 	Delta.assign(QO.Delta);
288 	prec = QO.prec;
289 	xprec = QO.xprec;
290 	bigfloat::set_precision(prec);
291 
292 	R.assign(QO.R);
293 	R_xbig.assign(QO.R_xbig); //MM
294 	R_xbig_prec = QO.R_xbig_prec;
295 	h.assign(QO.h);
296 	L.assign(QO.L);
297 	Cfunc.assign(QO.Cfunc);
298 	CL = QO.CL;
299 	gens = QO.gens;
300 	hfact.assign(QO.hfact);
301 	disc_fact.assign(QO.disc_fact);
302 	fact_base = QO.fact_base;
303 	contributors = QO.contributors;
304 	minima_qn = QO.minima_qn;
305 	fund_unit = QO.fund_unit;
306 	DET = QO.DET;
307 	AMAT = QO.AMAT;
308 
309 	subused = QO.subused;
310 	prin_list = QO.prin_list;
311 	FI = QO.FI;
312 	Q = QO.Q;
313 	nucomp_bound = QO.nucomp_bound;
314 	RHNF = QO.RHNF;
315 	U = QO.U;
316 	UINV = QO.UINV;
317 
318 	bigfloat::set_precision(oprec);
319 
320 	qo_l().set_last(this);
321 }
322 
323 
324 
325 //
326 // destuctor
327 //     - deletes this order from the list of current orders.
328 //
329 
~quadratic_order()330 quadratic_order::~quadratic_order()
331 {
332 	debug_handler("quadratic_order", "~quadratic_order");
333 
334 	TR.kill();
335 	qo_l().nullify(this);
336 }
337 
338 
339 
340 //
341 // quadratic_order::verbose()
342 //
343 // Task:
344 //      sets the verbosity of commands.  Currently, the following levels are
345 //      supported:
346 //         0 - nothing
347 //         1 - run times and some run-time data (data only for subexp algs.)
348 //         >1 - debug information for subexponential algorithms
349 //
350 
351 void
verbose(int state)352 quadratic_order::verbose(int state)
353 {
354 	debug_handler("quadratic_order", "verbose");
355 
356 	if (state <= 0) {
357 		info = 0;
358 	}
359 	else {
360 		info = state;
361 	}
362 
363 	//MM, seems to be unused, qo_info = info;
364 	qo_info = info;
365 }
366 
367 
368 
369 //
370 // quadratic_order::verification()
371 //
372 // Task:
373 //      sets the level of verifications.  Currently, the following levels are
374 //      supported:
375 //         0 - nothing
376 //         1-3 - varies levels of verification
377 //
378 
379 void
verification(int level)380 quadratic_order::verification(int level)
381 {
382 	debug_handler("quadratic_order", "verification");
383 
384 	if (level <= 0)
385 		do_verify = 0;
386 	else
387 		do_verify = level;
388 }
389 
390 
391 
392 //
393 // quadratic_order::assign(long)
394 //
395 // Task:
396 //    set to the quadratic order of discriminant D.  If D is not a quadratic
397 //    discriminant, 0 the quadratic order will be set to the empty order
398 //    (discriminant 0) and false will be returned.
399 //
400 
401 bool
assign(const long D)402 quadratic_order::assign(const long D)
403 {
404 	debug_handler("quadratic_order", "assign(long)");
405 
406 	return assign(bigint(D));
407 }
408 
409 
410 
411 //
412 // quadratic_order::assign(bigint)
413 //
414 // Task:
415 //    set to the quadratic order of discriminant D.  If D is not a quadratic
416 //    discriminant, 0 the quadratic order will be set to the empty order
417 //    (discriminant 0) and false will be returned.
418 //
419 
420 bool
assign(const bigint & D)421 quadratic_order::assign(const bigint & D)
422 {
423 	debug_handler("quadratic_order", "assign(bigint)");
424 
425 	long oprec = bigfloat::get_precision();
426 	bigfloat temp;
427 	bool is_disc;
428 
429 	is_disc = is_quadratic_discriminant(D);
430 
431 	if (is_disc)
432 		Delta.assign(D);
433 	else
434 		Delta.assign_zero();
435 
436 	prec = static_cast<long>(decimal_length(Delta) >> 1) + 10;
437 	xprec = static_cast<long>(std::ceil(static_cast<double>(prec)*std::log(10.0L) / std::log(2.0L)));
438 
439 	bigfloat::set_precision(prec);
440 
441 	CL.set_mode(EXPAND);
442 	gens.set_mode(EXPAND);
443 	fact_base.set_mode(EXPAND);
444 	contributors.set_mode(EXPAND);
445 	minima_qn.set_mode(EXPAND);
446 
447 	R.assign_zero();
448 	R_xbig.assign_zero(); // MM
449 	h.assign_zero();
450 	L.assign_zero();
451 	Cfunc.assign_zero();
452 	CL.reset();
453 	gens.reset();
454 	fact_base.reset();
455 	contributors.reset();
456 	TR.kill();
457 	minima_qn.reset();
458 	fund_unit.reset();
459 	DET = 1;
460 	AMAT.kill();
461 	AMAT.set_zero_element(0);
462 	AMAT.set_representation(matrix_flags::sparse_representation);
463 	AMAT.set_orientation(matrix_flags::row_oriented);
464 
465 	subused = false;
466 	hfact = rational_factorization();
467 	disc_fact = rational_factorization();
468 	RHNF.kill();
469 	RHNF.set_representation(matrix_flags::sparse_representation);
470 	U.kill();
471 	UINV.kill();
472 	U.resize(1, 1);
473 	UINV.resize(1, 1);
474 	prin_list.empty();
475 	FI.assign_zero();
476 	Q = 0;
477 
478 	if (Delta.is_le_zero())
479 		newton_root(nucomp_bound, abs(Delta >> 2), 4);
480 	else {
481 		nucomp_bound.assign_zero();
482 		min_log = log((sqrt(bigfloat(Delta-4)) + sqrt(bigfloat(Delta))) /
483 			      bigfloat(2.0)) / bigfloat(2.0);
484 	}
485 
486 	sieve.reset();
487 
488 	bigfloat::set_precision(oprec);
489 
490 	// set last-used quadratic_order
491 	qo_l().set_last(this);
492 
493 	return is_disc;
494 }
495 
496 
497 
498 //
499 // quadratic_order::assign(quadratic_order)
500 //
501 // Task:
502 //    set to a copy of QO.
503 //
504 
505 void
assign(const quadratic_order & QO)506 quadratic_order::assign(const quadratic_order & QO)
507 {
508 	debug_handler("quadratic_order", "assign(quadratic_order)");
509 
510 	long oprec = bigfloat::get_precision();
511 
512 	Delta.assign(QO.Delta);
513 	prec = QO.prec;
514 	xprec = QO.xprec;
515 	bigfloat::set_precision(prec);
516 
517 	R.assign(QO.R);
518 	R_xbig.assign(QO.R_xbig); //MM
519 	R_xbig_prec = QO.R_xbig_prec;
520 	h.assign(QO.h);
521 	L.assign(QO.L);
522 	Cfunc.assign(QO.Cfunc);
523 	CL = QO.CL;
524 	gens = QO.gens;
525 	fact_base = QO.fact_base;
526 	contributors = QO.contributors;
527 	minima_qn = QO.minima_qn;
528 	fund_unit = QO.fund_unit;
529 	DET = QO.DET;
530 	AMAT = QO.AMAT;
531 
532 	subused = QO.subused;
533 	hfact.assign(QO.hfact);
534 	disc_fact.assign(QO.disc_fact);
535 	prin_list = QO.prin_list;
536 	FI = QO.FI;
537 	Q = QO.Q;
538 	nucomp_bound = QO.nucomp_bound;
539 	RHNF = QO.RHNF;
540 	U = QO.U;
541 	UINV = QO.UINV;
542 	min_log = QO.min_log;
543 
544 	bigfloat::set_precision(oprec);
545 
546 	// set last-used quadratic_order
547 	qo_l().set_last(this);
548 }
549 
550 
551 
552 //
553 // operator =
554 //
555 // Task:
556 //    set to a copy of QO.
557 //
558 
operator =(const quadratic_order & QO)559 quadratic_order & quadratic_order::operator = (const quadratic_order & QO)
560 {
561 	debug_handler("quadratic_order", "operator = ");
562 
563 	this->assign(QO);
564 	return *this;
565 }
566 
567 
568 
569 //
570 // quadratic_order::discriminant()
571 //
572 // Task:
573 //      returns the discriminant
574 //
575 
576 bigint
discriminant() const577 quadratic_order::discriminant() const
578 {
579 	return Delta;
580 }
581 
582 
583 
584 //
585 // quadratic_order::nu_bound()
586 //
587 // Task:
588 //      returns the bound used for NUCOMP and NUDUPL
589 //
590 
591 bigint
nu_bound() const592 quadratic_order::nu_bound() const
593 {
594 	return nucomp_bound;
595 }
596 
597 
598 
599 //
600 // quadratic_order::is_zero()
601 //
602 // Task:
603 //      tests if the quadratic order is the empty order
604 //
605 
606 bool
is_zero() const607 quadratic_order::is_zero() const
608 {
609 	debug_handler("quadratic_order", "is_zero");
610 
611 	return Delta.is_zero();
612 }
613 
614 
615 
616 //
617 // quadratic_order::is_equal()
618 //
619 // Task:
620 //      tests if the quadratic orders are equal
621 //
622 
623 bool
is_equal(const quadratic_order & QO) const624 quadratic_order::is_equal(const quadratic_order & QO) const
625 {
626 	debug_handler("quadratic_order", "is_equal");
627 
628 	return !Delta.compare(QO.Delta);
629 }
630 
631 
632 
633 //
634 // quadratic_order::is_subset()
635 //
636 // Task:
637 //      tests if the quadratic order is a subset of QO2.
638 //
639 
640 bool
is_subset(quadratic_order & QO2)641 quadratic_order::is_subset(quadratic_order & QO2)
642 {
643 	debug_handler("quadratic_order", "is_subset");
644 
645 	return ((this->is_proper_subset(QO2)) || (this->is_equal(QO2)));
646 }
647 
648 
649 
650 //
651 // quadratic_order::is_proper_subset()
652 //
653 // Task:
654 //      tests if the quadratic order is a proper subset of QO2.
655 //
656 
657 bool
is_proper_subset(quadratic_order & QO2)658 quadratic_order::is_proper_subset(quadratic_order & QO2)
659 {
660 	debug_handler("quadratic_order", "is_proper_subset");
661 
662 	bigint c1, c2, m1, m2, r;
663 
664 	// set last-used quadratic_order
665 	qo_l().set_last(this);
666 
667 	c1.assign(conductor());
668 	c2.assign(QO2.conductor());
669 	divide(m1, Delta, c1*c1);
670 	divide(m2, QO2.Delta, c2*c2);
671 	remainder(r, c1, c2);
672 
673 	return ((m1 == m2) && (r.is_zero()) && (c1 != c2));
674 }
675 
676 
677 
678 //
679 // operator ==
680 //
681 // Task:
682 //      tests if QO1 and QO2 are equal (same discriminant)
683 //
684 
operator ==(const quadratic_order & QO1,const quadratic_order & QO2)685 bool operator == (const quadratic_order & QO1, const quadratic_order & QO2)
686 {
687 	debug_handler("quadratic_order", "operator == ");
688 
689 	return !QO1.Delta.compare(QO2.Delta);
690 }
691 
692 
693 
694 //
695 // operator !=
696 //
697 // Task:
698 //      tests if QO1 and QO2 are not equal
699 //
700 
operator !=(const quadratic_order & QO1,const quadratic_order & QO2)701 bool operator != (const quadratic_order & QO1, const quadratic_order & QO2)
702 {
703 	debug_handler("quadratic_order", "operator != ");
704 
705 	return QO1.Delta.compare(QO2.Delta);
706 }
707 
708 
709 
710 //
711 // operator <=
712 //
713 // Task:
714 //      tests if QO1 is a subset of QO2
715 //
716 
operator <=(quadratic_order & QO1,quadratic_order & QO2)717 bool operator <= (quadratic_order & QO1, quadratic_order & QO2)
718 {
719 	debug_handler("quadratic_order", "operator <= ");
720 
721 	return (QO1.is_subset(QO2));
722 }
723 
724 
725 
726 //
727 // operator <
728 //
729 // Task:
730 //      tests if QO1 is a proper subset of QO2
731 //
732 
operator <(quadratic_order & QO1,quadratic_order & QO2)733 bool operator < (quadratic_order & QO1, quadratic_order & QO2)
734 {
735 	debug_handler("quadratic_order", "operator < ");
736 
737 	return (QO1.is_proper_subset(QO2));
738 }
739 
740 
741 
742 //
743 // operator >=
744 //
745 // Task:
746 //      tests if QO2 is a subset of QO1
747 //
748 
operator >=(quadratic_order & QO1,quadratic_order & QO2)749 bool operator >= (quadratic_order & QO1, quadratic_order & QO2)
750 {
751 	debug_handler("quadratic_order", "operator >= ");
752 
753 	return (QO2.is_subset(QO1));
754 }
755 
756 
757 
758 //
759 // operator >
760 //
761 // Task:
762 //      tests if QO2 is a proper subset of QO1
763 //
764 
operator >(quadratic_order & QO1,quadratic_order & QO2)765 bool operator > (quadratic_order & QO1, quadratic_order & QO2)
766 {
767 	debug_handler("quadratic_order", "operator >");
768 
769 	return (QO2.is_proper_subset(QO1));
770 }
771 
772 
773 
774 //
775 // operator !
776 //
777 // Task:
778 //      tests if QO is the empty order
779 //
780 
operator !(const quadratic_order & QO)781 bool operator ! (const quadratic_order & QO)
782 {
783 	debug_handler("quadratic_order", "operator !");
784 
785 	return (QO.is_zero());
786 }
787 
788 
789 
790 //
791 // is_quadratic_discriminant(long)
792 //
793 // Task:
794 //     returns true if D is a quadratic discriminant (congruent to 0 or 1 mod
795 //     4 and not a perfect square).
796 //
797 
798 bool
is_quadratic_discriminant(const long D)799 is_quadratic_discriminant(const long D)
800 {
801 	debug_handler("quadratic_order", "is_quadratic_discriminant(long)");
802 
803 	return is_quadratic_discriminant(bigint(D));
804 }
805 
806 
807 
808 //
809 // is_quadratic_discriminant(bigint)
810 //
811 // Task:
812 //     returns true if D is a quadratic discriminant (congruent to 0 or 1 mod
813 //     4 and not a perfect square).
814 //
815 
816 bool
is_quadratic_discriminant(const bigint & D)817 is_quadratic_discriminant(const bigint & D)
818 {
819 	debug_handler("quadratic_order", "is_quadratic_discriminant(bigint)");
820 
821 	long k;
822 	bool is_disc;
823 
824 	is_disc = true;
825 
826 	// testing whether D is = 0,1 (mod 4)
827 	k = remainder(D, 4);
828 	if (k < 0)
829 		k += 4;
830 	if (k > 1)
831 		is_disc = false;
832 	else if (D.is_gt_zero()) {
833 		// testing whether D is a perfect square
834 		if (D.is_square())
835 			is_disc = false;
836 	}
837 
838 	return is_disc;
839 }
840 
841 
842 
843 //
844 // quadratic_order::is_imaginary()
845 //
846 // Task:
847 //      returns true if the quadratic order is imaginary (discriminant is less
848 //      than 0)
849 //
850 
851 bool
is_imaginary() const852 quadratic_order::is_imaginary() const
853 {
854 	debug_handler("quadratic_order", "is_imaginary");
855 
856 	return Delta.is_lt_zero();
857 }
858 
859 
860 
861 //
862 // quadratic_order::is_real()
863 //
864 // Task:
865 //      returns true if the quadratic order is real (discriminant is greater
866 //      than 0)
867 //
868 
869 bool
is_real() const870 quadratic_order::is_real() const
871 {
872 	debug_handler("quadratic_order", "is_real");
873 
874 	return Delta.is_gt_zero();
875 }
876 
877 
878 
879 //
880 // quadratic_order::is_R_computed()
881 //
882 // Task:
883 //      returns true if the regulator has already been computed.
884 //
885 
886 bool
is_R_computed() const887 quadratic_order::is_R_computed() const
888 {
889 	debug_handler("quadratic_order", "is_R_computed");
890 
891 	long oprec = bigfloat::get_precision();
892 	bool iscomp;
893 
894 	bigfloat::set_precision(prec);
895 	iscomp = (R.is_gt_zero()) || (Delta.is_zero());
896 	bigfloat::set_precision(oprec);
897 
898 	return iscomp;
899 }
900 
901 
902 
903 //
904 // quadratic_order::is_h_computed()
905 //
906 // Task:
907 //      returns true if the class number has already been computed.
908 //
909 
910 bool
is_h_computed() const911 quadratic_order::is_h_computed() const
912 {
913 	debug_handler("quadratic_order", "is_h_computed");
914 
915 	return (h.is_gt_zero()) || (Delta.is_zero());
916 }
917 
918 
919 
920 //
921 // quadratic_order::is_L_computed()
922 //
923 // Task:
924 //      returns true if the L(1) function has already been computed.
925 //
926 
927 bool
is_L_computed() const928 quadratic_order::is_L_computed() const
929 {
930 	debug_handler("quadratic_order", "is_L_computed");
931 
932 	long oprec = bigfloat::get_precision();
933 	bool iscomp;
934 
935 	bigfloat::set_precision(prec);
936 	iscomp = (L.is_gt_zero()) || (Delta.is_zero());
937 	bigfloat::set_precision(oprec);
938 
939 	return iscomp;
940 }
941 
942 
943 
944 //
945 // quadratic_order::is_C_computed()
946 //
947 // Task:
948 //      returns true if the C(Delta) function has already been computed.
949 //
950 
951 bool
is_C_computed() const952 quadratic_order::is_C_computed() const
953 {
954 	debug_handler("quadratic_order", "is_C_computed");
955 
956 	long oprec = bigfloat::get_precision();
957 	bool iscomp;
958 
959 	bigfloat::set_precision(prec);
960 	iscomp = (Cfunc.is_gt_zero()) || (Delta.is_zero());
961 	bigfloat::set_precision(oprec);
962 
963 	return iscomp;
964 }
965 
966 
967 
968 //
969 // quadratic_order::is_CL_computed()
970 //
971 // Task:
972 //      returns true if the class group has already been computed.
973 //
974 
975 bool
is_CL_computed() const976 quadratic_order::is_CL_computed() const
977 {
978 	debug_handler("quadratic_order", "is_CL_computed");
979 
980 	return (CL.size() > 0) || (Delta.is_zero());
981 }
982 
983 
984 
985 //
986 // quadratic_order::is_subexp_computed()
987 //
988 // Task:
989 //      returns true if the class group has already been computed, and it was
990 //      done with the subexponential algorithm.
991 //
992 
993 bool
is_subexp_computed() const994 quadratic_order::is_subexp_computed() const
995 {
996 	debug_handler("quadratic_order", "is_subexp_computed");
997 
998 	return (subused);
999 }
1000 
1001 
1002 
1003 //<MM>
1004 //
1005 // quadratic_order::is_fundamental_unit_computed()
1006 //
1007 // Task:
1008 //      returns true if the fundamental unit has already been computed.
1009 //
1010 
1011 bool quadratic_order
is_fundamental_unit_computed() const1012 ::is_fundamental_unit_computed() const
1013 {
1014 	debug_handler("quadratic_order", "is_fundamental_unit_computed");
1015 
1016 	return fund_unit.is_initialized();
1017 }
1018 
1019 
1020 
1021 //<MM>
1022 
1023 
1024 
1025 
1026 //
1027 // quadratic_order::bach_bound()
1028 //
1029 // Task:
1030 //      returns Bach's bound on the norms of the prime ideals which generate
1031 //      the class group.
1032 //
1033 
1034 int
bach_bound()1035 quadratic_order::bach_bound()
1036 {
1037 	debug_handler("quadratic_order", "bach_bound()");
1038 
1039 	return bach_bound(conductor());
1040 }
1041 
1042 
1043 
1044 //
1045 // quadratic_order::bach_bound(const bigint & f)
1046 //
1047 // Task:
1048 //      returns Bach's bound on the norms of the prime ideals which generate
1049 //      the class group, using Table 3 of Bach90.
1050 //
1051 
1052 int
bach_bound(const bigint & f)1053 quadratic_order::bach_bound(const bigint & f)
1054 {
1055 	debug_handler("quadratic_order", "bach_bound(bigint)");
1056 
1057 	bigfloat temp, C1, C2;
1058 	bigint D;
1059 	int BachBound, dec;
1060 	long oprec = bigfloat::get_precision();
1061 
1062 	bigfloat::set_precision(prec);
1063 
1064 	if (f.is_one()) {
1065 		// order is maximal
1066 		dec = decimal_length(Delta*Delta);
1067 
1068 		if (dec >= 1000) {
1069 			C1.assign(1.011);
1070 			C2.assign(14.236);
1071 		}
1072 		else if (dec >= 100) {
1073 			C1.assign(1.044);
1074 			C2.assign(2.336);
1075 		}
1076 		else if (dec >= 50) {
1077 			C1.assign(1.065);
1078 			C2.assign(0.439);
1079 		}
1080 		else if (dec >= 20) {
1081 			C1.assign(1.110);
1082 			C2.assign(-1.367);
1083 		}
1084 		else if (dec >= 10) {
1085 			C1.assign(1.166);
1086 			C2.assign(-2.350);
1087 		}
1088 		else if (dec >= 5) {
1089 			C1.assign(1.259);
1090 			C2.assign(-3.077);
1091 		}
1092 		else {
1093 			C1.assign(1.571);
1094 			C2.assign(-3.570);
1095 		}
1096 
1097 		temp = C1 * log(bigfloat(Delta*Delta)) + C2;
1098 		square(temp, temp);
1099 		floor(temp, temp);
1100 		temp.intify(BachBound);
1101 	}
1102 	else {
1103 		D = Delta*Delta / (f*f);
1104 		dec = decimal_length(D);
1105 
1106 		if (dec >= 1000) {
1107 			C1.assign(1.076);
1108 			C2.assign(4.206);
1109 		}
1110 		else if (dec >= 100) {
1111 			C1.assign(1.287);
1112 			C2.assign(-0.887);
1113 		}
1114 		else if (dec >= 50) {
1115 			C1.assign(1.420);
1116 			C2.assign(-1.701);
1117 		}
1118 		else if (dec >= 20) {
1119 			C1.assign(1.682);
1120 			C2.assign(-2.471);
1121 		}
1122 		else if (dec >= 10) {
1123 			C1.assign(1.974);
1124 			C2.assign(-2.885);
1125 		}
1126 		else if (dec >= 5) {
1127 			C1.assign(2.379);
1128 			C2.assign(-3.189);
1129 		}
1130 		else {
1131 			C1.assign(3.166);
1132 			C2.assign(-3.468);
1133 		}
1134 
1135 		temp = C1 * log(bigfloat(D)) + C2;
1136 		square(temp, temp);
1137 		floor(temp, temp);
1138 		temp.intify(BachBound);
1139 	}
1140 
1141 	bigfloat::set_precision(oprec);
1142 
1143 	return BachBound;
1144 }
1145 
1146 
1147 
1148 //
1149 // swap
1150 //
1151 // Task:
1152 //      swaps Q1 and Q2
1153 //
1154 
1155 void
swap(quadratic_order & Q1,quadratic_order & Q2)1156 swap(quadratic_order & Q1, quadratic_order & Q2)
1157 {
1158 	debug_handler("quadratic_order", "swap");
1159 
1160 	quadratic_order C;
1161 
1162 	C.assign(Q1);
1163 	Q1.assign(Q2);
1164 	Q2.assign(C);
1165 }
1166 
1167 
1168 
1169 //
1170 // quadratic_order::conductor()
1171 //
1172 // Task:
1173 //      computes the conductor of the order, i.e., the integer f such that
1174 //      Delta = f^2 Delta_0 for some fundamental discriminant Delta_0.
1175 //      The discriminant is factored if it has not been already.
1176 //
1177 
1178 bigint
conductor()1179 quadratic_order::conductor()
1180 {
1181 	debug_handler("quadratic_order", "conductor");
1182 
1183 	bigint cond, b, temp;
1184 	int e;
1185 	long m4;
1186 	lidia_size_t num_facts, i;
1187 
1188 	// set last-used quadratic_order
1189 	qo_l().set_last(this);
1190 
1191 	cond.assign_one();
1192 	factor_discriminant();
1193 	num_facts = disc_fact.no_of_comp();
1194 	for (i = 0; i < num_facts; ++i) {
1195 		b = disc_fact.base(i);
1196 		e = disc_fact.exponent(i);
1197 		if ((b > 0) && (e > 1)) {
1198 			if (b == 2) {
1199 				if ((e % 2) == 1)
1200 					shift_right(temp, Delta, e-1);
1201 				else
1202 					shift_right(temp, Delta, e);
1203 				remainder(m4, temp, 4);
1204 				if (m4 < 0)
1205 					m4 += 4;
1206 				if (m4 == 1)
1207 					e >>= 1;
1208 				else
1209 					e = (e-2) >> 1;
1210 			}
1211 			else
1212 				e >>= 1;
1213 
1214 			if (e > 0) {
1215 				power(temp, b, e);
1216 				multiply(cond, cond, temp);
1217 			}
1218 		}
1219 	}
1220 
1221 	return cond;
1222 }
1223 
1224 
1225 
1226 //
1227 // quadratic_order::is_maximal()
1228 //
1229 // Task:
1230 //      returns true if the quadratic order is maximal (conductor = 1).
1231 //
1232 
1233 bool
is_maximal()1234 quadratic_order::is_maximal()
1235 {
1236 	debug_handler("quadratic_order", "is_maximal");
1237 
1238 	return (conductor() == bigint(1));
1239 }
1240 
1241 
1242 
1243 //
1244 // quadratic_order::maximize()
1245 //
1246 // Task:
1247 //      maximizes the quadratic order
1248 //
1249 
1250 void
maximize()1251 quadratic_order::maximize()
1252 {
1253 	debug_handler("quadratic_order", "maximize()");
1254 
1255 	bigint cond;
1256 
1257 	// set last-used quadratic_order
1258 	qo_l().set_last(this);
1259 
1260 	square(cond, conductor());
1261 	assign(Delta / cond);
1262 }
1263 
1264 
1265 
1266 //
1267 // quadratic_order::maximize(quadratic_order)
1268 //
1269 // Task:
1270 //      returns the maximal quadratic order of which this quadratic order is
1271 //      a subset.
1272 //
1273 
1274 void
maximize(quadratic_order & max_ord)1275 quadratic_order::maximize(quadratic_order & max_ord)
1276 {
1277 	debug_handler("quadratic_order", "maximize()");
1278 
1279 	bigint cond;
1280 
1281 	// set last-used quadratic_order
1282 	qo_l().set_last(this);
1283 
1284 	square(cond, conductor());
1285 	max_ord.assign(Delta / cond);
1286 }
1287 
1288 
1289 
1290 //
1291 // kronecker
1292 //
1293 // Task:
1294 //      computes the kronecker symbol (D/p)
1295 //
1296 
1297 int
kronecker(const bigint & D,const bigint & p)1298 kronecker(const bigint & D, const bigint & p)
1299 {
1300 	debug_handler("quadratic_order", "kronecker");
1301 
1302 	bigint c1;
1303 	int kro;
1304 
1305 	// if p > 2, use jacobi symbol
1306 	if (p.compare(bigint(2))) {
1307 		remainder(c1, D, p);
1308 		if (c1.is_lt_zero())
1309 			add(c1, c1, p);
1310 		kro = jacobi(c1, p);
1311 	}
1312 	else {
1313 		remainder(c1, D, 8);
1314 		if (c1.is_lt_zero())
1315 			add(c1, c1, 8);
1316 		if (c1.is_one())
1317 			kro = 1;
1318 		else if (!c1.compare(bigint(5)))
1319 			kro = -1;
1320 		else
1321 			kro = 0;
1322 	}
1323 
1324 	return kro;
1325 }
1326 
1327 
1328 
1329 //
1330 // quadratic_order::generate_optimal_Q()
1331 //
1332 // Task:
1333 //      computes the optimal value of Q for computing h* such that
1334 //      h* < h < 2h*
1335 //
1336 
1337 long
generate_optimal_Q()1338 quadratic_order::generate_optimal_Q()
1339 {
1340 	debug_handler("quadratic_order", "generate_optimal_Q");
1341 
1342 	long OQ;
1343 	bigfloat A, l2;
1344 	long oprec = bigfloat::get_precision();
1345 
1346 	bigfloat::set_precision(prec);
1347 
1348 	// set last-used quadratic_order
1349 	qo_l().set_last(this);
1350 
1351 	if ((PL.get_lower_bound() != 2) || (PL.get_upper_bound() < 1000000))
1352 		PL.resize(2, 1000000);
1353 
1354 	l2 = log(sqrt(bigfloat(2.0)));
1355 	OQ = PL.get_first_prime(3);
1356 	A = estimate_L1_error(OQ);
1357 	while (A >= l2) {
1358 		OQ = PL.get_next_prime();
1359 		A = estimate_L1_error(OQ);
1360 	}
1361 
1362 	bigfloat::set_precision(oprec);
1363 
1364 	return OQ;
1365 }
1366 
1367 
1368 
1369 //
1370 // quadratic_order::get_optimal_Q()
1371 //
1372 // Task:
1373 //      returns a value of Q from a pre-computed table which will compute h*
1374 //      such that h* < h < 2h*.
1375 //
1376 
1377 long
get_optimal_Q()1378 quadratic_order::get_optimal_Q()
1379 {
1380 	debug_handler("quadratic_order", "get_optimal_Q");
1381 
1382 	long Dlog;
1383 	bigfloat temp;
1384 
1385 	// set last-used quadratic_order
1386 	qo_l().set_last(this);
1387 
1388 	temp = floor(log(bigfloat(abs(Delta))) / std::log(10.0));
1389 	temp.longify(Dlog);
1390 
1391 	return OQvals[Dlog / 5];
1392 }
1393 
1394 
1395 
1396 //
1397 // quadratic_order::generate_optimal_Q_cnum()
1398 //
1399 // Task:
1400 //      computes the optimal value of Q required to determine whether the
1401 //      approximation of h gleaned from the analytic class number formula is
1402 //      actually the class number, provided that h <=3.
1403 //
1404 
1405 long
generate_optimal_Q_cnum(const bigfloat & h2,const bigfloat & t)1406 quadratic_order::generate_optimal_Q_cnum(const bigfloat & h2,
1407                                          const bigfloat & t)
1408 {
1409 	debug_handler("quadratic_order", "generate_optimal_Q_cnum");
1410 
1411 	long OQ;
1412 	bigfloat A, val;
1413 	long oprec = bigfloat::get_precision();
1414 
1415 	bigfloat::set_precision(prec);
1416 
1417 	// set last-used quadratic_order
1418 	qo_l().set_last(this);
1419 
1420 	if ((PL.get_lower_bound() != 2) || (PL.get_upper_bound() < 1000000))
1421 		PL.resize(2, 1000000);
1422 
1423 	val = log((h2 + bigfloat(1.0)) / (h2 + t));
1424 	OQ = PL.get_first_prime(3);
1425 	A = estimate_L1_error(OQ);
1426 	while (A >= val) {
1427 		OQ = PL.get_next_prime();
1428 		A = estimate_L1_error(OQ);
1429 	}
1430 
1431 	bigfloat::set_precision(oprec);
1432 
1433 	return OQ;
1434 }
1435 
1436 
1437 
1438 //
1439 // quadratic_order::get_optimal_Q_cnum()
1440 //
1441 // Task:
1442 //      returns a value of Q from a pre-computed table which is sufficient
1443 //      to determine whether the approximation of h gleaned from the analytic
1444 //      class number formula is actually the class number, provided that h <=3.
1445 //
1446 
1447 long
get_optimal_Q_cnum()1448 quadratic_order::get_optimal_Q_cnum()
1449 {
1450 	debug_handler("quadratic_order", "get_optimal_Q_cnum");
1451 
1452 	long Dlog;
1453 	bigfloat temp;
1454 
1455 	// set last-used quadratic_order
1456 	qo_l().set_last(this);
1457 
1458 	temp = floor(log(bigfloat(abs(Delta))) / std::log(10.0));
1459 	temp.longify(Dlog);
1460 
1461 	return OQvals_cnum[Dlog / 5];
1462 }
1463 
1464 
1465 
1466 //
1467 // quadratic_order::generate_optimal_Q_cfunc()
1468 //
1469 // Task:
1470 //      computes the optimal value of Q required to compute C(Delta) to 8
1471 //      significant digits.
1472 //
1473 
1474 long
generate_optimal_Q_cfunc()1475 quadratic_order::generate_optimal_Q_cfunc()
1476 {
1477 	debug_handler("quadratic_order", "generate_optimal_Q_cfunc");
1478 
1479 	long OQ;
1480 	bigfloat B, sb, val, lDelta, lQ, temp1, temp2;
1481 	long oprec = bigfloat::get_precision();
1482 
1483 	bigfloat::set_precision(prec);
1484 
1485 	// set last-used quadratic_order
1486 	qo_l().set_last(this);
1487 
1488 	sqrt(val, bigfloat(1.0000002));
1489 	val.divide_by_2();
1490 	inc(val);
1491 	val.assign(log(val));
1492 
1493 	lDelta.assign(log(bigfloat(abs(Delta))));
1494 
1495 	OQ = 100;
1496 	lQ.assign(log(bigfloat(OQ)));
1497 
1498 	multiply(temp1, Pi(), lQ);
1499 	temp1.invert();
1500 	square(temp2, lQ);
1501 	divide(temp2, bigfloat(5.3), temp2);
1502 	add(temp1, temp1, temp2);
1503 	multiply(B, lDelta, temp1);
1504 
1505 	divide(temp1, bigfloat(4.0), lQ);
1506 	divide(temp2, bigfloat(1.0), Pi());
1507 	add(B, B, temp1);
1508 	add(B, B, temp2);
1509 
1510 	power(temp1, bigfloat(OQ), bigfloat(1.5));
1511 	multiply(temp1, temp1, 9);
1512 	multiply(temp2, lQ, 13);
1513 	add(temp2, temp2, 8);
1514 	divide(sb, temp2, temp1);
1515 	multiply(sb, sb, B);
1516 
1517 	square(temp1, bigfloat(OQ));
1518 	divide(temp1, bigfloat(2.0), temp1);
1519 	add(sb, sb, temp1);
1520 
1521 	while (sb >= val) {
1522 		if (OQ == 100)
1523 			OQ = 1000;
1524 		else if (OQ < 5000)
1525 			OQ += 1000;
1526 		else
1527 			OQ += 5000;
1528 
1529 		lQ.assign(log(bigfloat(OQ)));
1530 
1531 		multiply(temp1, Pi(), lQ);
1532 		temp1.invert();
1533 		square(temp2, lQ);
1534 		divide(temp2, bigfloat(5.3), temp2);
1535 		add(temp1, temp1, temp2);
1536 		multiply(B, lDelta, temp1);
1537 
1538 		divide(temp1, bigfloat(4.0), lQ);
1539 		divide(temp2, bigfloat(1.0), Pi());
1540 		add(B, B, temp1);
1541 		add(B, B, temp2);
1542 
1543 		power(temp1, bigfloat(OQ), bigfloat(1.5));
1544 		multiply(temp1, temp1, 9);
1545 		multiply(temp2, lQ, 13);
1546 		add(temp2, temp2, 8);
1547 		divide(sb, temp2, temp1);
1548 		multiply(sb, sb, B);
1549 
1550 		square(temp1, bigfloat(OQ));
1551 		divide(temp1, bigfloat(2.0), temp1);
1552 		add(sb, sb, temp1);
1553 	}
1554 
1555 	bigfloat::set_precision(oprec);
1556 
1557 	return OQ;
1558 }
1559 
1560 
1561 
1562 //
1563 // quadratic_order::get_optimal_Q_cfunc()
1564 //
1565 // Task:
1566 //      returns a value of Q from a pre-computed list which is sufficient to
1567 //      compute an approximation of C(Delta) accurate to 8 significant digits.
1568 //
1569 
1570 long
get_optimal_Q_cfunc()1571 quadratic_order::get_optimal_Q_cfunc()
1572 {
1573 	debug_handler("quadratic_order", "get_optimal_Q_cfunc");
1574 
1575 	long Dlog;
1576 	bigfloat temp;
1577 
1578 	// set last-used quadratic_order
1579 	qo_l().set_last(this);
1580 
1581 	temp = floor(log(bigfloat(abs(Delta))) / std::log(10.0));
1582 	temp.longify(Dlog);
1583 
1584 	return OQvals_cfunc[Dlog / 5];
1585 }
1586 
1587 
1588 
1589 //
1590 // quadratic_order::estimate_C
1591 //
1592 // Task:
1593 //      computes a truncated product approximation of C(D) using primes
1594 //      less than Q.
1595 //
1596 
1597 bigfloat
estimate_C(const long nQ)1598 quadratic_order::estimate_C(const long nQ)
1599 {
1600 	debug_handler("quadratic_order", "estimate_C");
1601 
1602 	int kron;
1603 	long P, QQ;
1604 	double E;
1605 
1606 	// set last-used quadratic_order
1607 	qo_l().set_last(this);
1608 
1609 	E = 1.0;
1610 
1611 	// generate sufficiently many primes
1612 	QQ = 1000000;
1613 	if (QQ > nQ)
1614 		QQ = nQ+10;
1615 	P = 3;
1616 
1617 	if ((QQ > static_cast<long>(PL.get_upper_bound())) || (PL.get_lower_bound() != 2))
1618 		PL.resize(2, QQ);
1619 	P = PL.get_first_prime(P);
1620 
1621 	// compute partial product for LQ <= p < QQ
1622 	while ((P< nQ) && (P > 0)) {
1623 		kron = kronecker(Delta, bigint(P));
1624 		E *= 1.0L - (static_cast<double>(kron) / (P - 1.0));
1625 
1626 		P = PL.get_next_prime();
1627 		if (P == 0) {
1628 			QQ = PL.get_upper_bound() + 1000000;
1629 			if (QQ > nQ)
1630 				QQ = nQ+10;
1631 			PL.resize(PL.get_upper_bound() + 1, QQ);
1632 			P = PL.get_first_prime();
1633 		}
1634 	}
1635 
1636 	return bigfloat(E);
1637 }
1638 
1639 
1640 
1641 //
1642 // quadratic_order::estimate_L
1643 //
1644 // Task:
1645 //      computes a truncated product approximation of L(s,X) using primes
1646 //      less than Q.
1647 //
1648 
1649 bigfloat
estimate_L(const int s,const long nQ)1650 quadratic_order::estimate_L(const int s, const long nQ)
1651 {
1652 	debug_handler("quadratic_order", "estimate_L");
1653 
1654 	int kron, i;
1655 	long m8, P, s2, QQ;
1656 	double E, Ps;
1657 	bigfloat estL;
1658 
1659 	// set last-used quadratic_order
1660 	qo_l().set_last(this);
1661 
1662 	// compute partial product for p = 2
1663 	m8 = remainder(Delta, 8);
1664 	if (m8 < 0)
1665 		m8 += 8;
1666 	if (m8 == 1)
1667 		kron = 1;
1668 	else if (m8 == 5)
1669 		kron = -1;
1670 	else
1671 		kron = 0;
1672 	s2 = 1;
1673 	s2 <<= s;
1674 	E = std::log(static_cast<double>(s2) / static_cast<double>(s2-kron));
1675 
1676 	// generate sufficiently many primes
1677 	QQ = 1000000;
1678 	if (QQ > nQ)
1679 		QQ = nQ+10;
1680 	P = 3;
1681 
1682 	if ((QQ > static_cast<long>(PL.get_upper_bound())) || (PL.get_lower_bound() != 2))
1683 		PL.resize(2, QQ);
1684 	P = PL.get_first_prime(P);
1685 
1686 	// compute partial product for LQ <= p < QQ
1687 	while ((P< nQ) && (P > 0)) {
1688 		kron = kronecker(Delta, bigint(P));
1689 		Ps = P;
1690 		for (i = 1; i < s; ++i)
1691 			Ps *= P;
1692 		E += std::log(Ps / (Ps - kron));
1693 
1694 		P = PL.get_next_prime();
1695 		if (P == 0) {
1696 			QQ = PL.get_upper_bound() + 1000000;
1697 			if (QQ > nQ)
1698 				QQ = nQ+10;
1699 			PL.resize(PL.get_upper_bound() + 1, QQ);
1700 			P = PL.get_first_prime();
1701 		}
1702 	}
1703 
1704 	return std::exp(E);
1705 }
1706 
1707 
1708 
1709 //
1710 // quadratic_order::estimate_L1
1711 //
1712 // Task:
1713 //      computes a truncated product approximation of L(1,X) using primes
1714 //      less than 2Q.  This function uses Bach's method of a weighted average
1715 //      of truncated Euler products.
1716 //
1717 
1718 bigfloat
estimate_L1(const long nQ)1719 quadratic_order::estimate_L1(const long nQ)
1720 {
1721 	debug_handler("quadratic_order", "estimate_L1");
1722 
1723 	int kron;
1724 	long Q2, m8, P, QQ;
1725 	double E, C, wt, s;
1726 	register long i;
1727 	bigfloat nFI;
1728 
1729 	// set last-used quadratic_order
1730 	qo_l().set_last(this);
1731 
1732 	// if nQ = Q, then the estimate has already been computed
1733 	if ((nQ == Q) && (!FI.is_zero())) {
1734 		nFI = FI;
1735 		return(nFI);
1736 	}
1737 
1738 	// compute partial product for p = 2
1739 	Q = nQ;
1740 	m8 = remainder(Delta, 8);
1741 	if (m8 < 0)
1742 		m8 += 8;
1743 	if (m8 == 1)
1744 		E = std::log(2.0);
1745 	else if (m8 == 5)
1746 		E = std::log(2.0 / 3.0);
1747 	else
1748 		E = 0.0;
1749 
1750 	// compute weights
1751 	Q2 = Q << 1;
1752 	C = 0.0;
1753 	for (i = Q; i <= Q2-1; ++i)
1754 		C += static_cast<double>(i) * std::log(static_cast<double>(i));
1755 
1756 	// generate list of primes
1757 	QQ = 1000000;
1758 	if (QQ > Q2)
1759 		QQ = Q2+10;
1760 	P = 3;
1761 
1762 	if ((QQ > static_cast<long>(PL.get_upper_bound())) || (PL.get_lower_bound() != 2))
1763 		PL.resize(2, QQ);
1764 	P = PL.get_first_prime(P);
1765 
1766 	// compute partial product for p < Q
1767 	while ((P< Q) && (P > 0)) {
1768 		kron = kronecker(Delta, bigint(P));
1769 		E += std::log(static_cast<double>(P) / (P - kron));
1770 
1771 		P = PL.get_next_prime();
1772 		if (P == 0) {
1773 			QQ = PL.get_upper_bound() + 1000000;
1774 			if (QQ > Q2)
1775 				QQ = Q2+10;
1776 			PL.resize(PL.get_upper_bound() + 1, QQ);
1777 			P = PL.get_first_prime();
1778 		}
1779 	}
1780 
1781 	// computed weighted partial products for Q < p < 2Q
1782 	wt = 1.0;
1783 	s = 0.0;
1784 	for (i = Q; i <= P; ++i)
1785 		s += static_cast<double>(i) * std::log(static_cast<double>(i));
1786 	wt -= s / C;
1787 
1788 	while ((P< Q2) && (P > 0)) {
1789 		kron = kronecker(Delta, bigint(P));
1790 		E += wt*std::log(static_cast<double>(P) / (P - kron));
1791 
1792 		i = P + 1;
1793 		P = PL.get_next_prime();
1794 		if (P == 0) {
1795 			QQ = PL.get_upper_bound() + 1000000;
1796 			if (QQ > Q2)
1797 				QQ = Q2+10;
1798 			PL.resize(PL.get_upper_bound() + 1, QQ);
1799 			P = PL.get_first_prime();
1800 		}
1801 		s = 0.0;
1802 		for (; i <= P; ++i)
1803 			s += static_cast<double>(i) * std::log(static_cast<double>(i));
1804 		wt -= s / C;
1805 	}
1806 
1807 	nFI = exp(bigfloat(E));
1808 
1809 	// set static value
1810 	FI = nFI;
1811 
1812 	return(nFI);
1813 }
1814 
1815 
1816 
1817 //
1818 // quadratic_order::estimate_L1_error
1819 //
1820 // Task:
1821 //      computes an upper bound of the error in the approximation of L(1,X)
1822 //      computed by the function estimate_L1
1823 //
1824 
1825 bigfloat
estimate_L1_error(const long nQ) const1826 quadratic_order::estimate_L1_error(const long nQ) const
1827 {
1828 	debug_handler("quadratic_order", "estimate_L1_error");
1829 
1830 	bigfloat A, lq;
1831 
1832 	lq.assign(log(bigfloat(nQ)));
1833 	if (Q >= 100000)
1834 		A = (6.338*log(bigfloat(abs(Delta))) + 17.031) / (lq * sqrt(bigfloat(nQ)));
1835 	else if (Q >= 50000)
1836 		A = (6.378*log(bigfloat(abs(Delta))) + 17.397) / (lq * sqrt(bigfloat(nQ)));
1837 	else if (Q >= 10000)
1838 		A = (6.510*log(bigfloat(abs(Delta))) + 18.606) / (lq * sqrt(bigfloat(nQ)));
1839 	else if (Q >= 5000)
1840 		A = (6.593*log(bigfloat(abs(Delta))) + 19.321) / (lq * sqrt(bigfloat(nQ)));
1841 	else if (Q >= 1000)
1842 		A = (6.897*log(bigfloat(abs(Delta))) + 21.528) / (lq * sqrt(bigfloat(nQ)));
1843 	else
1844 		A = (7.106*log(bigfloat(abs(Delta))) + 22.845) / (lq * sqrt(bigfloat(nQ)));
1845 
1846 	return A;
1847 }
1848 
1849 
1850 
1851 //
1852 // quadratic_order::Lfunction
1853 //
1854 // Task:
1855 //      returns an approximation of L(1,X) computed via the analytic class
1856 //      number formula.
1857 //
1858 
1859 bigfloat
Lfunction()1860 quadratic_order::Lfunction()
1861 {
1862 	debug_handler("quadratic_order", "Lfunction");
1863 
1864 	long oprec = bigfloat::get_precision();
1865 	bigfloat outL;
1866 
1867 	bigfloat::set_precision(prec);
1868 
1869 	// set last-used quadratic_order
1870 	qo_l().set_last(this);
1871 
1872 	if (!is_L_computed()) {
1873 		// L has not been computed - compute it
1874 
1875 		// compute class number, if needed
1876 		if (h.is_zero())
1877 			class_number();
1878 
1879 		if (is_imaginary()) {
1880 			// L = pi*h / sqrt(Delta)
1881 			multiply(L, Pi(), bigfloat(h));
1882 			divide(L, L, sqrt(bigfloat(-Delta)));
1883 			if (Delta == -4)  L.divide_by_2();
1884 			if (Delta == -3)  divide(L, L, bigfloat(3.0));
1885 		}
1886 		else {
1887 			// L = 2hR / sqrt(Delta)
1888 			multiply(L, R, bigfloat(h));
1889 			L.multiply_by_2();
1890 			divide(L, L, sqrt(bigfloat(Delta)));
1891 		}
1892 	}
1893 
1894 	outL.assign(L);
1895 	bigfloat::set_precision(oprec);
1896 
1897 	return outL;
1898 }
1899 
1900 
1901 
1902 //
1903 // quadratic_order::LDfunction
1904 //
1905 // Task:
1906 //      returns an approximation of the LD(1) function defined by Shanks
1907 //      (L(1,X) without the 2-factor) computed via the analytic class
1908 //      number formula.
1909 //
1910 
1911 bigfloat
LDfunction()1912 quadratic_order::LDfunction()
1913 {
1914 	debug_handler("quadratic_order", "LDfunction");
1915 
1916 	bigfloat LD;
1917 	long m8;
1918 	long oprec = bigfloat::get_precision();
1919 
1920 	bigfloat::set_precision(prec);
1921 
1922 	// set last-used quadratic_order
1923 	qo_l().set_last(this);
1924 
1925 	m8 = remainder(Delta, 8);
1926 	if (m8 < 0)
1927 		m8 += 8;
1928 	LD.assign(Lfunction());
1929 	if (m8 == 1)
1930 		LD.divide_by_2();
1931 	else if (m8 == 5)
1932 		multiply(LD, LD, bigfloat(1.5));
1933 
1934 	bigfloat::set_precision(oprec);
1935 
1936 	return LD;
1937 }
1938 
1939 
1940 
1941 //
1942 // quadratic_order::LLI
1943 //
1944 // Task:
1945 //      returns an approximation of the Lower Littlewood Index (defined by
1946 //      Shanks).
1947 //
1948 
1949 bigfloat
LLI()1950 quadratic_order::LLI()
1951 {
1952 	debug_handler("quadratic_order", "LLI");
1953 
1954 	bigfloat LLI, c;
1955 	long oprec = bigfloat::get_precision();
1956 
1957 	bigfloat::set_precision(prec);
1958 
1959 	// set last-used quadratic_order
1960 	qo_l().set_last(this);
1961 
1962 	if (Delta.is_zero())
1963 		LLI.assign_zero();
1964 	else {
1965 		c.assign(exp(Euler()) / (Pi() * Pi()));
1966 
1967 		if (Delta.is_odd())
1968 			c *= bigfloat(12.0);
1969 		else
1970 			c *= bigfloat(8.0);
1971 
1972 		LLI = Lfunction() * c * log(log(bigfloat(abs(Delta))));
1973 	}
1974 
1975 	bigfloat::set_precision(oprec);
1976 
1977 	return LLI;
1978 }
1979 
1980 
1981 
1982 //
1983 // quadratic_order::LLI_D
1984 //
1985 // Task:
1986 //      returns an approximation of the LLI_D value (defined by Shanks).
1987 //
1988 
1989 bigfloat
LLI_D()1990 quadratic_order::LLI_D()
1991 {
1992 	debug_handler("quadratic_order", "LLI_D");
1993 
1994 	bigfloat LLI, c;
1995 	long oprec = bigfloat::get_precision();
1996 
1997 	bigfloat::set_precision(prec);
1998 
1999 	// set last-used quadratic_order
2000 	qo_l().set_last(this);
2001 
2002 	if (Delta.is_zero())
2003 		LLI.assign_zero();
2004 	else {
2005 		c.assign(bigfloat(8.0) * exp(Euler()) / (Pi() * Pi()));
2006 		LLI = LDfunction() * c * log(log(bigfloat(abs(Delta << 2))));
2007 	}
2008 
2009 	bigfloat::set_precision(oprec);
2010 
2011 	return LLI;
2012 }
2013 
2014 
2015 
2016 //
2017 // quadratic_order::ULI
2018 //
2019 // Task:
2020 //      returns an approximation of the Upper Littlewood Index (defined by
2021 //      Shanks).
2022 //
2023 
2024 bigfloat
ULI()2025 quadratic_order::ULI()
2026 {
2027 	debug_handler("quadratic_order", "ULI");
2028 
2029 	bigfloat ULI, c;
2030 	long oprec = bigfloat::get_precision();
2031 
2032 	bigfloat::set_precision(prec);
2033 
2034 	// set last-used quadratic_order
2035 	qo_l().set_last(this);
2036 
2037 	if (Delta.is_zero())
2038 		ULI.assign_zero();
2039 	else {
2040 		c.assign(exp(Euler()));
2041 
2042 		if (Delta.is_odd())
2043 			c.multiply_by_2();
2044 
2045 		ULI = Lfunction() / (c * log(log(bigfloat(abs(Delta)))));
2046 	}
2047 
2048 	bigfloat::set_precision(oprec);
2049 
2050 	return ULI;
2051 }
2052 
2053 
2054 
2055 //
2056 // quadratic_order::ULI_D
2057 //
2058 // Task:
2059 //      returns an approximation of the ULI_D value (defined by Shanks).
2060 //
2061 
2062 bigfloat
ULI_D()2063 quadratic_order::ULI_D()
2064 {
2065 	debug_handler("quadratic_order", "ULI_D");
2066 
2067 	bigfloat ULI, c;
2068 	long oprec = bigfloat::get_precision();
2069 
2070 	bigfloat::set_precision(prec);
2071 
2072 	// set last-used quadratic_order
2073 	qo_l().set_last(this);
2074 
2075 	if (Delta.is_zero())
2076 		ULI.assign_zero();
2077 	else {
2078 		c.assign(exp(Euler()));
2079 		ULI = LDfunction() / (c * log(log(bigfloat(abs(Delta << 2)))));
2080 	}
2081 
2082 	bigfloat::set_precision(oprec);
2083 
2084 	return ULI;
2085 }
2086 
2087 
2088 
2089 //
2090 // quadratic_order::Cfunction
2091 //
2092 // Task:
2093 //      returns an approximation of C(Delta) accurate to 8 significant digits.
2094 //
2095 
2096 bigfloat
Cfunction()2097 quadratic_order::Cfunction()
2098 {
2099 	debug_handler("quadratic_order", "Cfunction");
2100 
2101 	bigfloat temp, E4, outC;
2102 	bigint PP;
2103 	int kron;
2104 	long Q2, m8, P, num, QQ;
2105 	double E, E2, E3, C, wt, dP, dP2, s;
2106 	register long i;
2107 	long oprec = bigfloat::get_precision();
2108 
2109 	bigfloat::set_precision(prec);
2110 
2111 	// set last-used quadratic_order
2112 	qo_l().set_last(this);
2113 
2114 	if (!is_C_computed()) {
2115 		Q2 = get_optimal_Q_cfunc();
2116 
2117 		if (h.is_zero()) {
2118 			Q = Q2 >> 1;
2119 
2120 			// compute partial product for p = 2
2121 			m8 = remainder(Delta, 8);
2122 			if (m8 < 0)
2123 				m8 += 8;
2124 			if (m8 == 1) {
2125 				E = std::log(2.0);
2126 				E2 = std::log(4.0 / 3.0);
2127 				Cfunc.assign(2.5);
2128 			}
2129 			else if (m8 == 5) {
2130 				E = std::log(2.0 / 3.0);
2131 				E2 = std::log(4.0 / 5.0);
2132 				Cfunc.assign(0.5);
2133 			}
2134 			else {
2135 				E = 0.0;
2136 				E2 = 0.0;
2137 				Cfunc.assign(15.0/16.0);
2138 			}
2139 			E3 = 0.0;
2140 
2141 			// compute weights
2142 			C = 0.0;
2143 			for (i = Q; i <= Q2-1; ++i)
2144 				C += static_cast<double>(i) * std::log(static_cast<double>(i));
2145 
2146 			// generate list of primes
2147 			QQ = 1000000;
2148 			if (QQ > Q2)
2149 				QQ = Q2+10;
2150 			P = 3;
2151 
2152 			if ((QQ > static_cast<long>(PL.get_upper_bound())) || (PL.get_lower_bound() != 2))
2153 				PL.resize(2, QQ);
2154 			P = PL.get_first_prime(P);
2155 
2156 			// compute partial product for p < Q
2157 			while ((P< Q) && (P > 0)) {
2158 				kron = kronecker(Delta, bigint(P));
2159 				E += std::log(static_cast<double>(P) / (P - kron));
2160 
2161 				dP = static_cast<double>(P);
2162 				dP2 = dP*dP;
2163 				E2 += std::log(dP2 / (dP2 - kron));
2164 
2165 				if (kron == 1) {
2166 					dP2 = (dP - 1.0) * (dP - 1.0);
2167 					dP2 *= dP;
2168 					E3 += std::log(1.0 - (2.0 / dP2));
2169 				}
2170 
2171 				P = PL.get_next_prime();
2172 				if (P == 0) {
2173 					QQ = PL.get_upper_bound() + 1000000;
2174 					if (QQ > Q2)
2175 						QQ = Q2+10;
2176 					PL.resize(PL.get_upper_bound() + 1, QQ);
2177 					P = PL.get_first_prime();
2178 				}
2179 			}
2180 
2181 			// computed weighted partial products for Q < p < 2Q
2182 			wt = 1.0;
2183 			s = 0.0;
2184 			for (i = Q; i <= P; ++i)
2185 				s += static_cast<double>(i) * std::log(static_cast<double>(i));
2186 			wt -= s / C;
2187 			while ((P< Q2) && (P > 0)) {
2188 				kron = kronecker(Delta, bigint(P));
2189 				E += wt * std::log(static_cast<double>(P) / (P - kron));
2190 
2191 				dP = static_cast<double>(P);
2192 				dP2 = dP*dP;
2193 				E2 += std::log(dP2 / (dP2 - kron));
2194 
2195 				if (kron == 1) {
2196 					dP2 = (dP - 1.0) * (dP - 1.0);
2197 					dP2 *= dP;
2198 					E3 += std::log(1.0 - (2.0 / dP2));
2199 				}
2200 
2201 				i = P + 1;
2202 				P = PL.get_next_prime();
2203 				if (P == 0) {
2204 					QQ = PL.get_upper_bound() + 1000000;
2205 					if (QQ > Q2)
2206 						QQ = Q2+10;
2207 					PL.resize(PL.get_upper_bound() + 1, QQ);
2208 					P = PL.get_first_prime();
2209 				}
2210 				s = 0.0;
2211 				for (; i <= P; ++i)
2212 					s += static_cast<double>(i) * std::log(static_cast<double>(i));
2213 				wt -= s / C;
2214 			}
2215 
2216 			FI = exp(bigfloat(E));
2217 			class_number();
2218 		}
2219 		else {
2220 			// compute partial product for p = 2
2221 			m8 = remainder(Delta, 8);
2222 			if (m8 < 0)
2223 				m8 += 8;
2224 			if (m8 == 1) {
2225 				E2 = std::log(4.0 / 3.0);
2226 				Cfunc.assign(2.5);
2227 			}
2228 			else if (m8 == 5) {
2229 				E2 = std::log(4.0 / 5.0);
2230 				Cfunc.assign(0.5);
2231 			}
2232 			else {
2233 				E2 = 0.0;
2234 				Cfunc.assign(15.0/16.0);
2235 			}
2236 			E3 = 0.0;
2237 
2238 			// generate list of primes
2239 			QQ = 1000000;
2240 			if (QQ > Q2)
2241 				QQ = Q2+10;
2242 			P = 3;
2243 
2244 			if ((QQ > static_cast<long>(PL.get_upper_bound())) || (PL.get_lower_bound() != 2))
2245 				PL.resize(2, QQ);
2246 			P = PL.get_first_prime(P);
2247 
2248 			// compute partial product for p < Q
2249 			while ((P< Q2) && (P > 0)) {
2250 				kron = kronecker(Delta, bigint(P));
2251 
2252 				dP = static_cast<double>(P);
2253 				dP2 = dP*dP;
2254 				E2 += std::log(dP2 / (dP2 - kron));
2255 
2256 				if (kron == 1) {
2257 					dP2 = (dP - 1.0) * (dP - 1.0);
2258 					dP2 *= dP;
2259 					E3 += std::log(1.0 - (2.0 / dP2));
2260 				}
2261 
2262 				P = PL.get_next_prime();
2263 				if (P == 0) {
2264 					QQ = PL.get_upper_bound() + 1000000;
2265 					if (QQ > Q2)
2266 						QQ = Q2+10;
2267 					PL.resize(PL.get_upper_bound() + 1, QQ);
2268 					P = PL.get_first_prime();
2269 				}
2270 			}
2271 		}
2272 
2273 		// compute C
2274 		if (is_imaginary()) {
2275 			multiply(Cfunc, Cfunc, sqrt(bigfloat(-Delta)));
2276 			power(temp, Pi(), 3);
2277 			divide(temp, temp, bigfloat(h));
2278 			divide(temp, temp, 90);
2279 			multiply(Cfunc, Cfunc, temp);
2280 
2281 			// multiply by the number of roots of unity
2282 			if (Delta == -4)  Cfunc.multiply_by_2();
2283 			if (Delta == -3)  multiply(Cfunc, Cfunc, bigfloat(3.0));
2284 		}
2285 		else {
2286 			multiply(Cfunc, Cfunc, sqrt(bigfloat(Delta)));
2287 			power(temp, Pi(), 4);
2288 			divide(temp, temp, bigfloat(h));
2289 			divide(temp, temp, R);
2290 			divide(temp, temp, 180);
2291 			multiply(Cfunc, Cfunc, temp);
2292 		}
2293 
2294 		factor_discriminant();
2295 		E4.assign_one();
2296 		num = disc_fact.no_of_comp();
2297 		for (i = 0; i < num; ++i) {
2298 			PP = disc_fact.base(static_cast<lidia_size_t>(i));
2299 			if (PP > 2) {
2300 				power(temp, bigfloat(PP), 4);
2301 				divide(temp, temp-bigfloat(1.0), temp);
2302 				E4 *= temp;
2303 			}
2304 		}
2305 
2306 		multiply(Cfunc, Cfunc, E4);
2307 		divide(Cfunc, Cfunc, exp(bigfloat(E2)));
2308 		multiply(Cfunc, Cfunc, exp(bigfloat(E3)));
2309 	}
2310 
2311 	outC.assign(Cfunc);
2312 	bigfloat::set_precision(oprec);
2313 
2314 	return outC;
2315 }
2316 
2317 
2318 
2319 //
2320 // quadratic_order::Cfunction_bigfloat
2321 //
2322 // Task:
2323 //      returns an approximation of C(Delta) accurate to 8 significant digits
2324 //      using bigfloats
2325 //
2326 
2327 bigfloat_int
Cfunction_bigfloat(int pmult)2328 quadratic_order::Cfunction_bigfloat(int pmult)
2329 {
2330 	debug_handler("quadratic_order", "Cfunction_bigfloat");
2331 
2332 	bigfloat_int temp, E4, outC, Cfunc_int;
2333 	bigint PP;
2334 	int kron;
2335 	long Q2, m8, P, num, QQ;
2336 	bigfloat_int E, E2, E3, C, wt, dP, dP2, s;
2337 	register long i;
2338 	long oprec = bigfloat::get_precision();
2339 
2340 	bigfloat::set_precision(prec*pmult);
2341 
2342 	// set last-used quadratic_order
2343 	qo_l().set_last(this);
2344 
2345 	if (!is_C_computed()) {
2346 		Q2 = get_optimal_Q_cfunc();
2347 
2348 		if (h.is_zero()) {
2349 			Q = Q2 >> 1;
2350 
2351 			// compute partial product for p = 2
2352 			m8 = remainder(Delta, 8);
2353 			if (m8 < 0)
2354 				m8 += 8;
2355 			if (m8 == 1) {
2356 				E = log(bigfloat_int(2.0));
2357 				E2.assign(bigint(4), bigint(3));
2358 				Cfunc_int.assign(2.5);
2359 			}
2360 			else if (m8 == 5) {
2361 				E = log(bigfloat_int(2.0) / bigfloat_int(3.0));
2362 				E2.assign(0.8);
2363 				Cfunc_int.assign(0.5);
2364 			}
2365 			else {
2366 				E.assign_zero();
2367 				E2.assign(bigfloat_int(1));
2368 				Cfunc_int.assign(0.9375);
2369 			}
2370 			E3.assign(bigfloat_int(1));
2371 
2372 			// compute weights
2373 			C.assign_zero();
2374 			for (i = Q; i <= Q2-1; ++i)
2375 				C += bigfloat_int(i) * log(bigfloat_int(i));
2376 
2377 			// generate list of primes
2378 			QQ = 1000000;
2379 			if (QQ > Q2)
2380 				QQ = Q2+10;
2381 			P = 3;
2382 
2383 			if ((QQ > static_cast<long>(PL.get_upper_bound())) || (PL.get_lower_bound() != 2))
2384 				PL.resize(2, QQ);
2385 			P = PL.get_first_prime(P);
2386 
2387 			// compute partial product for p < Q
2388 			while ((P< Q) && (P > 0)) {
2389 				kron = kronecker(Delta, bigint(P));
2390 				E += log(bigfloat_int(P) / bigfloat_int(P - kron));
2391 
2392 				dP.assign(P);
2393 				dP2 = dP*dP;
2394 				multiply(E2, E2, (dP2 / (dP2 - kron)));
2395 
2396 				if (kron == 1) {
2397 					dP2 = (dP - bigfloat_int(1)) * (dP - bigfloat_int(1));
2398 					dP2 *= dP;
2399 					multiply(E3, E3, (bigfloat_int(1) - (bigfloat_int(2) / dP2)));
2400 				}
2401 
2402 				P = PL.get_next_prime();
2403 				if (P == 0) {
2404 					QQ = PL.get_upper_bound() + 1000000;
2405 					if (QQ > Q2)
2406 						QQ = Q2+10;
2407 					PL.resize(PL.get_upper_bound() + 1, QQ);
2408 					P = PL.get_first_prime();
2409 				}
2410 			}
2411 
2412 			// computed weighted partial products for Q < p < 2Q
2413 			wt.assign_one();
2414 			s.assign_zero();
2415 			for (i = Q; i <= P; ++i)
2416 				s += bigfloat_int(i)*log(bigfloat_int(i));
2417 			wt -= s / C;
2418 			while ((P< Q2) && (P > 0)) {
2419 				kron = kronecker(Delta, bigint(P));
2420 				E += wt*log(bigfloat_int(P) / bigfloat_int(P - kron));
2421 
2422 				dP.assign(P);
2423 				dP2 = dP*dP;
2424 				multiply(E2, E2, (dP2 / (dP2 - kron)));
2425 
2426 				if (kron == 1) {
2427 					dP2 = (dP - bigfloat_int(1)) * (dP - bigfloat_int(1));
2428 					dP2 *= dP;
2429 					multiply(E3, E3, (bigfloat_int(1) - (bigfloat_int(2) / dP2)));
2430 				}
2431 
2432 				i = P + 1;
2433 				P = PL.get_next_prime();
2434 				if (P == 0) {
2435 					QQ = PL.get_upper_bound() + 1000000;
2436 					if (QQ > Q2)
2437 						QQ = Q2+10;
2438 					PL.resize(PL.get_upper_bound() + 1, QQ);
2439 					P = PL.get_first_prime();
2440 				}
2441 				s.assign_zero();
2442 				for (; i <= P; ++i)
2443 					s += bigfloat_int(i)*log(bigfloat_int(i));
2444 				wt -= s / C;
2445 			}
2446 
2447 			exp(bigfloat_int(E)).approx(FI);
2448 			class_number();
2449 		}
2450 		else {
2451 			// compute partial product for p = 2
2452 			m8 = remainder(Delta, 8);
2453 			if (m8 < 0)
2454 				m8 += 8;
2455 			if (m8 == 1) {
2456 				E2.assign(bigint(4), bigint(3));
2457 				Cfunc_int.assign(2.5);
2458 			}
2459 			else if (m8 == 5) {
2460 				E2.assign(bigfloat_int(0.8));
2461 				Cfunc_int.assign(0.5);
2462 			}
2463 			else {
2464 				E2.assign(bigfloat_int(1));
2465 				Cfunc_int.assign(0.9375);
2466 			}
2467 			E3.assign(bigfloat_int(1));
2468 
2469 			// generate list of primes
2470 			QQ = 1000000;
2471 			if (QQ > Q2)
2472 				QQ = Q2+10;
2473 			P = 3;
2474 
2475 			if ((QQ > static_cast<long>(PL.get_upper_bound())) || (PL.get_lower_bound() != 2))
2476 				PL.resize(2, QQ);
2477 			P = PL.get_first_prime(P);
2478 
2479 			// compute partial product for p < Q
2480 			while ((P< Q2) && (P > 0)) {
2481 				kron = kronecker(Delta, bigint(P));
2482 
2483 				dP.assign(P);
2484 				square(dP2, dP);
2485 				multiply(E2, E2, (dP2 / (dP2 - bigfloat_int(kron))));
2486 
2487 				if (kron == 1) {
2488 					dP2 = (dP - bigfloat_int(1)) * (dP - bigfloat_int(1));
2489 					dP2 *= dP;
2490 					multiply(E3, E3, (bigfloat_int(1) - (bigfloat_int(2) / dP2)));
2491 				}
2492 
2493 				P = PL.get_next_prime();
2494 				if (P == 0) {
2495 					QQ = PL.get_upper_bound() + 1000000;
2496 					if (QQ > Q2)
2497 						QQ = Q2+10;
2498 					PL.resize(PL.get_upper_bound() + 1, QQ);
2499 					P = PL.get_first_prime();
2500 				}
2501 			}
2502 		}
2503 
2504 		// compute C
2505 		if (is_imaginary()) {
2506 			multiply(Cfunc_int, Cfunc_int, sqrt(bigfloat_int(-Delta)));
2507 			temp = Pi();
2508 			power(temp, temp, 3);
2509 			divide(temp, temp, bigfloat_int(h));
2510 			divide(temp, temp, 90);
2511 			multiply(Cfunc_int, Cfunc_int, temp);
2512 
2513 			// multiply by the number of roots of unity
2514 			if (Delta == -4)  multiply(Cfunc_int, Cfunc_int, bigfloat_int(2));
2515 			if (Delta == -3)  multiply(Cfunc_int, Cfunc_int, bigfloat_int(3));
2516 		}
2517 		else {
2518 			multiply(Cfunc_int, Cfunc_int, sqrt(bigfloat_int(Delta)));
2519 			temp = Pi();
2520 			power(temp, temp, 4);
2521 			divide(temp, temp, bigfloat_int(h));
2522 			divide(temp, temp, R);
2523 			divide(temp, temp, 180);
2524 			multiply(Cfunc_int, Cfunc_int, temp);
2525 		}
2526 
2527 		factor_discriminant();
2528 		E4.assign(bigfloat_int(1));
2529 		num = disc_fact.no_of_comp();
2530 		for (i = 0; i < num; ++i) {
2531 			PP = disc_fact.base(static_cast<lidia_size_t>(i));
2532 			if (PP > 2) {
2533 				power(PP, PP, 4);
2534 				temp.assign(PP-1, PP);
2535 				E4 *= temp;
2536 			}
2537 		}
2538 		divide(Cfunc_int, Cfunc_int, E2);
2539 		multiply(Cfunc_int, Cfunc_int, E3);
2540 		multiply(Cfunc_int, Cfunc_int, E4);
2541 	}
2542 
2543 	Cfunc_int.approx(Cfunc);
2544 	outC.assign(Cfunc_int);
2545 	bigfloat::set_precision(oprec);
2546 
2547 	return outC;
2548 }
2549 
2550 
2551 
2552 //
2553 // quadratic_order::LLI
2554 //
2555 // Task:
2556 //      returns an approximation of the Lower Littlewood Index (defined by
2557 //      Shanks).
2558 //
2559 
2560 bigfloat
LLI(const bigint & nh)2561 quadratic_order::LLI(const bigint & nh)
2562 {
2563 	debug_handler("quadratic_order", "LLI(bigint)");
2564 
2565 	bigfloat LL;
2566 	bigint oh = h;
2567 
2568 	h = nh;
2569 	LL = LLI();
2570 	h = oh;
2571 
2572 	return LL;
2573 }
2574 
2575 
2576 
2577 //
2578 // quadratic_order::LLI
2579 //
2580 // Task:
2581 //      returns an approximation of the Lower Littlewood Index (defined by
2582 //      Shanks).
2583 //
2584 
2585 bigfloat
LLI(const bigint & nh,const bigfloat & nR)2586 quadratic_order::LLI(const bigint & nh, const bigfloat & nR)
2587 {
2588 	debug_handler("quadratic_order", "LLI(bigint, bigfloat)");
2589 
2590 	bigfloat LL, oR = R;
2591 	bigint oh = h;
2592 
2593 	h = nh;
2594 	R = nR;
2595 	LL = LLI();
2596 	h = oh;
2597 	R = oR;
2598 
2599 	return LL;
2600 }
2601 
2602 
2603 
2604 //
2605 // quadratic_order::ULI
2606 //
2607 // Task:
2608 //      returns an approximation of the Lower Littlewood Index (defined by
2609 //      Shanks).
2610 //
2611 
2612 bigfloat
ULI(const bigint & nh)2613 quadratic_order::ULI(const bigint & nh)
2614 {
2615 	debug_handler("quadratic_order", "ULI(bigint)");
2616 
2617 	bigfloat UL;
2618 	bigint oh = h;
2619 
2620 	h = nh;
2621 	UL = ULI();
2622 	h = oh;
2623 
2624 	return UL;
2625 }
2626 
2627 
2628 
2629 //
2630 // quadratic_order::ULI
2631 //
2632 // Task:
2633 //      returns an approximation of the Lower Littlewood Index (defined by
2634 //      Shanks).
2635 //
2636 
2637 bigfloat
ULI(const bigint & nh,const bigfloat & nR)2638 quadratic_order::ULI(const bigint & nh, const bigfloat & nR)
2639 {
2640 	debug_handler("quadratic_order", "ULI(bigint, bigfloat)");
2641 
2642 	bigfloat UL, oR = R;
2643 	bigint oh = h;
2644 
2645 	h = nh;
2646 	R = nR;
2647 	UL = ULI();
2648 	h = oh;
2649 	R = oR;
2650 
2651 	return UL;
2652 }
2653 
2654 
2655 
2656 //
2657 // quadratic_order::Cfunction
2658 //
2659 // Task:
2660 //      returns an approximation of C(Delta) accurate to 8 significant digits.
2661 //
2662 
2663 bigfloat
Cfunction(const bigint & nh,const rational_factorization & dfact,int pmult)2664 quadratic_order::Cfunction(const bigint & nh,
2665                            const rational_factorization & dfact,
2666                            int pmult)
2667 {
2668 	debug_handler("quadratic_order", "Cfunction(bigint)");
2669 
2670 	bigfloat CC;
2671 	bigfloat_int CC_int;
2672 	bigint oh = h;
2673 	rational_factorization ofact = disc_fact;
2674 
2675 	h = nh;
2676 	disc_fact = dfact;
2677 	if (pmult) {
2678 		CC_int = Cfunction_bigfloat(pmult);
2679 		CC_int.approx(CC);
2680 		std::cout << "-->C in " << CC_int << ", = " << CC << std::endl;
2681 	}
2682 	else
2683 		CC = Cfunction();
2684 	h = oh;
2685 	disc_fact = ofact;
2686 
2687 	return CC;
2688 }
2689 
2690 
2691 
2692 //
2693 // quadratic_order::Cfunction
2694 //
2695 // Task:
2696 //      returns an approximation of C(Delta) accurate to 8 significant digits.
2697 //
2698 
2699 bigfloat
Cfunction(const bigint & nh,const bigfloat & nR,const rational_factorization & dfact,int pmult)2700 quadratic_order::Cfunction(const bigint & nh,
2701                            const bigfloat & nR,
2702                            const rational_factorization & dfact,
2703                            int pmult)
2704 {
2705 	debug_handler("quadratic_order", "Cfunction(bigint, bigfloat)");
2706 
2707 	bigfloat CC, oR = R;
2708 	bigfloat_int CC_int;
2709 	bigint oh = h;
2710 	rational_factorization ofact = disc_fact;
2711 
2712 	h = nh;
2713 	R = nR;
2714 	disc_fact = dfact;
2715 	if (pmult) {
2716 		CC_int = Cfunction_bigfloat(pmult);
2717 		CC_int.approx(CC);
2718 		std::cout << "-->C in " << CC_int << ", = " << CC << std::endl;
2719 	}
2720 	else
2721 		CC = Cfunction();
2722 	h = oh;
2723 	R = oR;
2724 	disc_fact = ofact;
2725 
2726 	return CC;
2727 }
2728 
2729 
2730 
2731 //
2732 // quadratic_order::regulator
2733 //
2734 // Task:
2735 //      returns an approximation of the regulator.  The algorithm is chosen
2736 //      based on the size of the discriminant.
2737 //
2738 
2739 bigfloat
regulator()2740 quadratic_order::regulator()
2741 {
2742 	debug_handler("quadratic_order", "regulator");
2743 
2744 	int num;
2745 	bigfloat outR;
2746 	quadratic_order *QO = qi_class::get_current_order_ptr();
2747 	long oprec = bigfloat::get_precision();
2748 
2749 	bigfloat::set_precision(prec);
2750 	qi_class::set_current_order(*this);
2751 
2752 	// set last-used quadratic_order
2753 	qo_l().set_last(this);
2754 
2755 	if (!is_R_computed()) {
2756 		// R has not been computed yet - compute it
2757 
2758 		if (is_imaginary())
2759 			R.assign_one();
2760 
2761 		else {
2762 			num = decimal_length(Delta);
2763 			if (num < RBJTB)
2764 				regulator_BJT();
2765 			else if (num < RSHANKSB)
2766 				regulator_shanks();
2767 			else
2768 				class_group_siqs();
2769 		}
2770 	}
2771 
2772 	outR.assign(R);
2773 	bigfloat::set_precision(oprec);
2774 	qi_class::set_current_order(*QO);
2775 
2776 	return outR;
2777 }
2778 
2779 
2780 
2781 //
2782 // quadratic_order::regulator(long k)
2783 //
2784 // Task:
2785 //      returns an absolute k approximation to the regulator.
2786 //      The algorithm is chosen based on the size of the discriminant.
2787 //
2788 
2789 xbigfloat quadratic_order
regulator(long k)2790 ::regulator(long k)
2791 {
2792 	debug_handler("quadratic_order", "regulator(long k");
2793 
2794 	xbigfloat outR;
2795 
2796 	if (is_imaginary())
2797 		outR.assign_one();
2798 
2799 	else {
2800 		lidia_error_handler("quadratic_order::regulator(long)",
2801 				    "Not yet implemented.");
2802 
2803 		if (R_xbig.is_zero() || R_xbig_prec < k+1) {
2804 			this->fundamental_unit();
2805 			R_xbig = fund_unit.get_absolute_Ln_approximation(k);
2806 			R_xbig_prec = k;
2807 			outR = R_xbig;
2808 		}
2809 		else {
2810 			truncate(outR, R_xbig, R_xbig.exponent()+k+1);
2811 		}
2812 	}
2813 
2814 	return outR;
2815 }
2816 
2817 
2818 
2819 //
2820 // quadratic_order::class_number
2821 //
2822 // Task:
2823 //      returns the class number.  The algorithm is chosen based on the size
2824 //      of the discriminant.
2825 //
2826 
2827 bigint
class_number()2828 quadratic_order::class_number()
2829 {
2830 	debug_handler("quadratic_order", "class_number");
2831 
2832 	int num;
2833 
2834 	// set last-used quadratic_order
2835 	qo_l().set_last(this);
2836 
2837 	if (!is_h_computed()) {
2838 		quadratic_order *QO = qi_class::get_current_order_ptr();
2839 		long oprec = bigfloat::get_precision();
2840 
2841 		bigfloat::set_precision(prec);
2842 		qi_class::set_current_order(*this);
2843 
2844 		// h has not been computed yet - compute it
2845 		num = decimal_length(Delta);
2846 
2847 		if (is_real()) {
2848 			if (num < CNBJT_RB)
2849 				class_group_BJT();
2850 			if (num < CNSHANKS_RB)
2851 				class_number_shanks();
2852 			else
2853 				class_group_siqs();
2854 		}
2855 		else {
2856 			if (num < CNBJT_IB)
2857 				class_group_BJT();
2858 			if (num < CNSHANKS_IB)
2859 				class_number_shanks();
2860 			else
2861 				class_group_siqs();
2862 		}
2863 
2864 		bigfloat::set_precision(oprec);
2865 		qi_class::set_current_order(*QO);
2866 	}
2867 
2868 	return h;
2869 }
2870 
2871 
2872 
2873 //
2874 // quadratic_order::class_group
2875 //
2876 // Task:
2877 //      returns a vector containing the invariants of the class group.  The
2878 //      algorithm is chosen based on the size of the discriminant.
2879 //
2880 
2881 base_vector< bigint >
class_group()2882 quadratic_order::class_group()
2883 {
2884 	debug_handler("quadratic_order", "class_group");
2885 
2886 	int num;
2887 
2888 	// set last-used quadratic_order
2889 	qo_l().set_last(this);
2890 
2891 	if (!is_CL_computed()) {
2892 		quadratic_order *QO = qi_class::get_current_order_ptr();
2893 		long oprec = bigfloat::get_precision();
2894 
2895 		bigfloat::set_precision(prec);
2896 		qi_class::set_current_order(*this);
2897 
2898 		// CL has not been computed yet - compute it
2899 
2900 		if (h.is_gt_zero()) {
2901 			// class number is known
2902 //      class_group_h();
2903 			h.assign_zero();
2904 			class_group();
2905 		}
2906 		else {
2907 			num = decimal_length(Delta);
2908 
2909 			if (is_real()) {
2910 				if (num < CGBJT_RB)
2911 					class_group_BJT();
2912 				else if (num < CGRAND_RB)
2913 					class_group_randexp();
2914 				else
2915 					class_group_siqs();
2916 			}
2917 			else {
2918 				if (num < CGBJT_IB)
2919 					class_group_BJT();
2920 				else if (num < CGSHANKS_IB)
2921 					class_group_shanks();
2922 				else
2923 					class_group_siqs();
2924 			}
2925 		}
2926 
2927 		bigfloat::set_precision(oprec);
2928 		qi_class::set_current_order(*QO);
2929 	}
2930 
2931 	return CL;
2932 }
2933 
2934 
2935 
2936 //
2937 // quadratic_order::regulator(bool)
2938 //
2939 // Task:
2940 //      returns an approximation of the regulator.  If the parameter is true,
2941 //      then the regulator will be computed using the subexponential algorithm
2942 //      and recomputed if it was previously computed with a different
2943 //      algorithm.
2944 //
2945 
2946 bigfloat
regulator(bool usesub)2947 quadratic_order::regulator(bool usesub)
2948 {
2949 	debug_handler("quadratic_order", "regulator(bool)");
2950 
2951 	bigfloat outR;
2952 	quadratic_order *QO = qi_class::get_current_order_ptr();
2953 	long oprec = bigfloat::get_precision();
2954 
2955 	bigfloat::set_precision(prec);
2956 	qi_class::set_current_order(*this);
2957 
2958 	// set last-used quadratic_order
2959 	qo_l().set_last(this);
2960 
2961 	if (usesub) {
2962 		if ((!is_R_computed()) || ((is_R_computed()) && (!is_subexp_computed()))) {
2963 			assign(Delta);
2964 			class_group_siqs();
2965 		}
2966 	}
2967 	else
2968 		regulator();
2969 
2970 	outR.assign(R);
2971 	bigfloat::set_precision(oprec);
2972 	qi_class::set_current_order(*QO);
2973 
2974 	return outR;
2975 }
2976 
2977 
2978 
2979 //
2980 // quadratic_order::class_group(bool)
2981 //
2982 // Task:
2983 //      returns a vector containing the invariants of the class group.  If the
2984 //      parameter is true, then the class group will be computed using the
2985 //      subexponential algorithm, and recomputed if it was previously computed
2986 //      with a different algorithm.
2987 //
2988 
2989 base_vector< bigint >
class_group(bool usesub)2990 quadratic_order::class_group(bool usesub)
2991 {
2992 	debug_handler("quadratic_order", "class_group(bool)");
2993 
2994 	quadratic_order *QO = qi_class::get_current_order_ptr();
2995 	long oprec = bigfloat::get_precision();
2996 
2997 	bigfloat::set_precision(prec);
2998 	qi_class::set_current_order(*this);
2999 
3000 	// set last-used quadratic_order
3001 	qo_l().set_last(this);
3002 
3003 	if (usesub) {
3004 		if ((!is_CL_computed()) ||
3005 		    ((is_CL_computed()) && (!is_subexp_computed()))) {
3006 			assign(Delta);
3007 			class_group_siqs();
3008 		}
3009 	}
3010 	else
3011 		class_group();
3012 
3013 	bigfloat::set_precision(oprec);
3014 	qi_class::set_current_order(*QO);
3015 
3016 	return CL;
3017 }
3018 
3019 
3020 
3021 //
3022 // quadratic_order::regulator_BJT
3023 //
3024 // Task:
3025 //      computes an approximation of the regulator using a variation of
3026 //      baby-step giant-step due to Buchmann.
3027 //
3028 
3029 bigfloat
regulator_BJT()3030 quadratic_order::regulator_BJT()
3031 {
3032 	debug_handler("quadratic_order", "regulator_BJT");
3033 
3034 	qi_class_real A, UNIT;
3035 	timer t;
3036 	bigfloat outR;
3037 	quadratic_order *QO = qi_class::get_current_order_ptr();
3038 	long oprec = bigfloat::get_precision();
3039 
3040 	bigfloat::set_precision(prec);
3041 	qi_class::set_current_order(*this);
3042 
3043 	// set last-used quadratic_order
3044 	qo_l().set_last(this);
3045 
3046 	if (!is_R_computed()) {
3047 		if (is_imaginary())
3048 			R.assign_one();
3049 		else {
3050 			if (prin_list.no_of_elements() > 0) {
3051 				A = prin_list.last_entry();
3052 				if ((A.is_one()) && (!R.is_zero())) {
3053 					UNIT.assign_one();
3054 					A.assign(nearest(UNIT, A.get_distance()));
3055 					R.assign(A.get_distance());
3056 				}
3057 			}
3058 
3059 			if (R.is_zero()) {
3060 				if (info)
3061 					t.start_timer();
3062 
3063 				rBJT();
3064 
3065 				if (info) {
3066 					t.stop_timer();
3067 					std::cout << "\nelapsed CPU time (regulator_BJT) = ";
3068 					MyTime(t.user_time());
3069 					std::cout << "\n";
3070 				}
3071 			}
3072 
3073 			if (do_verify) {
3074 				if (verify_regulator()) {
3075 					if (info)
3076 						std::cout << "Regulator is correct." << std::endl;
3077 				}
3078 				else
3079 					std::cout << "Regulator is not correct!" << std::endl;
3080 			}
3081 		}
3082 	}
3083 
3084 	outR.assign(R);
3085 	bigfloat::set_precision(oprec);
3086 	qi_class::set_current_order(*QO);
3087 
3088 	return outR;
3089 }
3090 
3091 
3092 
3093 //
3094 // quadratic_order::regulator_shanks
3095 //
3096 // Task:
3097 //      computes an approximation of the regulator using the O(D^1/5) algorithm
3098 //      due to Mollin and Williams.
3099 //
3100 
3101 bigfloat
regulator_shanks()3102 quadratic_order::regulator_shanks()
3103 {
3104 	debug_handler("quadratic_order", "regulator_shanks");
3105 
3106 	qi_class_real A, UNIT;
3107 	timer t;
3108 	bigfloat outR;
3109 	quadratic_order *QO = qi_class::get_current_order_ptr();
3110 	long oprec = bigfloat::get_precision();
3111 
3112 	bigfloat::set_precision(prec);
3113 	qi_class::set_current_order(*this);
3114 
3115 	// set last-used quadratic_order
3116 	qo_l().set_last(this);
3117 
3118 	if (!is_R_computed()) {
3119 		if (is_imaginary())
3120 			R.assign_one();
3121 		else {
3122 			if (prin_list.no_of_elements() > 0) {
3123 				A = prin_list.last_entry();
3124 				if ((A.is_one()) && (!R.is_zero()))
3125 					UNIT.assign_one();
3126 				A.assign(nearest(UNIT, A.get_distance()));
3127 				R.assign(A.get_distance());
3128 			}
3129 
3130 			if (R.is_zero()) {
3131 				if (info)
3132 					t.start_timer();
3133 
3134 				rshanks();
3135 
3136 				if (info) {
3137 					t.stop_timer();
3138 					std::cout << "\nelapsed CPU time (regulator_shanks) = ";
3139 					MyTime(t.user_time());
3140 					std::cout << "\n";
3141 				}
3142 			}
3143 
3144 			if (do_verify) {
3145 				if (verify_regulator()) {
3146 					if (info)
3147 						std::cout << "Regulator is correct." << std::endl;
3148 				}
3149 				else
3150 					std::cout << "Regulator is not correct!" << std::endl;
3151 			}
3152 		}
3153 	}
3154 
3155 	outR.assign(R);
3156 	bigfloat::set_precision(oprec);
3157 	qi_class::set_current_order(*QO);
3158 
3159 	return outR;
3160 }
3161 
3162 
3163 
3164 //
3165 // quadratic_order::verify_regulator()
3166 //
3167 // Task:
3168 //      verifies whether the regulator is correct.
3169 //
3170 
3171 bool
verify_regulator()3172 quadratic_order::verify_regulator()
3173 {
3174 	debug_handler("quadratic_order", "verify_regulator");
3175 
3176 	bool is_reg;
3177 	timer t;
3178 	quadratic_order *QO = qi_class::get_current_order_ptr();
3179 
3180 	if (!is_real())  return true;
3181 
3182 	qi_class::set_current_order(*this);
3183 
3184 	// set last-used quadratic_order
3185 	qo_l().set_last(this);
3186 
3187 	if (!is_R_computed())
3188 		regulator();
3189 
3190 	if (info)
3191 		t.start_timer();
3192 
3193 
3194 	if (info > 1) {
3195 		std::cout << "\nVerifying regulator:" << std::endl;
3196 		std::cout << "  R = " << R << std::endl;
3197 	}
3198 
3199 	//<MM>
3200 	// Requires that R is an absolute 3-approximation
3201 	// to the regulator.
3202 	//
3203 
3204 	is_reg = could_be_regulator_multiple(R);
3205 
3206 	//</MM>
3207 
3208 
3209 	if (info) {
3210 		t.stop_timer();
3211 		std::cout << "\nelapsed CPU time (verify_regulator) = ";
3212 		MyTime(t.user_time());
3213 		std::cout << "\n";
3214 	}
3215 
3216 	qi_class::set_current_order(*QO);
3217 
3218 	return is_reg;
3219 }
3220 
3221 
3222 
3223 //
3224 // quadratic_order::class_number_shanks
3225 //
3226 // Task:
3227 //      computes the class number using the O(D^1/5) algorithm due to Shanks.
3228 //
3229 
3230 bigint
class_number_shanks()3231 quadratic_order::class_number_shanks()
3232 {
3233 	debug_handler("quadratic_order", "class_number_shanks");
3234 
3235 	timer t;
3236 
3237 	// set last-used quadratic_order
3238 	qo_l().set_last(this);
3239 
3240 	if (!is_h_computed()) {
3241 		quadratic_order *QO = qi_class::get_current_order_ptr();
3242 		long oprec = bigfloat::get_precision();
3243 
3244 		bigfloat::set_precision(prec);
3245 		qi_class::set_current_order(*this);
3246 
3247 		if (info)
3248 			t.start_timer();
3249 
3250 		if (is_imaginary())
3251 			cns_imag();
3252 		else
3253 			cns_real();
3254 
3255 		if (info) {
3256 			t.stop_timer();
3257 			std::cout << "\nelapsed CPU time (class_number_shanks) = ";
3258 			MyTime(t.user_time());
3259 			std::cout << "\n";
3260 		}
3261 
3262 		bigfloat::set_precision(oprec);
3263 		qi_class::set_current_order(*QO);
3264 	}
3265 
3266 	return h;
3267 }
3268 
3269 
3270 
3271 //
3272 // quadratic_order::class_group_BJT
3273 //
3274 // Task:
3275 //      computes the class group using a variation of baby-step giant-step
3276 //      due to Buchmann, Jacobson, and Teske.
3277 //
3278 
3279 base_vector< bigint >
class_group_BJT()3280 quadratic_order::class_group_BJT()
3281 {
3282 	debug_handler("quadratic_order", "class_group_BJT");
3283 
3284 	timer t;
3285 
3286 	// set last-used quadratic_order
3287 	qo_l().set_last(this);
3288 
3289 	if (!is_CL_computed()) {
3290 		quadratic_order *QO = qi_class::get_current_order_ptr();
3291 		long oprec = bigfloat::get_precision();
3292 
3293 		bigfloat::set_precision(prec);
3294 		qi_class::set_current_order(*this);
3295 
3296 		if (info)
3297 			t.start_timer();
3298 
3299 		if (is_imaginary())
3300 			cgBJT_imag();
3301 		else
3302 			cgBJT_real();
3303 
3304 		if (info) {
3305 			t.stop_timer();
3306 			std::cout << "\nelapsed CPU time (class_group_BJT) = ";
3307 			MyTime(t.user_time());
3308 			std::cout << "\n";
3309 		}
3310 
3311 		if (do_verify) {
3312 			if (verify_class_group()) {
3313 				if (info)
3314 					std::cout << "Class group is correct." << std::endl;
3315 			}
3316 			else
3317 				std::cout << "Class group is not correct!" << std::endl;
3318 		}
3319 
3320 		bigfloat::set_precision(oprec);
3321 		qi_class::set_current_order(*QO);
3322 	}
3323 
3324 	return CL;
3325 }
3326 
3327 
3328 
3329 //
3330 // quadratic_order::class_group_shanks
3331 //
3332 // Task:
3333 //      computes the class group using an unpublished variation of Shanks
3334 //      baby-step giant-step algorithm.  It is essentially the same as the
3335 //      class number algorith, except we continue computing the subgroup until
3336 //      we have the entire class group.
3337 //
3338 
3339 base_vector< bigint >
class_group_shanks()3340 quadratic_order::class_group_shanks()
3341 {
3342 	debug_handler("quadratic_order", "class_group_shanks");
3343 
3344 	bigint temp;
3345 	bigfloat d1, d2;
3346 	timer t;
3347 
3348 	// set last-used quadratic_order
3349 	qo_l().set_last(this);
3350 
3351 	if (!is_CL_computed()) {
3352 		quadratic_order *QO = qi_class::get_current_order_ptr();
3353 		long oprec = bigfloat::get_precision();
3354 
3355 		bigfloat::set_precision(prec);
3356 		qi_class::set_current_order(*this);
3357 
3358 		d1.assign(0);
3359 		d2.assign(0);
3360 
3361 		if (info)
3362 			t.start_timer();
3363 
3364 		if (is_imaginary())
3365 			cgs_imag(d1, d2, temp);
3366 		else
3367 			cgs_real(d1, temp);
3368 
3369 		if (info) {
3370 			t.stop_timer();
3371 			std::cout << "\nelapsed CPU time (class_group_shanks) = ";
3372 			MyTime(t.user_time());
3373 			std::cout << "\n";
3374 		}
3375 
3376 		if (do_verify) {
3377 			if (verify_class_group()) {
3378 				if (info)
3379 					std::cout << "Class group is correct." << std::endl;
3380 			}
3381 			else
3382 				std::cout << "Class group is not correct!" << std::endl;
3383 		}
3384 
3385 		bigfloat::set_precision(oprec);
3386 		qi_class::set_current_order(*QO);
3387 	}
3388 
3389 	return CL;
3390 }
3391 
3392 
3393 
3394 //
3395 // quadratic_order::class_group_h
3396 //
3397 // Task:
3398 //      computes the class group given the class number.
3399 //
3400 
3401 base_vector< bigint >
class_group_h()3402 quadratic_order::class_group_h()
3403 {
3404 	debug_handler("quadratic_order", "class_group_h");
3405 
3406 	timer t;
3407 
3408 	// set last-used quadratic_order
3409 	qo_l().set_last(this);
3410 
3411 	if (!is_CL_computed()) {
3412 		quadratic_order *QO = qi_class::get_current_order_ptr();
3413 		long oprec = bigfloat::get_precision();
3414 
3415 		bigfloat::set_precision(prec);
3416 		qi_class::set_current_order(*this);
3417 
3418 		if (info)
3419 			t.start_timer();
3420 
3421 		// compute h if necessary
3422 		if (h.is_zero());
3423 		class_number();
3424 
3425 		if (is_imaginary())
3426 			cgh_imag();
3427 		else
3428 			cgh_real();
3429 
3430 		if (info) {
3431 			t.stop_timer();
3432 			std::cout << "\nelapsed CPU time (class_group_h) = ";
3433 			MyTime(t.user_time());
3434 			std::cout << "\n";
3435 		}
3436 
3437 		if (do_verify) {
3438 			if (verify_class_group()) {
3439 				if (info)
3440 					std::cout << "Class group is correct." << std::endl;
3441 			}
3442 			else
3443 				std::cout << "Class group is not correct!" << std::endl;
3444 		}
3445 
3446 		bigfloat::set_precision(oprec);
3447 		qi_class::set_current_order(*QO);
3448 	}
3449 
3450 	return CL;
3451 }
3452 
3453 
3454 
3455 //
3456 // quadratic_order::class_group_randexp
3457 //
3458 // Task:
3459 //      computes the class group using a subexponential algorithm with
3460 //      random exponent vectors (Duellmann's variation)
3461 //
3462 
3463 base_vector< bigint >
class_group_randexp()3464 quadratic_order::class_group_randexp()
3465 {
3466 	debug_handler("quadratic_order", "class_group_randexp");
3467 
3468 	int decimal_digits;
3469 	timer t;
3470 
3471 	// set last-used quadratic_order
3472 	qo_l().set_last(this);
3473 
3474 	if (!is_CL_computed()) {
3475 		quadratic_order *QO = qi_class::get_current_order_ptr();
3476 		long oprec = bigfloat::get_precision();
3477 
3478 		bigfloat::set_precision(prec);
3479 		qi_class::set_current_order(*this);
3480 
3481 		if (info)
3482 			t.start_timer();
3483 
3484 		decimal_digits = decimal_length(Delta);
3485 		if (decimal_digits < 7)
3486 			class_group();
3487 		else {
3488 			subused = true;
3489 			cg_randexp();
3490 		}
3491 
3492 		if (info) {
3493 			t.stop_timer();
3494 			std::cout << "\nelapsed CPU time (class_group_randexp) = ";
3495 			MyTime(t.user_time());
3496 			std::cout << "\n";
3497 		}
3498 
3499 		if (do_verify) {
3500 			if (verify_regulator()) {
3501 				if (info)
3502 					std::cout << "Regulator is correct." << std::endl;
3503 			}
3504 			else
3505 				std::cout << "Regulator is not correct!" << std::endl;
3506 
3507 			if (verify_cg(false)) {
3508 				if (info)
3509 					std::cout << "Class group is correct." << std::endl;
3510 			}
3511 			else
3512 				std::cout << "Class group is not correct!" << std::endl;
3513 		}
3514 
3515 		bigfloat::set_precision(oprec);
3516 		qi_class::set_current_order(*QO);
3517 	}
3518 
3519 	return CL;
3520 }
3521 
3522 
3523 
3524 //
3525 // quadratic_order::class_group_mpqs
3526 //
3527 // Task:
3528 //      computes the class group using an algorithm based on the MPQS
3529 //      factoring algorithm due to Jacobson.
3530 //
3531 
3532 base_vector< bigint >
class_group_mpqs()3533 quadratic_order::class_group_mpqs()
3534 {
3535 	debug_handler("quadratic_order", "class_group_mpqs");
3536 
3537 	int decimal_digits;
3538 	timer t;
3539 
3540 	// set last-used quadratic_order
3541 	qo_l().set_last(this);
3542 
3543 	if (!is_CL_computed()) {
3544 		quadratic_order *QO = qi_class::get_current_order_ptr();
3545 		long oprec = bigfloat::get_precision();
3546 
3547 		bigfloat::set_precision(prec);
3548 		qi_class::set_current_order(*this);
3549 
3550 		if (info)
3551 			t.start_timer();
3552 
3553 		decimal_digits = decimal_length(Delta);
3554 		if (decimal_digits < 7)
3555 			class_group();
3556 		else {
3557 			subused = true;
3558 			cg_mpqs(false, false);
3559 		}
3560 
3561 		if (info) {
3562 			t.stop_timer();
3563 			std::cout << "\nelapsed CPU time (class_group_mpqs) = ";
3564 			MyTime(t.user_time());
3565 			std::cout << "\n";
3566 		}
3567 
3568 		if (do_verify) {
3569 			if (verify_regulator()) {
3570 				if (info)
3571 					std::cout << "Regulator is correct." << std::endl;
3572 			}
3573 			else
3574 				std::cout << "Regulator is not correct!" << std::endl;
3575 
3576 			if (verify_class_group()) {
3577 				if (info)
3578 					std::cout << "Class group is correct." << std::endl;
3579 			}
3580 			else
3581 				std::cout << "Class group is not correct!" << std::endl;
3582 		}
3583 
3584 		bigfloat::set_precision(oprec);
3585 		qi_class::set_current_order(*QO);
3586 	}
3587 
3588 	return CL;
3589 }
3590 
3591 
3592 
3593 //
3594 // quadratic_order::class_group_siqs
3595 //
3596 // Task:
3597 //      computes the class group using an algorithm based on the
3598 //      self-initializing MPQS factoring algorithm due to Jacobson.
3599 //
3600 
3601 base_vector< bigint >
class_group_siqs()3602 quadratic_order::class_group_siqs()
3603 {
3604 	debug_handler("quadratic_order", "class_group_siqs");
3605 
3606 	int decimal_digits;
3607 	timer t;
3608 
3609 	// set last-used quadratic_order
3610 	qo_l().set_last(this);
3611 
3612 	if (!is_CL_computed()) {
3613 		quadratic_order *QO = qi_class::get_current_order_ptr();
3614 		long oprec = bigfloat::get_precision();
3615 
3616 		bigfloat::set_precision(prec);
3617 		qi_class::set_current_order(*this);
3618 
3619 		if (info)
3620 			t.start_timer();
3621 
3622 		decimal_digits = decimal_length(Delta);
3623 		if (decimal_digits < 7)
3624 			class_group();
3625 		else {
3626 			subused = true;
3627 			cg_mpqs(true, false);
3628 		}
3629 
3630 		if (info) {
3631 			t.stop_timer();
3632 			std::cout << "\nelapsed CPU time (class_group_siqs) = ";
3633 			MyTime(t.user_time());
3634 			std::cout << "\n";
3635 		}
3636 
3637 		if (do_verify) {
3638 			if (verify_regulator()) {
3639 				if (info)
3640 					std::cout << "Regulator is correct." << std::endl;
3641 			}
3642 			else
3643 				std::cout << "Regulator is not correct!" << std::endl;
3644 
3645 			if (verify_class_group()) {
3646 				if (info)
3647 					std::cout << "Class group is correct." << std::endl;
3648 			}
3649 			else
3650 				std::cout << "Class group is not correct!" << std::endl;
3651 		}
3652 
3653 		bigfloat::set_precision(oprec);
3654 		qi_class::set_current_order(*QO);
3655 	}
3656 
3657 	return CL;
3658 }
3659 
3660 
3661 
3662 //
3663 //
3664 // quadratic_order::class_group_lp
3665 //
3666 // Task:
3667 //      computes the class group using large primes an algorithm based on the
3668 //      self-initializing MPQS factoring algorithm due to Jacobson.
3669 //
3670 
3671 base_vector< bigint >
class_group_lp()3672 quadratic_order::class_group_lp()
3673 {
3674 	debug_handler("quadratic_order", "class_group_lp");
3675 
3676 	int decimal_digits;
3677 	timer t;
3678 
3679 	// set last-used quadratic_order
3680 	qo_l().set_last(this);
3681 
3682 	if (!is_CL_computed()) {
3683 		quadratic_order *QO = qi_class::get_current_order_ptr();
3684 		long oprec = bigfloat::get_precision();
3685 
3686 		bigfloat::set_precision(prec);
3687 		qi_class::set_current_order(*this);
3688 
3689 		if (info)
3690 			t.start_timer();
3691 
3692 		decimal_digits = decimal_length(Delta);
3693 		if (decimal_digits < 7)
3694 			class_group();
3695 		else {
3696 			subused = true;
3697 			cg_mpqs(true, true);
3698 		}
3699 
3700 		if (info) {
3701 			t.stop_timer();
3702 			std::cout << "\nelapsed CPU time (class_group_lp) = ";
3703 			MyTime(t.user_time());
3704 			std::cout << "\n";
3705 		}
3706 
3707 		if (do_verify) {
3708 			if (verify_regulator()) {
3709 				if (info)
3710 					std::cout << "Regulator is correct." << std::endl;
3711 			}
3712 			else
3713 				std::cout << "Regulator is not correct!" << std::endl;
3714 
3715 			if (verify_class_group()) {
3716 				if (info)
3717 					std::cout << "Class group is correct." << std::endl;
3718 			}
3719 			else
3720 				std::cout << "Class group is not correct!" << std::endl;
3721 		}
3722 
3723 		bigfloat::set_precision(oprec);
3724 		qi_class::set_current_order(*QO);
3725 	}
3726 
3727 	return CL;
3728 }
3729 
3730 
3731 
3732 //
3733 // quadratic_order::count_primes_to_verify(int BachBound)
3734 //
3735 // Task:
3736 //      counts the number of primes not in the factor base but less than
3737 //      Bach's bound.
3738 //
3739 
3740 int
count_primes_to_verify(int BachBound)3741 quadratic_order::count_primes_to_verify(int BachBound)
3742 {
3743 	debug_handler("quadratic_order", "count_primes_to_verify(int)");
3744 
3745 	int nump = 0, p, size_FB;
3746 	qi_class P;
3747 	prime_list PList;
3748 	quadratic_order *QO = qi_class::get_current_order_ptr();
3749 
3750 	qi_class::set_current_order(*this);
3751 
3752 	if (fact_base.size() == 0) {
3753 		qo_read_parameters(decimal_length(Delta), size_FB);
3754 		qo_get_fb(size_FB);
3755 		fact_base[fact_base.size()-1].get_a().intify(p);
3756 		fact_base.reset();
3757 	}
3758 	else
3759 		fact_base[fact_base.size()-1].get_a().intify(p);
3760 
3761 	if (p == 2)
3762 		++p;
3763 	else
3764 		p += 2;
3765 	if (p <= BachBound) {
3766 		PList.resize(p, BachBound+10);
3767 		p = PList.get_first_prime();
3768 		while ((p <= BachBound) && (p > 0)) {
3769 			if (generate_prime_ideal(P, p))
3770 				++nump;
3771 			p = PList.get_next_prime();
3772 		}
3773 	}
3774 
3775 	qi_class::set_current_order(*QO);
3776 
3777 	return nump;
3778 }
3779 
3780 
3781 
3782 //
3783 // quadratic_order::verify_class_group(int level)
3784 //
3785 // Task:
3786 //      verifies whether the class group is correct using the given level
3787 //
3788 
3789 bool
verify_class_group(int level)3790 quadratic_order::verify_class_group(int level)
3791 {
3792 	debug_handler("quadratic_order", "verify_class_group(int)");
3793 
3794 	int templ;
3795 	bool temp = true;;
3796 
3797 	if (level) {
3798 		templ = do_verify;
3799 		do_verify = level;
3800 
3801 		if (sieve.is_initialized())
3802 			temp = verify_cg(true);
3803 		else
3804 			temp = verify_cg(false);
3805 
3806 		do_verify = templ;
3807 	}
3808 
3809 	return temp;
3810 }
3811 
3812 
3813 
3814 //
3815 // quadratic_order::verify_class_group()
3816 //
3817 // Task:
3818 //      verifies whether the class group is correct.
3819 //
3820 
3821 bool
verify_class_group()3822 quadratic_order::verify_class_group()
3823 {
3824 	debug_handler("quadratic_order", "verify_class_group");
3825 
3826 	if (sieve.is_initialized())
3827 		return verify_cg(true);
3828 	else
3829 		return verify_cg(false);
3830 }
3831 
3832 
3833 
3834 //
3835 // quadratic_order::verify_cg()
3836 //
3837 // Task:
3838 //      verifies whether the class group is correct.
3839 //
3840 
3841 bool
verify_cg(bool use_sieve)3842 quadratic_order::verify_cg(bool use_sieve)
3843 {
3844 	debug_handler("quadratic_order", "verify_cg");
3845 
3846 	lidia_size_t i, j, numFB, numCL, matcols, matrows;
3847 	int p, BachBound = 0, num, nump = 0, countp = 0;
3848 	bigint ord;
3849 	qi_class A, C, P;
3850 	math_vector< bigint > vec, rvec, nvec, pvec;
3851 	sort_vector< int > D_primes;
3852 	base_vector< qi_class > new_FB;
3853 	matrix< bigint > vmat;
3854 	bool OK;
3855 	timer t, tprime;
3856 	prime_list PList;
3857 	quadratic_number_standard q;
3858 	quadratic_order *QO = qi_class::get_current_order_ptr();
3859 	long oprec = bigfloat::get_precision();
3860 
3861 	// FOR TABLE GENERATION
3862 	if (special)
3863 		t.start_timer();
3864 
3865 	bigfloat::set_precision(prec);
3866 	qi_class::set_current_order(*this);
3867 
3868 	// set last-used quadratic_order
3869 	qo_l().set_last(this);
3870 
3871 	if (info)
3872 		tprime.start_timer();
3873 
3874 	if (info)
3875 		t.start_timer();
3876 
3877 	OK = true;
3878 
3879 	// verify primes less than Bach bound
3880 	numFB = fact_base.size();
3881 	numCL = CL.size();
3882 	if (numFB == 0) {
3883 		qo_read_parameters(decimal_length(Delta), numFB);
3884 		std::cout << "\nInput factor base size:  ";
3885 		std::cin >> numFB;
3886 		std::cout << std::endl;
3887 		qo_get_fb(numFB);
3888 	}
3889 
3890 
3891 	if (do_verify > 1) {
3892 		BachBound = bach_bound();
3893 
3894 #if 0
3895 		std::cout << "Input bach bound:  ";
3896 		std::cin >> BachBound;
3897 		std::cout << std::endl;
3898 #endif
3899 
3900 		nump = count_primes_to_verify(BachBound);
3901 
3902 		// FOR TABLE GENERATION
3903 		if (special)
3904 			std::cout << nump << "  " << std::flush;
3905 	}
3906 	else {
3907 		// FOR TABLE GENERATION
3908 		if (special)
3909 			std::cout << "0  " << std::flush;
3910 	}
3911 
3912 	if (info > 1) {
3913 		std::cout << "\nVerifying class group (level " << do_verify << "):" << std::endl;
3914 		std::cout << "  CL = " << CL << std::endl;
3915 		std::cout << "  generators = " << gens << std::endl;
3916 	}
3917 
3918 	num = decimal_length(Delta);
3919 	if (is_subexp_computed() || !is_CL_computed()) {
3920 		if (info > 1) {
3921 			if (do_verify == 2) {
3922 				std::cout << "  BachBound = " << BachBound << std::endl;
3923 				std::cout << "Representing each prime ideal less than Bach's bound over "
3924 					"the factor base:" << std::endl;
3925 				std::cout << nump << " primes to verify!" << std::endl;
3926 			}
3927 			else if (do_verify > 2) {
3928 				std::cout << "  BachBound = " << BachBound << std::endl;
3929 				std::cout << "Representing each prime ideal less than Bach's bound over "
3930 					"the generators:" << std::endl;
3931 				std::cout << nump+fact_base.size() << " primes to verify!" << std::endl;
3932 			}
3933 		}
3934 
3935 		pvec.set_mode(EXPAND);
3936 		nvec.set_mode(EXPAND);
3937 		vec.set_mode(EXPAND);
3938 		vec.set_size(numFB);
3939 		rvec.set_capacity(numCL);
3940 		rvec.set_size(numCL);
3941 
3942 		if (do_verify > 2) {
3943 			generators();
3944 
3945 			// verify primes in factor base
3946 			for (i = 0; i < numFB; ++i) {
3947 				for (j = 0; j < numFB; ++j)
3948 					vec[j].assign_zero();
3949 				vec[i].assign_one();
3950 				rvec = represent_over_generators(vec);
3951 
3952 				C.assign_one();
3953 				for (j = 0; j < numCL; ++j) {
3954 					power(A, gens[j], rvec[j]);
3955 					multiply(C, C, A);
3956 				}
3957 
3958 				if (fact_base[i].is_equivalent(C)) {
3959 					if (info > 2) {
3960 						std::cout << "  " << fact_base[i] << " ~ " << rvec << " ~ " << C << std::endl;
3961 					}
3962 				}
3963 				else {
3964 					OK = false;
3965 					if (info > 2) {
3966 						std::cout << "  ERROR:  " << fact_base[i] << " <>" << rvec << " <>";
3967 						std::cout << C << std::endl;
3968 					}
3969 				}
3970 			}
3971 		}
3972 
3973 		if (do_verify > 1) {
3974 			matcols = 0;
3975 			vmat.set_representation(matrix_flags::sparse_representation);
3976 			vmat.set_orientation(matrix_flags::column_oriented);
3977 
3978 			new_FB.set_mode(EXPAND);
3979 			new_FB = fact_base;
3980 
3981 			fact_base[numFB-1].get_a().intify(p);
3982 			p += 2;
3983 			if (p <= BachBound) {
3984 				PList.resize(p, BachBound+10);
3985 			}
3986 			p = PList.get_first_prime();
3987 			while ((p <= BachBound) && (p > 0)) {
3988 				if (generate_prime_ideal(P, p)) {
3989 					TR.touch_files();
3990 					++countp;
3991 
3992 					if ((info > 2) && (info <= 4)) {
3993 						if ((countp % 200) == 0) {
3994 							tprime.stop_timer();
3995 							std::cout << countp << "   ";
3996 							MyTime(tprime.user_time());
3997 							std::cout << std::endl;
3998 							tprime.cont_timer();
3999 						}
4000 					}
4001 
4002 					if (info > 4) {
4003 						std::cout << "P = " << P << std::flush;
4004 					}
4005 
4006 					if (use_sieve) {
4007 						if (do_verify > 2)
4008 							vec = represent_over_FB_sieve(P, q);
4009 						else
4010 							represent_over_FB_sieve(P);
4011 					}
4012 					else
4013 						vec = represent_over_FB(P, q);
4014 
4015 					if (info > 4) {
4016 						std::cout << " --- done!" << std::endl;
4017 					}
4018 
4019 					if (do_verify > 2) {
4020 						pvec = vec;
4021 
4022 						// compute representation over original factor base
4023 						for (j = vec.size()-1; j >= numFB; --j) {
4024 							if (!vec[j].is_zero()) {
4025 								nvec = (math_vector< bigint > ) vmat.get_column_vector(j-numFB);
4026 								multiply(nvec, nvec, vec[j]);
4027 								vec[j].assign_zero();
4028 								nvec.set_size(vec.size());
4029 								add(vec, vec, nvec);
4030 
4031 							}
4032 						}
4033 						vec.set_size(numFB);
4034 
4035 
4036 						// compute representation over generators
4037 						fact_base.set_size(numFB);
4038 
4039 						rvec = represent_over_generators(vec);
4040 
4041 						fact_base = new_FB;
4042 
4043 						C.assign_one();
4044 						for (j = 0; j < numCL; ++j) {
4045 							power(A, gens[j], rvec[j]);
4046 							multiply(C, C, A);
4047 						}
4048 
4049 						if (P.is_equivalent(C)) {
4050 							if (info > 2) {
4051 								std::cout << "  " << P << " ~ " << rvec << " ~ " << C << std::endl;
4052 							}
4053 						}
4054 						else {
4055 							OK = false;
4056 							if (info > 2) {
4057 								std::cout << "  ERROR:  " << P << " = " << rvec << " <>" << C;
4058 								std::cout << std::endl;
4059 							}
4060 						}
4061 					}
4062 
4063 					// add verified prime to factor base
4064 					if ((fact_base.size() < (numFB*4)) && (fact_base.size() < 8000)
4065 					    && (num > 28)) {
4066 						new_FB[new_FB.size()] = P;
4067 						fact_base[fact_base.size()] = P;
4068 
4069 						if (do_verify > 2) {
4070 							matrows = pvec.size();
4071 							vmat.resize(matrows, matcols+1);
4072 							vmat.sto_column_vector(pvec, matrows, matcols);
4073 							++matcols;
4074 						}
4075 					}
4076 				}
4077 				p = PList.get_next_prime();
4078 			}
4079 
4080 			fact_base.set_size(numFB);
4081 		}
4082 	}
4083 
4084 	if (do_verify > 0 && is_CL_computed()) {
4085 		generators();
4086 
4087 		// verify orders of generators
4088 		if (info > 1)
4089 			std::cout << "Checking orders of generators:" << std::endl;
4090 
4091 		for (i = 0; i < numCL; ++i) {
4092 			A.assign(gens[i]);
4093 			ord.assign(A.order_in_CL());
4094 			if (ord == CL[i]) {
4095 				if (info > 2) {
4096 					std::cout << "  ord(G[" << i << "]) = " << ord << " = " << CL[i] << std::endl;
4097 				}
4098 			}
4099 			else {
4100 				OK = false;
4101 				if (info > 2) {
4102 					std::cout << "  ERROR:  ord(G[" << i << "]) = " << ord << " <>" << CL[i];
4103 					std::cout << std::endl;
4104 				}
4105 			}
4106 		}
4107 	}
4108 
4109 	if (info) {
4110 		tprime.stop_timer();
4111 		std::cout << "\nelapsed CPU time (verify_class_group) = ";
4112 		MyTime(tprime.user_time());
4113 		std::cout << "\n";
4114 	}
4115 
4116 	bigfloat::set_precision(oprec);
4117 	qi_class::set_current_order(*QO);
4118 
4119 	// FOR TABLE GENERATION
4120 	if (special) {
4121 		t.stop_timer();
4122 		std::cout << t.user_time() << std::endl;
4123 	}
4124 
4125 	return OK;
4126 }
4127 
4128 
4129 
4130 //
4131 // quadratic_order::is_in_lattice(qi_class_real)
4132 //
4133 // Task:
4134 //      determines whether the factorization of A over the factorbase used
4135 //      by the subexponential algorithm is in the lattice of vectors
4136 //      representing the principal class.
4137 //
4138 
4139 bool
is_in_lattice(const qi_class_real & AA)4140 quadratic_order::is_in_lattice(const qi_class_real & AA)
4141 {
4142 	debug_handler("quadratic_order", "is_in_lattice(qi_class_real)");
4143 
4144 	math_vector< bigint > vec;
4145 	bool isin;
4146 	quadratic_number_standard q;
4147 	quadratic_order *QO = qi_class::get_current_order_ptr();
4148 
4149 	// set last-used quadratic_order
4150 	qo_l().set_last(this);
4151 
4152 	qi_class::set_current_order(*this);
4153 
4154 	vec = represent_over_FB_sieve(AA, q);
4155 	isin = is_in_lattice(vec);
4156 
4157 	qi_class::set_current_order(*QO);
4158 
4159 	return isin;
4160 }
4161 
4162 
4163 
4164 //
4165 // quadratic_order::is_in_lattice(math_vector)
4166 //
4167 // Task:
4168 //      determines whether the factorization of A over the factorbase used
4169 //      by the subexponential algorithm is in the lattice of vectors
4170 //      representing the principal class.
4171 //
4172 
4173 bool
is_in_lattice(math_vector<bigint> & vec)4174 quadratic_order::is_in_lattice(math_vector< bigint > & vec)
4175 {
4176 	debug_handler("quadratic_order", "is_in_lattice(math_vector)");
4177 
4178 	// set last-used quadratic_order
4179 	qo_l().set_last(this);
4180 
4181 	math_vector< bigint > x;
4182 	bool is_in_lat;
4183 
4184 	is_in_lat = RHNF.solve_hnf(vec, x);
4185 	return is_in_lat;
4186 }
4187 
4188 
4189 
4190 //
4191 // quadratic_order::is_in_lattice(qi_class_real, bigfloat)
4192 //
4193 // Task:
4194 //      determines whether the factorization of A over the factorbase used
4195 //      by the subexponential algorithm is in the lattice of vectors
4196 //      representing the principal class.  If so, the distance from (1) is
4197 //      also returned.
4198 //
4199 
4200 bool
is_in_lattice(const qi_class_real & AA,bigfloat & dist)4201 quadratic_order::is_in_lattice(const qi_class_real & AA, bigfloat & dist)
4202 {
4203 	debug_handler("quadratic_order", "is_in_lattice(qi_class_real, bigfloat)");
4204 
4205 	lidia_size_t i;
4206 	bigint max, DET2;
4207 	bigfloat ldist, rd2, a, b, rtemp1, rtemp2, maxlog;
4208 	qi_class_real UNIT, C, D;
4209 	math_vector< bigint > vec, x;
4210 	math_vector< bigfloat > logs;
4211 	bool in_lat;
4212 	quadratic_order *QO = qi_class::get_current_order_ptr();
4213 	long oprec = bigfloat::get_precision(), xprec2;
4214 	quadratic_number_standard q;
4215 
4216 	xbigfloat ldist_x;
4217 	timer t;
4218 
4219 	// set last-used quadratic_order
4220 	qo_l().set_last(this);
4221 
4222 	qi_class::set_current_order(*this);
4223 
4224 	bigfloat::set_precision(prec);
4225 
4226 	logs.set_capacity(minima_qn.size());
4227 	logs.set_size(minima_qn.size());
4228 
4229 	// FOR TABLE GENERATION
4230 	if (special2)
4231 		t.start_timer();
4232 
4233 	// find dependency over relation matrix
4234 	vec = represent_over_FB_sieve(AA, q);
4235 
4236 	in_lat = RHNF.solve_hnf(vec, x);
4237 
4238 	if (in_lat) {
4239 		if (info > 3) {
4240 			C.assign_one();
4241 			qi_class_real PP;
4242 			for (i = 0; i < vec.size(); ++i)
4243 				if (!vec[i].is_zero()) {
4244 					power_real(PP, static_cast<qi_class_real>(fact_base[i]), vec[i]);
4245 					C *= PP;
4246 				}
4247 			std::cout << "A = " << AA << std::endl;
4248 			std::cout << "C = " << C << std::endl;
4249 		}
4250 
4251 		DET2 = 1;
4252 		math_vector< bigint > w;
4253 		w.set_mode(EXPAND);
4254 
4255 		// compute vector
4256 		x = TR.multiply(vec);
4257 
4258 		// remove common factors with DET
4259 		if (DET != 1) {
4260 			bigint G = DET;
4261 			for (i = 0; i < x.size(); ++i)
4262 				G = gcd(G, x[i]);
4263 			if (info > 3) {
4264 				std::cout << "DET = " << DET << std::endl;
4265 				std::cout << "G = " << G << std::endl;
4266 				std::cout << "DET/G = " << DET/G << std::endl;
4267 			}
4268 			if (DET.is_lt_zero())
4269 				G.negate();
4270 			if (abs(G) > 1)
4271 				divide(x, x, G);
4272 			divide(DET2, DET, G);
4273 		}
4274 
4275 		matrix< bigint > B;
4276 		B.set_storage_mode(matrix_flags::sparse_representation);
4277 		B.resize(AMAT.get_no_of_rows(), AMAT.get_no_of_columns());
4278 		for (lidia_size_t ii1 = 0; ii1 < AMAT.get_no_of_rows(); ++ii1)
4279 			for (lidia_size_t ii2 = 0; ii2 < AMAT.get_no_of_columns(); ++ii2)
4280 				if (AMAT.member(ii1, ii2))
4281 					B.sto(ii1, ii2, bigint(AMAT.member(ii1, ii2)));
4282 
4283 		multiply(w, B, x);
4284 
4285 		if (info > 3) {
4286 			std::cout << "DET2 = " << DET2 << std::endl;
4287 			std::cout << "Ax = " << w << std::endl;
4288 			std::cout << "vec = " << vec << std::endl;
4289 			if (DET2 != 1)
4290 				std::cout << "Ax/DET2 = " << w / DET2 << std::endl;
4291 		}
4292 
4293 		if ((w/DET2) != vec) {
4294 			subtract(w, w, vec);
4295 			if (info > 3)
4296 				std::cout << "Ax - vec = " << w << std::endl;
4297 
4298 			if (!RHNF.solve_hnf(w, w)) {
4299 				std::cout << "ERROR!  Couldn't solve Hx = w!!!" << std::endl;
4300 				std::cout << "w = " << w << std::endl;
4301 			}
4302 
4303 			w = TR.multiply(w);
4304 			subtract(x, x, w);
4305 
4306 			multiply(w, B, x);
4307 			if (info > 3) {
4308 				std::cout << "Ax = " << w << std::endl;
4309 				std::cout << "vec = " << vec << std::endl;
4310 				if (DET2 != 1)
4311 					std::cout << "Ax/DET2 = " << w / DET2 << std::endl;
4312 			}
4313 		}
4314 
4315 		// FOR TABLE GENERATION
4316 		if (special2) {
4317 			t.stop_timer();
4318 			std::cout << t.user_time() << std::endl;
4319 			t.start_timer();
4320 		}
4321 
4322 
4323 		// Markus' Method
4324 		bigint q;
4325 		base_vector< xbigfloat > lbw, lbr;
4326 		base_vector< long > plbw, plbr;
4327 		long br, lowregbound, e_tmp;
4328 		xbigfloat Rx, ldist_x, LL;
4329 		math_vector< bigint > fud;
4330 		bigfloat tmp;
4331 
4332 		if (info > 3) {
4333 			max = x[0];
4334 			for (i = 1; i < x.size(); ++i)
4335 				if (x[i] > max)
4336 					max = x[i];
4337 			std::cout << "MAX in x:  " << decimal_length(max) << " decimal digits" << std::endl;
4338 		}
4339 
4340 
4341 		//<MM>
4342 		//
4343 		// NOTE: The reduction modulo the regulator is a hack that must be
4344 		//       removed by using a not yet existing function
4345 		//       quadratic_number_power_product::short_principal_ideal_generator.
4346 		//
4347 		//       !! HERE, we assume that the basis of "fund_unit" is "minima_qn" !!
4348 		//
4349 		//</MM>
4350 
4351 //    const math_vector<bigint> & fund_unit_exp = fund_unit.get_exponents();
4352 		math_vector< bigint > fund_unit_exp = fund_unit.get_exponents();
4353 		fund_unit_exp.set_size(minima_qn.size());
4354 
4355 		if (DET2 == 1) {
4356 			Rx.assign(R_xbig * fu_mult);
4357 			br = Rx.b_value();
4358 
4359 			// compute distance modulo R
4360 			reduce_modulo_regulator(q, minima_qn, x, lbw, plbw, minima_qn, fund_unit_exp, lbr, plbr, br);
4361 
4362 			// compute new vector
4363 			if (info > 3) {
4364 				std::cout << "q = " << decimal_length(q) << " decimal digits" << std::endl;
4365 			}
4366 			x = x - q*fund_unit_exp;
4367 		}
4368 		else {
4369 			Rx.assign(R_xbig * fu_mult * DET2);
4370 			br = Rx.b_value();
4371 
4372 			fud = DET2*fund_unit_exp;
4373 
4374 			// compute distance modulo R
4375 			reduce_modulo_regulator(q, minima_qn, x, lbw, plbw, minima_qn, fud, lbr, plbr, br);
4376 
4377 			// compute new vector
4378 			if (info > 3) {
4379 				std::cout << "q = " << decimal_length(q) << " decimal digits" << std::endl;
4380 			}
4381 			x = x - q*fud;
4382 		}
4383 
4384 		if (info > 3) {
4385 			max = x[0];
4386 			for (i = 1; i < x.size(); ++i)
4387 				if (x[i] > max)
4388 					max = x[i];
4389 			std::cout << "MAX in x:  " << decimal_length(max) << " decimal digits" << std::endl;
4390 		}
4391 
4392 
4393 		lowregbound = lower_regulator_bound(LL, Delta, h, 0);
4394 
4395 		if (DET2 == 1) {
4396 			xprec2 = static_cast<long>(std::ceil(static_cast<double>(bigfloat::get_precision()) *
4397 							     std::log(10.0L) / std::log(2.0L)));
4398 
4399 			relative_Ln_approximation (ldist_x, minima_qn, x, xprec2, lowregbound);
4400 		}
4401 		else {
4402 			max = abs(x[0]);
4403 			for (i = 1; i < x.size(); ++i)
4404 				if (abs(x[i]) > max)
4405 					max.assign(abs(x[i]));
4406 
4407 			xprec2 = static_cast<long>(std::ceil(static_cast<double>(4*decimal_length(max) +
4408 										 5 + decimal_length(DET2)) *
4409 							     std::log(10.0L) / std::log(2.0L)));
4410 
4411 			relative_Ln_approximation (ldist_x, minima_qn, x, xprec2, lowregbound);
4412 
4413 			std::cout << ldist_x << std::endl;
4414 			ldist.assign(ldist_x.get_mantissa());
4415 			e_tmp = ldist_x.get_exponent()-ldist_x.get_mantissa().bit_length();
4416 			if (e_tmp > 0)
4417 				shift_left(ldist, ldist, e_tmp);
4418 			else
4419 				shift_right(ldist, ldist, -e_tmp);
4420 
4421 			std::cout << "ldist = " << ldist << std::endl;
4422 
4423 			divide(ldist_x, ldist_x, xbigfloat(DET2), xprec2);
4424 
4425 			std::cout << ldist_x << std::endl;
4426 		}
4427 
4428 
4429 		// store in bigfloat
4430 		ldist.assign(ldist_x.get_mantissa());
4431 		e_tmp = ldist_x.get_exponent()-ldist_x.get_mantissa().bit_length();
4432 		if (e_tmp > 0)
4433 			shift_left(ldist, ldist, e_tmp);
4434 		else
4435 			shift_right(ldist, ldist, -e_tmp);
4436 
4437 		if (info > 3) {
4438 			std::cout << "ldist = " << ldist << std::endl;
4439 		}
4440 
4441 		// final reduction mod R
4442 		while (ldist < 0)
4443 			add(ldist, ldist, R);
4444 		while (ldist > R)
4445 			subtract(ldist, ldist, R);
4446 
4447 		if (info > 3)
4448 			std::cout << "ldist = " << ldist << std::endl;
4449 
4450 		// FOR TABLE GENERATION
4451 		if (special2) {
4452 			t.stop_timer();
4453 			std::cout << t.user_time() << std::endl;
4454 		}
4455 
4456 		// convert to standard formulation
4457 		UNIT.assign_one();
4458 		dist = AA.convert_distance(UNIT, ldist);
4459 	}
4460 	else {
4461 		// FOR TABLE GENERATION
4462 		if (special2) {
4463 			t.stop_timer();
4464 			std::cout << t.user_time() << std::endl;
4465 			std::cout << "0" << std::endl;
4466 		}
4467 
4468 		dist.assign(R);
4469 	}
4470 
4471 	qi_class::set_current_order(*QO);
4472 	bigfloat::set_precision(oprec);
4473 
4474 	return in_lat;
4475 }
4476 
4477 
4478 
4479 //
4480 // quadratic_order::factor_h
4481 //
4482 // Task:
4483 //      computes the class number and factors it.
4484 //
4485 
4486 rational_factorization
factor_h()4487 quadratic_order::factor_h()
4488 {
4489 	debug_handler("quadratic_order", "factor_h");
4490 
4491 	// set last-used quadratic_order
4492 	qo_l().set_last(this);
4493 
4494 	// compute class number, if needed
4495 	if (!is_h_computed())
4496 		class_number();
4497 
4498 	// factor h, if necessary
4499 	if (h.is_gt_zero()) {
4500 		if (hfact != h) {
4501 			hfact.assign(h);
4502 			hfact.factor();
4503 		}
4504 	}
4505 
4506 	return hfact;
4507 }
4508 
4509 
4510 
4511 //
4512 // quadratic_order::exponent
4513 //
4514 // Task:
4515 //      computes the exponent of the class group.
4516 //
4517 
4518 bigint
exponent()4519 quadratic_order::exponent()
4520 {
4521 	debug_handler("quadratic_order", "exponent");
4522 
4523 	bigint expo;
4524 
4525 	// set last-used quadratic_order
4526 	qo_l().set_last(this);
4527 
4528 	// compute class group, if needed
4529 	if (!is_CL_computed())
4530 		class_group();
4531 
4532 	// return largest non-cyclic factor
4533 	if (!is_CL_computed())
4534 		expo.assign_zero();
4535 	else
4536 		expo.assign(CL[CL.size()-1]);
4537 
4538 	return expo;
4539 }
4540 
4541 
4542 
4543 //
4544 // quadratic_order::p_rank
4545 //
4546 // Task:
4547 //      computes the p_rank of the class group.
4548 //
4549 
4550 int
p_rank(const bigint & p)4551 quadratic_order::p_rank(const bigint & p)
4552 {
4553 	debug_handler("quadratic_order", "p_rank");
4554 
4555 	int prank;
4556 	lidia_size_t i, rank;
4557 	bigint rem;
4558 
4559 	// set last-used quadratic_order
4560 	qo_l().set_last(this);
4561 
4562 	// compute class group, if needed
4563 	if (!is_CL_computed())
4564 		class_group();
4565 
4566 	if (h.is_zero())
4567 		prank = 0;
4568 	else {
4569 		// compute p-rank
4570 		rank = CL.size();
4571 		i = 0;
4572 		remainder(rem, CL[i], p);
4573 		while ((!rem.is_zero()) && (i < rank)) {
4574 			++i;
4575 			if (i < rank)
4576 				remainder(rem, CL[i], p);
4577 		}
4578 		prank = rank - i;
4579 	}
4580 
4581 	return prank;
4582 }
4583 
4584 
4585 
4586 //
4587 // quadratic_order::generators
4588 //
4589 // Task:
4590 //      computes a set of generators of the class group.
4591 //
4592 
4593 base_vector< qi_class >
generators()4594 quadratic_order::generators()
4595 {
4596 	debug_handler("quadratic_order", "generators");
4597 
4598 	qi_class G, A;
4599 	lidia_size_t i, j, numCL, numFB;
4600 	quadratic_order *QO = qi_class::get_current_order_ptr();
4601 
4602 	// set last-used quadratic_order
4603 	qo_l().set_last(this);
4604 
4605 	qi_class::set_current_order(*this);
4606 
4607 	// compute class group, if needed
4608 	if (!is_CL_computed())
4609 		class_group();
4610 
4611 	if ((gens.size() == 0) && (!Delta.is_zero())) {
4612 		if (h.is_one()) {
4613 			G.assign_one();
4614 			gens[0] = G;
4615 		}
4616 		else if (fact_base.size() == 1)
4617 			gens[0].assign(fact_base[0]);
4618 		else {
4619 			// reduce transformation matrix elements and compute inverse
4620 			set_transformations();
4621 
4622 			// columns of UINV are the exponents vectors for the factor base elements
4623 			// which correspond to each generator
4624 
4625 			numCL = CL.size();
4626 			numFB = contributors.size();
4627 			for (i = 0; i < numCL; ++i) {
4628 				G.assign_one();
4629 				for (j = 0; j < numFB; ++j) {
4630 					power(A, fact_base[contributors[j]], UINV.member(j, i));
4631 					G *= A;
4632 				}
4633 				gens[i] = G;
4634 			}
4635 		}
4636 	}
4637 
4638 	qi_class::set_current_order(*QO);
4639 
4640 	return gens;
4641 }
4642 
4643 
4644 
4645 //
4646 // quadratic_order::factor_discriminant
4647 //
4648 // Task:
4649 //      factors the discriminant by computing the class group and utiliziting
4650 //      its 2-Sylow subgroup.
4651 //
4652 
4653 rational_factorization
factor_discriminant()4654 quadratic_order::factor_discriminant()
4655 {
4656 	debug_handler("quadratic_order", "factor_discriminant");
4657 
4658 	timer t;
4659 	quadratic_order *QO = qi_class::get_current_order_ptr();
4660 	long oprec = bigfloat::get_precision();
4661 
4662 	bigfloat::set_precision(prec);
4663 	qi_class::set_current_order(*this);
4664 
4665 	// set last-used quadratic_order
4666 	qo_l().set_last(this);
4667 
4668 	if (info)
4669 		t.start_timer();
4670 
4671 	if (!Delta.is_zero()) {
4672 		if (is_imaginary())
4673 			factor_imag();
4674 		else
4675 			factor_real();
4676 	}
4677 
4678 	if (info) {
4679 		t.stop_timer();
4680 		std::cout << "\nelapsed CPU time (factor) = ";
4681 		MyTime(t.user_time());
4682 		std::cout << "\n";
4683 	}
4684 
4685 	bigfloat::set_precision(oprec);
4686 	qi_class::set_current_order(*QO);
4687 
4688 	return disc_fact;
4689 }
4690 
4691 
4692 
4693 //<MM>
4694 
4695 //
4696 // ::fundamental_unit
4697 //
4698 // Returns the fundamental unit of the order.
4699 //
4700 //
4701 
4702 const quadratic_number_power_product & quadratic_order
fundamental_unit()4703 ::fundamental_unit ()
4704 {
4705 	lidia_error_handler("quadratic_order::fundamental_unit()",
4706 			    "Not yet implemented.");
4707 
4708 	if (!is_fundamental_unit_computed()) {
4709 
4710 		// fundamental unit has not been computed yet - compute it
4711 
4712 		if (is_imaginary())
4713 			fund_unit.assign_one(*this);
4714 
4715 		else {
4716 
4717 			int num;
4718 			quadratic_order *QO = qi_class::get_current_order_ptr();
4719 			long oprec = bigfloat::get_precision();
4720 
4721 			bigfloat::set_precision(prec);
4722 			qi_class::set_current_order(*this);
4723 
4724 			// set last-used quadratic_order
4725 			qo_l().set_last(this);
4726 
4727 			num = decimal_length(Delta);
4728 			if (num < RBJTB)
4729 				regulator_BJT();
4730 			else if (num < RSHANKSB)
4731 				regulator_shanks();
4732 			else
4733 				class_group_siqs();
4734 
4735 			bigfloat::set_precision(oprec);
4736 			qi_class::set_current_order(*QO);
4737 
4738 			// The result of class_group_siqs should be
4739 			// a compact representation of the fundamental
4740 			// unit or an absolute 4 approximation to the
4741 			// regulator. Same applies for the regulator
4742 			// functions.
4743 			//
4744 			if (!is_fundamental_unit_computed()) {
4745 				fund_unit.assign(*this, R, 1);
4746 			}
4747 		}
4748 	}
4749 
4750 	return fund_unit;
4751 }
4752 
4753 
4754 
4755 //
4756 // ::could_be_regulator_multiple
4757 //
4758 // If the function returns false, then l is not an absolute
4759 // 3 approximation to the regulator.
4760 //
4761 bool quadratic_order
could_be_regulator_multiple(const xbigfloat & l)4762 ::could_be_regulator_multiple (const xbigfloat & l)
4763 {
4764 	debug_handler ("quadratic_order",
4765 		       "could_be_regulator_multiple()");
4766 
4767 	// Determine the accuracy, 2^{-k+2} < ln2.
4768 	//
4769 	long k = 3;
4770 
4771 	// Determine minimum in O close to l
4772 	//
4773 	quadratic_ideal J, K;
4774 	quadratic_number_power_product beta;
4775 	xbigfloat b;
4776 
4777 	K.assign(*this);
4778 	J.assign(*this);
4779 	J.order_close(beta, b, l, k+1);
4780 
4781 	if (info > 1)
4782 		std::cout << "  nearest from (1) = " << J << std::endl;
4783 
4784 	// Found O again ?
4785 	//
4786 	if (J != K) {
4787 		quadratic_ideal A1, A2;
4788 		A1 = J; A1.rho();
4789 		A2 = J; A2.inverse_rho();
4790 		if (K != A1 && K != A2)
4791 			return false;
4792 	}
4793 	return true;
4794 }
4795 
4796 
4797 
4798 //</MM>
4799 
4800 
4801 
4802 
4803 //
4804 // operator >>
4805 //
4806 // Task:
4807 //      inputs a quadratic_order from the std::istream in.
4808 //
4809 
operator >>(std::istream & in,quadratic_order & QO)4810 std::istream & operator >> (std::istream & in, quadratic_order & QO)
4811 {
4812 	debug_handler("quadratic_order", "operator >>");
4813 
4814 	bigint D;
4815 
4816 	in >> D;
4817 	QO.assign(D);
4818 
4819 	return in;
4820 }
4821 
4822 
4823 
4824 //
4825 // operator <<
4826 //
4827 // Task:
4828 //      outputs a quadratic_order to the std::ostream out.
4829 //
4830 
operator <<(std::ostream & out,const quadratic_order & QO)4831 std::ostream & operator << (std::ostream & out, const quadratic_order & QO)
4832 {
4833 	debug_handler("quadratic_order", "operator << ");
4834 
4835 	int num;
4836 	long oprec = bigfloat::get_precision();
4837 
4838 	num = decimal_length(QO.Delta);
4839 	bigfloat::set_precision(QO.prec);
4840 
4841 	out << "Quadratic Order:" << std::endl;
4842 	out << "   Delta = " << QO.Delta << " (" << num << ")" << std::endl;
4843 	if (!QO.Delta.is_zero()) {
4844 		if (QO.disc_fact.no_of_comp() > 0)
4845 			out << " = " << QO.disc_fact << std::endl;
4846 		if (QO.is_R_computed())
4847 			out << "   R = " << QO.R << std::endl;
4848 		if (QO.is_h_computed())
4849 			out << "   h = " << QO.h << std::endl;
4850 		if (QO.is_CL_computed())
4851 			out << "   CL = " << QO.CL << std::endl;
4852 		if (QO.gens.size() > 0)
4853 			out << "   generators = " << QO.gens << std::endl;
4854 		if (QO.is_L_computed())
4855 			out << "   L(1, X) = " << QO.L << std::endl;
4856 		if (QO.is_C_computed()) {
4857 			bigfloat::set_precision(8);
4858 			out << "   C(Delta) = " << QO.Cfunc << std::endl;
4859 		}
4860 	}
4861 
4862 	bigfloat::set_precision(oprec);
4863 
4864 	return out;
4865 }
4866 
4867 
4868 
4869 #ifdef LIDIA_NAMESPACE
4870 }	// end of namespace LiDIA
4871 #endif
4872