Skip to content

Commit 0c6b414

Browse files
ibiremeibireme
committed
Use python to generate fp table, remove gmp/mpfr.
1 parent 14f6f49 commit 0c6b414

File tree

4 files changed

+164
-197
lines changed

4 files changed

+164
-197
lines changed

‎CMakeLists.txt

Lines changed: 0 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -524,17 +524,7 @@ if(YYJSON_BUILD_MISC)
524524
endif()
525525

526526
# make tables
527-
find_package(GMP REQUIRED)
528-
find_package(MPFR REQUIRED)
529527
add_executable(make_tables "misc/make_tables.c")
530-
target_include_directories(make_tables PRIVATE
531-
${GMP_INCLUDE_DIR}
532-
${MPFR_INCLUDES}
533-
)
534-
target_link_libraries(make_tables PRIVATE
535-
${GMP_LIBRARIES}
536-
${MPFR_LIBRARIES}
537-
)
538528
if(XCODE)
539529
set_default_xcode_property(make_tables)
540530
endif()

‎cmake/FindGMP.cmake

Lines changed: 0 additions & 24 deletions
This file was deleted.

‎cmake/FindMPFR.cmake

Lines changed: 0 additions & 83 deletions
This file was deleted.

‎misc/make_tables.c

Lines changed: 164 additions & 80 deletions
Original file line numberDiff line numberDiff line change
@@ -13,9 +13,6 @@
1313
#include <string.h>
1414
#include <math.h>
1515

16-
#include <gmp.h>
17-
#include <mpfr.h>
18-
1916
typedef float f32;
2017
typedef double f64;
2118
typedef int8_t i8;
@@ -30,82 +27,170 @@ typedef size_t usize;
3027

3128
/*----------------------------------------------------------------------------*/
3229

