Skip to content

Commit e36f2b9

Browse files
fredrik-johanssonfredrik-johansson
authored
Merge pull request #2639 from edgarcosta/fix/fmpz-sizeinbase
fix(fmpz): return exact size from fmpz_sizeinbase
2 parents 70aeeba + 9edf7a3 commit e36f2b9

8 files changed

Lines changed: 181 additions & 43 deletions

File tree

‎src/fmpz/sizeinbase.c‎

Lines changed: 41 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -9,16 +9,55 @@
99
(at your option) any later version. See <https://www.gnu.org/licenses/>.
1010
*/
1111

12+
#include <math.h>
1213
#include <gmp.h>
1314
#include "long_extras.h"
15+
#include "ulong_extras.h"
1416
#include "fmpz.h"
1517

1618
size_t fmpz_sizeinbase(const fmpz_t f, int b)
1719
{
1820
fmpz d = *f;
21+
mpz_srcptr z;
22+
long e;
23+
double s, lg;
24+
slong lo, hi;
1925

2026
if (!COEFF_IS_MPZ(d))
2127
return z_sizeinbase(d, b);
22-
else
23-
return mpz_sizeinbase(COEFF_TO_PTR(d), b);
28+
29+
z = COEFF_TO_PTR(d);
30+
31+
/* Power-of-2 bases: mpz_sizeinbase is exact. */
32+
if ((b & (b - 1)) == 0)
33+
return mpz_sizeinbase(z, b);
34+
35+
if (b == 10 && mpz_size(z) == 1)
36+
return (size_t) n_nonzero_sizeinbase10(mpz_getlimbn(z, 0));
37+
38+
/* Otherwise approximate log_b(|z|) in double precision. If its
39+
floor is robust to a tiny relative perturbation, return floor + 1
40+
in O(1). When |z| is near a power of b the check is inconclusive
41+
and we compare |f| against b^hi directly. */
42+
s = mpz_get_d_2exp(&e, z);
43+
lg = (log(fabs(s)) + (double) e * 0.69314718055994530942)
44+
* ((b == 10) ? 0.43429448190325176 : 1.0 / log((double) b));
45+
46+
lo = (slong) (lg * (1.0 - 1e-12));
47+
hi = (slong) (lg * (1.0 + 1e-12));
48+
49+
if (lo == hi)
50+
return (size_t) lo + 1;
51+
52+
{
53+
fmpz_t pow;
54+
int cmp;
55+
56+
fmpz_init(pow);
57+
fmpz_ui_pow_ui(pow, (ulong) b, (ulong) hi);
58+
cmp = fmpz_cmpabs(f, pow);
59+
fmpz_clear(pow);
60+
61+
return (size_t) ((cmp >= 0) ? hi + 1 : hi);
62+
}
2463
}

‎src/fmpz/test/t-sizeinbase.c‎

Lines changed: 34 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -23,8 +23,9 @@ TEST_FUNCTION_START(fmpz_sizeinbase, state)
2323
{
2424
fmpz_t a;
2525
mpz_t b;
26+
char * str;
2627
int base;
27-
size_t r1, r2;
28+
size_t r1, exact;
2829

2930
fmpz_init(a);
3031
mpz_init(b);
@@ -34,15 +35,21 @@ TEST_FUNCTION_START(fmpz_sizeinbase, state)
3435
fmpz_get_mpz(b, a);
3536

3637
r1 = fmpz_sizeinbase(a, base);
37-
r2 = mpz_sizeinbase(b, base);
38-
result = (r1 == r2 || r1 + 1 == r2);
38+
39+
/* exact size = length of |b| written in base b */
40+
str = mpz_get_str(NULL, base, b);
41+
exact = strlen(str) - (str[0] == '-');
42+
flint_free(str);
43+
44+
result = (r1 == exact);
3945

4046
if (!result)
4147
{
4248
flint_printf("FAIL:\n");
4349
gmp_printf("b = %Zd\n", b);
4450
flint_printf("base = %d\n", base);
45-
flint_printf("r1 = %wu\n, r2 = %wu\n", (ulong) r1, (ulong) r2);
51+
flint_printf("r1 = %wu, exact = %wu\n",
52+
(ulong) r1, (ulong) exact);
4653
fflush(stdout);
4754
flint_abort();
4855
}
@@ -51,5 +58,28 @@ TEST_FUNCTION_START(fmpz_sizeinbase, state)
5158
mpz_clear(b);
5259
}
5360

61+
/* 10^19 - 1 has 19 decimal digits, not 20. */
62+
{
63+
fmpz_t a;
64+
size_t r;
65+
66+
fmpz_init(a);
67+
fmpz_set_ui(a, 10);
68+
fmpz_pow_ui(a, a, 19);
69+
fmpz_sub_ui(a, a, 1);
70+
71+
r = fmpz_sizeinbase(a, 10);
72+
if (r != 19)
73+
{
74+
flint_printf("FAIL (10^19 - 1):\n");
75+
flint_printf("fmpz_sizeinbase(10^19 - 1, 10) = %wu, expected 19\n",
76+
(ulong) r);
77+
fflush(stdout);
78+
flint_abort();
79+
}
80+
81+
fmpz_clear(a);
82+
}
83+
5484
TEST_FUNCTION_END(state);
5585
}

