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 : Thomas Pfahler (TPf)
14 // Changes : See CVS log
15 //
16 //==============================================================================================
17
18
19 #ifdef HAVE_CONFIG_H
20 # include "config.h"
21 #endif
22 #include "LiDIA/Fp_polynomial.h"
23
24
25
26 #ifdef LIDIA_NAMESPACE
27 namespace LiDIA {
28 #endif
29
30
31
32 crossover_class Fp_polynomial::crossovers;
33
34
spline(const int x[],const int y[],int n,double yp1,double ypn,double y2[])35 void spline(const int x[], const int y[], int n, double yp1, double ypn,
36 double y2[])
37 //see Numerical Recipes in C, ISBN 0-521-43108-5, p. 113ff
38 {
39 int i, k;
40 double p, qn, sig, un;
41 double *u = new double[n];
42 memory_handler(u, "Fp_polynomial", "spline::Error in memory allocation");
43
44 if (yp1 > 0.99e30)
45 y2[0] = u[0] = 0.0;
46 else {
47 y2[0] = -0.5;
48 u[0] = (3.0/(x[1]-x[0])) * static_cast<double>((y[1]-y[0]) / static_cast<double>(x[1]-x[0]) - yp1);
49 }
50
51 for (i = 1; i < n-1; i++) {
52 sig = static_cast<double>(x[i]-x[i-1]) / static_cast<double>(x[i+1]-x[i-1]);
53 p = sig*y2[i-1]+2.0;
54 y2[i] = (sig-1.0)/p;
55 u[i] = static_cast<double>(y[i+1]-y[i]) / static_cast<double>(x[i+1]-x[i])
56 - static_cast<double>(y[i]-y[i-1]) / static_cast<double>(x[i]-x[i-1]);
57 u[i] = (6.0 * u[i] / static_cast<double>(x[i+1]-x[i-1]) - sig * u[i-1]) / p;
58 }
59 if (ypn > 0.99e30)
60 qn = un = 0.0;
61 else {
62 qn = 0.5;
63 un = (3.0 / static_cast<double>(x[n-1]-x[n-2]))
64 * (ypn - static_cast<double>(y[n-1]-y[n-2]) / static_cast<double>(x[n-1]-x[n-2]));
65 }
66 y2[n-1] = (un-qn*u[n-2])/(qn*y2[n-2]+1.0);
67 for (k = n-2; k >= 0; k--)
68 y2[k] = y2[k]*y2[k+1]+u[k];
69 delete[] u;
70 }
71
72
73
splint(const int xa[],const int ya[],const double y2a[],int n,double x)74 double splint(const int xa[], const int ya[], const double y2a[], int n,
75 double x)
76 //see Numerical Recipes in C, ISBN 0-521-43108-5, p. 113ff
77 {
78 int klo, khi, k;
79 double h, b, a, y;
80
81 klo = 0;
82 khi = n-1;
83 while (khi-klo > 1) {
84 k = (khi+klo) >> 1;
85 if (xa[k] > x) khi = k;
86 else klo = k;
87 }
88 h = xa[khi]-xa[klo];
89 if (h == 0.0) {
90 lidia_error_handler("crossover_class",
91 "splint(...)::identical x-values for interpolation)");
92 return 0.0;
93 }
94 a = static_cast<double>(xa[khi]-x) / h;
95 b = static_cast<double>(x-xa[klo]) / h;
96 y = a*ya[klo] + b*ya[khi] +
97 ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h) / 6.0;
98 return y;
99 }
100
101
102
crossover_class()103 crossover_class::crossover_class()
104 {
105 //
106 #include "LiDIA/Fp_pol_crossover.h"
107 //
108 init(x_val, fftmul_val, fftdiv_val, fftrem_val, inv_val, gcd_val, xgcd_val);
109 halfgcd = halfgcd_val;
110 log2_newton = log2_newton_val;
111 }
112
113
114
init(const int x_val[CROV_NUM_VALUES],const int fftmul_val[CROV_NUM_VALUES],const int fftdiv_val[CROV_NUM_VALUES],const int fftrem_val[CROV_NUM_VALUES],const int inv_val[CROV_NUM_VALUES],const int gcd_val[CROV_NUM_VALUES],const int xgcd_val[CROV_NUM_VALUES])115 void crossover_class::init(const int x_val[CROV_NUM_VALUES],
116 const int fftmul_val[CROV_NUM_VALUES],
117 const int fftdiv_val[CROV_NUM_VALUES],
118 const int fftrem_val[CROV_NUM_VALUES],
119 const int inv_val[CROV_NUM_VALUES],
120 const int gcd_val[CROV_NUM_VALUES],
121 const int xgcd_val[CROV_NUM_VALUES])
122 {
123 int i;
124 for (i = 0; i < CROV_NUM_VALUES; i++) {
125 x[i] = x_val[i];
126 y_fftmul[i] = fftmul_val[i];
127 y_fftdiv[i] = fftdiv_val[i];
128 y_fftrem[i] = fftrem_val[i];
129 y_inv[i] = inv_val[i];
130 y_gcd[i] = gcd_val[i];
131 }
132 spline(x, y_fftmul, CROV_NUM_VALUES, 0.0, 0.0, y2_fftmul);
133 spline(x, y_fftdiv, CROV_NUM_VALUES, 0.0, 0.0, y2_fftdiv);
134 spline(x, y_fftrem, CROV_NUM_VALUES, 0.0, 0.0, y2_fftrem);
135 spline(x, y_inv, CROV_NUM_VALUES, 0.0, 0.0, y2_inv);
136 spline(x, y_gcd, CROV_NUM_VALUES, 0.0, 0.0, y2_gcd);
137 spline(x, y_xgcd, CROV_NUM_VALUES, 0.0, 0.0, y2_xgcd);
138 }
139
140
141
fftmul_crossover(const bigint & modulus) const142 int crossover_class::fftmul_crossover(const bigint &modulus) const
143 {
144 int bl = modulus.bit_length(), r;
145 if (bl >= x[CROV_NUM_VALUES - 1])
146 r = y_fftmul[CROV_NUM_VALUES - 1];
147 else
148 r = int(splint(x, y_fftmul, y2_fftmul, CROV_NUM_VALUES, bl));
149 return comparator< int >::max(r, 1);
150 }
151
152
153
fftdiv_crossover(const bigint & modulus) const154 int crossover_class::fftdiv_crossover(const bigint &modulus) const
155 {
156 int bl = modulus.bit_length(), r;
157 if (bl >= x[CROV_NUM_VALUES - 1])
158 r = y_fftdiv[CROV_NUM_VALUES - 1];
159 else
160 r = int(splint(x, y_fftdiv, y2_fftdiv, CROV_NUM_VALUES, bl));
161 return comparator< int >::max(r, 1);
162 }
163
164
165
fftrem_crossover(const bigint & modulus) const166 int crossover_class::fftrem_crossover(const bigint &modulus) const
167 {
168 int bl = modulus.bit_length(), r;
169 if (bl >= x[CROV_NUM_VALUES - 1])
170 r = y_fftrem[CROV_NUM_VALUES - 1];
171 else
172 r = int(splint(x, y_fftrem, y2_fftrem, CROV_NUM_VALUES, bl));
173 return comparator< int >::max(r, 1);
174 }
175
176
177
inv_crossover(const bigint & modulus) const178 int crossover_class::inv_crossover(const bigint &modulus) const
179 {
180 int bl = modulus.bit_length(), r;
181 if (bl >= x[CROV_NUM_VALUES - 1])
182 r = y_inv[CROV_NUM_VALUES - 1];
183 else
184 r = int(splint(x, y_inv, y2_inv, CROV_NUM_VALUES, bl));
185 return comparator< int >::max(r, 1);
186 }
187
188
189
halfgcd_crossover(const bigint & modulus) const190 int crossover_class::halfgcd_crossover(const bigint &modulus) const
191 {
192 (void)modulus;
193 return halfgcd;
194 }
195
196
197
gcd_crossover(const bigint & modulus) const198 int crossover_class::gcd_crossover(const bigint &modulus) const
199 {
200 int bl = modulus.bit_length(), r;
201 if (bl >= x[CROV_NUM_VALUES - 1])
202 r = y_gcd[CROV_NUM_VALUES - 1];
203 else
204 r = int(splint(x, y_gcd, y2_gcd, CROV_NUM_VALUES, bl));
205 return comparator< int >::max(r, 1);
206 }
207
208
209
xgcd_crossover(const bigint & modulus) const210 int crossover_class::xgcd_crossover(const bigint &modulus) const
211 {
212 int bl = modulus.bit_length(), r;
213 if (bl >= x[CROV_NUM_VALUES - 1])
214 r = y_xgcd[CROV_NUM_VALUES - 1];
215 else
216 r = int(splint(x, y_xgcd, y2_xgcd, CROV_NUM_VALUES, bl));
217 return comparator< int >::max(r, 1);
218 }
219
220
221
log2_newton_crossover(const bigint & modulus) const222 int crossover_class::log2_newton_crossover(const bigint &modulus) const
223 {
224 (void)modulus;
225 return log2_newton;
226 }
227
228
229
230 #ifdef LIDIA_NAMESPACE
231 } // end of namespace LiDIA
232 #endif
233