1 /*
2     Copyright (C) 2015 William Hart
3     Copyright (C) 2015 Fredrik Johansson
4     Copyright (C) 2015 Kushagra Singh
5 
6     This file is part of FLINT.
7 
8     FLINT is free software: you can redistribute it and/or modify it under
9     the terms of the GNU Lesser General Public License (LGPL) as published
10     by the Free Software Foundation; either version 2.1 of the License, or
11     (at your option) any later version.  See <http://www.gnu.org/licenses/>.
12 */
13 
14 #include <gmp.h>
15 #define ulong ulongxx /* interferes with system includes */
16 #include <math.h>
17 #undef ulong
18 #define ulong mp_limb_t
19 #include "flint.h"
20 #include "ulong_extras.h"
21 #include "longlong.h"
22 
23 
24 /* this table contains the value of UWORD_MAX / n, for n in range [1, FLINT_BITS] */
25 
26 static const mp_limb_t mul_factor[] = {
27 #ifdef FLINT64
28                         UWORD(0), UWORD_MAX,        UWORD(9223372036854775807),
29                         UWORD(6148914691236517205), UWORD(4611686018427387903),
30                         UWORD(3689348814741910323), UWORD(3074457345618258602),
31                         UWORD(2635249153387078802), UWORD(2305843009213693951),
32                         UWORD(2049638230412172401), UWORD(1844674407370955161),
33                         UWORD(1676976733973595601), UWORD(1537228672809129301),
34                         UWORD(1418980313362273201), UWORD(1317624576693539401),
35                         UWORD(1229782938247303441), UWORD(1152921504606846975),
36                         UWORD(1085102592571150095), UWORD(1024819115206086200),
37                         UWORD(970881267037344821), UWORD(922337203685477580),
38                         UWORD(878416384462359600), UWORD(838488366986797800),
39                         UWORD(802032351030850070), UWORD(768614336404564650),
40                         UWORD(737869762948382064), UWORD(709490156681136600),
41                         UWORD(683212743470724133), UWORD(658812288346769700),
42                         UWORD(636094623231363848), UWORD(614891469123651720),
43                         UWORD(595056260442243600), UWORD(576460752303423487),
44                         UWORD(558992244657865200), UWORD(542551296285575047),
45                         UWORD(527049830677415760), UWORD(512409557603043100),
46                         UWORD(498560650640798692), UWORD(485440633518672410),
47                         UWORD(472993437787424400), UWORD(461168601842738790),
48                         UWORD(449920587163647600), UWORD(439208192231179800),
49                         UWORD(428994048225803525), UWORD(419244183493398900),
50                         UWORD(409927646082434480), UWORD(401016175515425035),
51                         UWORD(392483916461905353), UWORD(384307168202282325),
52                         UWORD(376464164769582686), UWORD(368934881474191032),
53                         UWORD(361700864190383365), UWORD(354745078340568300),
54                         UWORD(348051774975651917), UWORD(341606371735362066),
55                         UWORD(335395346794719120), UWORD(329406144173384850),
56                         UWORD(323627089012448273), UWORD(318047311615681924),
57                         UWORD(312656679215416129), UWORD(307445734561825860),
58                         UWORD(302405640552615600), UWORD(297528130221121800),
59                         UWORD(292805461487453200), UWORD(288230376151711743)
60 #else
61                         UWORD(0), UWORD_MAX, UWORD(2147483647), UWORD(1431655765),
62                         UWORD(1073741823), UWORD(858993459), UWORD(715827882),
63                         UWORD(613566756), UWORD(536870911), UWORD(477218588),
64                         UWORD(429496729), UWORD(390451572), UWORD(357913941),
65                         UWORD(330382099), UWORD(306783378), UWORD(286331153),
66                         UWORD(268435455), UWORD(252645135), UWORD(238609294),
67                         UWORD(226050910), UWORD(214748364), UWORD(204522252),
68                         UWORD(195225786), UWORD(186737708), UWORD(178956970),
69                         UWORD(171798691), UWORD(165191049), UWORD(159072862),
70                         UWORD(153391689), UWORD(148102320), UWORD(143165576),
71                         UWORD(138547332), UWORD(134217727),
72 #endif
73                                                             };
74                         /* this table consists of 65 values in case of FLINT64,
75                            otherwise 33 */
76 
77 /* function to get a good approximation of the cube root */
78 /* Algorithm for this approximation is mentioned in this article */
79 /* http://en.wikipedia.org/wiki/Fast_inverse_square_root */
80 /* Intead of the inverse square root, we calculate the nth root */
81 
82 mp_limb_t
n_root_estimate(double a,int n)83 n_root_estimate(double a, int n)
84 {
85     typedef union {
86         slong      uword_val;
87 #ifdef FLINT64
88         double     double_val;
89 #else
90         float      double_val;
91 #endif
92     } uni;
93 
94     uni alias;
95     ulong i, hi, lo, s;
96     alias.double_val = a;
97 
98 #ifdef FLINT64
99     s = ((1 << 10) - 1);
100     s <<= 52;
101 #else
102     s = ((1 << 7) - 1);
103     s <<= 23;
104 #endif
105 
106     i = alias.uword_val;
107     i -= s;
108     umul_ppmm(hi, lo, i, mul_factor[n]);
109     i = hi;
110     i += s;
111     alias.uword_val = i;
112     return (mp_limb_t)alias.double_val;
113 }
114