Line data Source code
1 : /* Auxiliary functions to evaluate Lucas sequences.
2 :
3 : Copyright 2002, 2003, 2005, 2006, 2008, 2011, 2012, 2015
4 : Paul Zimmermann, Alexander Kruppa, Dave Newman.
5 :
6 : This file is part of the ECM Library.
7 :
8 : The ECM Library is free software; you can redistribute it and/or modify
9 : it under the terms of the GNU Lesser General Public License as published by
10 : the Free Software Foundation; either version 3 of the License, or (at your
11 : option) any later version.
12 :
13 : The ECM Library is distributed in the hope that it will be useful, but
14 : WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 : or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 : License for more details.
17 :
18 : You should have received a copy of the GNU Lesser General Public License
19 : along with the ECM Library; see the file COPYING.LIB. If not, see
20 : http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 : 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22 :
23 : /* References:
24 :
25 : A p+1 Method of Factoring, H. C. Williams, Mathematics of Computation,
26 : volume 39, number 159, pages 225-234, 1982.
27 :
28 : Evaluating recurrences of form X_{m+n} = f(X_m, X_n, X_{m-n}) via
29 : Lucas chains, Peter L. Montgomery, December 1983, revised January 1992. */
30 :
31 : #include "ecm-impl.h"
32 :
33 : /* P <- V_2(Q) */
34 : static void
35 534411622 : pp1_duplicate (mpres_t P, mpres_t Q, mpmod_t n)
36 : {
37 534411622 : mpres_sqr (P, Q, n);
38 534411622 : mpres_sub_ui (P, P, 2, n);
39 534411622 : }
40 :
41 : /* P <- V_{m+n} where Q = V_m, R = V_n, S = V_{m-n}.
42 : t is an auxiliary variable.
43 : Warning: P may equal Q, R or S.
44 : */
45 : static void
46 4338659693 : pp1_add3 (mpres_t P, mpres_t Q, mpres_t R, mpres_t S, mpmod_t n, mpres_t t)
47 : {
48 4338659693 : mpres_mul (t, Q, R, n);
49 4338659693 : mpres_sub (P, t, S, n);
50 4338659693 : }
51 :
52 : /* computes V_k(P) from P=A and puts the result in P=A. Assumes k>2.
53 : Uses auxiliary variables t, B, C, T, T2.
54 : */
55 : void
56 115424819 : pp1_mul_prac (mpres_t A, ecm_uint k, mpmod_t n, mpres_t t, mpres_t B,
57 : mpres_t C, mpres_t T, mpres_t T2)
58 : {
59 : ecm_uint d, e, r;
60 : static double val = 0.61803398874989485; /* 1/(golden ratio) */
61 :
62 : /* Note: we used to use several (4) values of "val", but:
63 : (1) the code to estimate the best value was buggy;
64 : (2) even after fixing the bug, the overhead to choose the
65 : best value was larger than the corresponding gain (for a c155
66 : and B1=10^7). */
67 :
68 115424819 : d = k;
69 115424819 : r = (ecm_uint) ((double) d * val + 0.5);
70 :
71 : /* first iteration always begins by Condition 3, then a swap */
72 115424819 : d = k - r;
73 115424819 : e = 2 * r - k;
74 115424819 : mpres_set (B, A, n); /* B=A */
75 115424819 : mpres_set (C, A, n); /* C=A */
76 115424819 : pp1_duplicate (A, A, n); /* A = 2*A */
77 4035458155 : while (d != e)
78 : {
79 3920033336 : if (d < e)
80 : {
81 2780657906 : r = d;
82 2780657906 : d = e;
83 2780657906 : e = r;
84 2780657906 : mpres_swap (A, B, n);
85 : }
86 : /* do the first line of Table 4 whose condition qualifies */
87 3920033336 : if (d - e <= e / 4 && ((d + e) % 3) == 0)
88 : { /* condition 1 */
89 84096583 : d = (2 * d - e) / 3;
90 84096583 : e = (e - d) / 2;
91 84096583 : pp1_add3 (T, A, B, C, n, t); /* T = f(A,B,C) */
92 84096583 : pp1_add3 (T2, T, A, B, n, t); /* T2 = f(T,A,B) */
93 84096583 : pp1_add3 (B, B, T, A, n, t); /* B = f(B,T,A) */
94 84096583 : mpres_swap (A, T2, n); /* swap A and T2 */
95 : }
96 3835936753 : else if (d - e <= e / 4 && (d - e) % 6 == 0)
97 : { /* condition 2 */
98 17664628 : d = (d - e) / 2;
99 17664628 : pp1_add3 (B, A, B, C, n, t); /* B = f(A,B,C) */
100 17664628 : pp1_duplicate (A, A, n); /* A = 2*A */
101 : }
102 3818272125 : else if ((d + 3) / 4 <= e) /* <==> (d <= 4 * e) */
103 : { /* condition 3 */
104 3416949950 : d -= e;
105 3416949950 : pp1_add3 (C, B, A, C, n, t); /* C = f(B,A,C) */
106 3416949950 : mpres_swap (B, C, n);
107 : }
108 401322175 : else if ((d + e) % 2 == 0)
109 : { /* condition 4 */
110 172809930 : d = (d - e) / 2;
111 172809930 : pp1_add3 (B, B, A, C, n, t); /* B = f(B,A,C) */
112 172809930 : pp1_duplicate (A, A, n); /* A = 2*A */
113 : }
114 : /* d+e is now odd */
115 228512245 : else if (d % 2 == 0)
116 : { /* condition 5 */
117 146946383 : d /= 2;
118 146946383 : pp1_add3 (C, C, A, B, n, t); /* C = f(C,A,B) */
119 146946383 : pp1_duplicate (A, A, n); /* A = 2*A */
120 : }
121 : /* d is odd, e even */
122 81565862 : else if (d % 3 == 0)
123 : { /* condition 6 */
124 31579726 : d = d / 3 - e;
125 31579726 : pp1_duplicate (T, A, n); /* T = 2*A */
126 31579726 : pp1_add3 (T2, A, B, C, n, t); /* T2 = f(A,B,C) */
127 31579726 : pp1_add3 (A, T, A, A, n, t); /* A = f(T,A,A) */
128 31579726 : pp1_add3 (C, T, T2, C, n, t); /* C = f(T,T2,C) */
129 31579726 : mpres_swap (B, C, n);
130 : }
131 49986136 : else if ((d + e) % 3 == 0) /* d+e <= val[i]*k < k < 2^32 */
132 : { /* condition 7 */
133 30124979 : d = (d - 2 * e) / 3;
134 30124979 : pp1_add3 (T, A, B, C, n, t); /* T1 = f(A,B,C) */
135 30124979 : pp1_add3 (B, T, A, B, n, t); /* B = f(T1,A,B) */
136 30124979 : pp1_duplicate (T, A, n);
137 30124979 : pp1_add3 (A, A, T, A, n, t); /* A = 3*A */
138 : }
139 19861157 : else if ((d - e) % 3 == 0)
140 : { /* condition 8: never happens? */
141 5799481 : d = (d - e) / 3;
142 5799481 : pp1_add3 (T, A, B, C, n, t); /* T1 = f(A,B,C) */
143 5799481 : pp1_add3 (C, C, A, B, n, t); /* C = f(A,C,B) */
144 5799481 : mpres_swap (B, T, n); /* swap B and T */
145 5799481 : pp1_duplicate (T, A, n);
146 5799481 : pp1_add3 (A, A, T, A, n, t); /* A = 3*A */
147 : }
148 : else /* necessarily e is even */
149 : { /* condition 9: never happens? */
150 14061676 : e /= 2;
151 14061676 : pp1_add3 (C, C, B, A, n, t); /* C = f(C,B,A) */
152 14061676 : pp1_duplicate (B, B, n); /* B = 2*B */
153 : }
154 : }
155 :
156 115424819 : pp1_add3 (A, A, B, C, n, t);
157 :
158 : ASSERT(d == 1);
159 115424819 : }
|