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