33-
void make_pow10_sig_table(void) {
34-
static const int DEF_PREC = 5000;
35-
static const int BUF_LEN = 2000;
36-
char buf[BUF_LEN];
37-
38-
mpfr_t sigMax, sigMin, half;
39-
mpfr_t pow10, pow2, div, sub;
40-
mpfr_inits2(DEF_PREC, sigMax, sigMin, half, NULL);
41-
mpfr_inits2(DEF_PREC, pow10, pow2, div, sub, NULL);
42-
43-
mpfr_set_ui(sigMax, 0xFFFFFFFFFFFFFFFFULL, MPFR_RNDN);
44-
mpfr_set_ui(sigMin, 0x8000000000000000ULL, MPFR_RNDN);
45-
mpfr_set_d(half, 0.5, MPFR_RNDN);
46-
47-
int e10min = -343, e10max = 324, e10step = 1;
48-
49-
printf("#define POW10_SIG_TABLE_MIN_EXP %d\n", e10min);
50-
printf("#define POW10_SIG_TABLE_MAX_EXP %d\n", e10max);
51-
printf("#define POW10_SIG_TABLE_MIN_EXACT_EXP %d\n", 0);
52-
printf("#define POW10_SIG_TABLE_MAX_EXACT_EXP %d\n", 55);
53-
printf("static const u64 pow10_sig_table[] = {\n");
54-
55-
for (int e10 = e10min; e10 <= e10max; e10 += e10step) {
56-
mpfr_set_d(pow10, 10, MPFR_RNDN);
57-
mpfr_pow_si(pow10, pow10, e10, MPFR_RNDN); // pow10 = 10^e10
58-
59-
// 10^e10 = 2^e2
60-
// e2 = floor(log2(pow(10, e10)))
61-
// e2 = floor(log2(10) * e10)
62-
int e2 = (int)floor(log2(10) * e10) - 64 + 1;
63-
mpfr_set_d(pow2, 2, MPFR_RNDN);
64-
mpfr_pow_si(pow2, pow2, e2, MPFR_RNDN); // pow2 = 2^e2
65-
mpfr_div(div, pow10, pow2, MPFR_RNDN); // div = pow10 / pow2;
66-
if (mpfr_cmp(div, sigMin) < 0 || mpfr_cmp(div, sigMax) > 0) {
67-
printf("err!\n"); // make sure the highest bit is 1 (normalized)
68-
}
69-
70-
mpfr_set_d(pow2, 2, MPFR_RNDN);
71-
mpfr_pow_si(pow2, pow2, e2, MPFR_RNDN); // pow2 = 2^e2
72-
mpfr_div(div, pow10, pow2, MPFR_RNDN); // div = pow10 / pow2;
73-
74-
mpfr_snprintf(buf, BUF_LEN, "%.1000Rg", div);
75-
u64 val = strtoull(buf, NULL, 0);
76-
mpfr_sub_ui(sub, div, val, MPFR_RNDN); // sub = div - (uint64_t)div
77-
int cmp = mpfr_cmp(sub, half);
78-
if (cmp == 0) printf("err!\n"); // avoid round to even
79-
if (cmp > 0 && val == UINT64_MAX) printf("err!\n"); // avoid round up overflow
80-
81-
printf(" ");
82-
printf("U64(0x%.8X, 0x%.8X),", (u32)(val >> 32), (u32)val);
83-
84-
mpfr_set_d(pow2, 2, MPFR_RNDN);
85-
mpfr_pow_si(pow2, pow2, 64, MPFR_RNDN); // pow2 = 2^64
86-
mpfr_mul(sub, sub, pow2, MPFR_RNDN); // sub *= 2^64
87-
88-
mpfr_snprintf(buf, BUF_LEN, "%.1000Rg", sub);
89-
u64 val2 = strtoull(buf, NULL, 0);
90-
mpfr_sub_ui(sub, sub, val2, MPFR_RNDN); // sub -= (uint64_t)sub
91-
int cmp2 = mpfr_cmp(sub, half);
92-
if (cmp2 == 0) printf("err!\n"); // avoid round to even
93-
if ((cmp > 0) && (val2 < ((u64)1) << 63)) printf("err!\n"); // avoid round up overflow
94-
bool is_exact = mpfr_cmp_ui(sub, 0) == 0;
95-
96-
printf(" ");
97-
printf("U64(0x%.8X, 0x%.8X)", (u32)(val2 >> 32), (u32)val2);
98-
printf("%c", e10 < e10max ? ',' : ' ');
99-
printf(" /* %s 10^%d */", is_exact ? "==" : "~=", e10);
100-
printf("\n");
101-
}
102-
103-
printf("};\n");
104-
printf("\n");
105-
106-
mpfr_clears(sigMax, sigMin, half, NULL);
107-
mpfr_clears(pow10, pow2, div, sub, NULL);
108-
}
30+
//Generate fp table with C (gmp/mpfr):
31+
//#include <gmp.h>
32+
//#include <mpfr.h>
33+
//void make_pow10_sig_table(void) {
34+
// static const int DEF_PREC = 5000;
35+
// static const int BUF_LEN = 2000;
36+
// char buf[BUF_LEN];
37+
//
38+
// mpfr_t sigMax, sigMin, half;
39+
// mpfr_t pow10, pow2, div, sub;
40+
// mpfr_inits2(DEF_PREC, sigMax, sigMin, half, NULL);
41+
// mpfr_inits2(DEF_PREC, pow10, pow2, div, sub, NULL);
42+
//
43+
// mpfr_set_ui(sigMax, 0xFFFFFFFFFFFFFFFFULL, MPFR_RNDN);
44+
// mpfr_set_ui(sigMin, 0x8000000000000000ULL, MPFR_RNDN);
45+
// mpfr_set_d(half, 0.5, MPFR_RNDN);
46+
//
47+
// int e10min = -343, e10max = 324, e10step = 1;
48+
//
49+
// printf("#define POW10_SIG_TABLE_MIN_EXP %d\n", e10min);
50+
// printf("#define POW10_SIG_TABLE_MAX_EXP %d\n", e10max);
51+
// printf("#define POW10_SIG_TABLE_MIN_EXACT_EXP %d\n", 0);
52+
// printf("#define POW10_SIG_TABLE_MAX_EXACT_EXP %d\n", 55);
53+
// printf("static const u64 pow10_sig_table[] = {\n");
54+
//
55+
// for (int e10 = e10min; e10 <= e10max; e10 += e10step) {
56+
// mpfr_set_d(pow10, 10, MPFR_RNDN);
57+
// mpfr_pow_si(pow10, pow10, e10, MPFR_RNDN); // pow10 = 10^e10
58+
//
59+
// // 10^e10 = 2^e2
60+
// // e2 = floor(log2(pow(10, e10)))
61+
// // e2 = floor(log2(10) * e10)
62+
// int e2 = (int)floor(log2(10) * e10) - 64 + 1;
63+
// mpfr_set_d(pow2, 2, MPFR_RNDN);
64+
// mpfr_pow_si(pow2, pow2, e2, MPFR_RNDN); // pow2 = 2^e2
65+
// mpfr_div(div, pow10, pow2, MPFR_RNDN); // div = pow10 / pow2;
66+
// if (mpfr_cmp(div, sigMin) < 0 || mpfr_cmp(div, sigMax) > 0) {
67+
// printf("err!\n"); // make sure the highest bit is 1 (normalized)
68+
// }
69+
//
70+
// mpfr_set_d(pow2, 2, MPFR_RNDN);
71+
// mpfr_pow_si(pow2, pow2, e2, MPFR_RNDN); // pow2 = 2^e2
72+
// mpfr_div(div, pow10, pow2, MPFR_RNDN); // div = pow10 / pow2;
73+
//
74+
// mpfr_snprintf(buf, BUF_LEN, "%.1000Rg", div);
75+
// u64 val = strtoull(buf, NULL, 0);
76+
// mpfr_sub_ui(sub, div, val, MPFR_RNDN); // sub = div - (uint64_t)div
77+
// int cmp = mpfr_cmp(sub, half);
78+
// if (cmp == 0) printf("err!\n"); // avoid round to even
79+
// if (cmp > 0 && val == UINT64_MAX) printf("err!\n"); // avoid round up overflow
80+
//
81+
// printf(" ");
82+
// printf("U64(0x%.8X, 0x%.8X),", (u32)(val >> 32), (u32)val);
83+
//
84+
// mpfr_set_d(pow2, 2, MPFR_RNDN);
85+
// mpfr_pow_si(pow2, pow2, 64, MPFR_RNDN); // pow2 = 2^64
86+
// mpfr_mul(sub, sub, pow2, MPFR_RNDN); // sub *= 2^64
87+
//
88+
// mpfr_snprintf(buf, BUF_LEN, "%.1000Rg", sub);
89+
// u64 val2 = strtoull(buf, NULL, 0);
90+
// mpfr_sub_ui(sub, sub, val2, MPFR_RNDN); // sub -= (uint64_t)sub
91+
// int cmp2 = mpfr_cmp(sub, half);
92+
// if (cmp2 == 0) printf("err!\n"); // avoid round to even
93+
// if ((cmp > 0) && (val2 < ((u64)1) << 63)) printf("err!\n"); // avoid round up overflow
94+
// bool is_exact = mpfr_cmp_ui(sub, 0) == 0;
95+
//
96+
// printf(" ");
97+
// printf("U64(0x%.8X, 0x%.8X)", (u32)(val2 >> 32), (u32)val2);
98+
// printf("%c", e10 < e10max ? ',' : ' ');
99+
// printf(" /* %s 10^%d */", is_exact ? "==" : "~=", e10);
100+
// printf("\n");
101+
// }
102+
//
103+
// printf("};\n");
104+
// printf("\n");
105+
//
106+
// mpfr_clears(sigMax, sigMin, half, NULL);
107+
// mpfr_clears(pow10, pow2, div, sub, NULL);
108+
//}
109+
110+
111+
112+
//Generate fp table with Python:
113+
//import decimal
114+
//from decimal import Decimal
115+
//
116+
//POW10_SIG_TABLE_MIN_EXP = -343
117+
//POW10_SIG_TABLE_MAX_EXP = 324
118+
//POW10_SIG_TABLE_MIN_EXACT_EXP = 0
119+
//POW10_SIG_TABLE_MAX_EXACT_EXP = 55
120+
//
121+
//
122+
//def calc_pow10_u128(p: int) -> int:
123+
// """
124+
// Calculate the power of 10 and return the high 128 bits of the result.
125+
// """
126+
//
127+
// # Calculate 10^p with high precision
128+
// decimal.getcontext().prec = 5000
129+
// sig = Decimal(10) ** p
130+
//
131+
// # Normalize the sig to range [0.5,1)
132+
// while sig < 1:
133+
// sig *= 2
134+
// while sig >= 1:
135+
// sig /= 2
136+
//
137+
// # Calculate the highest 128 bits of the sig
138+
// all = sig * (2 ** 128)
139+
// top = int(all)
140+
// return top
141+
//
142+
//
143+
//def is_pow10_u128_exact(p: int) -> bool:
144+
// """
145+
// Check if calc_pow10_u128(p) is exact value.
146+
// """
147+
//
148+
// # Calculate 10^p with high precision
149+
// decimal.getcontext().prec = 5000
150+
// sig = Decimal(10) ** p
151+
//
152+
// # Normalize the sig to range [0.5,1)
153+
// while sig < 1:
154+
// sig *= 2
155+
// while sig >= 1:
156+
// sig /= 2
157+
//
158+
// # Calculate the highest 128 bits of the sig
159+
// all = sig * (2 ** 128)
160+
// top = int(all)
161+
// return top == all
162+
//
163+
//
164+
//def print_pow10_u128_table():
165+
// """
166+
// Print the power of 10 table for yy_strtod() and yy_dtoa().
167+
// """
168+
// print(f"#define POW10_SIG_TABLE_MIN_EXP {POW10_SIG_TABLE_MIN_EXP}")
169+
// print(f"#define POW10_SIG_TABLE_MAX_EXP {POW10_SIG_TABLE_MAX_EXP}")
170+
// print(f"#define POW10_SIG_TABLE_MIN_EXACT_EXP {POW10_SIG_TABLE_MIN_EXACT_EXP}")
171+
// print(f"#define POW10_SIG_TABLE_MAX_EXACT_EXP {POW10_SIG_TABLE_MAX_EXACT_EXP}")
172+
// print("static const u64 pow10_sig_table[] = {")
173+
// for p in range(POW10_SIG_TABLE_MIN_EXP, POW10_SIG_TABLE_MAX_EXP + 1):
174+
// is_exact = is_pow10_u128_exact(p)
175+
// assert is_exact == (p in range(POW10_SIG_TABLE_MIN_EXACT_EXP, POW10_SIG_TABLE_MAX_EXACT_EXP + 1))
176+
//
177+
// c = calc_pow10_u128(p)
178+
// s = f"{c:X}"
179+
// line = f" U64(0x{s[0:8]}, 0x{s[8:16]}), U64(0x{s[16:24]}, 0x{s[24:32]})"
180+
// if is_exact:
181+
// line += f", /* == 10^{p} */"
182+
// elif p == POW10_SIG_TABLE_MAX_EXP:
183+
// line += f" /* ~= 10^{p} */"
184+
// else:
185+
// line += f", /* ~= 10^{p} */"
186+
// print(line)
187+
// print("};")
188+
//
189+
//
190+
//if __name__ == "__main__":
191+
// print_pow10_u128_table()
192+
193+
109194

110195
/*----------------------------------------------------------------------------*/
111196

@@ -471,7 +556,6 @@ static void make_esc_single_char_table(void) {
471556
}
472557

473558
int main(void) {
474-
make_pow10_sig_table();
475559
make_dec_trailing_zero_table();
476560
make_char_table();
477561
make_digit_table();

0 commit comments

Comments
 (0)