hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-04 17:01:40 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

mag methods for powers, factorials, bernoulli numbers (work in progress)

Browse source
This commit is contained in:
Fredrik Johansson 2014-06-03 17:43:30 +02:00
parent 0c375f88ae
commit bf914e438c
4 changed files with 792 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

11
mag.h
View file

@ -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

87
mag/bernoulli_div_fac_ui.c Normal file
View file

@ -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 */
}
}

593
mag/fac_ui.c Normal file
View file

@ -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;
}
}

101
mag/pow_ui.c Normal file
View file

@ -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);
}
}