From bf914e438c1c25fcc2d099bf668375f149b6e40e Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Tue, 3 Jun 2014 17:43:30 +0200 Subject: [PATCH] mag methods for powers, factorials, bernoulli numbers (work in progress) --- mag.h | 11 + mag/bernoulli_div_fac_ui.c | 87 ++++++ mag/fac_ui.c | 593 +++++++++++++++++++++++++++++++++++++ mag/pow_ui.c | 101 +++++++ 4 files changed, 792 insertions(+) create mode 100644 mag/bernoulli_div_fac_ui.c create mode 100644 mag/fac_ui.c create mode 100644 mag/pow_ui.c diff --git a/mag.h b/mag.h index b194fd20..fa2c7e31 100644 --- a/mag.h +++ b/mag.h @@ -639,6 +639,17 @@ _mag_vec_clear(mag_ptr v, long n) double mag_d_log_upper_bound(double x); double mag_d_log_lower_bound(double x); +/* TODO: document/test */ +void mag_pow_ui(mag_t z, const mag_t x, ulong e); +void mag_pow_ui_lower(mag_t z, const mag_t x, ulong e); + +/* TODO: document/test */ +void mag_fac_ui(mag_t z, ulong n); +void mag_rfac_ui(mag_t z, ulong n); + +/* TODO: document/test */ +void mag_bernoulli_div_fac_ui(mag_t z, ulong n); + #ifdef __cplusplus } #endif diff --git a/mag/bernoulli_div_fac_ui.c b/mag/bernoulli_div_fac_ui.c new file mode 100644 index 00000000..37212d6f --- /dev/null +++ b/mag/bernoulli_div_fac_ui.c @@ -0,0 +1,87 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2014 Fredrik Johansson + +******************************************************************************/ + +#include "mag.h" + +const int mag_bernoulli_div_fac_ui_tab[32] = +{ + 1, 536870912, + -3, 715827883, + -9, 763549742, + -14, 581752185, + -20, 930803495, + -25, 752164440, + -30, 609225645, + -36, 987456906, + -41, 800365631, + -46, 648744507, + -52, 1051701970, + -57, 852476898, + -62, 690991623, + -67, 560096688, + -73, 907994540, + -78, 735992650, +}; + +/* computes upper bound for B_n / n! */ +void +mag_bernoulli_div_fac_ui(mag_t z, ulong n) +{ + if (n % 2 == 1) + { + if (n == 1) + { + mag_one(z); + mag_mul_2exp_si(z, z, -1); + } + else + { + mag_zero(z); + } + } + else if (n < 32) + { + _fmpz_demote(MAG_EXPREF(z)); + MAG_EXP(z) = mag_bernoulli_div_fac_ui_tab[n]; + MAG_MAN(z) = mag_bernoulli_div_fac_ui_tab[n + 1]; + } + else + { + /* upper bound for 1/(2pi) */ + _fmpz_demote(MAG_EXPREF(z)); + MAG_EXP(z) = -2; + MAG_MAN(z) = 683565276; + + /* 2 (1/(2pi))^n */ + mag_pow_ui(z, z, n); + mag_mul_2exp_si(z, z, 1); + + /* we should multiply by an upper bound for zeta(n), + but the error in the upper bound for 1/(2pi) + used above is already larger than this factor + when n >= 32, so there is no need */ + } +} + diff --git a/mag/fac_ui.c b/mag/fac_ui.c new file mode 100644 index 00000000..ed7133dc --- /dev/null +++ b/mag/fac_ui.c @@ -0,0 +1,593 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2014 Fredrik Johansson + +******************************************************************************/ + +#include "mag.h" + +#define MAG_FAC_TABLE_NUM 256 + +const int mag_fac_tab[MAG_FAC_TABLE_NUM * 2] = +{ + 1, 536870912, + 1, 536870912, + 2, 536870912, + 3, 805306368, + 5, 805306368, + 7, 1006632960, + 10, 754974720, + 13, 660602880, + 16, 660602880, + 19, 743178240, + 22, 928972800, + 26, 638668800, + 29, 958003200, + 33, 778377600, + 37, 681080400, + 41, 638512875, + 45, 638512875, + 49, 678419930, + 53, 763222421, + 57, 906326625, + 62, 566454141, + 66, 743471060, + 70, 1022272707, + 75, 734758508, + 80, 551068881, + 84, 861045127, + 89, 699599166, + 94, 590286796, + 98, 1033001893, + 103, 936157966, + 108, 877648093, + 113, 850221590, + 118, 850221590, + 123, 876791015, + 128, 931590453, + 133, 1018927058, + 139, 573146470, + 144, 662700606, + 149, 786956970, + 154, 959103807, + 160, 599439879, + 165, 768032345, + 170, 1008042453, + 176, 677278523, + 181, 931257969, + 187, 654790760, + 192, 941261717, + 198, 691239074, + 203, 1036858610, + 209, 793844874, + 215, 620191308, + 220, 988429896, + 226, 803099291, + 232, 665066600, + 238, 561149944, + 243, 964476466, + 249, 843916908, + 255, 751613496, + 261, 681149731, + 267, 627934908, + 273, 588688976, + 279, 561094181, + 285, 543559988, + 290, 1070133725, + 296, 1070133725, + 303, 543427282, + 309, 560409385, + 315, 586678575, + 321, 623345986, + 327, 672044891, + 333, 735049099, + 339, 815445095, + 345, 917375731, + 351, 1046381693, + 358, 604939417, + 364, 708913379, + 370, 841834637, + 376, 1012832298, + 383, 617194682, + 389, 761849685, + 395, 952312106, + 402, 602635005, + 408, 772126100, + 414, 1001351036, + 421, 657136617, + 427, 872759570, + 434, 586385336, + 440, 797117566, + 447, 548018327, + 453, 762087986, + 459, 1071686230, + 466, 761901929, + 473, 547617012, + 479, 795755970, + 486, 584383290, + 492, 867443946, + 499, 650582960, + 505, 986039798, + 512, 754936721, + 519, 583896370, + 525, 912338078, + 532, 719891765, + 539, 573663750, + 545, 923240098, + 552, 750132579, + 559, 615343132, + 565, 1019162061, + 572, 851955786, + 579, 718837694, + 586, 612135224, + 592, 1052107416, + 599, 912374400, + 606, 798327600, + 613, 704773585, + 620, 627688974, + 627, 563939312, + 633, 1022140003, + 640, 934299847, + 647, 861307671, + 654, 800746976, + 661, 750700290, + 668, 709646368, + 675, 676381694, + 682, 649960534, + 689, 629649268, + 696, 614891863, + 703, 605284178, + 710, 600555395, + 717, 600555395, + 724, 605247234, + 731, 614704222, + 738, 629111352, + 745, 648771082, + 752, 674113702, + 759, 705712782, + 766, 744306450, + 773, 790825603, + 780, 846430528, + 787, 912557913, + 794, 990980859, + 802, 541942657, + 809, 596983708, + 816, 662278801, + 823, 739889598, + 830, 832375798, + 837, 942925709, + 845, 537762319, + 852, 617586413, + 859, 714084290, + 866, 831238743, + 873, 974107902, + 881, 574571458, + 888, 682303606, + 895, 815566029, + 902, 981227879, + 910, 594102818, + 917, 724062809, + 924, 888108289, + 932, 548129335, + 939, 680879408, + 946, 851099260, + 953, 1070523287, + 961, 677440518, + 968, 862678159, + 976, 552653196, + 983, 712404510, + 990, 923899599, + 998, 602700129, + 1005, 791043920, + 1012, 1044425175, + 1020, 693563593, + 1027, 926557612, + 1035, 622530896, + 1042, 841389414, + 1050, 571881867, + 1057, 781869740, + 1065, 537535447, + 1072, 743310735, + 1079, 1033666490, + 1087, 722758991, + 1094, 1016379831, + 1102, 718612303, + 1109, 1021776868, + 1117, 730410808, + 1124, 1049965536, + 1132, 758764157, + 1140, 551289583, + 1147, 805399625, + 1155, 591465350, + 1162, 873335555, + 1170, 648178733, + 1177, 967204202, + 1185, 725403152, + 1193, 546885970, + 1200, 828874048, + 1208, 631368904, + 1215, 966783634, + 1223, 743970218, + 1231, 575414466, + 1238, 894589677, + 1246, 698898186, + 1254, 548744279, + 1261, 865987064, + 1269, 686700680, + 1277, 547214604, + 1284, 876398390, + 1292, 705226829, + 1300, 570242007, + 1307, 926643260, + 1315, 756517349, + 1323, 620580638, + 1330, 1022988396, + 1338, 847162265, + 1346, 704865479, + 1354, 589223486, + 1361, 989711324, + 1369, 835068930, + 1377, 707851397, + 1385, 602779706, + 1392, 1031318402, + 1400, 886289252, + 1408, 765116894, + 1416, 663499806, + 1424, 577970534, + 1431, 1011448435, + 1439, 888968351, + 1447, 784792372, + 1455, 695890111, + 1463, 619777131, + 1471, 554410011, + 1478, 996205489, + 1486, 898919797, + 1494, 814646066, + 1502, 741455208, + 1510, 677736401, + 1518, 622140837, + 1526, 573536084, + 1533, 1061937906, + 1541, 987270397, + 1549, 921709472, + 1557, 864102630, + 1565, 813471617, + 1573, 768984888, + 1581, 729934874, + 1589, 695719177, + 1597, 665824993, + 1605, 639816204, + 1613, 617322666, + 1621, 598031333, + 1629, 581678914, + 1637, 568045814, + 1645, 556951169, + 1653, 548248807, + 1661, 541824017, + 1669, 537591016, + 1676, 1070982102, +}; + +const int mag_rfac_tab[MAG_FAC_TABLE_NUM * 2] = +{ + 1, 536870912, + 1, 536870912, + 0, 536870912, + -2, 715827883, + -4, 715827883, + -6, 572662307, + -9, 763549742, + -12, 872628277, + -15, 872628277, + -18, 775669579, + -21, 620535663, + -25, 902597328, + -28, 601731552, + -32, 740592680, + -36, 846391634, + -40, 902817743, + -44, 902817743, + -48, 849710817, + -52, 755298504, + -56, 636040846, + -61, 1017665353, + -65, 775364078, + -69, 563901148, + -74, 784558119, + -79, 1046077491, + -83, 669489595, + -88, 823987193, + -93, 976577414, + -97, 558044237, + -102, 615772951, + -107, 656824481, + -112, 678012368, + -117, 678012368, + -122, 657466538, + -127, 618792036, + -132, 565752719, + -138, 1005782611, + -143, 869866042, + -148, 732518772, + -153, 601041044, + -159, 961665670, + -164, 750568328, + -169, 571861583, + -175, 851142821, + -180, 619012961, + -186, 880373989, + -191, 612434079, + -197, 833952789, + -202, 555968526, + -208, 726162972, + -214, 929488605, + -219, 583208536, + -225, 717795122, + -231, 866771468, + -237, 1027284702, + -242, 597692918, + -248, 683077620, + -254, 766964346, + -260, 846305485, + -266, 918026288, + -272, 979228041, + -278, 1027386797, + -284, 1060528307, + -289, 538681045, + -295, 538681045, + -302, 1060787288, + -308, 1028642219, + -314, 982583612, + -320, 924784576, + -326, 857771201, + -332, 784247955, + -338, 706927734, + -344, 628380208, + -350, 550908676, + -357, 952923115, + -363, 813161058, + -369, 684767207, + -375, 569157159, + -382, 934001491, + -388, 756659436, + -394, 605327549, + -401, 956566990, + -407, 746588871, + -413, 575682985, + -420, 877231215, + -426, 660503503, + -433, 983074981, + -439, 723181595, + -446, 1051900502, + -452, 756422833, + -458, 537900681, + -465, 756607552, + -472, 1052671376, + -478, 724419012, + -485, 986442909, + -491, 664551013, + -498, 886068017, + -504, 584622197, + -511, 763588175, + -518, 987265519, + -524, 631849933, + -531, 800760310, + -538, 1004875684, + -544, 624388774, + -551, 768478492, + -558, 936811875, + -564, 565622264, + -571, 676632242, + -578, 801934508, + -585, 941721257, + -591, 547910550, + -598, 631824778, + -605, 722085461, + -612, 817937513, + -619, 918385980, + -626, 1022203525, + -632, 563974359, + -639, 616997589, + -646, 669285520, + -653, 719903753, + -660, 767897336, + -667, 812321149, + -674, 852271370, + -681, 886916547, + -688, 915526759, + -695, 937499401, + -702, 952380344, + -709, 959879401, + -716, 959879401, + -723, 952438476, + -730, 937785576, + -737, 916309571, + -744, 888542614, + -751, 855138756, + -758, 816848961, + -765, 774493830, + -772, 728935369, + -779, 681049104, + -786, 631697720, + -793, 581707253, + -801, 1063693262, + -808, 965622252, + -815, 870420058, + -822, 779117255, + -829, 692548671, + -836, 611353310, + -844, 1071961968, + -851, 933409060, + -858, 807272701, + -865, 693496012, + -872, 591783263, + -880, 1003288181, + -887, 844874258, + -894, 706822909, + -901, 587489171, + -909, 970304695, + -916, 796147442, + -923, 649088360, + -931, 1051687470, + -938, 846641485, + -945, 677313188, + -952, 538485019, + -960, 850939290, + -967, 668222264, + -975, 1043078656, + -982, 809176170, + -989, 623943071, + -997, 956463629, + -1004, 728734194, + -1011, 551940691, + -1019, 831157747, + -1026, 622153167, + -1034, 925995411, + -1041, 685129553, + -1049, 1008006698, + -1056, 737284900, + -1064, 1072414399, + -1071, 775531317, + -1078, 557685442, + -1086, 797583648, + -1093, 567170594, + -1101, 802186034, + -1108, 564174793, + -1116, 789228126, + -1123, 549028262, + -1131, 759736405, + -1139, 1045658708, + -1146, 715744998, + -1154, 974631486, + -1161, 660067885, + -1169, 889354624, + -1176, 596007287, + -1184, 794676383, + -1192, 1054078518, + -1199, 695474487, + -1207, 913033172, + -1214, 596266561, + -1222, 774843856, + -1230, 1001818319, + -1237, 644385653, + -1245, 824813635, + -1253, 1050508908, + -1260, 665669011, + -1268, 839464369, + -1276, 1053445483, + -1283, 657761082, + -1291, 817411830, + -1299, 1010905451, + -1306, 622095663, + -1314, 761992773, + -1322, 928905476, + -1329, 563506640, + -1337, 680460848, + -1345, 817830878, + -1353, 978339742, + -1360, 582453428, + -1368, 690315174, + -1376, 814381034, + -1384, 956337361, + -1391, 558955170, + -1399, 650420562, + -1407, 753428343, + -1415, 868818269, + -1423, 997387789, + -1430, 569935879, + -1438, 648460378, + -1446, 734539189, + -1454, 828378998, + -1462, 930109752, + -1470, 1039773347, + -1477, 578656472, + -1485, 641281631, + -1493, 707621110, + -1501, 777472121, + -1509, 850567790, + -1517, 926575975, + -1525, 1005099363, + -1532, 542838475, + -1540, 583893485, + -1548, 625425658, + -1556, 667120702, + -1564, 708642737, + -1572, 749638598, + -1580, 789742721, + -1588, 828582526, + -1596, 865784191, + -1604, 900978670, + -1612, 933807852, + -1620, 963930686, + -1628, 991029139, + -1636, 1014813839, + -1644, 1035029254, + -1652, 1051458290, + -1660, 1063926174, + -1668, 1072303546, + -1675, 538254329, +}; + +void +mag_fac_ui(mag_t z, ulong n) +{ + if (n < MAG_FAC_TABLE_NUM) + { + _fmpz_demote(MAG_EXPREF(z)); + MAG_EXP(z) = mag_fac_tab[n * 2]; + MAG_MAN(z) = mag_fac_tab[n * 2 + 1]; + } + else + { + double x = n; + + x = ceil((((x+0.5)*log(x) - x) * 1.4426950408889634074 + 2) * 1.0000001); + + /* x + 1 could round down for huge x, but this doesn't matter + as long as the result was perturbed up above */ + fmpz_set_d(MAG_EXPREF(z), x + 1); + MAG_MAN(z) = MAG_ONE_HALF; + } +} + +void +mag_rfac_ui(mag_t z, ulong n) +{ + if (n < MAG_FAC_TABLE_NUM) + { + _fmpz_demote(MAG_EXPREF(z)); + MAG_EXP(z) = mag_rfac_tab[n * 2]; + MAG_MAN(z) = mag_rfac_tab[n * 2 + 1]; + } + else + { + double x = n; + + x = ceil((((x+0.5)*log(x) - x) * 1.4426950408889634074) * -0.9999999); + + /* x + 1 could round down for huge x, but this doesn't matter + as long as the value was perturbed up above */ + fmpz_set_d(MAG_EXPREF(z), x + 1); + MAG_MAN(z) = MAG_ONE_HALF; + } +} + diff --git a/mag/pow_ui.c b/mag/pow_ui.c new file mode 100644 index 00000000..96095986 --- /dev/null +++ b/mag/pow_ui.c @@ -0,0 +1,101 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2014 Fredrik Johansson + +******************************************************************************/ + +#include "mag.h" + +void +mag_pow_ui(mag_t z, const mag_t x, ulong e) +{ + if (mag_is_inf(z)) + { + mag_inf(z); + } + else if (e <= 2) + { + if (e == 0) + mag_one(z); + else if (e == 1) + mag_set(z, x); + else + mag_mul(z, x, x); + } + else + { + mag_t y; + int i, bits; + + mag_init_set(y, x); + + bits = FLINT_BIT_COUNT(e); + + for (i = bits - 2; i >= 0; i--) + { + mag_mul(y, y, y); + if (e & (1UL << i)) + mag_mul(y, y, x); + } + + mag_swap(z, y); + mag_clear(y); + } +} + +void +mag_pow_ui_lower(mag_t z, const mag_t x, ulong e) +{ + if (e <= 2) + { + if (e == 0) + mag_one(z); + else if (e == 1) + mag_set(z, x); + else + mag_mul_lower(z, x, x); + } + else if (mag_is_inf(z)) + { + mag_inf(z); + } + else + { + mag_t y; + int i, bits; + + mag_init_set(y, x); + + bits = FLINT_BIT_COUNT(e); + + for (i = bits - 2; i >= 0; i--) + { + mag_mul_lower(y, y, y); + if (e & (1UL << i)) + mag_mul_lower(y, y, x); + } + + mag_swap(z, y); + mag_clear(y); + } +} +