‎src/long_extras/sizeinbase.c‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@
1111

1212
#include <limits.h>
1313
#include "long_extras.h"
14+
#include "ulong_extras.h"
1415

1516
size_t z_sizeinbase(slong n, int b)
1617
{
@@ -21,6 +22,9 @@ size_t z_sizeinbase(slong n, int b)
2122
return 1;
2223
}
2324

25+
if (b == 10)
26+
return (size_t) n_nonzero_sizeinbase10((n < 0) ? -(ulong) n : (ulong) n);
27+
2428
if (n <= 0)
2529
{
2630
if (n > WORD_MIN)

‎src/radix/get_str.c‎

Lines changed: 1 addition & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -14,34 +14,6 @@
1414
#include "mpn_extras.h"
1515
#include "radix.h"
1616

17-
static const uint8_t n_sizeinbase10_tab1[64] = {
18-
1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 7, 8,
19-
8, 8, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 13,
20-
14, 14, 14, 15, 15, 15, 16, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19,
21-
19, 19, 20,
22-
};
23-
24-
static const ulong n_sizeinbase10_tab2[FLINT_BITS] = {
25-
1, 1, 1, 10, 10, 10, 100, 100, 100, 1000, 1000, 1000, 1000, 10000, 10000,
26-
10000, 100000, 100000, 100000, 1000000, 1000000, 1000000, 1000000, 10000000,
27-
10000000, 10000000, 100000000, 100000000, 100000000, 1000000000,
28-
1000000000, 1000000000,
29-
#if FLINT_BITS == 64
30-
1000000000, UWORD(10000000000), UWORD(10000000000), UWORD(10000000000),
31-
UWORD(100000000000), UWORD(100000000000), UWORD(100000000000),
32-
UWORD(1000000000000), UWORD(1000000000000), UWORD(1000000000000),
33-
UWORD(1000000000000), UWORD(10000000000000), UWORD(10000000000000),
34-
UWORD(10000000000000), UWORD(100000000000000), UWORD(100000000000000),
35-
UWORD(100000000000000), UWORD(1000000000000000), UWORD(1000000000000000),
36-
UWORD(1000000000000000), UWORD(1000000000000000), UWORD(10000000000000000),
37-
UWORD(10000000000000000), UWORD(10000000000000000),
38-
UWORD(100000000000000000), UWORD(100000000000000000),
39-
UWORD(100000000000000000), UWORD(1000000000000000000),
40-
UWORD(1000000000000000000), UWORD(1000000000000000000),
41-
UWORD(1000000000000000000), UWORD(10000000000000000000),
42-
#endif
43-
};
44-
4517
static const char dec_to_str_tab[200] =
4618
"000102030405060708091011121314151617181920212223242526272829"
4719
"303132333435363738394041424344454647484950515253545556575859"
@@ -103,14 +75,6 @@ n_get_str_nd(char * s, ulong x, int d)
10375
}
10476
}
10577

