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