1 /*
2 Copyright (C) 2014 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 /* See verify_taylor.py for code to generate tables and
15 proof of correctness */
16
17 #define TMP_ALLOC_LIMBS(size) TMP_ALLOC((size) * sizeof(mp_limb_t))
18
19 #define ODD_RECIPROCAL_TAB_SIZE 256
20
21 #if FLINT_BITS == 64
22
23 const mp_limb_t odd_reciprocal_tab_numer[ODD_RECIPROCAL_TAB_SIZE] = {
24 UWORD(13835020108241056725), UWORD(4611673369413685575),
25 UWORD(2767004021648211345), UWORD(1976431444034436675),
26 UWORD(1537224456471228525), UWORD(1257729100749186975),
27 UWORD(1064232316018542825), UWORD(922334673882737115),
28 UWORD(813824712249473925), UWORD(728158953065318775),
29 UWORD(658810481344812225), UWORD(601522613401785075),
30 UWORD(553400804329642269), UWORD(512408152157076175),
31 UWORD(477069658904864025), UWORD(446290971233582475),
32 UWORD(419243033583062325), UWORD(395286288806887335),
33 UWORD(373919462384893425), UWORD(354744105339514275),
34 UWORD(337439514835147725), UWORD(321744653680024575),
35 UWORD(307444891294245705), UWORD(294362129962575675),
36 UWORD(43378025254434585), UWORD(41676926224848915),
37 UWORD(40104212027684805), UWORD(38645877044859903),
38 UWORD(37289881359075345), UWORD(36025817584191435),
39 UWORD(34844643237168765), UWORD(33738464086782455),
40 UWORD(32700357499496841), UWORD(31724227424884995),
41 UWORD(30804684600975285), UWORD(29936947006581615),
42 UWORD(18078521432946975), UWORD(17596427528068389),
43 UWORD(17139377462404275), UWORD(16705469172216825),
44 UWORD(16292988451915175), UWORD(15900386320543725),
45 UWORD(15526259583589755), UWORD(15169334075921025),
46 UWORD(14828450164102575), UWORD(14502550160495925),
47 UWORD(14190667361345475), UWORD(34099971069969189),
48 UWORD(33396878882959515), UWORD(32722194461081545),
49 UWORD(32074230214327455), UWORD(31451429627641485),
50 UWORD(30852354777591171), UWORD(30275675249038065),
51 UWORD(29720158271991495), UWORD(29184659924748405),
52 UWORD(28668117271213035), UWORD(70414222249896525),
53 UWORD(69210560331094875), UWORD(68047357636454625),
54 UWORD(66922607923455375), UWORD(65834435436895125),
55 UWORD(64781084469904803), UWORD(63760909911323625),
56 UWORD(62772368672388375), UWORD(61814011898764125),
57 UWORD(60884477885248875), UWORD(59982485620282225),
58 UWORD(22849011931010265), UWORD(22520249169412995),
59 UWORD(22200813010981605), UWORD(21890312129709135),
60 UWORD(21588376789989009), UWORD(21294657377880315),
61 UWORD(21008823050660445), UWORD(20730560493698055),
62 UWORD(20459572774826185), UWORD(11592603279646389),
63 UWORD(11444926804746435), UWORD(11300965461290505),
64 UWORD(11160580797175095), UWORD(11023641155491965),
65 UWORD(10890021262698123), UWORD(10759601846378385),
66 UWORD(10632269280149055), UWORD(10507915253480645),
67 UWORD(3812050066792275), UWORD(3768483780314649),
68 UWORD(3725902042683975), UWORD(3684271852262925),
69 UWORD(3643561666050075), UWORD(3603741319973025),
70 UWORD(3564781954351695), UWORD(3526655944144725),
71 UWORD(3489336833624675), UWORD(56769492397408245),
72 UWORD(56181207502098315), UWORD(55604989989256281),
73 UWORD(55040472324390735), UWORD(54487301748266205),
74 UWORD(53945139541815795), UWORD(53413660334507265),
75 UWORD(52892551453194999), UWORD(52381512308719685),
76 UWORD(7543627748528235), UWORD(7472124167973465),
77 UWORD(7401963377663855), UWORD(7333107904383261),
78 UWORD(7265521656416595), UWORD(7199169860467585),
79 UWORD(7134019002001815), UWORD(7070036768800005),
80 UWORD(827619730851221), UWORD(820327927055175),
81 UWORD(813163491011025), UWORD(806123114465475),
82 UWORD(799203602753325), UWORD(792401869963935),
83 UWORD(785714934352425), UWORD(779139913981275),
84 UWORD(1317621181161285), UWORD(1306776562386295),
85 UWORD(1296108998611713), UWORD(1285614188906355),
86 UWORD(1275287970521565), UWORD(1265126313385935),
87 UWORD(1255125314861145), UWORD(1245281194744587),
88 UWORD(30567046059275355), UWORD(30331007093566665),
89 UWORD(30098585583271135), UWORD(29869699000888845),
90 UWORD(29644267310316099), UWORD(29422212873534705),
91 UWORD(29203460361463815), UWORD(28987936668759285),
92 UWORD(438859621842675), UWORD(435667915502001),
93 UWORD(432522298783575), UWORD(429421780512725),
94 UWORD(426365397733275), UWORD(423352214710425),
95 UWORD(420381321975615), UWORD(417451835411325),
96 UWORD(22634811266652585), UWORD(22479245553479715),
97 UWORD(22325803604309205), UWORD(22174442223941007),
98 UWORD(22025119380682145), UWORD(21877794167433435),
99 UWORD(21732426764327565), UWORD(21588978402846855),
100 UWORD(19609694213279247), UWORD(19481943762378405),
101 UWORD(19355847038997315), UWORD(19231372138424985),
102 UWORD(19108487971406295), UWORD(18987164238254509),
103 UWORD(18867371403943755), UWORD(18749080674138465),
104 UWORD(9266890993958325), UWORD(9209510863964775),
105 UWORD(9152836950955761), UWORD(9096856296821475),
106 UWORD(9041556258542925), UWORD(8986924498672575),
107 UWORD(8932948976158025), UWORD(8879617937494395),
108 UWORD(552205768818165), UWORD(548947917674695),
109 UWORD(545728281793905), UWORD(542546192687235),
110 UWORD(539400997367309), UWORD(536292057901215),
111 UWORD(533218750979145), UWORD(17071036886667995),
112 UWORD(16974317130936165), UWORD(16878687175268919),
113 UWORD(16784128703698785), UWORD(16690623808413555),
114 UWORD(16598154978450045), UWORD(16506705088761615),
115 UWORD(16416257389645113), UWORD(912010168967625),
116 UWORD(907067024420375), UWORD(902177175232125),
117 UWORD(897339764104875), UWORD(892553952029649),
118 UWORD(887818917801375), UWORD(883133857549125),
119 UWORD(29941527972993945), UWORD(29785175346503115),
120 UWORD(29630447162884917), UWORD(29477318236978535),
121 UWORD(29325763901569905), UWORD(29175759994144995),
122 UWORD(29027282844047565), UWORD(28880309260027071),
123 UWORD(1449018065980485), UWORD(1441754817529455),
124 UWORD(1434564020434545), UWORD(1427444596015515),
125 UWORD(1420395486899389), UWORD(1413415656496935),
126 UWORD(1406504088494505), UWORD(592554335551395),
127 UWORD(589684823030565), UWORD(586842968461743),
128 UWORD(584028373888785), UWORD(581240648953755),
129 UWORD(578479410716445), UWORD(575744283479015),
130 UWORD(433241218467057), UWORD(431211985593675),
131 UWORD(429201673306525), UWORD(427210018209975),
132 UWORD(425236761774825), UWORD(423281650226435),
133 UWORD(421344434435925), UWORD(2623079648743785),
134 UWORD(2611183595915015), UWORD(2599394956655805),
135 UWORD(2587712282693307), UWORD(2576134151674545),
136 UWORD(2564659166589135), UWORD(2553285955207365),
137 UWORD(210505107582315), UWORD(209579810406129),
138 UWORD(208662612111135), UWORD(207753406829605),
139 UWORD(206852090530995), UWORD(205958560982265),
140 UWORD(205072717709223), UWORD(3779907202434825),
141 UWORD(3763788195174975), UWORD(3747806079696525),
142 UWORD(3731959119528675), UWORD(3716245607446449),
143 UWORD(3700663864857575), UWORD(3685212241204725),
144 UWORD(4501257379467765), UWORD(4482618632554855),
145 UWORD(4464133607265969), UWORD(4445800409700195),
146 UWORD(4427617176940685), UWORD(4409582076423615),
147 UWORD(4391693305322505), UWORD(355589942551207),
148 UWORD(354158997108345), UWORD(352739522170035),
149 UWORD(351331380364965), UWORD(349934436506655),
150 UWORD(348548557550193), UWORD(347173612549995),
151 UWORD(6291002587483845), UWORD(6266380268159055),
152 };
153
154 const mp_limb_t odd_reciprocal_tab_denom[ODD_RECIPROCAL_TAB_SIZE] = {
155 UWORD(13835020108241056725), UWORD(13835020108241056725),
156 UWORD(13835020108241056725), UWORD(13835020108241056725),
157 UWORD(13835020108241056725), UWORD(13835020108241056725),
158 UWORD(13835020108241056725), UWORD(13835020108241056725),
159 UWORD(13835020108241056725), UWORD(13835020108241056725),
160 UWORD(13835020108241056725), UWORD(13835020108241056725),
161 UWORD(13835020108241056725), UWORD(13835020108241056725),
162 UWORD(13835020108241056725), UWORD(13835020108241056725),
163 UWORD(13835020108241056725), UWORD(13835020108241056725),
164 UWORD(13835020108241056725), UWORD(13835020108241056725),
165 UWORD(13835020108241056725), UWORD(13835020108241056725),
166 UWORD(13835020108241056725), UWORD(13835020108241056725),
167 UWORD(2125523237467294665), UWORD(2125523237467294665),
168 UWORD(2125523237467294665), UWORD(2125523237467294665),
169 UWORD(2125523237467294665), UWORD(2125523237467294665),
170 UWORD(2125523237467294665), UWORD(2125523237467294665),
171 UWORD(2125523237467294665), UWORD(2125523237467294665),
172 UWORD(2125523237467294665), UWORD(2125523237467294665),
173 UWORD(1319732064605129175), UWORD(1319732064605129175),
174 UWORD(1319732064605129175), UWORD(1319732064605129175),
175 UWORD(1319732064605129175), UWORD(1319732064605129175),
176 UWORD(1319732064605129175), UWORD(1319732064605129175),
177 UWORD(1319732064605129175), UWORD(1319732064605129175),
178 UWORD(1319732064605129175), UWORD(3239497251647072955),
179 UWORD(3239497251647072955), UWORD(3239497251647072955),
180 UWORD(3239497251647072955), UWORD(3239497251647072955),
181 UWORD(3239497251647072955), UWORD(3239497251647072955),
182 UWORD(3239497251647072955), UWORD(3239497251647072955),
183 UWORD(3239497251647072955), UWORD(8097635558738100375),
184 UWORD(8097635558738100375), UWORD(8097635558738100375),
185 UWORD(8097635558738100375), UWORD(8097635558738100375),
186 UWORD(8097635558738100375), UWORD(8097635558738100375),
187 UWORD(8097635558738100375), UWORD(8097635558738100375),
188 UWORD(8097635558738100375), UWORD(8097635558738100375),
189 UWORD(3130314634548406305), UWORD(3130314634548406305),
190 UWORD(3130314634548406305), UWORD(3130314634548406305),
191 UWORD(3130314634548406305), UWORD(3130314634548406305),
192 UWORD(3130314634548406305), UWORD(3130314634548406305),
193 UWORD(3130314634548406305), UWORD(1796853508345190295),
194 UWORD(1796853508345190295), UWORD(1796853508345190295),
195 UWORD(1796853508345190295), UWORD(1796853508345190295),
196 UWORD(1796853508345190295), UWORD(1796853508345190295),
197 UWORD(1796853508345190295), UWORD(1796853508345190295),
198 UWORD(659484661555063575), UWORD(659484661555063575),
199 UWORD(659484661555063575), UWORD(659484661555063575),
200 UWORD(659484661555063575), UWORD(659484661555063575),
201 UWORD(659484661555063575), UWORD(659484661555063575),
202 UWORD(659484661555063575), UWORD(10842973047904974795),
203 UWORD(10842973047904974795), UWORD(10842973047904974795),
204 UWORD(10842973047904974795), UWORD(10842973047904974795),
205 UWORD(10842973047904974795), UWORD(10842973047904974795),
206 UWORD(10842973047904974795), UWORD(10842973047904974795),
207 UWORD(1576618199442401115), UWORD(1576618199442401115),
208 UWORD(1576618199442401115), UWORD(1576618199442401115),
209 UWORD(1576618199442401115), UWORD(1576618199442401115),
210 UWORD(1576618199442401115), UWORD(1576618199442401115),
211 UWORD(186214439441524725), UWORD(186214439441524725),
212 UWORD(186214439441524725), UWORD(186214439441524725),
213 UWORD(186214439441524725), UWORD(186214439441524725),
214 UWORD(186214439441524725), UWORD(186214439441524725),
215 UWORD(317546704659869685), UWORD(317546704659869685),
216 UWORD(317546704659869685), UWORD(317546704659869685),
217 UWORD(317546704659869685), UWORD(317546704659869685),
218 UWORD(317546704659869685), UWORD(317546704659869685),
219 UWORD(7855730837233766235), UWORD(7855730837233766235),
220 UWORD(7855730837233766235), UWORD(7855730837233766235),
221 UWORD(7855730837233766235), UWORD(7855730837233766235),
222 UWORD(7855730837233766235), UWORD(7855730837233766235),
223 UWORD(119808676763050275), UWORD(119808676763050275),
224 UWORD(119808676763050275), UWORD(119808676763050275),
225 UWORD(119808676763050275), UWORD(119808676763050275),
226 UWORD(119808676763050275), UWORD(119808676763050275),
227 UWORD(6541460456062597065), UWORD(6541460456062597065),
228 UWORD(6541460456062597065), UWORD(6541460456062597065),
229 UWORD(6541460456062597065), UWORD(6541460456062597065),
230 UWORD(6541460456062597065), UWORD(6541460456062597065),
231 UWORD(5980956735050170335), UWORD(5980956735050170335),
232 UWORD(5980956735050170335), UWORD(5980956735050170335),
233 UWORD(5980956735050170335), UWORD(5980956735050170335),
234 UWORD(5980956735050170335), UWORD(5980956735050170335),
235 UWORD(2974672009060622325), UWORD(2974672009060622325),
236 UWORD(2974672009060622325), UWORD(2974672009060622325),
237 UWORD(2974672009060622325), UWORD(2974672009060622325),
238 UWORD(2974672009060622325), UWORD(2974672009060622325),
239 UWORD(186093344091721605), UWORD(186093344091721605),
240 UWORD(186093344091721605), UWORD(186093344091721605),
241 UWORD(186093344091721605), UWORD(186093344091721605),
242 UWORD(186093344091721605), UWORD(5991933947220466245),
243 UWORD(5991933947220466245), UWORD(5991933947220466245),
244 UWORD(5991933947220466245), UWORD(5991933947220466245),
245 UWORD(5991933947220466245), UWORD(5991933947220466245),
246 UWORD(5991933947220466245), UWORD(334707732011118375),
247 UWORD(334707732011118375), UWORD(334707732011118375),
248 UWORD(334707732011118375), UWORD(334707732011118375),
249 UWORD(334707732011118375), UWORD(334707732011118375),
250 UWORD(11407722157710693045), UWORD(11407722157710693045),
251 UWORD(11407722157710693045), UWORD(11407722157710693045),
252 UWORD(11407722157710693045), UWORD(11407722157710693045),
253 UWORD(11407722157710693045), UWORD(11407722157710693045),
254 UWORD(575260172194252545), UWORD(575260172194252545),
255 UWORD(575260172194252545), UWORD(575260172194252545),
256 UWORD(575260172194252545), UWORD(575260172194252545),
257 UWORD(575260172194252545), UWORD(243539831911623345),
258 UWORD(243539831911623345), UWORD(243539831911623345),
259 UWORD(243539831911623345), UWORD(243539831911623345),
260 UWORD(243539831911623345), UWORD(243539831911623345),
261 UWORD(184127517848499225), UWORD(184127517848499225),
262 UWORD(184127517848499225), UWORD(184127517848499225),
263 UWORD(184127517848499225), UWORD(184127517848499225),
264 UWORD(184127517848499225), UWORD(1151531965798521615),
265 UWORD(1151531965798521615), UWORD(1151531965798521615),
266 UWORD(1151531965798521615), UWORD(1151531965798521615),
267 UWORD(1151531965798521615), UWORD(1151531965798521615),
268 UWORD(95358813734788695), UWORD(95358813734788695),
269 UWORD(95358813734788695), UWORD(95358813734788695),
270 UWORD(95358813734788695), UWORD(95358813734788695),
271 UWORD(95358813734788695), UWORD(1765216663537063275),
272 UWORD(1765216663537063275), UWORD(1765216663537063275),
273 UWORD(1765216663537063275), UWORD(1765216663537063275),
274 UWORD(1765216663537063275), UWORD(1765216663537063275),
275 UWORD(2165104799523994965), UWORD(2165104799523994965),
276 UWORD(2165104799523994965), UWORD(2165104799523994965),
277 UWORD(2165104799523994965), UWORD(2165104799523994965),
278 UWORD(2165104799523994965), UWORD(176017021562847465),
279 UWORD(176017021562847465), UWORD(176017021562847465),
280 UWORD(176017021562847465), UWORD(176017021562847465),
281 UWORD(176017021562847465), UWORD(176017021562847465),
282 UWORD(3202120317029277105), UWORD(3202120317029277105),
283 };
284
285 #else
286
287 const mp_limb_t odd_reciprocal_tab_numer[ODD_RECIPROCAL_TAB_SIZE] = {
288 UWORD(1673196525), UWORD(557732175),
289 UWORD(334639305), UWORD(239028075),
290 UWORD(185910725), UWORD(152108775),
291 UWORD(128707425), UWORD(111546435),
292 UWORD(98423325), UWORD(88062975),
293 UWORD(79676025), UWORD(72747675),
294 UWORD(66927861), UWORD(12806255),
295 UWORD(11923065), UWORD(11153835),
296 UWORD(10477845), UWORD(9879111),
297 UWORD(9345105), UWORD(60902835),
298 UWORD(57931965), UWORD(55237455),
299 UWORD(52782457), UWORD(50536395),
300 UWORD(48473685), UWORD(3267715),
301 UWORD(3144405), UWORD(3030063),
302 UWORD(2923745), UWORD(2824635),
303 UWORD(6310395), UWORD(6110065),
304 UWORD(5922063), UWORD(5745285),
305 UWORD(5578755), UWORD(33304425),
306 UWORD(32391975), UWORD(31528189),
307 UWORD(30709275), UWORD(29931825),
308 UWORD(18208955), UWORD(17770185),
309 UWORD(17352063), UWORD(16953165),
310 UWORD(16572195), UWORD(28280835),
311 UWORD(27672645), UWORD(27090063),
312 UWORD(26531505), UWORD(25995515),
313 UWORD(1157205), UWORD(1134735),
314 UWORD(1113121), UWORD(1092315),
315 UWORD(1442445), UWORD(1416455),
316 UWORD(1391385), UWORD(1367187),
317 UWORD(590359), UWORD(580437),
318 UWORD(570843), UWORD(561561),
319 UWORD(2146173), UWORD(2112375),
320 UWORD(2079625), UWORD(2047875),
321 UWORD(2570805), UWORD(2532719),
322 UWORD(2495745), UWORD(2459835),
323 UWORD(1016015), UWORD(1001805),
324 UWORD(987987), UWORD(974545),
325 UWORD(3580965), UWORD(3533535),
326 UWORD(3487345), UWORD(3442347),
327 UWORD(4172637), UWORD(4120151),
328 UWORD(4068969), UWORD(4019043),
329 UWORD(1608711), UWORD(1589445),
330 UWORD(1570635), UWORD(1552265),
331 UWORD(5544525), UWORD(5481159),
332 UWORD(5419225), UWORD(5358675),
333 UWORD(6330885), UWORD(6261695),
334 UWORD(6194001), UWORD(6127755),
335 UWORD(2396095), UWORD(2371005),
336 UWORD(2346435), UWORD(2322369),
337 UWORD(8119797), UWORD(8038191),
338 UWORD(7958209), UWORD(7879803),
339 UWORD(9128493), UWORD(9040295),
340 UWORD(8953785), UWORD(8868915),
341 UWORD(3405815), UWORD(3374133),
342 UWORD(3343035), UWORD(3312505),
343 UWORD(11389725), UWORD(11287575),
344 UWORD(11187241), UWORD(11088675),
345 UWORD(12648405), UWORD(12538895),
346 UWORD(12431265), UWORD(12325467),
347 UWORD(4665519), UWORD(4626477),
348 UWORD(4588083), UWORD(4550321),
349 UWORD(15437253), UWORD(15312255),
350 UWORD(15189265), UWORD(15068235),
351 UWORD(16973565), UWORD(16840439),
352 UWORD(16709385), UWORD(16580355),
353 UWORD(6202855), UWORD(6155685),
354 UWORD(6109227), UWORD(6063465),
355 UWORD(73983), UWORD(73437),
356 UWORD(72899), UWORD(77283),
357 UWORD(76725), UWORD(76175),
358 UWORD(80655), UWORD(80085),
359 UWORD(79523), UWORD(84099),
360 UWORD(83517), UWORD(82943),
361 UWORD(87615), UWORD(87021),
362 UWORD(86435), UWORD(91203),
363 UWORD(90597), UWORD(89999),
364 UWORD(94863), UWORD(94245),
365 UWORD(93635), UWORD(98595),
366 UWORD(97965), UWORD(97343),
367 UWORD(102399), UWORD(101757),
368 UWORD(101123), UWORD(106275),
369 UWORD(105621), UWORD(104975),
370 UWORD(110223), UWORD(109557),
371 UWORD(108899), UWORD(114243),
372 UWORD(113565), UWORD(112895),
373 UWORD(118335), UWORD(117645),
374 UWORD(116963), UWORD(122499),
375 UWORD(121797), UWORD(121103),
376 UWORD(126735), UWORD(126021),
377 UWORD(125315), UWORD(131043),
378 UWORD(130317), UWORD(129599),
379 UWORD(135423), UWORD(134685),
380 UWORD(133955), UWORD(139875),
381 UWORD(139125), UWORD(138383),
382 UWORD(144399), UWORD(143637),
383 UWORD(142883), UWORD(148995),
384 UWORD(148221), UWORD(147455),
385 UWORD(153663), UWORD(152877),
386 UWORD(152099), UWORD(158403),
387 UWORD(157605), UWORD(156815),
388 UWORD(163215), UWORD(162405),
389 UWORD(161603), UWORD(168099),
390 UWORD(167277), UWORD(166463),
391 UWORD(173055), UWORD(172221),
392 UWORD(171395), UWORD(178083),
393 UWORD(177237), UWORD(176399),
394 UWORD(183183), UWORD(182325),
395 UWORD(181475), UWORD(188355),
396 UWORD(187485), UWORD(186623),
397 UWORD(193599), UWORD(192717),
398 UWORD(191843), UWORD(198915),
399 UWORD(198021), UWORD(197135),
400 UWORD(204303), UWORD(203397),
401 UWORD(202499), UWORD(209763),
402 UWORD(208845), UWORD(207935),
403 UWORD(215295), UWORD(214365),
404 UWORD(213443), UWORD(220899),
405 UWORD(219957), UWORD(219023),
406 UWORD(226575), UWORD(225621),
407 UWORD(224675), UWORD(232323),
408 UWORD(231357), UWORD(230399),
409 UWORD(238143), UWORD(237165),
410 UWORD(236195), UWORD(244035),
411 UWORD(243045), UWORD(242063),
412 UWORD(249999), UWORD(248997),
413 UWORD(248003), UWORD(256035),
414 UWORD(255021), UWORD(254015),
415 UWORD(262143), UWORD(261117),
416 };
417
418 const mp_limb_t odd_reciprocal_tab_denom[ODD_RECIPROCAL_TAB_SIZE] = {
419 UWORD(1673196525), UWORD(1673196525),
420 UWORD(1673196525), UWORD(1673196525),
421 UWORD(1673196525), UWORD(1673196525),
422 UWORD(1673196525), UWORD(1673196525),
423 UWORD(1673196525), UWORD(1673196525),
424 UWORD(1673196525), UWORD(1673196525),
425 UWORD(1673196525), UWORD(345768885),
426 UWORD(345768885), UWORD(345768885),
427 UWORD(345768885), UWORD(345768885),
428 UWORD(345768885), UWORD(2375210565),
429 UWORD(2375210565), UWORD(2375210565),
430 UWORD(2375210565), UWORD(2375210565),
431 UWORD(2375210565), UWORD(166653465),
432 UWORD(166653465), UWORD(166653465),
433 UWORD(166653465), UWORD(166653465),
434 UWORD(384934095), UWORD(384934095),
435 UWORD(384934095), UWORD(384934095),
436 UWORD(384934095), UWORD(2364614175),
437 UWORD(2364614175), UWORD(2364614175),
438 UWORD(2364614175), UWORD(2364614175),
439 UWORD(1474925355), UWORD(1474925355),
440 UWORD(1474925355), UWORD(1474925355),
441 UWORD(1474925355), UWORD(2573555985),
442 UWORD(2573555985), UWORD(2573555985),
443 UWORD(2573555985), UWORD(2573555985),
444 UWORD(116877705), UWORD(116877705),
445 UWORD(116877705), UWORD(116877705),
446 UWORD(157226505), UWORD(157226505),
447 UWORD(157226505), UWORD(157226505),
448 UWORD(69072003), UWORD(69072003),
449 UWORD(69072003), UWORD(69072003),
450 UWORD(268271625), UWORD(268271625),
451 UWORD(268271625), UWORD(268271625),
452 UWORD(341917065), UWORD(341917065),
453 UWORD(341917065), UWORD(341917065),
454 UWORD(143258115), UWORD(143258115),
455 UWORD(143258115), UWORD(143258115),
456 UWORD(533563785), UWORD(533563785),
457 UWORD(533563785), UWORD(533563785),
458 UWORD(655104009), UWORD(655104009),
459 UWORD(655104009), UWORD(655104009),
460 UWORD(265437315), UWORD(265437315),
461 UWORD(265437315), UWORD(265437315),
462 UWORD(959202825), UWORD(959202825),
463 UWORD(959202825), UWORD(959202825),
464 UWORD(1145890185), UWORD(1145890185),
465 UWORD(1145890185), UWORD(1145890185),
466 UWORD(452861955), UWORD(452861955),
467 UWORD(452861955), UWORD(452861955),
468 UWORD(1599600009), UWORD(1599600009),
469 UWORD(1599600009), UWORD(1599600009),
470 UWORD(1871341065), UWORD(1871341065),
471 UWORD(1871341065), UWORD(1871341065),
472 UWORD(725438595), UWORD(725438595),
473 UWORD(725438595), UWORD(725438595),
474 UWORD(2517129225), UWORD(2517129225),
475 UWORD(2517129225), UWORD(2517129225),
476 UWORD(2896484745), UWORD(2896484745),
477 UWORD(2896484745), UWORD(2896484745),
478 UWORD(1105728003), UWORD(1105728003),
479 UWORD(1105728003), UWORD(1105728003),
480 UWORD(3782126985), UWORD(3782126985),
481 UWORD(3782126985), UWORD(3782126985),
482 UWORD(4294311945), UWORD(4294311945),
483 UWORD(4294311945), UWORD(4294311945),
484 UWORD(1618945155), UWORD(1618945155),
485 UWORD(1618945155), UWORD(1618945155),
486 UWORD(19901427), UWORD(19901427),
487 UWORD(19901427), UWORD(21252825),
488 UWORD(21252825), UWORD(21252825),
489 UWORD(22664055), UWORD(22664055),
490 UWORD(22664055), UWORD(24136413),
491 UWORD(24136413), UWORD(24136413),
492 UWORD(25671195), UWORD(25671195),
493 UWORD(25671195), UWORD(27269697),
494 UWORD(27269697), UWORD(27269697),
495 UWORD(28933215), UWORD(28933215),
496 UWORD(28933215), UWORD(30663045),
497 UWORD(30663045), UWORD(30663045),
498 UWORD(32460483), UWORD(32460483),
499 UWORD(32460483), UWORD(34326825),
500 UWORD(34326825), UWORD(34326825),
501 UWORD(36263367), UWORD(36263367),
502 UWORD(36263367), UWORD(38271405),
503 UWORD(38271405), UWORD(38271405),
504 UWORD(40352235), UWORD(40352235),
505 UWORD(40352235), UWORD(42507153),
506 UWORD(42507153), UWORD(42507153),
507 UWORD(44737455), UWORD(44737455),
508 UWORD(44737455), UWORD(47044437),
509 UWORD(47044437), UWORD(47044437),
510 UWORD(49429395), UWORD(49429395),
511 UWORD(49429395), UWORD(51893625),
512 UWORD(51893625), UWORD(51893625),
513 UWORD(54438423), UWORD(54438423),
514 UWORD(54438423), UWORD(57065085),
515 UWORD(57065085), UWORD(57065085),
516 UWORD(59774907), UWORD(59774907),
517 UWORD(59774907), UWORD(62569185),
518 UWORD(62569185), UWORD(62569185),
519 UWORD(65449215), UWORD(65449215),
520 UWORD(65449215), UWORD(68416293),
521 UWORD(68416293), UWORD(68416293),
522 UWORD(71471715), UWORD(71471715),
523 UWORD(71471715), UWORD(74616777),
524 UWORD(74616777), UWORD(74616777),
525 UWORD(77852775), UWORD(77852775),
526 UWORD(77852775), UWORD(81181005),
527 UWORD(81181005), UWORD(81181005),
528 UWORD(84602763), UWORD(84602763),
529 UWORD(84602763), UWORD(88119345),
530 UWORD(88119345), UWORD(88119345),
531 UWORD(91732047), UWORD(91732047),
532 UWORD(91732047), UWORD(95442165),
533 UWORD(95442165), UWORD(95442165),
534 UWORD(99250995), UWORD(99250995),
535 UWORD(99250995), UWORD(103159833),
536 UWORD(103159833), UWORD(103159833),
537 UWORD(107169975), UWORD(107169975),
538 UWORD(107169975), UWORD(111282717),
539 UWORD(111282717), UWORD(111282717),
540 UWORD(115499355), UWORD(115499355),
541 UWORD(115499355), UWORD(119821185),
542 UWORD(119821185), UWORD(119821185),
543 UWORD(124249503), UWORD(124249503),
544 UWORD(124249503), UWORD(128785605),
545 UWORD(128785605), UWORD(128785605),
546 UWORD(133430787), UWORD(133430787),
547 };
548
549 #endif
550
_arb_atan_taylor_rs(mp_ptr y,mp_limb_t * error,mp_srcptr x,mp_size_t xn,ulong N,int alternating)551 void _arb_atan_taylor_rs(mp_ptr y, mp_limb_t * error,
552 mp_srcptr x, mp_size_t xn, ulong N, int alternating)
553 {
554 mp_ptr s, t, xpow;
555 mp_limb_t new_denom, old_denom, c;
556 slong power, k, m;
557
558 TMP_INIT;
559 TMP_START;
560
561 if (N >= ODD_RECIPROCAL_TAB_SIZE)
562 {
563 flint_printf("_arb_atan_taylor_rs: N too large!\n");
564 flint_abort();
565 }
566
567 if (N <= 2)
568 {
569 if (N == 0)
570 {
571 flint_mpn_zero(y, xn);
572 error[0] = 0;
573 }
574 else if (N == 1)
575 {
576 flint_mpn_copyi(y, x, xn);
577 error[0] = 0;
578 }
579 else
580 {
581 t = TMP_ALLOC_LIMBS(3 * xn);
582
583 /* x * (1 - x^2 / 3) */
584
585 /* higher index ---> */
586 /* t = | | x^2 (lo) | x^2 (hi) | */
587 mpn_sqr(t + xn, x, xn);
588
589 /* t = | x^3 (lo) | x^3 (hi) | x^2 (hi) | */
590 mpn_mul_n(t, t + 2 * xn, x, xn);
591
592 /* y = x - x^3 / 3 */
593 mpn_divrem_1(t, 0, t + xn, xn, 3);
594
595 if (alternating)
596 mpn_sub_n(y, x, t, xn);
597 else
598 mpn_add_n(y, x, t, xn);
599
600 error[0] = 3;
601 }
602 }
603 else
604 {
605 /* Choose m ~= sqrt(num_terms) (m must be even, >= 2) */
606 m = 2;
607 while (m * m < N)
608 m += 2;
609
610 /* todo: merge allocations */
611 xpow = TMP_ALLOC_LIMBS((m + 1) * xn);
612 s = TMP_ALLOC_LIMBS(xn + 2);
613 t = TMP_ALLOC_LIMBS(2 * xn + 2); /* todo: 1 limb too much? */
614
615 /* higher index ---> */
616 /* | ---xn--- | */
617 /* xpow = | <temp> | (x^2)^m | (x^2)^(m-1) | ... | (x^2)^2 | x^2 | */
618
619 #define XPOW_WRITE(__k) (xpow + (m - (__k)) * xn)
620 #define XPOW_READ(__k) (xpow + (m - (__k) + 1) * xn)
621
622 mpn_sqr(XPOW_WRITE(1), x, xn);
623 mpn_sqr(XPOW_WRITE(2), XPOW_READ(1), xn);
624
625 for (k = 4; k <= m; k += 2)
626 {
627 mpn_mul_n(XPOW_WRITE(k - 1), XPOW_READ(k / 2), XPOW_READ(k / 2 - 1), xn);
628 mpn_sqr(XPOW_WRITE(k), XPOW_READ(k / 2), xn);
629 }
630
631 flint_mpn_zero(s, xn + 1);
632
633 /* todo: skip one nonscalar multiplication (use x^m)
634 when starting on x^0 */
635 power = (N - 1) % m;
636
637 for (k = N - 1; k >= 0; k--)
638 {
639 c = odd_reciprocal_tab_numer[k];
640 new_denom = odd_reciprocal_tab_denom[k];
641 old_denom = odd_reciprocal_tab_denom[k+1];
642
643 /* change denominators */
644 if (new_denom != old_denom && k < N - 1)
645 {
646 /* hack when s is negative: add 1 to get a positive number */
647 if (alternating && (k % 2 == 0))
648 s[xn] += old_denom;
649
650 /* multiply by new denominator */
651 s[xn + 1] = mpn_mul_1(s, s, xn + 1, new_denom);
652 /* divide by old denominator */
653 mpn_divrem_1(s, 0, s, xn + 2, old_denom);
654
655 if (s[xn + 1] != 0)
656 {
657 flint_printf("bad division!\n");
658 flint_abort();
659 }
660
661 /* subtract 1 */
662 if (alternating && (k % 2 == 0))
663 s[xn] -= new_denom;
664 }
665
666 if (power == 0)
667 {
668 /* sub/add c * x^0 -- only top limb is affected */
669 if (alternating & k)
670 s[xn] -= c;
671 else
672 s[xn] += c;
673
674 /* Outer polynomial evaluation: multiply by (x^2)^m */
675 if (k != 0)
676 {
677 mpn_mul(t, s, xn + 1, XPOW_READ(m), xn);
678 flint_mpn_copyi(s, t + xn, xn + 1);
679 }
680
681 power = m - 1;
682 }
683 else
684 {
685 if (alternating & k)
686 s[xn] -= mpn_submul_1(s, XPOW_READ(power), xn, c);
687 else
688 s[xn] += mpn_addmul_1(s, XPOW_READ(power), xn, c);
689
690 power--;
691 }
692 }
693
694 /* finally divide by denominator and multiply by x */
695 mpn_divrem_1(s, 0, s, xn + 1, odd_reciprocal_tab_denom[0]);
696 mpn_mul(t, s, xn + 1, x, xn);
697 flint_mpn_copyi(y, t + xn, xn);
698
699 /* error bound (ulp) */
700 error[0] = 2;
701 }
702
703 TMP_END;
704 }
705
706