1 /*
2 Copyright (C) 2012 Fredrik Johansson
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "arb.h"
13
14 /*
15 Fibonacci numbers using binary powering.
16
17 D. Takahashi, "A fast algorithm for computing large Fibonacci numbers",
18 Information Processing Letters 75 (2000) 243–246
19 */
20
arb_fib_fmpz(arb_t f,const fmpz_t n,slong prec)21 void arb_fib_fmpz(arb_t f, const fmpz_t n, slong prec)
22 {
23 arb_t t, u;
24 slong wp, sign, i;
25
26 if (fmpz_sgn(n) < 0)
27 {
28 fmpz_t m;
29 fmpz_init(m);
30 fmpz_neg(m, n);
31 arb_fib_fmpz(f, m, prec);
32 if (fmpz_is_even(m))
33 arb_neg(f, f);
34 fmpz_clear(m);
35 return;
36 }
37
38 if (fmpz_cmp_ui(n, 4) <= 0)
39 {
40 ulong x = fmpz_get_ui(n);
41 arb_set_ui(f, x - (x > 1));
42 return;
43 }
44
45 wp = ARF_PREC_ADD(prec, 3 * fmpz_bits(n));
46
47 arb_init(u);
48 arb_init(t);
49 arb_set_ui(f, UWORD(1));
50 arb_set_ui(u, UWORD(1));
51 sign = -1;
52
53 for (i = fmpz_flog_ui(n, UWORD(2)) - 1; i > 0; i--)
54 {
55 arb_mul(t, f, f, wp);
56 arb_add(f, f, u, wp);
57 arb_mul_2exp_si(f, f, -1);
58 arb_mul(f, f, f, wp);
59 arb_mul_2exp_si(f, f, 1);
60 arb_submul_ui(f, t, 3, wp);
61 arb_sub_si(f, f, 2 * sign, wp);
62 arb_mul_ui(u, t, 5, wp);
63 arb_add_si(u, u, 2 * sign, wp);
64 sign = 1;
65
66 if (fmpz_tstbit(n, i))
67 {
68 arb_set(t, f);
69 arb_add(f, f, u, wp);
70 arb_mul_2exp_si(f, f, -1);
71 arb_mul_2exp_si(t, t, 1);
72 arb_add(u, f, t, wp);
73 sign = -1;
74 }
75 }
76
77 if (fmpz_tstbit(n, 0))
78 {
79 arb_add(f, f, u, wp);
80 arb_mul_2exp_si(f, f, -1);
81 arb_mul(f, f, u, wp);
82 arb_sub_si(f, f, sign, prec);
83 }
84 else
85 {
86 arb_mul(f, f, u, prec);
87 }
88
89 arb_clear(u);
90 arb_clear(t);
91 }
92
arb_fib_ui(arb_t f,ulong n,slong prec)93 void arb_fib_ui(arb_t f, ulong n, slong prec)
94 {
95 fmpz_t t;
96 fmpz_init_set_ui(t, n);
97 arb_fib_fmpz(f, t, prec);
98 fmpz_clear(t);
99 }
100