Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mpn/strongfibo.c
Line
Count
Source (jump to first uncovered line)
1
/* mpn_fib2m -- calculate Fibonacci numbers, modulo m.
2
3
Contributed to the GNU project by Marco Bodrato.
4
5
   THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY.  THEY'RE ALMOST
6
   CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
7
   FUTURE GNU MP RELEASES.
8
9
Copyright 2001, 2002, 2005, 2009, 2018, 2022 Free Software Foundation, Inc.
10
11
This file is part of the GNU MP Library.
12
13
The GNU MP Library is free software; you can redistribute it and/or modify
14
it under the terms of either:
15
16
  * the GNU Lesser General Public License as published by the Free
17
    Software Foundation; either version 3 of the License, or (at your
18
    option) any later version.
19
20
or
21
22
  * the GNU General Public License as published by the Free Software
23
    Foundation; either version 2 of the License, or (at your option) any
24
    later version.
25
26
or both in parallel, as here.
27
28
The GNU MP Library is distributed in the hope that it will be useful, but
29
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
30
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
31
for more details.
32
33
You should have received copies of the GNU General Public License and the
34
GNU Lesser General Public License along with the GNU MP Library.  If not,
35
see https://www.gnu.org/licenses/.  */
36
37
#include <stdio.h>
38
#include "gmp-impl.h"
39
40
41
#if ! HAVE_NATIVE_mpn_rsblsh1_n && ! HAVE_NATIVE_mpn_sublsh1_n
42
/* Stores |{ap,n}-{bp,n}| in {rp,n},
43
   returns the sign of {ap,n}-{bp,n}. */
44
static int
45
abs_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
46
{
47
  mp_limb_t  x, y;
48
  while (--n >= 0)
49
    {
50
      x = ap[n];
51
      y = bp[n];
52
      if (x != y)
53
        {
54
          ++n;
55
          if (x > y)
56
            {
57
              ASSERT_NOCARRY (mpn_sub_n (rp, ap, bp, n));
58
              return 1;
59
            }
60
          else
61
            {
62
              ASSERT_NOCARRY (mpn_sub_n (rp, bp, ap, n));
63
              return -1;
64
            }
65
        }
66
      rp[n] = 0;
67
    }
68
  return 0;
69
}
70
#endif
71
72
/* Computes at most count terms of the sequence needed by the
73
   Lucas-Lehmer-Riesel test, indexing backward:
74
   L_i = L_{i+1}^2 - 2
75
76
   The sequence is computed modulo M = {mp, mn}.
77
   The starting point is given in L_{count+1} = {lp, mn}.
78
   The scratch pointed by sp, needs a space of at least 3 * mn + 1 limbs.
79
80
   Returns the index i>0 if L_i = 0 (mod M) is found within the
81
   computed count terms of the sequence.  Otherwise it returns zero.
82
83
   Note: (+/-2)^2-2=2, (+/-1)^2-2=-1, 0^2-2=-2
84
 */
