Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mpn/fib2m.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, based on the previous
4
fib2_ui.c file.
5
6
   THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY.  THEY'RE ALMOST
7
   CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
8
   FUTURE GNU MP RELEASES.
9
10
Copyright 2001, 2002, 2005, 2009, 2018 Free Software Foundation, Inc.
11
12
This file is part of the GNU MP Library.
13
14
The GNU MP Library is free software; you can redistribute it and/or modify
15
it under the terms of either:
16
17
  * the GNU Lesser General Public License as published by the Free
18
    Software Foundation; either version 3 of the License, or (at your
19
    option) any later version.
20
21
or
22
23
  * the GNU General Public License as published by the Free Software
24
    Foundation; either version 2 of the License, or (at your option) any
25
    later version.
26
27
or both in parallel, as here.
28
29
The GNU MP Library is distributed in the hope that it will be useful, but
30
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
31
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
32
for more details.
33
34
You should have received copies of the GNU General Public License and the
35
GNU Lesser General Public License along with the GNU MP Library.  If not,
36
see https://www.gnu.org/licenses/.  */
37
38
#include <stdio.h>
39
#include "gmp-impl.h"
40
#include "longlong.h"
41
42
43
/* Stores |{ap,n}-{bp,n}| in {rp,n},
44
   returns the sign of {ap,n}-{bp,n}. */
