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 : Andreas Mueller (AM), Volker Mueller (VM)
14 // Changes : See CVS log
15 //
16 //==============================================================================================
17
18
19 #ifdef HAVE_CONFIG_H
20 # include "config.h"
21 #endif
22 #include "LiDIA/rational_factorization.h"
23
24
25
26 #ifdef LIDIA_NAMESPACE
27 namespace LiDIA {
28 #endif
29
30
31
32 void rational_factorization::
trialdiv(lidia_size_t index,unsigned int lower_bound,unsigned int upper_bound,ecm_primes & prim)33 trialdiv(lidia_size_t index, unsigned int lower_bound,
34 unsigned int upper_bound, ecm_primes &prim)
35 {
36 bigint N(factors[index].single_base);
37 int lexp = factors[index].single_exponent;
38
39 rf_single_factor fact;
40
41 lidia_size_t n = no_of_comp();
42 long p, q;
43 int count = 0;
44
45 bigint TMP;
46
47 if (info) {
48 std::cout << "index " << index << ": trialdivision up to " << upper_bound << "\n";
49 std::cout.flush();
50 }
51
52 prim.resetprimes(lower_bound);
53
54 if (lower_bound < 2) {
55 while (N.is_even()) {
56 count++;
57 N.divide_by_2();
58 }
59
60 if (count) {
61 fact.single_base = bigint(2);
62 fact.single_exponent = count * lexp;
63 fact.factor_state = rf_single_factor::prime;
64 factors[n++] = fact;
65
66 if (info)
67 if (count > 1)
68 std::cout << "factor: 2 ^ " << count << "\n";
69 else
70 std::cout << "factor: 2 \n";
71
72 count = 0;
73 }
74 }
75
76 while ((static_cast<unsigned int>(p = prim.getprimes()) < upper_bound) &&
77 !N.is_one() && p != 1) {
78 div_rem(TMP, q, N, p);
79
80 while (!q) {
81 N.assign(TMP);
82 count++;
83 div_rem(TMP, q, N, p);
84 }
85
86 if (count) {
87 fact.single_base.assign(p);
88 fact.single_exponent = count * lexp;
89 fact.factor_state = rf_single_factor::prime;
90
91 factors[n++] = fact;
92
93 if (info)
94 if (count > 1)
95 std::cout << "factor: " << p << " ^ " << count << "\n";
96 else
97 std::cout << "factor: " << p << " \n";
98
99 count = 0;
100 }
101 }
102
103 if (N.is_one()) {
104 n--;
105 factors[index] = factors[n];
106 factors.remove_from(n);
107 }
108 else {
109 fact.single_base.assign(N);
110 fact.single_exponent = lexp;
111 fact.factor_state = rf_single_factor::dont_know;
112
113 factors[index] = fact;
114 }
115 }
116
117
118
119 void rational_factorization::
trialdiv_comp(lidia_size_t index,unsigned int upper_bound,unsigned int lower_bound)120 trialdiv_comp(lidia_size_t index, unsigned int upper_bound,
121 unsigned int lower_bound)
122 {
123 if ((index< 0) || (index >= no_of_comp())) {
124 lidia_error_handler("rational_factorization", "trialdiv_comp::index out of range");
125 return;
126 }
127
128 if (rf_single_factor::state(factors[index].factor_state) == rf_single_factor::prime) {
129 if (info)
130 std::cout << "prime number " << factors[index].single_base << "\n";
131 return;
132 }
133
134 if (lower_bound > upper_bound) {
135 lidia_warning_handler("rational_factorization", "trialdiv_comp::lower_bound > upper_bound");
136 return;
137 }
138
139 ecm_primes prim(lower_bound, upper_bound+200, 200000);
140
141 trialdiv(index, lower_bound, upper_bound, prim);
142
143 compose();
144 }
145
146
147
148 void rational_factorization::
trialdiv(unsigned int upper_bound,unsigned int lower_bound)149 trialdiv(unsigned int upper_bound, unsigned int lower_bound)
150 {
151 if (lower_bound > upper_bound) {
152 lidia_warning_handler("rational_factorization", "trialdiv::lower_bound > upper_bound");
153 return;
154 }
155
156 ecm_primes prim(lower_bound, upper_bound+200, 200000);
157
158 lidia_size_t index, len = no_of_comp();
159
160 for (index = 0; index < len; index++) {
161 if (rf_single_factor::state(factors[index].factor_state) == rf_single_factor::prime) {
162 if (info)
163 std::cout << "index " << index << " prime number " << factors[index].single_base << "\n";
164 }
165 else
166 trialdiv(index, lower_bound, upper_bound, prim);
167 }
168
169 compose();
170 }
171
172
173
trialdiv(const bigint & N,unsigned int upper_bound,unsigned int lower_bound)174 rational_factorization trialdiv(const bigint &N, unsigned int upper_bound,
175 unsigned int lower_bound)
176 {
177 rational_factorization f(N);
178
179 f.trialdiv(upper_bound, lower_bound);
180
181 return f;
182 }
183
184
185
186 #ifdef LIDIA_NAMESPACE
187 } // end of namespace LiDIA
188 #endif
189