LCOV - code coverage report
Current view: top level - ecm - lucas.c (source / functions) Hit Total Coverage
Test: unnamed Lines: 70 70 100.0 %
Date: 2022-03-21 11:19:20 Functions: 3 3 100.0 %

          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 : }

Generated by: LCOV version 1.14