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