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 FACTORIAL_TAB_SIZE 288
20 
21 #if FLINT_BITS == 64
22 
23 ARB_DLL const mp_limb_t factorial_tab_numer[FACTORIAL_TAB_SIZE] = {
24     UWORD(2432902008176640000),
25     UWORD(2432902008176640000),
26     UWORD(1216451004088320000),
27     UWORD(405483668029440000),
28     UWORD(101370917007360000),
29     UWORD(20274183401472000),
30     UWORD(3379030566912000),
31     UWORD(482718652416000),
32     UWORD(60339831552000),
33     UWORD(6704425728000),
34     UWORD(670442572800),
35     UWORD(60949324800),
36     UWORD(5079110400),
37     UWORD(390700800),
38     UWORD(27907200),
39     UWORD(1860480),
40     UWORD(116280),
41     UWORD(6840),
42     UWORD(380),
43     UWORD(20),
44     UWORD(1),
45     UWORD(169958063987712000),
46     UWORD(7725366544896000),
47     UWORD(335885501952000),
48     UWORD(13995229248000),
49     UWORD(559809169920),
50     UWORD(21531121920),
51     UWORD(797448960),
52     UWORD(28480320),
53     UWORD(982080),
54     UWORD(32736),
55     UWORD(1056),
56     UWORD(33),
57     UWORD(1),
58     UWORD(405179306820288000),
59     UWORD(11576551623436800),
60     UWORD(321570878428800),
61     UWORD(8691104822400),
62     UWORD(228713284800),
63     UWORD(5864443200),
64     UWORD(146611080),
65     UWORD(3575880),
66     UWORD(85140),
67     UWORD(1980),
68     UWORD(45),
69     UWORD(1),
70     UWORD(129210868410624000),
71     UWORD(2749167412992000),
72     UWORD(57274321104000),
73     UWORD(1168863696000),
74     UWORD(23377273920),
75     UWORD(458377920),
76     UWORD(8814960),
77     UWORD(166320),
78     UWORD(3080),
79     UWORD(56),
80     UWORD(1),
81     UWORD(13431688016947200),
82     UWORD(231580827878400),
83     UWORD(3925098777600),
84     UWORD(65418312960),
85     UWORD(1072431360),
86     UWORD(17297280),
87     UWORD(274560),
88     UWORD(4290),
89     UWORD(66),
90     UWORD(1),
91     UWORD(51698307350592000),
92     UWORD(760269225744000),
93     UWORD(11018394576000),
94     UWORD(157405636800),
95     UWORD(2216980800),
96     UWORD(30791400),
97     UWORD(421800),
98     UWORD(5700),
99     UWORD(76),
100     UWORD(1),
101     UWORD(166872705665702400),
102     UWORD(2139393662380800),
103     UWORD(27080932435200),
104     UWORD(338511655440),
105     UWORD(4179156240),
106     UWORD(50965320),
107     UWORD(614040),
108     UWORD(7310),
109     UWORD(86),
110     UWORD(1),
111     UWORD(4900933563926400),
112     UWORD(55692426862800),
113     UWORD(625757605200),
114     UWORD(6952862280),
115     UWORD(76405080),
116     UWORD(830490),
117     UWORD(8930),
118     UWORD(95),
119     UWORD(1),
120     UWORD(10385445095625600),
121     UWORD(107066444284800),
122     UWORD(1092514737600),
123     UWORD(11035502400),
124     UWORD(110355024),
125     UWORD(1092624),
126     UWORD(10712),
127     UWORD(104),
128     UWORD(1),
129     UWORD(20632508204394240),
130     UWORD(194646303815040),
131     UWORD(1819124334720),
132     UWORD(16843743840),
133     UWORD(154529760),
134     UWORD(1404816),
135     UWORD(12656),
136     UWORD(113),
137     UWORD(1),
138     UWORD(38823716505974400),
139     UWORD(337597534834560),
140     UWORD(2910323576160),
141     UWORD(24874560480),
142     UWORD(210801360),
143     UWORD(1771440),
144     UWORD(14762),
145     UWORD(122),
146     UWORD(1),
147     UWORD(69746158460160000),
148     UWORD(562469019840000),
149     UWORD(4499752158720),
150     UWORD(35712318720),
151     UWORD(281199360),
152     UWORD(2196870),
153     UWORD(17030),
154     UWORD(131),
155     UWORD(1),
156     UWORD(120384786168259200),
157     UWORD(905148768182400),
158     UWORD(6754841553600),
159     UWORD(50035863360),
160     UWORD(367910760),
161     UWORD(2685480),
162     UWORD(19460),
163     UWORD(140),
164     UWORD(1),
165     UWORD(1346749373249280),
166     UWORD(9484150515840),
167     UWORD(66322730880),
168     UWORD(460574520),
169     UWORD(3176376),
170     UWORD(21756),
171     UWORD(148),
172     UWORD(1),
173     UWORD(1961463552048000),
174     UWORD(13076423680320),
175     UWORD(86598832320),
176     UWORD(569729160),
177     UWORD(3723720),
178     UWORD(24180),
179     UWORD(156),
180     UWORD(1),
181     UWORD(2802505908948480),
182     UWORD(17737379170560),
183     UWORD(111555843840),
184     UWORD(697224024),
185     UWORD(4330584),
186     UWORD(26732),
187     UWORD(164),
188     UWORD(1),
189     UWORD(3935446603320960),
190     UWORD(23707509658560),
191     UWORD(141961135680),
192     UWORD(845006760),
193     UWORD(5000040),
194     UWORD(29412),
195     UWORD(172),
196     UWORD(1),
197     UWORD(5440248896544000),
198     UWORD(31265798256000),
199     UWORD(178661704320),
200     UWORD(1015123320),
201     UWORD(5735160),
202     UWORD(32220),
203     UWORD(180),
204     UWORD(1),
205     UWORD(7413519413139840),
206     UWORD(40733623149120),
207     UWORD(222588104640),
208     UWORD(1209717960),
209     UWORD(6539016),
210     UWORD(35156),
211     UWORD(188),
212     UWORD(1),
213     UWORD(9970981685683200),
214     UWORD(52478850977280),
215     UWORD(274758382080),
216     UWORD(1431033240),
217     UWORD(7414680),
218     UWORD(38220),
219     UWORD(196),
220     UWORD(1),
221     UWORD(13250183553129600),
222     UWORD(66920118955200),
223     UWORD(336282004800),
224     UWORD(1681410024),
225     UWORD(8365224),
226     UWORD(41412),
227     UWORD(204),
228     UWORD(1),
229     UWORD(17413448982209280),
230     UWORD(84531305738880),
231     UWORD(408363795840),
232     UWORD(1963287480),
233     UWORD(9393720),
234     UWORD(44732),
235     UWORD(212),
236     UWORD(1),
237     UWORD(22651084881532800),
238     UWORD(105846191035200),
239     UWORD(492307865280),
240     UWORD(2279203080),
241     UWORD(10503240),
242     UWORD(48180),
243     UWORD(220),
244     UWORD(1),
245     UWORD(29184853478054400),
246     UWORD(131463303955200),
247     UWORD(589521542400),
248     UWORD(2631792600),
249     UWORD(11696856),
250     UWORD(51756),
251     UWORD(228),
252     UWORD(1),
253     UWORD(37271720825539200),
254     UWORD(162050960111040),
255     UWORD(701519307840),
256     UWORD(3023790120),
257     UWORD(12977640),
258     UWORD(55460),
259     UWORD(236),
260     UWORD(1),
261     UWORD(47207892014680320),
262     UWORD(198352487456640),
263     UWORD(829926725760),
264     UWORD(3458028024),
265     UWORD(14348664),
266     UWORD(59292),
267     UWORD(244),
268     UWORD(1),
269     UWORD(59333143654512000),
270     UWORD(241191640872000),
271     UWORD(976484376000),
272     UWORD(3937437000),
273     UWORD(15813000),
274     UWORD(63252),
275     UWORD(252),
276     UWORD(1),
277     UWORD(284751785364480),
278     UWORD(1121070021120),
279     UWORD(4396353024),
280     UWORD(17173254),
281     UWORD(66822),
282     UWORD(259),
283     UWORD(1),
284     UWORD(334679402201760),
285     UWORD(1282296560160),
286     UWORD(4894261680),
287     UWORD(18609360),
288     UWORD(70490),
289     UWORD(266),
290     UWORD(1),
291     UWORD(391698294099840),
292     UWORD(1461560798880),
293     UWORD(5433311520),
294     UWORD(20123376),
295     UWORD(74256),
296     UWORD(273),
297     UWORD(1),
298     UWORD(456592291848000),
299     UWORD(1660335606720),
300     UWORD(6015708720),
301     UWORD(21717360),
302     UWORD(78120),
303     UWORD(280),
304     UWORD(1),
305     UWORD(530208485286480),
306     UWORD(1880171933640),
307     UWORD(6643717080),
308     UWORD(23393370),
309     UWORD(82082),
310     UWORD(287),
311     UWORD(1),
312 };
313 
314 ARB_DLL const mp_limb_t factorial_tab_denom[FACTORIAL_TAB_SIZE] = {
315     UWORD(2432902008176640000),
316     UWORD(2432902008176640000),
317     UWORD(2432902008176640000),
318     UWORD(2432902008176640000),
319     UWORD(2432902008176640000),
320     UWORD(2432902008176640000),
321     UWORD(2432902008176640000),
322     UWORD(2432902008176640000),
323     UWORD(2432902008176640000),
324     UWORD(2432902008176640000),
325     UWORD(2432902008176640000),
326     UWORD(2432902008176640000),
327     UWORD(2432902008176640000),
328     UWORD(2432902008176640000),
329     UWORD(2432902008176640000),
330     UWORD(2432902008176640000),
331     UWORD(2432902008176640000),
332     UWORD(2432902008176640000),
333     UWORD(2432902008176640000),
334     UWORD(2432902008176640000),
335     UWORD(2432902008176640000),
336     UWORD(3569119343741952000),
337     UWORD(3569119343741952000),
338     UWORD(3569119343741952000),
339     UWORD(3569119343741952000),
340     UWORD(3569119343741952000),
341     UWORD(3569119343741952000),
342     UWORD(3569119343741952000),
343     UWORD(3569119343741952000),
344     UWORD(3569119343741952000),
345     UWORD(3569119343741952000),
346     UWORD(3569119343741952000),
347     UWORD(3569119343741952000),
348     UWORD(3569119343741952000),
349     UWORD(13776096431889792000),
350     UWORD(13776096431889792000),
351     UWORD(13776096431889792000),
352     UWORD(13776096431889792000),
353     UWORD(13776096431889792000),
354     UWORD(13776096431889792000),
355     UWORD(13776096431889792000),
356     UWORD(13776096431889792000),
357     UWORD(13776096431889792000),
358     UWORD(13776096431889792000),
359     UWORD(13776096431889792000),
360     UWORD(13776096431889792000),
361     UWORD(5943699946888704000),
362     UWORD(5943699946888704000),
363     UWORD(5943699946888704000),
364     UWORD(5943699946888704000),
365     UWORD(5943699946888704000),
366     UWORD(5943699946888704000),
367     UWORD(5943699946888704000),
368     UWORD(5943699946888704000),
369     UWORD(5943699946888704000),
370     UWORD(5943699946888704000),
371     UWORD(5943699946888704000),
372     UWORD(765606216965990400),
373     UWORD(765606216965990400),
374     UWORD(765606216965990400),
375     UWORD(765606216965990400),
376     UWORD(765606216965990400),
377     UWORD(765606216965990400),
378     UWORD(765606216965990400),
379     UWORD(765606216965990400),
380     UWORD(765606216965990400),
381     UWORD(765606216965990400),
382     UWORD(3463786592489664000),
383     UWORD(3463786592489664000),
384     UWORD(3463786592489664000),
385     UWORD(3463786592489664000),
386     UWORD(3463786592489664000),
387     UWORD(3463786592489664000),
388     UWORD(3463786592489664000),
389     UWORD(3463786592489664000),
390     UWORD(3463786592489664000),
391     UWORD(3463786592489664000),
392     UWORD(12849198336259084800),
393     UWORD(12849198336259084800),
394     UWORD(12849198336259084800),
395     UWORD(12849198336259084800),
396     UWORD(12849198336259084800),
397     UWORD(12849198336259084800),
398     UWORD(12849198336259084800),
399     UWORD(12849198336259084800),
400     UWORD(12849198336259084800),
401     UWORD(12849198336259084800),
402     UWORD(426381220061596800),
403     UWORD(426381220061596800),
404     UWORD(426381220061596800),
405     UWORD(426381220061596800),
406     UWORD(426381220061596800),
407     UWORD(426381220061596800),
408     UWORD(426381220061596800),
409     UWORD(426381220061596800),
410     UWORD(426381220061596800),
411     UWORD(997002729180057600),
412     UWORD(997002729180057600),
413     UWORD(997002729180057600),
414     UWORD(997002729180057600),
415     UWORD(997002729180057600),
416     UWORD(997002729180057600),
417     UWORD(997002729180057600),
418     UWORD(997002729180057600),
419     UWORD(997002729180057600),
420     UWORD(2166413361461395200),
421     UWORD(2166413361461395200),
422     UWORD(2166413361461395200),
423     UWORD(2166413361461395200),
424     UWORD(2166413361461395200),
425     UWORD(2166413361461395200),
426     UWORD(2166413361461395200),
427     UWORD(2166413361461395200),
428     UWORD(2166413361461395200),
429     UWORD(4425903681681081600),
430     UWORD(4425903681681081600),
431     UWORD(4425903681681081600),
432     UWORD(4425903681681081600),
433     UWORD(4425903681681081600),
434     UWORD(4425903681681081600),
435     UWORD(4425903681681081600),
436     UWORD(4425903681681081600),
437     UWORD(4425903681681081600),
438     UWORD(8578777490599680000),
439     UWORD(8578777490599680000),
440     UWORD(8578777490599680000),
441     UWORD(8578777490599680000),
442     UWORD(8578777490599680000),
443     UWORD(8578777490599680000),
444     UWORD(8578777490599680000),
445     UWORD(8578777490599680000),
446     UWORD(8578777490599680000),
447     UWORD(15890791774210214400),
448     UWORD(15890791774210214400),
449     UWORD(15890791774210214400),
450     UWORD(15890791774210214400),
451     UWORD(15890791774210214400),
452     UWORD(15890791774210214400),
453     UWORD(15890791774210214400),
454     UWORD(15890791774210214400),
455     UWORD(15890791774210214400),
456     UWORD(189891661628148480),
457     UWORD(189891661628148480),
458     UWORD(189891661628148480),
459     UWORD(189891661628148480),
460     UWORD(189891661628148480),
461     UWORD(189891661628148480),
462     UWORD(189891661628148480),
463     UWORD(189891661628148480),
464     UWORD(292258069255152000),
465     UWORD(292258069255152000),
466     UWORD(292258069255152000),
467     UWORD(292258069255152000),
468     UWORD(292258069255152000),
469     UWORD(292258069255152000),
470     UWORD(292258069255152000),
471     UWORD(292258069255152000),
472     UWORD(439993427704911360),
473     UWORD(439993427704911360),
474     UWORD(439993427704911360),
475     UWORD(439993427704911360),
476     UWORD(439993427704911360),
477     UWORD(439993427704911360),
478     UWORD(439993427704911360),
479     UWORD(439993427704911360),
480     UWORD(649348689547958400),
481     UWORD(649348689547958400),
482     UWORD(649348689547958400),
483     UWORD(649348689547958400),
484     UWORD(649348689547958400),
485     UWORD(649348689547958400),
486     UWORD(649348689547958400),
487     UWORD(649348689547958400),
488     UWORD(941163059102112000),
489     UWORD(941163059102112000),
490     UWORD(941163059102112000),
491     UWORD(941163059102112000),
492     UWORD(941163059102112000),
493     UWORD(941163059102112000),
494     UWORD(941163059102112000),
495     UWORD(941163059102112000),
496     UWORD(1341847013778311040),
497     UWORD(1341847013778311040),
498     UWORD(1341847013778311040),
499     UWORD(1341847013778311040),
500     UWORD(1341847013778311040),
501     UWORD(1341847013778311040),
502     UWORD(1341847013778311040),
503     UWORD(1341847013778311040),
504     UWORD(1884515538594124800),
505     UWORD(1884515538594124800),
506     UWORD(1884515538594124800),
507     UWORD(1884515538594124800),
508     UWORD(1884515538594124800),
509     UWORD(1884515538594124800),
510     UWORD(1884515538594124800),
511     UWORD(1884515538594124800),
512     UWORD(2610286159966531200),
513     UWORD(2610286159966531200),
514     UWORD(2610286159966531200),
515     UWORD(2610286159966531200),
516     UWORD(2610286159966531200),
517     UWORD(2610286159966531200),
518     UWORD(2610286159966531200),
519     UWORD(2610286159966531200),
520     UWORD(3569757041352902400),
521     UWORD(3569757041352902400),
522     UWORD(3569757041352902400),
523     UWORD(3569757041352902400),
524     UWORD(3569757041352902400),
525     UWORD(3569757041352902400),
526     UWORD(3569757041352902400),
527     UWORD(3569757041352902400),
528     UWORD(4824681079766486400),
529     UWORD(4824681079766486400),
530     UWORD(4824681079766486400),
531     UWORD(4824681079766486400),
532     UWORD(4824681079766486400),
533     UWORD(4824681079766486400),
534     UWORD(4824681079766486400),
535     UWORD(4824681079766486400),
536     UWORD(6449852618650022400),
537     UWORD(6449852618650022400),
538     UWORD(6449852618650022400),
539     UWORD(6449852618650022400),
540     UWORD(6449852618650022400),
541     UWORD(6449852618650022400),
542     UWORD(6449852618650022400),
543     UWORD(6449852618650022400),
544     UWORD(8535224069048476800),
545     UWORD(8535224069048476800),
546     UWORD(8535224069048476800),
547     UWORD(8535224069048476800),
548     UWORD(8535224069048476800),
549     UWORD(8535224069048476800),
550     UWORD(8535224069048476800),
551     UWORD(8535224069048476800),
552     UWORD(11188270407479235840),
553     UWORD(11188270407479235840),
554     UWORD(11188270407479235840),
555     UWORD(11188270407479235840),
556     UWORD(11188270407479235840),
557     UWORD(11188270407479235840),
558     UWORD(11188270407479235840),
559     UWORD(11188270407479235840),
560     UWORD(14536620195355440000),
561     UWORD(14536620195355440000),
562     UWORD(14536620195355440000),
563     UWORD(14536620195355440000),
564     UWORD(14536620195355440000),
565     UWORD(14536620195355440000),
566     UWORD(14536620195355440000),
567     UWORD(14536620195355440000),
568     UWORD(72042201697213440),
569     UWORD(72042201697213440),
570     UWORD(72042201697213440),
571     UWORD(72042201697213440),
572     UWORD(72042201697213440),
573     UWORD(72042201697213440),
574     UWORD(72042201697213440),
575     UWORD(87016644572457600),
576     UWORD(87016644572457600),
577     UWORD(87016644572457600),
578     UWORD(87016644572457600),
579     UWORD(87016644572457600),
580     UWORD(87016644572457600),
581     UWORD(87016644572457600),
582     UWORD(104583444524657280),
583     UWORD(104583444524657280),
584     UWORD(104583444524657280),
585     UWORD(104583444524657280),
586     UWORD(104583444524657280),
587     UWORD(104583444524657280),
588     UWORD(104583444524657280),
589     UWORD(125106287966352000),
590     UWORD(125106287966352000),
591     UWORD(125106287966352000),
592     UWORD(125106287966352000),
593     UWORD(125106287966352000),
594     UWORD(125106287966352000),
595     UWORD(125106287966352000),
596     UWORD(148988584365500880),
597     UWORD(148988584365500880),
598     UWORD(148988584365500880),
599     UWORD(148988584365500880),
600     UWORD(148988584365500880),
601     UWORD(148988584365500880),
602     UWORD(148988584365500880),
603 };
604 
605 #else
606 
607 ARB_DLL const mp_limb_t factorial_tab_numer[FACTORIAL_TAB_SIZE] = {
608     UWORD(479001600),
609     UWORD(479001600),
610     UWORD(239500800),
611     UWORD(79833600),
612     UWORD(19958400),
613     UWORD(3991680),
614     UWORD(665280),
615     UWORD(95040),
616     UWORD(11880),
617     UWORD(1320),
618     UWORD(132),
619     UWORD(12),
620     UWORD(1),
621     UWORD(19535040),
622     UWORD(1395360),
623     UWORD(93024),
624     UWORD(5814),
625     UWORD(342),
626     UWORD(19),
627     UWORD(1),
628     UWORD(165765600),
629     UWORD(7893600),
630     UWORD(358800),
631     UWORD(15600),
632     UWORD(650),
633     UWORD(26),
634     UWORD(1),
635     UWORD(24165120),
636     UWORD(863040),
637     UWORD(29760),
638     UWORD(992),
639     UWORD(32),
640     UWORD(1),
641     UWORD(60233040),
642     UWORD(1771560),
643     UWORD(50616),
644     UWORD(1406),
645     UWORD(38),
646     UWORD(1),
647     UWORD(2961840),
648     UWORD(74046),
649     UWORD(1806),
650     UWORD(43),
651     UWORD(1),
652     UWORD(4669920),
653     UWORD(103776),
654     UWORD(2256),
655     UWORD(48),
656     UWORD(1),
657     UWORD(7027800),
658     UWORD(140556),
659     UWORD(2756),
660     UWORD(53),
661     UWORD(1),
662     UWORD(10182480),
663     UWORD(185136),
664     UWORD(3306),
665     UWORD(58),
666     UWORD(1),
667     UWORD(14295960),
668     UWORD(238266),
669     UWORD(3906),
670     UWORD(63),
671     UWORD(1),
672     UWORD(19545240),
673     UWORD(300696),
674     UWORD(4556),
675     UWORD(68),
676     UWORD(1),
677     UWORD(26122320),
678     UWORD(373176),
679     UWORD(5256),
680     UWORD(73),
681     UWORD(1),
682     UWORD(34234200),
683     UWORD(456456),
684     UWORD(6006),
685     UWORD(78),
686     UWORD(1),
687     UWORD(44102880),
688     UWORD(551286),
689     UWORD(6806),
690     UWORD(83),
691     UWORD(1),
692     UWORD(635970),
693     UWORD(7482),
694     UWORD(87),
695     UWORD(1),
696     UWORD(728910),
697     UWORD(8190),
698     UWORD(91),
699     UWORD(1),
700     UWORD(830490),
701     UWORD(8930),
702     UWORD(95),
703     UWORD(1),
704     UWORD(941094),
705     UWORD(9702),
706     UWORD(99),
707     UWORD(1),
708     UWORD(1061106),
709     UWORD(10506),
710     UWORD(103),
711     UWORD(1),
712     UWORD(1190910),
713     UWORD(11342),
714     UWORD(107),
715     UWORD(1),
716     UWORD(1330890),
717     UWORD(12210),
718     UWORD(111),
719     UWORD(1),
720     UWORD(1481430),
721     UWORD(13110),
722     UWORD(115),
723     UWORD(1),
724     UWORD(1642914),
725     UWORD(14042),
726     UWORD(119),
727     UWORD(1),
728     UWORD(1815726),
729     UWORD(15006),
730     UWORD(123),
731     UWORD(1),
732     UWORD(2000250),
733     UWORD(16002),
734     UWORD(127),
735     UWORD(1),
736     UWORD(2196870),
737     UWORD(17030),
738     UWORD(131),
739     UWORD(1),
740     UWORD(2405970),
741     UWORD(18090),
742     UWORD(135),
743     UWORD(1),
744     UWORD(2627934),
745     UWORD(19182),
746     UWORD(139),
747     UWORD(1),
748     UWORD(2863146),
749     UWORD(20306),
750     UWORD(143),
751     UWORD(1),
752     UWORD(3111990),
753     UWORD(21462),
754     UWORD(147),
755     UWORD(1),
756     UWORD(3374850),
757     UWORD(22650),
758     UWORD(151),
759     UWORD(1),
760     UWORD(3652110),
761     UWORD(23870),
762     UWORD(155),
763     UWORD(1),
764     UWORD(3944154),
765     UWORD(25122),
766     UWORD(159),
767     UWORD(1),
768     UWORD(4251366),
769     UWORD(26406),
770     UWORD(163),
771     UWORD(1),
772     UWORD(4574130),
773     UWORD(27722),
774     UWORD(167),
775     UWORD(1),
776     UWORD(4912830),
777     UWORD(29070),
778     UWORD(171),
779     UWORD(1),
780     UWORD(5267850),
781     UWORD(30450),
782     UWORD(175),
783     UWORD(1),
784     UWORD(5639574),
785     UWORD(31862),
786     UWORD(179),
787     UWORD(1),
788     UWORD(6028386),
789     UWORD(33306),
790     UWORD(183),
791     UWORD(1),
792     UWORD(6434670),
793     UWORD(34782),
794     UWORD(187),
795     UWORD(1),
796     UWORD(6858810),
797     UWORD(36290),
798     UWORD(191),
799     UWORD(1),
800     UWORD(7301190),
801     UWORD(37830),
802     UWORD(195),
803     UWORD(1),
804     UWORD(7762194),
805     UWORD(39402),
806     UWORD(199),
807     UWORD(1),
808     UWORD(8242206),
809     UWORD(41006),
810     UWORD(203),
811     UWORD(1),
812     UWORD(8741610),
813     UWORD(42642),
814     UWORD(207),
815     UWORD(1),
816     UWORD(9260790),
817     UWORD(44310),
818     UWORD(211),
819     UWORD(1),
820     UWORD(9800130),
821     UWORD(46010),
822     UWORD(215),
823     UWORD(1),
824     UWORD(10360014),
825     UWORD(47742),
826     UWORD(219),
827     UWORD(1),
828     UWORD(10940826),
829     UWORD(49506),
830     UWORD(223),
831     UWORD(1),
832     UWORD(11542950),
833     UWORD(51302),
834     UWORD(227),
835     UWORD(1),
836     UWORD(12166770),
837     UWORD(53130),
838     UWORD(231),
839     UWORD(1),
840     UWORD(12812670),
841     UWORD(54990),
842     UWORD(235),
843     UWORD(1),
844     UWORD(13481034),
845     UWORD(56882),
846     UWORD(239),
847     UWORD(1),
848     UWORD(14172246),
849     UWORD(58806),
850     UWORD(243),
851     UWORD(1),
852     UWORD(14886690),
853     UWORD(60762),
854     UWORD(247),
855     UWORD(1),
856     UWORD(15624750),
857     UWORD(62750),
858     UWORD(251),
859     UWORD(1),
860     UWORD(16386810),
861     UWORD(64770),
862     UWORD(255),
863     UWORD(1),
864     UWORD(66306),
865     UWORD(258),
866     UWORD(1),
867     UWORD(67860),
868     UWORD(261),
869     UWORD(1),
870     UWORD(69432),
871     UWORD(264),
872     UWORD(1),
873     UWORD(71022),
874     UWORD(267),
875     UWORD(1),
876     UWORD(72630),
877     UWORD(270),
878     UWORD(1),
879     UWORD(74256),
880     UWORD(273),
881     UWORD(1),
882     UWORD(75900),
883     UWORD(276),
884     UWORD(1),
885     UWORD(77562),
886     UWORD(279),
887     UWORD(1),
888     UWORD(79242),
889     UWORD(282),
890     UWORD(1),
891     UWORD(80940),
892     UWORD(285),
893     UWORD(1),
894     UWORD(82656),
895     UWORD(288),
896 };
897 
898 ARB_DLL const mp_limb_t factorial_tab_denom[FACTORIAL_TAB_SIZE] = {
899     UWORD(479001600),
900     UWORD(479001600),
901     UWORD(479001600),
902     UWORD(479001600),
903     UWORD(479001600),
904     UWORD(479001600),
905     UWORD(479001600),
906     UWORD(479001600),
907     UWORD(479001600),
908     UWORD(479001600),
909     UWORD(479001600),
910     UWORD(479001600),
911     UWORD(479001600),
912     UWORD(253955520),
913     UWORD(253955520),
914     UWORD(253955520),
915     UWORD(253955520),
916     UWORD(253955520),
917     UWORD(253955520),
918     UWORD(253955520),
919     UWORD(3315312000),
920     UWORD(3315312000),
921     UWORD(3315312000),
922     UWORD(3315312000),
923     UWORD(3315312000),
924     UWORD(3315312000),
925     UWORD(3315312000),
926     UWORD(652458240),
927     UWORD(652458240),
928     UWORD(652458240),
929     UWORD(652458240),
930     UWORD(652458240),
931     UWORD(652458240),
932     UWORD(1987690320),
933     UWORD(1987690320),
934     UWORD(1987690320),
935     UWORD(1987690320),
936     UWORD(1987690320),
937     UWORD(1987690320),
938     UWORD(115511760),
939     UWORD(115511760),
940     UWORD(115511760),
941     UWORD(115511760),
942     UWORD(115511760),
943     UWORD(205476480),
944     UWORD(205476480),
945     UWORD(205476480),
946     UWORD(205476480),
947     UWORD(205476480),
948     UWORD(344362200),
949     UWORD(344362200),
950     UWORD(344362200),
951     UWORD(344362200),
952     UWORD(344362200),
953     UWORD(549853920),
954     UWORD(549853920),
955     UWORD(549853920),
956     UWORD(549853920),
957     UWORD(549853920),
958     UWORD(843461640),
959     UWORD(843461640),
960     UWORD(843461640),
961     UWORD(843461640),
962     UWORD(843461640),
963     UWORD(1250895360),
964     UWORD(1250895360),
965     UWORD(1250895360),
966     UWORD(1250895360),
967     UWORD(1250895360),
968     UWORD(1802440080),
969     UWORD(1802440080),
970     UWORD(1802440080),
971     UWORD(1802440080),
972     UWORD(1802440080),
973     UWORD(2533330800),
974     UWORD(2533330800),
975     UWORD(2533330800),
976     UWORD(2533330800),
977     UWORD(2533330800),
978     UWORD(3484127520),
979     UWORD(3484127520),
980     UWORD(3484127520),
981     UWORD(3484127520),
982     UWORD(3484127520),
983     UWORD(53421480),
984     UWORD(53421480),
985     UWORD(53421480),
986     UWORD(53421480),
987     UWORD(64144080),
988     UWORD(64144080),
989     UWORD(64144080),
990     UWORD(64144080),
991     UWORD(76405080),
992     UWORD(76405080),
993     UWORD(76405080),
994     UWORD(76405080),
995     UWORD(90345024),
996     UWORD(90345024),
997     UWORD(90345024),
998     UWORD(90345024),
999     UWORD(106110600),
1000     UWORD(106110600),
1001     UWORD(106110600),
1002     UWORD(106110600),
1003     UWORD(123854640),
1004     UWORD(123854640),
1005     UWORD(123854640),
1006     UWORD(123854640),
1007     UWORD(143736120),
1008     UWORD(143736120),
1009     UWORD(143736120),
1010     UWORD(143736120),
1011     UWORD(165920160),
1012     UWORD(165920160),
1013     UWORD(165920160),
1014     UWORD(165920160),
1015     UWORD(190578024),
1016     UWORD(190578024),
1017     UWORD(190578024),
1018     UWORD(190578024),
1019     UWORD(217887120),
1020     UWORD(217887120),
1021     UWORD(217887120),
1022     UWORD(217887120),
1023     UWORD(248031000),
1024     UWORD(248031000),
1025     UWORD(248031000),
1026     UWORD(248031000),
1027     UWORD(281199360),
1028     UWORD(281199360),
1029     UWORD(281199360),
1030     UWORD(281199360),
1031     UWORD(317588040),
1032     UWORD(317588040),
1033     UWORD(317588040),
1034     UWORD(317588040),
1035     UWORD(357399024),
1036     UWORD(357399024),
1037     UWORD(357399024),
1038     UWORD(357399024),
1039     UWORD(400840440),
1040     UWORD(400840440),
1041     UWORD(400840440),
1042     UWORD(400840440),
1043     UWORD(448126560),
1044     UWORD(448126560),
1045     UWORD(448126560),
1046     UWORD(448126560),
1047     UWORD(499477800),
1048     UWORD(499477800),
1049     UWORD(499477800),
1050     UWORD(499477800),
1051     UWORD(555120720),
1052     UWORD(555120720),
1053     UWORD(555120720),
1054     UWORD(555120720),
1055     UWORD(615288024),
1056     UWORD(615288024),
1057     UWORD(615288024),
1058     UWORD(615288024),
1059     UWORD(680218560),
1060     UWORD(680218560),
1061     UWORD(680218560),
1062     UWORD(680218560),
1063     UWORD(750157320),
1064     UWORD(750157320),
1065     UWORD(750157320),
1066     UWORD(750157320),
1067     UWORD(825355440),
1068     UWORD(825355440),
1069     UWORD(825355440),
1070     UWORD(825355440),
1071     UWORD(906070200),
1072     UWORD(906070200),
1073     UWORD(906070200),
1074     UWORD(906070200),
1075     UWORD(992565024),
1076     UWORD(992565024),
1077     UWORD(992565024),
1078     UWORD(992565024),
1079     UWORD(1085109480),
1080     UWORD(1085109480),
1081     UWORD(1085109480),
1082     UWORD(1085109480),
1083     UWORD(1183979280),
1084     UWORD(1183979280),
1085     UWORD(1183979280),
1086     UWORD(1183979280),
1087     UWORD(1289456280),
1088     UWORD(1289456280),
1089     UWORD(1289456280),
1090     UWORD(1289456280),
1091     UWORD(1401828480),
1092     UWORD(1401828480),
1093     UWORD(1401828480),
1094     UWORD(1401828480),
1095     UWORD(1521390024),
1096     UWORD(1521390024),
1097     UWORD(1521390024),
1098     UWORD(1521390024),
1099     UWORD(1648441200),
1100     UWORD(1648441200),
1101     UWORD(1648441200),
1102     UWORD(1648441200),
1103     UWORD(1783288440),
1104     UWORD(1783288440),
1105     UWORD(1783288440),
1106     UWORD(1783288440),
1107     UWORD(1926244320),
1108     UWORD(1926244320),
1109     UWORD(1926244320),
1110     UWORD(1926244320),
1111     UWORD(2077627560),
1112     UWORD(2077627560),
1113     UWORD(2077627560),
1114     UWORD(2077627560),
1115     UWORD(2237763024),
1116     UWORD(2237763024),
1117     UWORD(2237763024),
1118     UWORD(2237763024),
1119     UWORD(2406981720),
1120     UWORD(2406981720),
1121     UWORD(2406981720),
1122     UWORD(2406981720),
1123     UWORD(2585620800),
1124     UWORD(2585620800),
1125     UWORD(2585620800),
1126     UWORD(2585620800),
1127     UWORD(2774023560),
1128     UWORD(2774023560),
1129     UWORD(2774023560),
1130     UWORD(2774023560),
1131     UWORD(2972539440),
1132     UWORD(2972539440),
1133     UWORD(2972539440),
1134     UWORD(2972539440),
1135     UWORD(3181524024),
1136     UWORD(3181524024),
1137     UWORD(3181524024),
1138     UWORD(3181524024),
1139     UWORD(3401339040),
1140     UWORD(3401339040),
1141     UWORD(3401339040),
1142     UWORD(3401339040),
1143     UWORD(3632352360),
1144     UWORD(3632352360),
1145     UWORD(3632352360),
1146     UWORD(3632352360),
1147     UWORD(3874938000),
1148     UWORD(3874938000),
1149     UWORD(3874938000),
1150     UWORD(3874938000),
1151     UWORD(4129476120),
1152     UWORD(4129476120),
1153     UWORD(4129476120),
1154     UWORD(4129476120),
1155     UWORD(16974336),
1156     UWORD(16974336),
1157     UWORD(16974336),
1158     UWORD(17575740),
1159     UWORD(17575740),
1160     UWORD(17575740),
1161     UWORD(18191184),
1162     UWORD(18191184),
1163     UWORD(18191184),
1164     UWORD(18820830),
1165     UWORD(18820830),
1166     UWORD(18820830),
1167     UWORD(19464840),
1168     UWORD(19464840),
1169     UWORD(19464840),
1170     UWORD(20123376),
1171     UWORD(20123376),
1172     UWORD(20123376),
1173     UWORD(20796600),
1174     UWORD(20796600),
1175     UWORD(20796600),
1176     UWORD(21484674),
1177     UWORD(21484674),
1178     UWORD(21484674),
1179     UWORD(22187760),
1180     UWORD(22187760),
1181     UWORD(22187760),
1182     UWORD(22906020),
1183     UWORD(22906020),
1184     UWORD(22906020),
1185     UWORD(23639616),
1186     UWORD(23639616),
1187 };
1188 
1189 #endif
1190 
_arb_exp_taylor_rs(mp_ptr y,mp_limb_t * error,mp_srcptr x,mp_size_t xn,ulong N)1191 void _arb_exp_taylor_rs(mp_ptr y, mp_limb_t * error,
1192     mp_srcptr x, mp_size_t xn, ulong N)
1193 {
1194     mp_ptr s, t, xpow;
1195     mp_limb_t new_denom, old_denom, c;
1196     slong power, k, m;
1197 
1198     TMP_INIT;
1199     TMP_START;
1200 
1201     if (N >= FACTORIAL_TAB_SIZE - 1)
1202     {
1203         flint_printf("_arb_exp_taylor_rs: N too large!\n");
1204         flint_abort();
1205     }
1206 
1207     if (N <= 3)
1208     {
1209         if (N <= 1)
1210         {
1211             flint_mpn_zero(y, xn);
1212             y[xn] = N;
1213             error[0] = 0;
1214         }
1215         else if (N == 2)
1216         {
1217             flint_mpn_copyi(y, x, xn);
1218             y[xn] = 1;
1219             error[0] = 0;
1220         }
1221         else
1222         {
1223             /* 1 + x + x^2 / 2 */
1224             t = TMP_ALLOC_LIMBS(2 * xn);
1225 
1226             mpn_sqr(t, x, xn);
1227             mpn_rshift(t + xn, t + xn, xn, 1);
1228             y[xn] = mpn_add_n(y, x, t + xn, xn) + 1;
1229 
1230             error[0] = 2;
1231         }
1232     }
1233     else
1234     {
1235         /* Choose m ~= sqrt(num_terms) (m must be even, >= 2) */
1236         /* TODO: drop evenness assumption since we don't have sign issues here? */
1237         /* TODO: then just need to fix power construction below... */
1238         m = 2;
1239         while (m * m < N)
1240             m += 2;
1241 
1242         /* todo: merge allocations */
1243         xpow = TMP_ALLOC_LIMBS((m + 1) * xn);
1244         s = TMP_ALLOC_LIMBS(xn + 2);
1245         t = TMP_ALLOC_LIMBS(2 * xn + 2);     /* todo: 1 limb too much? */
1246 
1247         /* higher index ---> */
1248         /*        | ---xn--- | */
1249         /* xpow = |  <temp>  | x^m | x^(m-1) | ... | x^2 | x | */
1250 
1251 #define XPOW_WRITE(__k) (xpow + (m - (__k)) * xn)
1252 #define XPOW_READ(__k) (xpow + (m - (__k) + 1) * xn)
1253 
1254         flint_mpn_copyi(XPOW_READ(1), x, xn);
1255         mpn_sqr(XPOW_WRITE(2), XPOW_READ(1), xn);
1256 
1257         for (k = 4; k <= m; k += 2)
1258         {
1259             mpn_mul_n(XPOW_WRITE(k - 1), XPOW_READ(k / 2), XPOW_READ(k / 2 - 1), xn);
1260             mpn_sqr(XPOW_WRITE(k), XPOW_READ(k / 2), xn);
1261         }
1262 
1263         flint_mpn_zero(s, xn + 1);
1264 
1265         /* todo: skip one nonscalar multiplication (use x^m)
1266            when starting on x^0 */
1267         power = (N - 1) % m;
1268 
1269         for (k = N - 1; k >= 0; k--)
1270         {
1271             c = factorial_tab_numer[k];
1272             new_denom = factorial_tab_denom[k];
1273             old_denom = factorial_tab_denom[k+1];
1274 
1275             /* change denominators */
1276             if (new_denom != old_denom && k < N - 1)
1277             {
1278                 mpn_divrem_1(s, 0, s, xn + 1, old_denom);
1279             }
1280 
1281             if (power == 0)
1282             {
1283                 /* add c * x^0 -- only top limb is affected */
1284                 s[xn] += c;
1285 
1286                 /* Outer polynomial evaluation: multiply by x^m */
1287                 if (k != 0)
1288                 {
1289                     mpn_mul(t, s, xn + 1, XPOW_READ(m), xn);
1290                     flint_mpn_copyi(s, t + xn, xn + 1);
1291                 }
1292 
1293                 power = m - 1;
1294             }
1295             else
1296             {
1297                 s[xn] += mpn_addmul_1(s, XPOW_READ(power), xn, c);
1298 
1299                 power--;
1300             }
1301         }
1302 
1303         /* finally divide by denominator */
1304         mpn_divrem_1(y, 0, s, xn + 1, factorial_tab_denom[0]);
1305 
1306         /* error bound (ulp) */
1307         error[0] = 2;
1308     }
1309 
1310     TMP_END;
1311 }
1312 
1313