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