Coverage Report

Created: 2025-03-09 06:52

/src/gmp-6.2.1/mpn/fib2_ui.c
Line
Count
Source
1
/* mpn_fib2_ui -- calculate Fibonacci numbers.
2
3
   THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY.  THEY'RE ALMOST
4
   CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
5
   FUTURE GNU MP RELEASES.
6
7
Copyright 2001, 2002, 2005, 2009, 2018 Free Software Foundation, Inc.
8
9
This file is part of the GNU MP Library.
10
11
The GNU MP Library is free software; you can redistribute it and/or modify
12
it under the terms of either:
13
14
  * the GNU Lesser General Public License as published by the Free
15
    Software Foundation; either version 3 of the License, or (at your
16
    option) any later version.
17
18
or
19
20
  * the GNU General Public License as published by the Free Software
21
    Foundation; either version 2 of the License, or (at your option) any
22
    later version.
23
24
or both in parallel, as here.
25
26
The GNU MP Library is distributed in the hope that it will be useful, but
27
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
28
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
29
for more details.
30
31
You should have received copies of the GNU General Public License and the
32
GNU Lesser General Public License along with the GNU MP Library.  If not,
33
see https://www.gnu.org/licenses/.  */
34
35
#include <stdio.h>
36
#include "gmp-impl.h"
37
38
/* change this to "#define TRACE(x) x" for diagnostics */
39
#define TRACE(x)
40
41
42
/* Store F[n] at fp and F[n-1] at f1p.  fp and f1p should have room for
43
   MPN_FIB2_SIZE(n) limbs.
44
45
   The return value is the actual number of limbs stored, this will be at
46
   least 1.  fp[size-1] will be non-zero, except when n==0, in which case
47
   fp[0] is 0 and f1p[0] is 1.  f1p[size-1] can be zero, since F[n-1]<F[n]
48
   (for n>0).
49
50
   Notes: F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k.
51
52
   In F[2k+1] with k even, +2 is applied to 4*F[k]^2 just by ORing into the
53
   low limb.
54
55
   In F[2k+1] with k odd, -2 is applied to F[k-1]^2 just by ORing into the
56
   low limb.
57
*/
58
59
mp_size_t
60
mpn_fib2_ui (mp_ptr fp, mp_ptr f1p, unsigned long int n)
61
702
{
62
702
  mp_size_t      size;
63
702
  unsigned long  nfirst, mask;
64
65
702
  TRACE (printf ("mpn_fib2_ui n=%lu\n", n));
66
67
702
  ASSERT (! MPN_OVERLAP_P (fp, MPN_FIB2_SIZE(n), f1p, MPN_FIB2_SIZE(n)));
68
69
  /* Take a starting pair from the table. */
70
702
  mask = 1;
71
2.08k
  for (nfirst = n; nfirst > FIB_TABLE_LIMIT; nfirst /= 2)
72
1.38k
    mask <<= 1;
73
702
  TRACE (printf ("nfirst=%lu mask=0x%lX\n", nfirst, mask));
74
75
702
  f1p[0] = FIB_TABLE ((int) nfirst - 1);
76
702
  fp[0]  = FIB_TABLE (nfirst);
77
702
  size = 1;
78
79
  /* Skip to the end if the table lookup gives the final answer. */
80
702
  if (mask != 1)
81
503
    {
82
503
      mp_size_t  alloc;
83
503
      mp_ptr        xp;
84
503
      TMP_DECL;
85
86
503
      TMP_MARK;
87
503
      alloc = MPN_FIB2_SIZE (n);
88
503
      xp = TMP_ALLOC_LIMBS (alloc);
89
90
503
      do
91
1.38k
  {
92
    /* Here fp==F[k] and f1p==F[k-1], with k being the bits of n from
93
       n&mask upwards.
94
95
       The next bit of n is n&(mask>>1) and we'll double to the pair
96
       fp==F[2k],f1p==F[2k-1] or fp==F[2k+1],f1p==F[2k], according as
97
       that bit is 0 or 1 respectively.  */
98
99
1.38k
    TRACE (printf ("k=%lu mask=0x%lX size=%ld alloc=%ld\n",
100
1.38k
       n >> refmpn_count_trailing_zeros(mask),
101
1.38k
       mask, size, alloc);
102
1.38k
     mpn_trace ("fp ", fp, size);
103
1.38k
     mpn_trace ("f1p", f1p, size));
104
105
    /* fp normalized, f1p at most one high zero */
106
1.38k
    ASSERT (fp[size-1] != 0);
107
1.38k
    ASSERT (f1p[size-1] != 0 || f1p[size-2] != 0);
108
109
    /* f1p[size-1] might be zero, but this occurs rarely, so it's not
110
       worth bothering checking for it */
111
1.38k
    ASSERT (alloc >= 2*size);
112
1.38k
    mpn_sqr (xp, fp,  size);
113
1.38k
    mpn_sqr (fp, f1p, size);
114
1.38k
    size *= 2;
115
116
    /* Shrink if possible.  Since fp was normalized there'll be at
117
       most one high zero on xp (and if there is then there's one on
118
       yp too).  */
119
1.38k
    ASSERT (xp[size-1] != 0 || fp[size-1] == 0);
120
1.38k
    size -= (xp[size-1] == 0);
121
1.38k
    ASSERT (xp[size-1] != 0);  /* only one xp high zero */
122
123
    /* Calculate F[2k-1] = F[k]^2 + F[k-1]^2. */
124
1.38k
    f1p[size] = mpn_add_n (f1p, xp, fp, size);
125
126
    /* Calculate F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k.
127
       n&mask is the low bit of our implied k.  */
128
129
1.38k
    ASSERT ((fp[0] & 2) == 0);
130
    /* fp is F[k-1]^2 == 0 or 1 mod 4, like all squares. */
131
1.38k
    fp[0] |= (n & mask ? 2 : 0);      /* possible -2 */
132
1.38k
#if HAVE_NATIVE_mpn_rsblsh2_n
133
1.38k
    fp[size] = mpn_rsblsh2_n (fp, fp, xp, size);
134
1.38k
    MPN_INCR_U(fp, size + 1, (n & mask ? 0 : 2)); /* possible +2 */
135
#else
136
    {
137
      mp_limb_t  c;
138
139
      c = mpn_lshift (xp, xp, size, 2);
140
      xp[0] |= (n & mask ? 0 : 2);  /* possible +2 */
141
      c -= mpn_sub_n (fp, xp, fp, size);
142
      fp[size] = c;
143
    }
144
#endif
145
1.38k
    ASSERT (alloc >= size+1);
146
1.38k
    size += (fp[size] != 0);
147
148
    /* now n&mask is the new bit of n being considered */
149
1.38k
    mask >>= 1;
150
151
    /* Calculate F[2k] = F[2k+1] - F[2k-1], replacing the unwanted one of
152
       F[2k+1] and F[2k-1].  */
153
1.38k
    if (n & mask)
154
1.38k
      ASSERT_NOCARRY (mpn_sub_n (f1p, fp, f1p, size));
155
718
    else {
156
718
      ASSERT_NOCARRY (mpn_sub_n ( fp, fp, f1p, size));
157
158
      /* Can have a high zero after replacing F[2k+1] with F[2k].
159
         f1p will have a high zero if fp does. */
160
718
      ASSERT (fp[size-1] != 0 || f1p[size-1] == 0);
161
718
      size -= (fp[size-1] == 0);
162
718
    }
163
1.38k
  }
164
1.38k
      while (mask != 1);
165
166
503
      TMP_FREE;
167
503
    }
168
169
702
  TRACE (printf ("done size=%ld\n", size);
170
702
   mpn_trace ("fp ", fp, size);
171
702
   mpn_trace ("f1p", f1p, size));
172
173
702
  return size;
174
702
}