106-
static slong
107-
n_nonzero_sizeinbase10(ulong n)
108-
{
109-
FLINT_ASSERT(n != 0);
110-
int b = FLINT_BIT_COUNT(n) - 1;
111-
return n_sizeinbase10_tab1[b] - (n < n_sizeinbase10_tab2[b]);
112-
}
113-
11478
static char * _radix_decimal_get_str(char * res, nn_srcptr t, slong decimal_limbs, int negative, slong digits_per_limb)
11579
{
11680
slong i;
@@ -159,7 +123,7 @@ char * radix_get_str_decimal(char * res, nn_srcptr x, slong n, int negative, con
159123
{
160124
return _radix_decimal_get_str(res, x, n, negative, radix->exp);
161125
}
162-
else if (n == 1 && LIMB_RADIX(radix) < n_sizeinbase10_tab2[FLINT_BITS - 1]) /* todo: could work even for 10/20-digit input */
126+
else if (n == 1 && n_nonzero_sizeinbase10(LIMB_RADIX(radix)) <= ((FLINT_BITS == 64) ? 19 : 9)) /* todo: could work even for 10/20-digit input */
163127
{
164128
return _radix_decimal_get_str(res, x, 1, negative, ((FLINT_BITS == 64) ? 19 : 9));
165129
}

‎src/ulong_extras.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -180,6 +180,7 @@ int n_is_perfect_power235(ulong n);
180180
int n_is_perfect_power(ulong * root, ulong n);
181181

182182
int n_sizeinbase(ulong n, int base);
183+
slong n_nonzero_sizeinbase10(ulong n);
183184

184185
ulong n_flog(ulong n, ulong b);
185186
ulong n_clog(ulong n, ulong b);

‎src/ulong_extras/nonzero_sizeinbase10.c‎

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,50 @@
1+
/*
2+
Copyright (C) 2026 Fredrik Johansson
3+
4+
This file is part of FLINT.
5+
6+
FLINT 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 3 of the License, or
9+
(at your option) any later version. See <https://www.gnu.org/licenses/>.
10+
*/
11+
12+
#include <stdint.h>
13+
#include "ulong_extras.h"
14+
15+
static const uint8_t n_sizeinbase10_tab1[64] = {
16+
1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 7, 8,
17+
8, 8, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 13,
18+
14, 14, 14, 15, 15, 15, 16, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19,
19+
19, 19, 20,
20+
};
21+
22+
static const ulong n_sizeinbase10_tab2[FLINT_BITS] = {
23+
1, 1, 1, 10, 10, 10, 100, 100, 100, 1000, 1000, 1000, 1000, 10000, 10000,
24+
10000, 100000, 100000, 100000, 1000000, 1000000, 1000000, 1000000, 10000000,
25+
10000000, 10000000, 100000000, 100000000, 100000000, 1000000000,
26+
1000000000, 1000000000,
27+
#if FLINT_BITS == 64
28+
1000000000, UWORD(10000000000), UWORD(10000000000), UWORD(10000000000),
29+
UWORD(100000000000), UWORD(100000000000), UWORD(100000000000),
30+
UWORD(1000000000000), UWORD(1000000000000), UWORD(1000000000000),
31+
UWORD(1000000000000), UWORD(10000000000000), UWORD(10000000000000),
32+
UWORD(10000000000000), UWORD(100000000000000), UWORD(100000000000000),
33+
UWORD(100000000000000), UWORD(1000000000000000), UWORD(1000000000000000),
34+
UWORD(1000000000000000), UWORD(1000000000000000), UWORD(10000000000000000),
35+
UWORD(10000000000000000), UWORD(10000000000000000),
36+
UWORD(100000000000000000), UWORD(100000000000000000),
37+
UWORD(100000000000000000), UWORD(1000000000000000000),
38+
UWORD(1000000000000000000), UWORD(1000000000000000000),
39+
UWORD(1000000000000000000), UWORD(10000000000000000000),
40+
#endif
41+
};
42+
43+
slong
44+
n_nonzero_sizeinbase10(ulong n)
45+
{
46+
int b;
47+
FLINT_ASSERT(n != 0);
48+
b = FLINT_BIT_COUNT(n) - 1;
49+
return n_sizeinbase10_tab1[b] - (n < n_sizeinbase10_tab2[b]);
50+
}

‎src/ulong_extras/test/main.c‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -78,6 +78,7 @@
7878
#include "t-mulmod_shoup.c"
7979
#include "t-mulmod_and_precomp_shoup.c"
8080
#include "t-nextprime.c"
81+
#include "t-nonzero_sizeinbase10.c"
8182
#include "t-nth_prime_bounds.c"
8283
#include "t-urandint.c"
8384
#include "t-pow.c"
@@ -179,6 +180,7 @@ test_struct tests[] =
179180
TEST_FUNCTION(n_mulmod_shoup),
180181
TEST_FUNCTION(n_mulmod_and_precomp_shoup),
181182
TEST_FUNCTION(n_nextprime),
183+
TEST_FUNCTION(n_nonzero_sizeinbase10),
182184
TEST_FUNCTION(n_nth_prime_bounds),
183185
TEST_FUNCTION(n_urandint),
184186
TEST_FUNCTION(n_pow),

‎src/ulong_extras/test/t-nonzero_sizeinbase10.c‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
/*
2+
Copyright (C) 2026 Fredrik Johansson
3+
4+
This file is part of FLINT.
5+
6+
FLINT 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 3 of the License, or
9+
(at your option) any later version. See <https://www.gnu.org/licenses/>.
10+
*/
11+
12+
#include <string.h>
13+
#include "test_helpers.h"
14+
#include "gmpcompat.h"
15+
#include "ulong_extras.h"
16+
17+
TEST_FUNCTION_START(n_nonzero_sizeinbase10, state)
18+
{
19+
ulong n;
20+
slong rep, size1, size2;
21+
mpz_t t;
22+
char * str;
23+
24+
mpz_init(t);
25+
str = flint_malloc((FLINT_BITS + 1) * sizeof(char));
26+
27+
for (rep = 0; rep < 10000 * flint_test_multiplier(); rep++)
28+
{
29+
n = n_randtest_not_zero(state);
30+
31+
size1 = n_nonzero_sizeinbase10(n);
32+
33+
flint_mpz_set_ui(t, n);
34+
mpz_get_str(str, 10, t);
35+
size2 = strlen(str);
36+
37+
if (size1 != size2)
38+
TEST_FUNCTION_FAIL(
39+
"n = %wu\n"
40+
"n_nonzero_sizeinbase10: %wd, strlen: %wd\n",
41+
n, size1, size2);
42+
}
43+
44+
flint_free(str);
45+
mpz_clear(t);
46+
47+
TEST_FUNCTION_END(state);
48+
}

0 commit comments

Comments
 (0)