85
86
static mp_bitcnt_t
87
mpn_llriter (mp_ptr lp, mp_srcptr mp, mp_size_t mn, mp_bitcnt_t count, mp_ptr sp)
88
258
{
89
258
  do
90
473
    {
91
473
      mpn_sqr (sp, lp, mn);
92
473
      mpn_tdiv_qr (sp + 2 * mn, lp, 0, sp, 2 * mn, mp, mn);
93
473
      if (lp[0] < 5)
94
258
  {
95
    /* If L^2 % M < 5, |L^2 % M - 2| <= 2 */
96
258
    if (mn == 1 || mpn_zero_p (lp + 1, mn - 1))
97
258
      return (lp[0] == 2) ? count : 0;
98
0
    else
99
0
      MPN_DECR_U (lp, mn, 2);
100
258
  }
101
215
      else
102
215
  lp[0] -= 2;
103
473
    } while (--count != 0);
104
0
  return 0;
105
258
}
106
107
/* Store the Lucas' number L[n] at lp (maybe), computed modulo m.  lp
108
   and scratch should have room for mn*2+1 limbs.
109
110
   Returns the size of L[n] normally.
111
112
   If F[n] is zero modulo m, or L[n] is, returns 0 and lp is
113
   undefined.
114
*/
115
116
static mp_size_t
117
mpn_lucm (mp_ptr lp, mp_srcptr np, mp_size_t nn, mp_srcptr mp, mp_size_t mn, mp_ptr scratch)
118
1.04k
{
119
1.04k
  int   neg;
120
1.04k
  mp_limb_t cy;
121
122
1.04k
  ASSERT (! MPN_OVERLAP_P (lp, MAX(2*mn+1,5), scratch, MAX(2*mn+1,5)));
123
1.04k
  ASSERT (nn > 0);
124
125
1.04k
  neg = mpn_fib2m (lp, scratch, np, nn, mp, mn);
126
127
  /* F[n] = +/-{lp, mn}, F[n-1] = +/-{scratch, mn} */
128
1.04k
  if (mpn_zero_p (lp, mn))
129
536
    return 0;
130
131
508
  if (neg) /* One sign is opposite, use sub instead of add. */
132
75
    {
133
75
#if HAVE_NATIVE_mpn_rsblsh1_n || HAVE_NATIVE_mpn_sublsh1_n
134
75
#if HAVE_NATIVE_mpn_rsblsh1_n
135
75
      cy = mpn_rsblsh1_n (lp, lp, scratch, mn); /* L[n] = +/-(2F[n-1]-(-F[n])) */
136
#else
137
      cy = mpn_sublsh1_n (lp, lp, scratch, mn); /* L[n] = -/+(F[n]-(-2F[n-1])) */
138
      if (cy != 0)
139
  cy = mpn_add_n (lp, lp, mp, mn) - cy;
140
#endif
141
75
      if (cy > 1)
142
18
  cy += mpn_add_n (lp, lp, mp, mn);
143
#else
144
      cy = mpn_lshift (scratch, scratch, mn, 1); /* 2F[n-1] */
145
      if (UNLIKELY (cy))
146
  cy -= mpn_sub_n (lp, scratch, lp, mn); /* L[n] = +/-(2F[n-1]-(-F[n])) */
147
      else
148
  abs_sub_n (lp, lp, scratch, mn);
149
#endif
150
75
      ASSERT (cy <= 1);
151
75
    }
152
433
  else
153
433
    {
154
433
#if HAVE_NATIVE_mpn_addlsh1_n
155
433
      cy = mpn_addlsh1_n (lp, lp, scratch, mn); /* L[n] = +/-(2F[n-1]+F[n])) */
156
#else
157
      cy = mpn_lshift (scratch, scratch, mn, 1);
158
      cy+= mpn_add_n (lp, lp, scratch, mn);
159
#endif
160
433
      ASSERT (cy <= 2);
161
433
    }
162
975
  while (cy || mpn_cmp (lp, mp, mn) >= 0)
163
467
    cy -= mpn_sub_n (lp, lp, mp, mn);
164
508
  MPN_NORMALIZE (lp, mn);
165
508
  return mn;
166
508
}
167
168
int
169
mpn_strongfibo (mp_srcptr mp, mp_size_t mn, mp_ptr scratch)
170
1.04k
{
171
1.04k
  mp_ptr  lp, sp;
172
1.04k
  mp_size_t en;
173
1.04k
  mp_bitcnt_t b0;
174
1.04k
  TMP_DECL;
175
176
1.04k
#if GMP_NUMB_BITS % 4 == 0
177
1.04k
  b0 = mpn_scan0 (mp, 0);
178
#else
179
  {
180
    mpz_t m = MPZ_ROINIT_N(mp, mn);
181
    b0 = mpz_scan0 (m, 0);
182
  }
183
  if (UNLIKELY (b0 == mn * GMP_NUMB_BITS))
184
    {
185
      en = 1;
186
      scratch [0] = 1;
187
    }
188
  else
189
#endif
190
1.04k
    {
191
1.04k
      int cnt = b0 % GMP_NUMB_BITS;
192
1.04k
      en = b0 / GMP_NUMB_BITS;
193
1.04k
      if (LIKELY (cnt != 0))
194
1.04k
  mpn_rshift (scratch, mp + en, mn - en, cnt);
195
0
      else
196
0
  MPN_COPY (scratch, mp + en, mn - en);
197
1.04k
      en = mn - en;
198
1.04k
      scratch [0] |= 1;
199
1.04k
      en -= scratch [en - 1] == 0;
200
1.04k
    }
201
1.04k
  TMP_MARK;
202
203
1.04k
  lp = TMP_ALLOC_LIMBS (4 * mn + 6);
204
1.04k
  sp = lp + 2 * mn + 3;
205
1.04k
  en = mpn_lucm (sp, scratch, en, mp, mn, lp);
206
1.04k
  if (en != 0 && LIKELY (--b0 != 0))
207
508
    {
208
508
      mpn_sqr (lp, sp, en);
209
508
      lp [0] |= 2; /* V^2 + 2 */
210
508
      if (LIKELY (2 * en >= mn))
211
508
  mpn_tdiv_qr (sp, lp, 0, lp, 2 * en, mp, mn);
212
0
      else
213
0
  MPN_ZERO (lp + 2 * en, mn - 2 * en);
214
508
      if (! mpn_zero_p (lp, mn) && LIKELY (--b0 != 0))
215
258
  b0 = mpn_llriter (lp, mp, mn, b0, lp + mn + 1);
216
508
    }
217
1.04k
  TMP_FREE;
218
1.04k
  return (b0 != 0);
219
1.04k
}