1/* mpn_fib2m -- calculate Fibonacci numbers, modulo m.
2
3Contributed to the GNU project by Marco Bodrato, based on the previous
4fib2_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
10Copyright 2001, 2002, 2005, 2009, 2018 Free Software Foundation, Inc.
11
12This file is part of the GNU MP Library.
13
14The GNU MP Library is free software; you can redistribute it and/or modify
15it 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
21or
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
27or both in parallel, as here.
28
29The GNU MP Library is distributed in the hope that it will be useful, but
30WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
31or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
32for more details.
33
34You should have received copies of the GNU General Public License and the
35GNU Lesser General Public License along with the GNU MP Library.  If not,
36see 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}. */
45static int
46abs_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
47{
48  mp_limb_t  x, y;
49  while (--n >= 0)
50    {
51      x = ap[n];
52      y = bp[n];
53      if (x != y)
54        {
55          ++n;
56          if (x > y)
57            {
58              ASSERT_NOCARRY (mpn_sub_n (rp, ap, bp, n));
59              return 1;
60            }
61          else
62            {
63              ASSERT_NOCARRY (mpn_sub_n (rp, bp, ap, n));
64              return -1;
65            }
66        }
67      rp[n] = 0;
68    }
69  return 0;
70}
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
91int
92mpn_fib2m (mp_ptr fp, mp_ptr f1p, mp_srcptr np, mp_size_t nn, mp_srcptr mp, mp_size_t mn)
93{
94  unsigned long	nfirst;
95  mp_limb_t	nh;
96  mp_bitcnt_t	nbi;
97  mp_size_t	sn, fn;
98  int		fcnt, ncnt;
99
100  ASSERT (! MPN_OVERLAP_P (fp, MAX(2*mn+1,5), f1p, MAX(2*mn+1,5)));
101  ASSERT (nn > 0 && np[nn - 1] != 0);
102
103  /* Estimate the maximal n such that fibonacci(n) fits in mn limbs. */
104#if GMP_NUMB_BITS % 16 == 0
105  if (UNLIKELY (ULONG_MAX / (23 * (GMP_NUMB_BITS / 16)) <= mn))
106    nfirst = ULONG_MAX;
107  else
108    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  sn = nn - 1;
127  nh = np[sn];
128  count_leading_zeros (ncnt, nh);
129  count_leading_zeros (fcnt, nfirst);
130
131  if (fcnt >= ncnt)
132    {
133      ncnt = fcnt - ncnt;
134      nh >>= ncnt;
135    }
136  else if (sn > 0)
137    {
138      ncnt -= fcnt;
139      nh <<= ncnt;
140      ncnt = GMP_NUMB_BITS - ncnt;
141      --sn;
142      nh |= np[sn] >> ncnt;
143    }
144  else
145    ncnt = 0;
146
147  nbi = sn * GMP_NUMB_BITS + ncnt;
148  if (nh > nfirst)
149    {
150      nh >>= 1;
151      ++nbi;
152    }
153
154  ASSERT (nh <= nfirst);
155  /* Take a starting pair from mpn_fib2_ui. */
156  fn = mpn_fib2_ui (fp, f1p, nh);
157  MPN_ZERO (fp + fn, mn - fn);
158  MPN_ZERO (f1p + fn, mn - fn);
159
160  if (nbi == 0)
161    {
162      if (fn == mn)
163	{
164	  mp_limb_t qp[2];
165	  mpn_tdiv_qr (qp, fp, 0, fp, fn, mp, mn);
166	  mpn_tdiv_qr (qp, f1p, 0, f1p, fn, mp, mn);
167	}
168
169      return 0;
170    }
171  else
172    {
173      mp_ptr	tp;
174      unsigned	pb = nh & 1;
175      int	neg;
176      TMP_DECL;
177
178      TMP_MARK;
179
180      tp = TMP_ALLOC_LIMBS (2 * mn + (mn < 2));
181
182      do
183	{
184	  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	  mpn_sqr (tp, fp,  mn);
193	  mpn_sqr (fp, f1p, mn);
194
195	  /* Calculate F[2k-1] = F[k]^2 + F[k-1]^2. */
196	  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	  ASSERT ((fp[0] & 2) == 0);
203	  ASSERT (pb == (pb & 1));
204	  ASSERT ((fp[0] + (pb ? 2 : 0)) == (fp[0] | (pb << 1)));
205	  fp[0] |= pb << 1;		/* possible -2 */
206#if HAVE_NATIVE_mpn_rsblsh2_n
207	  fp[2 * mn] = 1 + mpn_rsblsh2_n (fp, fp, tp, 2 * mn);
208	  MPN_INCR_U(fp, 2 * mn + 1, (1 ^ pb) << 1);	/* possible +2 */
209	  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	  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	  --nbi;
228	  pb = (np [nbi / GMP_NUMB_BITS] >> (nbi % GMP_NUMB_BITS)) & 1;
229	  rp = pb ? f1p : fp;
230	  if (neg)
231	    {
232	      /* Calculate -(F[2k+1] - F[2k-1]) */
233	      rp[2 * mn] = f1p[2 * mn] + 1 - mpn_sub_n (rp, f1p, fp, 2 * mn);
234	      neg = ! pb;
235	      if (pb) /* fp not overwritten, negate it. */
236		fp [2 * mn] = 1 ^ mpn_neg (fp, fp, 2 * mn);
237	    }
238	  else
239	    {
240	      neg = abs_sub_n (rp, fp, f1p, 2 * mn + 1) < 0;
241	    }
242
243	  mpn_tdiv_qr (tp, fp, 0, fp, 2 * mn + 1, mp, mn);
244	  mpn_tdiv_qr (tp, f1p, 0, f1p, 2 * mn + 1, mp, mn);
245	}
246      while (nbi != 0);
247
248      TMP_FREE;
249
250      return neg;
251    }
252}
253