45
static int
46
abs_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
47
203k
{
48
203k
  mp_limb_t  x, y;
49
530k
  while (--n >= 0)
50
530k
    {
51
530k
      x = ap[n];
52
530k
      y = bp[n];
53
530k
      if (x != y)
54
203k
        {
55
203k
          ++n;
56
203k
          if (x > y)
57
160k
            {
58
160k
              ASSERT_NOCARRY (mpn_sub_n (rp, ap, bp, n));
59
160k
              return 1;
60
160k
            }
61
42.8k
          else
62
42.8k
            {
63
42.8k
              ASSERT_NOCARRY (mpn_sub_n (rp, bp, ap, n));
64
42.8k
              return -1;
65
42.8k
            }
66
203k
        }
67
327k
      rp[n] = 0;
68
327k
    }
69
0
  return 0;
70
203k
}
71
72
/* Store F[n] at fp and F[n-1] at f1p.  Both are computed modulo m.
73
   fp and f1p should have room for mn*2+1 limbs.
74
75
   The sign of one or both the values may be flipped (n-F, instead of F),
76
   the return value is 0 (zero) if the signs are coherent (both positive
77
   or both negative) and 1 (one) otherwise.
78
79
   Notes:
80
81
   In F[2k+1] with k even, +2 is applied to 4*F[k]^2 just by ORing into the
82
   low limb.
83
84
   In F[2k+1] with k odd, -2 is applied to F[k-1]^2 just by ORing into the
85
   low limb.
86
87
   TODO: Should {tp, 2 * mn} be passed as a scratch pointer?
88
   Should the call to mpn_fib2_ui() obtain (up to) 2*mn limbs?
89
*/
90
91
int
92
mpn_fib2m (mp_ptr fp, mp_ptr f1p, mp_srcptr np, mp_size_t nn, mp_srcptr mp, mp_size_t mn)
93
1.04k
{
94
1.04k
  unsigned long nfirst;
95
1.04k
  mp_limb_t nh;
96
1.04k
  mp_bitcnt_t nbi;
97
1.04k
  mp_size_t sn, fn;
98
1.04k
  int   fcnt, ncnt;
99
100
1.04k
  ASSERT (! MPN_OVERLAP_P (fp, MAX(2*mn+1,5), f1p, MAX(2*mn+1,5)));
101
1.04k
  ASSERT (nn > 0 && np[nn - 1] != 0);
102
103
  /* Estimate the maximal n such that fibonacci(n) fits in mn limbs. */
104
1.04k
#if GMP_NUMB_BITS % 16 == 0
105
1.04k
  if (UNLIKELY (ULONG_MAX / (23 * (GMP_NUMB_BITS / 16)) <= mn))
106
0
    nfirst = ULONG_MAX;
107
1.04k
  else
108
1.04k
    nfirst = mn * (23 * (GMP_NUMB_BITS / 16));
109
#else
110
  {
111
    mp_bitcnt_t mbi;
112
    mbi = (mp_bitcnt_t) mn * GMP_NUMB_BITS;
113
114
    if (UNLIKELY (ULONG_MAX / 23 < mbi))
115
      {
116
  if (UNLIKELY (ULONG_MAX / 23 * 16 <= mbi))
117
    nfirst = ULONG_MAX;
118
  else
119
    nfirst = mbi / 16 * 23;
120
      }
121
    else
122
      nfirst = mbi * 23 / 16;
123
  }
124
#endif
125
126
1.04k
  sn = nn - 1;
127
1.04k
  nh = np[sn];
128
1.04k
  count_leading_zeros (ncnt, nh);
129
1.04k
  count_leading_zeros (fcnt, nfirst);
130
131
1.04k
  if (fcnt >= ncnt)
132
848
    {
133
848
      ncnt = fcnt - ncnt;
134
848
      nh >>= ncnt;
135
848
    }
136
196
  else if (sn > 0)
137
196
    {
138
196
      ncnt -= fcnt;
139
196
      nh <<= ncnt;
140
196
      ncnt = GMP_NUMB_BITS - ncnt;
141
196
      --sn;
142
196
      nh |= np[sn] >> ncnt;
143
196
    }
144
0
  else
145
0
    ncnt = 0;
146
147
1.04k
  nbi = sn * GMP_NUMB_BITS + ncnt;
148
1.04k
  if (nh > nfirst)
149
621
    {
150
621
      nh >>= 1;
151
621
      ++nbi;
152
621
    }
153
154
1.04k
  ASSERT (nh <= nfirst);
155
  /* Take a starting pair from mpn_fib2_ui. */
156
1.04k
  fn = mpn_fib2_ui (fp, f1p, nh);
157
1.04k
  MPN_ZERO (fp + fn, mn - fn);
158
1.04k
  MPN_ZERO (f1p + fn, mn - fn);
159
160
1.04k
  if (nbi == 0)
161
0
    {
162
0
      if (fn == mn)
163
0
  {
164
0
    mp_limb_t qp[2];
165
0
    mpn_tdiv_qr (qp, fp, 0, fp, fn, mp, mn);
166
0
    mpn_tdiv_qr (qp, f1p, 0, f1p, fn, mp, mn);
167
0
  }
168
169
0
      return 0;
170
0
    }
171
1.04k
  else
172
1.04k
    {
173
1.04k
      mp_ptr  tp;
174
1.04k
      unsigned  pb = nh & 1;
175
1.04k
      int neg;
176
1.04k
      TMP_DECL;
177
178
1.04k
      TMP_MARK;
179
180
1.04k
      tp = TMP_ALLOC_LIMBS (2 * mn + (mn < 2));
181
182
1.04k
      do
183
271k
  {
184
271k
    mp_ptr  rp;
185
    /* Here fp==F[k] and f1p==F[k-1], with k being the bits of n from
186
       nbi upwards.
187
188
       Based on the next bit of n, we'll double to the pair
189
       fp==F[2k],f1p==F[2k-1] or fp==F[2k+1],f1p==F[2k], according as
190
       that bit is 0 or 1 respectively.  */
191
192
271k
    mpn_sqr (tp, fp,  mn);
193
271k
    mpn_sqr (fp, f1p, mn);
194
195
    /* Calculate F[2k-1] = F[k]^2 + F[k-1]^2. */
196
271k
    f1p[2 * mn] = mpn_add_n (f1p, tp, fp, 2 * mn);
197
198
    /* Calculate F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k.
199
       pb is the low bit of our implied k.  */
200
201
    /* fp is F[k-1]^2 == 0 or 1 mod 4, like all squares. */
202
271k
    ASSERT ((fp[0] & 2) == 0);
203
271k
    ASSERT (pb == (pb & 1));
204
271k
    ASSERT ((fp[0] + (pb ? 2 : 0)) == (fp[0] | (pb << 1)));
205
271k
    fp[0] |= pb << 1;   /* possible -2 */
206
271k
#if HAVE_NATIVE_mpn_rsblsh2_n
207
271k
    fp[2 * mn] = 1 + mpn_rsblsh2_n (fp, fp, tp, 2 * mn);
208
271k
    MPN_INCR_U(fp, 2 * mn + 1, (1 ^ pb) << 1);  /* possible +2 */
209
271k
    fp[2 * mn] = (fp[2 * mn] - 1) & GMP_NUMB_MAX;
210
#else
211
    {
212
      mp_limb_t  c;
213
214
      c = mpn_lshift (tp, tp, 2 * mn, 2);
215
      tp[0] |= (1 ^ pb) << 1; /* possible +2 */
216
      c -= mpn_sub_n (fp, tp, fp, 2 * mn);
217
      fp[2 * mn] = c & GMP_NUMB_MAX;
218
    }
219
#endif
220
271k
    neg = fp[2 * mn] == GMP_NUMB_MAX;
221
222
    /* Calculate F[2k-1] = F[k]^2 + F[k-1]^2 */
223
    /* Calculate F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k */
224
225
    /* Calculate F[2k] = F[2k+1] - F[2k-1], replacing the unwanted one of
226
       F[2k+1] and F[2k-1].  */
227
271k
    --nbi;
228
271k
    pb = (np [nbi / GMP_NUMB_BITS] >> (nbi % GMP_NUMB_BITS)) & 1;
229
271k
    rp = pb ? f1p : fp;
230
271k
    if (neg)
231
67.2k
      {
232
        /* Calculate -(F[2k+1] - F[2k-1]) */
233
67.2k
        rp[2 * mn] = f1p[2 * mn] + 1 - mpn_sub_n (rp, f1p, fp, 2 * mn);
234
67.2k
        neg = ! pb;
235
67.2k
        if (pb) /* fp not overwritten, negate it. */
236
33.7k
    fp [2 * mn] = 1 ^ mpn_neg (fp, fp, 2 * mn);
237
67.2k
      }
238
203k
    else
239
203k
      {
240
203k
        neg = abs_sub_n (rp, fp, f1p, 2 * mn + 1) < 0;
241
203k
      }
242
243
271k
    mpn_tdiv_qr (tp, fp, 0, fp, 2 * mn + 1, mp, mn);
244
271k
    mpn_tdiv_qr (tp, f1p, 0, f1p, 2 * mn + 1, mp, mn);
245
271k
  }
246
271k
      while (nbi != 0);
247
248
1.04k
      TMP_FREE;
249
250
1.04k
      return neg;
251
1.04k
    }
252
1.04k
}