Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

Arbitrary precision integer binary logarithm (in Prolog)

163 views
Skip to first unread message

Julio P. Di Egidio

unread,
Feb 7, 2017, 1:46:08 AM2/7/17
to
[Cross-posted to comp.theory, sci.math, comp.lang.prolog. Please set
follow-ups as needed.]

Hello all,

I have put together an algorithm (plus reference implementation in Prolog)
for the arbitrary precision integer binary logarithm with remainder: I
needed this myself but could not find anything on the web, so I thought I'd
share.

Three Prolog predicates are in fact offered, which are progressive
refinements of the naive approach that shifts the input to the right by one
bit at a time:

ilog2_naive(+X, -Y, -R) :
Naive implementation: the input is shifted to the right by 1 bit at a time
repeatedly until the result of shifting is zero. The total number of
repetitions corresponds to the ceiling of the real logarithm. The
complexity is O(log2(X^3).

ilog2_n(+K, +X, -Y, -R) :
Faster implementation that only uses elementary integer arithmetic: the
input is shifted to the right by 2^M bits at a time, with M set equal to K
initially, and progressively decremented down to 1 while the input gets
smaller. The complexity is O(log2(X^(1/2)).

ilog2_f(+K, +X, -Y, -R) :
Even faster implementation that also uses FPU-based arithmetic: the input is
shifted to the right by 2^M bits at a time, with M set equal to K initially,
and progressively decremented down to 32 while the input gets smaller.
Concludes by using FPU-based arithmetic. The complexity is O(log2(X^(1/8)).

I'd appreciate feedback overall, and in particular on the following two
issues:

- The proposed algorithm is simple but should be pretty decent, or at least
I cannot guess any obvious better approach, given that the size of the input
is unbounded. Although, if anybody is aware of some better approach, please
let me know.

- It would also be great if somebody could check and in fact fix my
estimates for the algorithmic complexities. All I have done is some visual
inspection, and I am far from being an expert in that matter anyway...

The code was developed and can be run in SWI-Prolog, but should be
straightforward to port to other languages because, besides the overall
simplicity, the approach is in fact purely functional. It just assumes
availability of some big integers in arbitrary precision.

The source code, which includes full documentation, is available here:
<https://gist.github.com/jp-diegidio/eb2d6c1ae432c0e3c75d9f1e5bba602c>

An HTML extract of the documentation, for easier reading, is available here:
<http://julio.diegidio.name/Projects/Nan.Numerics.Gist/Prolog/ilog2.html>

Thank you,

Julio

burs...@gmail.com

unread,
Feb 7, 2017, 6:13:53 AM2/7/17
to
Here is an O(3) algorithm in number of big num ops,
or if you want the time is O(n) if you measure
it the big num ops itself will take, where
n is the bit size of the big num ops arguments.

ilog2(X, Y, R) :-
Y is msb(X), R is X-(1 << Y).

Here is some example runs:

?- ilog2(1000, Y, R).
Y = 9,
R = 488.

?- ilog2(1024, Y, R).
Y = 10,
R = 0.

The internal implementation of msb is verify fast,
for big nums that are binary represented internally.

msb(+IntExpr)
http://www.swi-prolog.org/pldoc/man?function=msb/1

That big nums are internally binary represented is
nowadays a quasi standard, since conversion to
decimals might incure some additional cost, which
is outweighted by the faster ops and the compacter

storage use. So the O(n) here is only a result of
the copying for R.

j4n bur53

unread,
Feb 7, 2017, 6:29:21 AM2/7/17
to
Here is how msb/1 is implemented in Java, its called
bitLength() there:

private static int bitLength(int[] val, int len) {
if (len == 0)
return 0;
return ((len - 1) << 5) + bitLengthForInt(val[0]);
}

static int bitLengthForInt(int n) {
return 32 - Integer.numberOfLeadingZeros(n);
}

public static int numberOfLeadingZeros(int i) {
// HD, Figure 5-6
if (i == 0)
return 32;
int n = 1;
if (i >>> 16 == 0) { n += 16; i <<= 16; }
if (i >>> 24 == 0) { n += 8; i <<= 8; }
if (i >>> 28 == 0) { n += 4; i <<= 4; }
if (i >>> 30 == 0) { n += 2; i <<= 2; }
n -= i >>> 31;
return n;
}

This is from the source of OpenJDK 1.8:
http://grepcode.com/file/repository.grepcode.com/java/root/jdk/openjdk/8u40-b25/java/math/BigInteger.java

(Somehow they stopped caching bitlength, numberOfLeadingZeros
when CPU compiled is just some register ops, nevertheless I
wouldn't call bitLength() in a loop too often)

burs...@gmail.com schrieb:

burs...@gmail.com

unread,
Feb 7, 2017, 7:01:45 AM2/7/17
to
BTW: bitlength is one off msb. Recently I used, when
porting from Jekejeke to SWI-Prolog:

bitlength(X, Y) :-
Y is msb(X)+1.

Kaz Kylheku

unread,
Feb 7, 2017, 10:16:34 AM2/7/17
to
["Followup-To:" header set to comp.theory.]
On 2017-02-07, Julio P. Di Egidio <ju...@diegidio.name> wrote:
> [Cross-posted to comp.theory, sci.math, comp.lang.prolog. Please set
> follow-ups as needed.]
>
> Hello all,
>
> I have put together an algorithm (plus reference implementation in Prolog)
> for the arbitrary precision integer binary logarithm with remainder: I
> needed this myself but could not find anything on the web, so I thought I'd
> share.

Pragmatically speaking, you probably won't run into an implementation of
arbitrary precision ("bignum") integers in which you cannot obtain an
estimate of this value from the implementation.

E.g. if a bignum consists of seventeen limbs of sixty-four bits,
then it consists of at least 16 x 64 bits, and you can apply a fast
algorithm on just the most significant 64 bit word to determine the
precise value.

Below is some C code for that which I hadded to a bignum library
called MPI (written by M. Fromberger).

We implement a binary search with bitmasks, rather than shifting,
which is practical due to the fixed size. The binary search tree
is expanded into open coded cases, so there is no loop, and all
the masks are literal constants:

static int s_highest_bit(mp_digit n)
{
#if MP_DIGIT_SIZE == 8
if (n & 0xFFFFFFFF00000000) {
if (n & 0xFFFF000000000000) {
if (n & 0xFF00000000000000) {
if (n & 0xF000000000000000) {
if (n & 0xC000000000000000)
return (n & 0x8000000000000000) ? 64 : 63;
else
return (n & 0x2000000000000000) ? 62 : 61;
} else {
if (n & 0x0C00000000000000)
return (n & 0x0800000000000000) ? 60 : 59;
else
return (n & 0x0200000000000000) ? 58 : 57;
}
} else {
if (n & 0x00F0000000000000) {
if (n & 0x00C0000000000000)
return (n & 0x0080000000000000) ? 56 : 55;
else
return (n & 0x0020000000000000) ? 54 : 53;
} else {
if (n & 0x000C000000000000)
return (n & 0x0008000000000000) ? 52 : 51;
else
return (n & 0x0002000000000000) ? 50 : 49;
}
}
} else {
if (n & 0x0000FF0000000000) {
if (n & 0x0000F00000000000) {
if (n & 0x0000C00000000000)
return (n & 0x0000800000000000) ? 48 : 47;
else
return (n & 0x0000200000000000) ? 46 : 45;
} else {
if (n & 0x00000C0000000000)
return (n & 0x0000080000000000) ? 44 : 43;
else
return (n & 0x0000020000000000) ? 42 : 41;
}
} else {
if (n & 0x000000F000000000) {
if (n & 0x000000C000000000)
return (n & 0x0000008000000000) ? 40 : 39;
else
return (n & 0x0000002000000000) ? 38 : 37;
} else {
if (n & 0x0000000C00000000)
return (n & 0x0000000800000000) ? 36 : 35;
else
return (n & 0x0000000200000000) ? 34 : 33;
}
}
}
} else {
if (n & 0x00000000FFFF0000) {
if (n & 0x00000000FF000000) {
if (n & 0x00000000F0000000) {
if (n & 0x00000000C0000000)
return (n & 0x0000000080000000) ? 32 : 31;
else
return (n & 0x0000000020000000) ? 30 : 29;
} else {
if (n & 0x000000000C000000)
return (n & 0x0000000008000000) ? 28 : 27;
else
return (n & 0x0000000002000000) ? 26 : 25;
}
} else {
if (n & 0x0000000000F00000) {
if (n & 0x0000000000C00000)
return (n & 0x0000000000800000) ? 24 : 23;
else
return (n & 0x0000000000200000) ? 22 : 21;
} else {
if (n & 0x00000000000C0000)
return (n & 0x0000000000080000) ? 20 : 19;
else
return (n & 0x0000000000020000) ? 18 : 17;
}
}
} else {
if (n & 0x000000000000FF00) {
if (n & 0x000000000000F000) {
if (n & 0x000000000000C000)
return (n & 0x0000000000008000) ? 16 : 15;
else
return (n & 0x0000000000002000) ? 14 : 13;
} else {
if (n & 0x0000000000000C00)
return (n & 0x0000000000000800) ? 12 : 11;
else
return (n & 0x0000000000000200) ? 10 : 9;
}
} else {
if (n & 0x00000000000000F0) {
if (n & 0x00000000000000C0)
return (n & 0x0000000000000080) ? 8 : 7;
else
return (n & 0x0000000000000020) ? 6 : 5;
} else {
if (n & 0x000000000000000C)
return (n & 0x0000000000000008) ? 4 : 3;
else
return (n & 0x0000000000000002) ? 2 : (n ? 1 : 0);
}
}
}
}
#elif MP_DIGIT_SIZE == 4
if (n & 0xFFFF0000) {
if (n & 0xFF000000) {
if (n & 0xF0000000) {
if (n & 0xC0000000)
return (n & 0x80000000) ? 32 : 31;
else
return (n & 0x20000000) ? 30 : 29;
} else {
if (n & 0x0C000000)
return (n & 0x08000000) ? 28 : 27;
else
return (n & 0x02000000) ? 26 : 25;
}
} else {
if (n & 0x00F00000) {
if (n & 0x00C00000)
return (n & 0x00800000) ? 24 : 23;
else
return (n & 0x00200000) ? 22 : 21;
} else {
if (n & 0x000C0000)
return (n & 0x00080000) ? 20 : 19;
else
return (n & 0x00020000) ? 18 : 17;
}
}
} else {
if (n & 0x0000FF00) {
if (n & 0x0000F000) {
if (n & 0x0000C000)
return (n & 0x00008000) ? 16 : 15;
else
return (n & 0x00002000) ? 14 : 13;
} else {
if (n & 0x00000C00)
return (n & 0x00000800) ? 12 : 11;
else
return (n & 0x00000200) ? 10 : 9;
}
} else {
if (n & 0x000000F0) {
if (n & 0x000000C0)
return (n & 0x00000080) ? 8 : 7;
else
return (n & 0x00000020) ? 6 : 5;
} else {
if (n & 0x0000000C)
return (n & 0x00000008) ? 4 : 3;
else
return (n & 0x00000002) ? 2 : (n ? 1 : 0);
}
}
}
#elif MP_DIGIT_SIZE == 2
if (n & 0xFF00) {
if (n & 0xF000) {
if (n & 0xC000)
return (n & 0x8000) ? 16 : 15;
else
return (n & 0x2000) ? 14 : 13;
} else {
if (n & 0x0C00)
return (n & 0x0800) ? 12 : 11;
else
return (n & 0x0200) ? 10 : 9;
}
} else {
if (n & 0x00F0) {
if (n & 0x00C0)
return (n & 0x0080) ? 8 : 7;
else
return (n & 0x0020) ? 6 : 5;
} else {
if (n & 0x000C)
return (n & 0x0008) ? 4 : 3;
else
return (n & 0x0002) ? 2 : (n ? 1 : 0);
}
}
#elif MP_DIGIT_SIZE == 1
if (n & 0xF0) {
if (n & 0xC0)
return (n & 0x80) ? 8 : 7;
else
return (n & 0x20) ? 6 : 5;
} else {
if (n & 0x0C)
return (n & 0x08) ? 4 : 3;
else
return (n & 0x02) ? 2 : (n ? 1 : 0);
}
#else
#error fixme: unsupported MP_DIGIT_SIZE
#endif
/* notreached */
abort();
}

With this, the function for the bits in the whole bignum is
just:

int s_highest_bit_mp(mp_int *a)
{
int nd1 = USED(a) - 1;
return s_highest_bit(DIGIT(a, nd1)) + nd1 * MP_DIGIT_BIT;
}

Ben Bacarisse

unread,
Feb 7, 2017, 12:11:54 PM2/7/17
to
"Julio P. Di Egidio" <ju...@diegidio.name> writes:

> I have put together an algorithm (plus reference implementation in
> Prolog) for the arbitrary precision integer binary logarithm with
> remainder: I needed this myself but could not find anything on the
> web, so I thought I'd share.

Thanks for sharing. I have a few remarks...

> Three Prolog predicates are in fact offered, which are progressive
> refinements of the naive approach that shifts the input to the right
> by one bit at a time:
>
> ilog2_naive(+X, -Y, -R) :
> Naive implementation: the input is shifted to the right by 1 bit at a
> time repeatedly until the result of shifting is zero. The total number
> of repetitions corresponds to the ceiling of the real logarithm. The
> complexity is O(log2(X^3).

All of you complexities are of the form O(log(x^C)). These are all
O(log x). If you want to capture the differences between the
algorithms, you'll need something more precise than big-O. Also, it
really helps if you say what is being counted. Here, I think you are
counting arithmetic operations (or maybe just shifts?). The actual time
complexity might depend on how those are implemented.

Have you tried a binary search? Double log until (1<<log) > x and then
search by halving the interval to find the log. It may not pay off, of
course but seems a reasonable algorithm.

This may have been done for other reasons, but it's a shame if it's
really needed. Internally Prolog will probably know the log2 of a
number! I don't think Haskell provides this function either. Seems
like an omission. Oh well...

<snip>
--
Ben.

Julio P. Di Egidio

unread,
Feb 7, 2017, 3:23:10 PM2/7/17
to
"Kaz Kylheku" wrote in message news:201702070...@kylheku.com...
> On 2017-02-07, Julio P. Di Egidio <ju...@diegidio.name> wrote:
<snip>
> > I have put together an algorithm (plus reference implementation in
> > Prolog)
> > for the arbitrary precision integer binary logarithm with remainder: I
> > needed this myself but could not find anything on the web, so I thought
> > I'd
> > share.
>
> Pragmatically speaking, you probably won't run into an implementation of
> arbitrary precision ("bignum") integers in which you cannot obtain an
> estimate of this value from the implementation.

Yes, I'd agree on that, which may very well be the reason why all bit
fiddling hacks are in fixed precision, and why my
enterprise here is almost pointless... Although, at least at a first
cursory inspection (benchmarking against the msb/1 arith
function in SWI, itself on top of GMP), my approach seem to fair pretty well
nevertheless, which I find interesting: I may very
well be wrong but I have a hunch that my algorithm is somehow optimal for
the general case. This was one of my couple of
questions anyway...

> Below is some C code for that which I hadded to a bignum library
> called MPI (written by M. Fromberger).
>
> We implement a binary search with bitmasks, rather than shifting,
> which is practical due to the fixed size.

But I am talking arbitrary precision: fixed size is easy.

Julio

Julio P. Di Egidio

unread,
Feb 7, 2017, 3:32:59 PM2/7/17
to
"Ben Bacarisse" wrote in message news:87lgtia...@bsb.me.uk...

> All of you complexities are of the form O(log(x^C)). These are all
> O(log x). If you want to capture the differences between the
> algorithms, you'll need something more precise than big-O.

Yes, I do want to capture the differences; and no, I don’t know exactly what
to use instead of big O, otherwise I'd not be asking.

Julio

Julio P. Di Egidio

unread,
Feb 7, 2017, 6:11:25 PM2/7/17
to
"Ben Bacarisse" wrote in message news:87lgtia...@bsb.me.uk...
<snip>
> All of you complexities are of the form O(log(x^C)). These are all
> O(log x). If you want to capture the differences between the
> algorithms, you'll need something more precise than big-O. Also, it
> really helps if you say what is being counted. Here, I think you are
> counting arithmetic operations (or maybe just shifts?). The actual time
> complexity might depend on how those are implemented.

That is a part I most need help with. I am indeed counting the arithmetic
operations, not finer than that, and that is all there is beside the control
flow (just note that in some places unification in Prolog is hiding the
assignment of return values). Anyway, below is a clean version of the
algorithm for easier reading: it is only for ilog2_n/4, the others should
then be straightforward.

Please be aware that I have also found a significant error in my
documentation: wherever I say "the input is shifted to the right by 2^M bits
at a time, with M set equal to K initially, and progressively decremented",
I should rather be saying "the input is shifted to the right by 2^2^N bits
at a time, with N set equal to K initially, and progressively decremented".

----------------------------------------
ilog2_n(in k:int>=0, in x:int>=0, out y:int, out r:int>=0)
if (x == 0)
y = -1
r = 0
else
m = 1 << K % M = 2^K
ilog2_n__do(m, x, 0, y)
r = x - (1 << y) % R = X - 2^Y

ilog2_n__do(in m:int>0, in x0:int>=0, in y0:int>=0, out y:int>=0)
x1 = x0 >> M % X1 = X0 // 2^M
if (x1 =\= 0)
y1 = y0 + m % Y1 = Y0 + M
ilog2_n__do(m, x1, y1, y)
else if (m =\= 1)
m1 = m >> 1 % M1 = M // 2
ilog2_n__do(m1, x0, y0, y)
else
y = y0
----------------------------------------

> Have you tried a binary search? Double log until (1<<log) > x and then
> search by halving the interval to find the log. It may not pay off, of
> course but seems a reasonable algorithm.

The naive approach is even simpler than that: "[t]he total number of
repetitions [already] corresponds to the ceiling of the real logarithm".
The refinement is to do it in larger chunks instead of one bit at a time, in
fact in chunks as large as possible with the shifted value decreasing. The
approach I am trying essentially reduces out chunks of 2^2^k for decreasing
k, where k initial can be fine tuned to the input, i.e. for anything about
the input distribution that is known. What an optimal value of k is when no
information at all is available I suppose could be interesting to the worst
case analysis...

> Internally Prolog will probably know the log2 of a number!

SWI-Prolog provides arbitrary precision big integers on top of GMP and in
fact it does provide an integer binary log (called msb).

Julio

Julio P. Di Egidio

unread,
Feb 7, 2017, 6:17:53 PM2/7/17
to
"Julio P. Di Egidio" wrote in message news:o7dk3h$4cr$1...@dont-email.me...

> > Have you tried a binary search? Double log until (1<<log) > x and then
> > search by halving the interval to find the log. It may not pay off, of
> > course but seems a reasonable algorithm.
>
> The naive approach is even simpler than that

Ah, sorry, now I see what you meant: that's a good idea, I will try it...

Julio

Peter Percival

unread,
Feb 7, 2017, 6:30:14 PM2/7/17
to
Julio P. Di Egidio wrote:
> [...]
> Please be aware that I have also found a significant error in my
> documentation: wherever I say "the input is shifted to the right by 2^M
> bits at a time, with M set equal to K initially, and progressively
> decremented", I should rather be saying "the input is shifted to the
> right by 2^2^N bits at a time, with N set equal to K initially, and
> progressively decremented".

Has that been updated in GitHub?


--
Do, as a concession to my poor wits, Lord Darlington, just explain
to me what you really mean.
I think I had better not, Duchess. Nowadays to be intelligible is
to be found out. -- Oscar Wilde, Lady Windermere's Fan

Ben Bacarisse

unread,
Feb 7, 2017, 7:25:28 PM2/7/17
to
"Julio P. Di Egidio" <ju...@diegidio.name> writes:

Ah, I was not sure you were asking about that. I think your best bet is
simple to try to find a function that is an upper bound on the number of
arithmetic operations involved. You'd end with statements like "the
number of operations is no more than 3*log2(x)".

--
Ben.

James Waldby

unread,
Feb 8, 2017, 3:33:55 PM2/8/17
to
On Tue, 07 Feb 2017 21:23:03 +0100, Julio P. Di Egidio wrote:

> "Kaz Kylheku" wrote in message news:201702070...@kylheku.com...
...
>> Below is some C code for that which I added to a bignum library
>> called MPI (written by M. Fromberger).
>>
>> We implement a binary search with bitmasks, rather than shifting,
>> which is practical due to the fixed size.
>
> But I am talking arbitrary precision: fixed size is easy.

The numbers treated by Kaz's code are arbitrary precision.

The "fixed size" referred to is the fixed size of units (words,
dwords, etc) used to contain numbers. Each number occupies an
arbitrary number of units.

--
jiw

Kaz Kylheku

unread,
Feb 8, 2017, 7:57:42 PM2/8/17
to
Exactly, so that if we have the number

0x07FFFFFF 09DFD96C AAABDC99 ... 00A071D7
-------- ------------------------------
27 bits + 17 words x 32 bits = 544 bits wide

We need to solve the problem only for the highest order limb of the
bignum integer, then add that value to the bit width of the array below
that, trivially obtained.

Not saying that it's not a valid theoretical problem, or that
it couldn't arise in practice.

(Say we are working in the "user space" of some language whose "kernel"
is not accessible for extension (closed source, or not easy to upstream
patches, or, in the modern, packaged OSS world, to have them
disseminated everywhere where we care our code to run in the short term,
etc)).

Ben's algorithm seems like an excellent start: keep doubling some intial
guess until it exceeds the value, then binary search.


--
TXR Programming Lanuage: http://nongnu.org/txr
Music DIY Mailing List: http://www.kylheku.com/diy
ADA MP-1 Mailing List: http://www.kylheku.com/mp1

Ben Bacarisse

unread,
Feb 8, 2017, 10:05:27 PM2/8/17
to
Kaz Kylheku <336-98...@kylheku.com> writes:
<snip>
> Ben's algorithm seems like an excellent start: keep doubling some intial
> guess until it exceeds the value, then binary search.

My Prolog is not up to writing a comparative test, but here's the binary
search and the naive algorithm in Haskell. The more complex algorithm
is not much slower for small numbers but very much better for large
ones.

import System.Environment

ilog2 :: Integer -> Int
ilog2 x = ilog_search (bounds 2 (0, 1))
where
bounds pow (lo, hi)
| pow >= x = (lo, hi)
| otherwise = bounds (pow*pow) (hi, hi*2)
ilog_search (lo, hi)
| hi <= lo = lo
| 2^mid >= x = ilog_search (lo, mid)
| otherwise = ilog_search ((mid+1), hi)
where mid = (lo + hi) `div` 2

ilog2_n :: Integer -> Int
ilog2_n x = ilog2_n' 0 1
where ilog2_n' lg pow | pow >= x = lg
| otherwise = ilog2_n' (lg+1) (pow*2)

main = do
[n', f'] <- getArgs
let n = read n'
-- decide which function to call
let f = if f' == "n" then ilog2_n else ilog2
-- use somewhat arbitrary list of large numbers to test
print $ sum [f (3^x) | x <- [n..2*n]]


--
Ben.

Julio P. Di Egidio

unread,
Feb 8, 2017, 11:54:05 PM2/8/17
to
"Ben Bacarisse" wrote in message news:87tw85a...@bsb.me.uk...
That's great, thank you.

Let me try and ask you a further one: would you count the arithmetic
operations only or would you also count the branching operations? -- As I
see it, the loops can be completely unwinded, but not so for the
conditionals, so I am not clear what the right approach should be: of
course, relative to some canonical/technically correct way to do an
algorithmic analysis here. I also understand I have not made a clear
distinction between the algorithm and the implementation: again, any advise
as to how to straighten things up (in this case, of course I am not asking
you to give a full course in CS) would be much appreciated.

Julio

Julio P. Di Egidio

unread,
Feb 9, 2017, 12:19:22 AM2/9/17
to
"Kaz Kylheku" wrote in message news:201702081...@kylheku.com...
> On 2017-02-08, James Waldby <n...@valid.invalid> wrote:
> > On Tue, 07 Feb 2017 21:23:03 +0100, Julio P. Di Egidio wrote:
> >> "Kaz Kylheku" wrote in message news:201702070...@kylheku.com...
> > ...
> >>> Below is some C code for that which I added to a bignum library
> >>> called MPI (written by M. Fromberger).
> >>>
> >>> We implement a binary search with bitmasks, rather than shifting,
> >>> which is practical due to the fixed size.
> >>
> >> But I am talking arbitrary precision: fixed size is easy.
> >
> > The numbers treated by Kaz's code are arbitrary precision.
> >
> > The "fixed size" referred to is the fixed size of units (words,
> > dwords, etc) used to contain numbers. Each number occupies an
> > arbitrary number of units.
>
> Exactly, so that if we have the number
>
> 0x07FFFFFF 09DFD96C AAABDC99 ... 00A071D7
> -------- ------------------------------
> 27 bits + 17 words x 32 bits = 544 bits wide
>
> We need to solve the problem only for the highest order limb of the
> bignum integer, then add that value to the bit width of the array below
> that, trivially obtained.

But we do not have that information. Namely, assume here we do not have any
access or even any info about the internal structure of the numbers.

> Not saying that it's not a valid theoretical problem, or that
> it couldn't arise in practice.

I do have started this mostly for an "academic" interest (I have provided a
reference implementation, but the point was an algorithm), which I will
elaborate a little bit upon below: in the normal case I think your point is
indeed valid, that it simply should not arise!

> (Say we are working in the "user space" of some language whose "kernel"
> is not accessible for extension (closed source, or not easy to upstream
> patches, or, in the modern, packaged OSS world, to have them
> disseminated everywhere where we care our code to run in the short term,
> etc)).

I am finding just two aspects might be of interest from what I have
presented:

The first has to do with mathematical foundations: in that context, e.g. the
"natural numbers" as the prototypical most primitive set of labels
("indexes") are given with no intrinsic structure, just some order/successor
relation...

Another point I am finding of interest is that my implementation (and by
incorporating Ben's suggestion, I am pretty sure this will improve) seems to
perform almost as well as corresponding the built-in: this I am finding at
least puzzling, maybe a hint to some kind of optimality (to be sought
for)...

> Ben's algorithm seems like an excellent start: keep doubling some intial
> guess until it exceeds the value, then binary search.

Ben's idea is good because it might solve the problem of finding initial k:
we are all already doing "binary search", the point is the size of the mask.

Julio

Julio P. Di Egidio

unread,
Feb 9, 2017, 12:55:31 AM2/9/17
to
"Julio P. Di Egidio" wrote in message news:o7gsi0$pdn$1...@dont-email.me...

> the loops can be completely unwinded, but not so for the conditionals, so
> I am not clear what the right approach should be

Silly me, of course the conditionals can be unwinded just as well, and this
I suppose answers the question. But any further advice still welcome.

Julio

Ben Bacarisse

unread,
Feb 9, 2017, 6:12:08 AM2/9/17
to
I don't follow this, but then I'm not entirely sure I know what the
algorithms (other than the naive one) are. The fact that you've put
binary search in scare quotes does not help because I'm using in the
normal sense of homing in on a result by halving the possible range.
Are you using it in some other way? (Obviously, because this is
log2(x), linear changes in the log cover exponential ranges for the
parameter x, but that's not a binary search.)

Just to be clear, here's an example. To find ilog2(1000):

Get first of l = 1, 2, 4, 8, 16, ... such that 2^l >= 1000.
This gives 16 so 8 < ilog2(1000) <= 16.

Do binary search in 8..16:
Mid point is 12. 2^12 >= 1000 so search in 8..12
Mid point is 10. 2^10 >= 1000 so search in 8..10
Mid point is 9. 2^9 < 1000 so search in 10..10
Range is now exhausted so return 10.

It's a shame I don't know Prolog well enough to implement the idea for a
direct comparison. I've written som Prolog but so long ago it would
take me a while to get used to inside-out style again. I've posted
Haskell code for both the binary search method and the naive search, but
you may not know Haskell well enough to implement your improved methods
in that.

--
Ben.

Ben Bacarisse

unread,
Feb 9, 2017, 6:17:40 AM2/9/17
to
"Julio P. Di Egidio" <ju...@diegidio.name> writes:

> "Ben Bacarisse" wrote in message news:87tw85a...@bsb.me.uk...
>> "Julio P. Di Egidio" <ju...@diegidio.name> writes:
>> > "Ben Bacarisse" wrote in message news:87lgtia...@bsb.me.uk...
>> >
>> >> All of you complexities are of the form O(log(x^C)). These are all
>> >> O(log x). If you want to capture the differences between the
>> >> algorithms, you'll need something more precise than big-O.
>> >
>> > Yes, I do want to capture the differences; and no, I don’t know
>> > exactly what to use instead of big O, otherwise I'd not be asking.
>>
>> Ah, I was not sure you were asking about that. I think your best bet is
>> simple to try to find a function that is an upper bound on the number of
>> arithmetic operations involved. You'd end with statements like "the
>> number of operations is no more than 3*log2(x)".
>
> That's great, thank you.
>
> Let me try and ask you a further one: would you count the arithmetic
> operations only or would you also count the branching operations?

I'd just count arithmetic.

> I also understand I have not made a clear
> distinction between the algorithm and the implementation: again, any
> advise as to how to straighten things up (in this case, of course I am
> not asking you to give a full course in CS) would be much appreciated.

That's a hard one. Once you leave the broad-brush analysis of big-O
it's not always clear if you are analysing an implementation or
something more abstract. The problem crops up even with big-O notation
but then it's often much easier to see that any variations in
implementation won't alter the asymptotic complexity. In the end, it
doesn't matter much as you are interested in fast integer logs. What
you want is a fast implementation not a fast algorithm, so analyse that.

--
Ben.

burs...@gmail.com

unread,
Feb 9, 2017, 6:58:44 AM2/9/17
to
Hi,

I am just considering using ilog2 for some range reduction
in exact reals. For n = floor(log2(x)) and x>=1, we have
the nice property 2^n+1 > x >= 2^n.

But I would go with some builtin such as msb/1(*) or
bitlength/1 which are O(1) namely constant speed. On
the other hand a Prolog realization and not a low level

realization could serve as a specification, and help
in reasoning about programs. But then if used as a
specification, it wouldn't matter so much if it

were unefficient, at least if we are able to switch
to msb/1 for the realization. Another problem would
be if we would need ilog3(x),

then msb/1 doesn't offer a direct path, I wonder
how one would go about for an efficient floor(log3(x))
in big nums. Ideas presented here could be used.

(*)SWI-Prolog Function msb/1
http://www.swi-prolog.org/pldoc/man?function=msb/1

Bye

Am Donnerstag, 9. Februar 2017 12:17:40 UTC+1 schrieb Ben Bacarisse:

Julio P. Di Egidio

unread,
Feb 10, 2017, 1:21:31 PM2/10/17
to
On 09/02/2017 12:12, Ben Bacarisse wrote:
> "Julio P. Di Egidio" <ju...@diegidio.name> writes:
>> "Kaz Kylheku" wrote in message news:201702081...@kylheku.com...
<snip>
>>> Ben's algorithm seems like an excellent start: keep doubling some intial
>>> guess until it exceeds the value, then binary search.
>>
>> Ben's idea is good because it might solve the problem of finding
>> initial k: we are all already doing "binary search", the point is the
>> size of the mask.
>
> I don't follow this, but then I'm not entirely sure I know what the
> algorithms (other than the naive one) are.

At some point I have replied to you with the algorithm in pseudo-code:
but I call it ilog2_n, you have called your naive one ilog2_n.

> The fact that you've put
> binary search in scare quotes does not help because I'm using in the
> normal sense of homing in on a result by halving the possible range.
> Are you using it in some other way? (Obviously, because this is
> log2(x), linear changes in the log cover exponential ranges for the
> parameter x, but that's not a binary search.)

Those were not scary quotes: all the approaches, beside the naive one,
use a form or other of binary search, in fact in this case a two-level
one: but even the simple masking tables in fixed precision are an
unwinding of a binary search, i.e. one needn't literally have lo, hi and
mid: in fact the point is to optimise down to the fewest and simplest
operations, and, for example, while your idea is in fact good, your code
is terribly inefficient... :)

> It's a shame I don't know Prolog well enough

Never mind, as an implementation problem the answer just is use the
built-in, then the question about mathematical foundations I'll leave
for another day... False alarm and case closed, as far as I can see.

Thanks again for your help,

Julio

Ben Bacarisse

unread,
Feb 10, 2017, 5:19:12 PM2/10/17
to
"Julio P. Di Egidio" <ju...@diegidio.name> writes:

> On 09/02/2017 12:12, Ben Bacarisse wrote:
>> "Julio P. Di Egidio" <ju...@diegidio.name> writes:
>>> "Kaz Kylheku" wrote in message news:201702081...@kylheku.com...
> <snip>
>>>> Ben's algorithm seems like an excellent start: keep doubling some intial
>>>> guess until it exceeds the value, then binary search.
>>>
>>> Ben's idea is good because it might solve the problem of finding
>>> initial k: we are all already doing "binary search", the point is the
>>> size of the mask.
>>
>> I don't follow this, but then I'm not entirely sure I know what the
>> algorithms (other than the naive one) are.
>
> At some point I have replied to you with the algorithm in pseudo-code:
> but I call it ilog2_n, you have called your naive one ilog2_n.

Yes, I said, I understand that one. It's the other with the same big-O
complexity that was not sure about.

>> The fact that you've put
>> binary search in scare quotes does not help because I'm using in the
>> normal sense of homing in on a result by halving the possible range.
>> Are you using it in some other way? (Obviously, because this is
>> log2(x), linear changes in the log cover exponential ranges for the
>> parameter x, but that's not a binary search.)
>
> Those were not scary quotes: all the approaches, beside the naive one,
> use a form or other of binary search,

That seems unlikely unless your analysis was wrong. They all had the
same completiy and that should not be the case if all but one use a
binary search of some sort (unless they don't use it to any effect but
that would be even less likely).

> in fact in this case a two-level
> one: but even the simple masking tables in fixed precision are an
> unwinding of a binary search, i.e. one needn't literally have lo, hi
> and mid: in fact the point is to optimise down to the fewest and
> simplest operations, and, for example, while your idea is in fact
> good, your code is terribly inefficient... :)

Is it? I can't compare it with anything of yours at the moment. Can
you tell me the problem you see with it? Obviously it's fast compared
to the Haskell naive version but presumably you see that as inefficient
too.

It uses 2^x rather than a bit shift (as your code does) but that's just
to avoid about using another module. It's much faster (though with the
same complexity measure in terms of arithmetic operation count) when a
bit shoft is used, but it was written for clarity.

>> It's a shame I don't know Prolog well enough
>
> Never mind, as an implementation problem the answer just is use the
> built-in, then the question about mathematical foundations I'll leave
> for another day... False alarm and case closed, as far as I can see.

Oh! It would have been better had you said that when using a built-in
operation was first suggested. I was hoping to compare
implementations. Oh well...

--
Ben.

Julio P. Di Egidio

unread,
Feb 11, 2017, 3:04:15 AM2/11/17
to
On 10/02/2017 23:19, Ben Bacarisse wrote:
> "Julio P. Di Egidio" <ju...@diegidio.name> writes:
>> On 09/02/2017 12:12, Ben Bacarisse wrote:
>>> "Julio P. Di Egidio" <ju...@diegidio.name> writes:
>>>> "Kaz Kylheku" wrote in message news:201702081...@kylheku.com...
>> <snip>
>>>>> Ben's algorithm seems like an excellent start: keep doubling some intial
>>>>> guess until it exceeds the value, then binary search.
>>>>
>>>> Ben's idea is good because it might solve the problem of finding
>>>> initial k: we are all already doing "binary search", the point is the
>>>> size of the mask.
>>>
>>> I don't follow this, but then I'm not entirely sure I know what the
>>> algorithms (other than the naive one) are.
>>
>> At some point I have replied to you with the algorithm in pseudo-code:
>> but I call it ilog2_n, you have called your naive one ilog2_n.
>
> Yes, I said, I understand that one. It's the other with the same big-O
> complexity that was not sure about.

I am annoyed by your tone. Yes, you do not understand, which I must
think is because you cannot be bothered to read carefully and have to
play the smart ass instead... But I have given names to the things I
have presented, it is not my fault if you cannot be bothered to follow.

>>> The fact that you've put
>>> binary search in scare quotes does not help because I'm using in the
>>> normal sense of homing in on a result by halving the possible range.
>>> Are you using it in some other way? (Obviously, because this is
>>> log2(x), linear changes in the log cover exponential ranges for the
>>> parameter x, but that's not a binary search.)
>>
>> Those were not scary quotes: all the approaches, beside the naive one,
>> use a form or other of binary search,
>
> That seems unlikely unless your analysis was wrong. They all had the
> same completiy and that should not be the case if all but one use a
> binary search of some sort (unless they don't use it to any effect but
> that would be even less likely).
>
>> in fact in this case a two-level
>> one: but even the simple masking tables in fixed precision are an
>> unwinding of a binary search, i.e. one needn't literally have lo, hi
>> and mid: in fact the point is to optimise down to the fewest and
>> simplest operations, and, for example, while your idea is in fact
>> good, your code is terribly inefficient... :)
>
> Is it? I can't compare it with anything of yours at the moment.

Again, it's not my fault if you cannot compare: I am not proficient in
Haskell but I can read your code.

> Can
> you tell me the problem you see with it? Obviously it's fast compared
> to the Haskell naive version but presumably you see that as inefficient
> too.
>
> It uses 2^x rather than a bit shift (as your code does) but that's just
> to avoid about using another module. It's much faster (though with the
> same complexity measure in terms of arithmetic operation count) when a
> bit shoft is used, but it was written for clarity.

It's not just that: the main thing is you bring around the exponents and
redo the powering at every iteration, which is extremely inefficient not
only because of taking all those powers but also and primarily because
of the possibly large size of the operands which means creating large
numbers at every operation/assignment. That is the reason why, not only
I rather bring around the power itself, not just the exponent, I also
strictly reduce the input along the way, which is not terribly efficient
with small numbers, but is way better with the bigger ones.

>>> It's a shame I don't know Prolog well enough
>>
>> Never mind, as an implementation problem the answer just is use the
>> built-in, then the question about mathematical foundations I'll leave
>> for another day... False alarm and case closed, as far as I can see.
>
> Oh! It would have been better had you said that when using a built-in
> operation was first suggested. I was hoping to compare
> implementations. Oh well...

Well, I just do not see what interesting there is here left to discover,
do you? But if you clarify what you would like to test or compare, I'd
be glad to do the coding effort...

Julio

burs...@gmail.com

unread,
Feb 11, 2017, 7:54:43 AM2/11/17
to
How would one solve this sub problem, the result
of rlog2 is not anymore a natural number >=0,
but an integer, and the input is a positive rational.

Mathematical specification, floor rounds
towards -infinity, it is not same as trunc
which rounds towards zero:

rlog2(x) = floor(log2(x))

Example:
rlog2(1000 rdiv 1) = 9
rlog2(1 rdiv 3) = -2

How would one produce some code, assuming ilog2/1
is already available. The input is a rational as
a compound rdiv(P,Q).

Am Samstag, 11. Februar 2017 09:04:15 UTC+1 schrieb Julio Di Egidio:

Ben Bacarisse

unread,
Feb 11, 2017, 8:24:56 AM2/11/17
to
"Julio P. Di Egidio" <ju...@diegidio.name> writes:
> On 10/02/2017 23:19, Ben Bacarisse wrote:
<snip>
>> Yes, I said, I understand that one. It's the other with the same big-O
>> complexity that was not sure about.
>
> I am annoyed by your tone.

Sorry about that. Can you get past it and talk about the technicalities?

>> "Julio P. Di Egidio" <ju...@diegidio.name> writes:
<snip>
>>> Never mind, as an implementation problem the answer just is use the
>>> built-in, then the question about mathematical foundations I'll leave
>>> for another day... False alarm and case closed, as far as I can see.
>>
>> Oh! It would have been better had you said that when using a built-in
>> operation was first suggested. I was hoping to compare
>> implementations. Oh well...
>
> Well, I just do not see what interesting there is here left to
> discover, do you? But if you clarify what you would like to test or
> compare, I'd be glad to do the coding effort...

Yes. I think an efficient algorithm for ilog2 is an interesting problem
because you don't always have access to the underlying "big num"
operations. I'd like to see a few algorithms implemented in a way that
allows their performance can be measured. I'd like to know if binary
search is indeed worth while. One would have to agree some base set of
operations (do we assume that bit shifts and tests possible, or do we
need to do multiplications, divisions and exp operations?), but all the
algorithms will probably show the same relative performance given the
same agreed set of basic operations.

I think doing a binary search *on the exponent* will be faster (for big
numbers) than what I understand you are doing. I think the reason you
say you are doing a binary search is because a linear change in the
exponent corresponds to an exponential change in the power, but I am
proposing a binary search on the exponent itself -- halving the exponent
range at each step.

I, too, would be happy to do more coding, but I could not understand the
details of your other algorithms (which is obviously my fault -- if I
knew Prolog well they would be unambiguous). They seem to do block
shifts first and then single-step shifts to home in on the log. That's
going to have the same complexity as the naive algorithm (which I did
understand) but, as I say, that may not be what they do.

Do we have a programming language in common? I can do almost any
conventional or functional language (for which a Linux implementation is
available). Prolog is just about the only gap I have! Is it the only
programming language you work in?

--
Ben.

Julio P. Di Egidio

unread,
Feb 11, 2017, 9:51:19 AM2/11/17
to
On 11/02/2017 14:24, Ben Bacarisse wrote:

> Do we have a programming language in common? I can do almost any
> conventional or functional language (for which a Linux implementation is
> available).

What about JavaScript? I am not inclined in trying Haskell because one
thing is reading it, another is writing efficient code.

> Prolog is just about the only gap I have! Is it the only
> programming language you work in?

I have stopped counting how many languages I have worked with...
Anyway, beside JS, functional is the domain I have explored less.

Julio

Ben Bacarisse

unread,
Feb 11, 2017, 12:22:06 PM2/11/17
to
"Julio P. Di Egidio" <ju...@diegidio.name> writes:

> On 11/02/2017 14:24, Ben Bacarisse wrote:
>
>> Do we have a programming language in common? I can do almost any
>> conventional or functional language (for which a Linux implementation is
>> available).
>
> What about JavaScript? I am not inclined in trying Haskell because
> one thing is reading it, another is writing efficient code.

JS would be fine though rather messy since JS does not have arbitrary
precision arithmetic built in. Which library should I use? I suggest
node.js's bignumber because it's easy to install, but it is rather slow.
Mind you, that may be an advantage since it will highlight the algorithm
rather than the arithmetic.

>> Prolog is just about the only gap I have! Is it the only
>> programming language you work in?
>
> I have stopped counting how many languages I have worked with...
> Anyway, beside JS, functional is the domain I have explored less.

Do you know any that have built-in arbitrary precision arithmetic such
as Python? The algorithms are clearer when you can write

b = 98765432123456789
print(b * 2)

rather than

var b = new BigNumber('98765432123456789');
console.log(b.times(2).toString());

Lots of the the others are at least partially functional (Scheme, ML,
Lisp). Even the Unix bc command's language might be less messy than JS
with a big number library. But if the only common language is JS,
that's fine.

--
Ben.

Julio P. Di Egidio

unread,
Feb 11, 2017, 1:59:53 PM2/11/17
to
On 11/02/2017 18:22, Ben Bacarisse wrote:
> "Julio P. Di Egidio" <ju...@diegidio.name> writes:
<snip>
>> What about JavaScript? I am not inclined in trying Haskell because
>> one thing is reading it, another is writing efficient code.
>
> JS would be fine though rather messy since JS does not have arbitrary
> precision arithmetic built in. Which library should I use? I suggest
> node.js's bignumber because it's easy to install, but it is rather slow.

You are right, JS won't do: bignums are too slow for anything but toy
applications, I have even tried in the past...

> Mind you, that may be an advantage since it will highlight the algorithm
> rather than the arithmetic.

In the very name of generality, maybe we should not find a common
language, as long as we are at least able to read each other's code and
pseudo-code.

I'd rather propose the following problem:

1. Let's write and compare our best algorithms.

2. I'll still use Prolog, and some pseudo-code if needed. You are
welcome to keep using Haskell, at least I can read and run Haskell.

3. Let's instrument our code to count the arithmetic operations instead
of an execution time (something as simple as incrementing a counter
after each operation). Then we can benchmark our implementations over a
common set of test cases (which we can decide later, something simple)
by comparing the operation counts after each invocation.

4. The interface to implement is `ilog2(in x, out y, out c), for x > 0`
where x is the input, y is the result, and c is the total count of
*elementary* arithmetic operations counted (e.g. `m = (h-l)/2` is
equivalent to `t = h-l, m = t/2` so counts as 2 operations). Of course,
the code being shared, we cannot easily cheat on `c`... :)

How does it sound?

Julio

Ben Bacarisse

unread,
Feb 11, 2017, 4:13:20 PM2/11/17
to
"Julio P. Di Egidio" <ju...@diegidio.name> writes:

> On 11/02/2017 18:22, Ben Bacarisse wrote:
>> "Julio P. Di Egidio" <ju...@diegidio.name> writes:
> <snip>
>>> What about JavaScript? I am not inclined in trying Haskell because
>>> one thing is reading it, another is writing efficient code.
>>
>> JS would be fine though rather messy since JS does not have arbitrary
>> precision arithmetic built in. Which library should I use? I suggest
>> node.js's bignumber because it's easy to install, but it is rather slow.
>
> You are right, JS won't do: bignums are too slow for anything but toy
> applications, I have even tried in the past...

The fact that it's slow won't matter much but you don't want to go use
it, OK.

>> Mind you, that may be an advantage since it will highlight the algorithm
>> rather than the arithmetic.
>
> In the very name of generality, maybe we should not find a common
> language, as long as we are at least able to read each other's code
> and pseudo-code.
>
> I'd rather propose the following problem:
>
> 1. Let's write and compare our best algorithms.
>
> 2. I'll still use Prolog, and some pseudo-code if needed. You are
> welcome to keep using Haskell, at least I can read and run Haskell.
>
> 3. Let's instrument our code to count the arithmetic operations
> instead of an execution time (something as simple as incrementing a
> counter after each operation). Then we can benchmark our
> implementations over a common set of test cases (which we can decide
> later, something simple) by comparing the operation counts after each
> invocation.
>
> 4. The interface to implement is `ilog2(in x, out y, out c), for x > 0`
> where x is the input, y is the result, and c is the total count of
> *elementary* arithmetic operations counted (e.g. `m = (h-l)/2` is
> equivalent to `t = h-l, m = t/2` so counts as 2 operations). Of
> course, the code being shared, we cannot easily cheat on `c`... :)
>
> How does it sound?

That would work, but I think a common language would be better. Can't
we find one? Python, Perl, C, C++, Go, ML, Scheme, Common Lisp... that
all have bindings to a big-num library.

Anyway, to speed thing along, here is the fastest Haskell version I have
with counted operations. It returns a pair: the log and the count.

import Data.Bits
import System.Environment

-- bit n (from Data.Bits) is 1 << n in other languages

ilog2 :: Integer -> (Int, Int)
ilog2 x = ilog_search (bounds (0, 1, 0))
where
bounds (lo, hi, ops)
| bit hi >= x = (lo, hi, ops+1)
| otherwise = bounds (hi, hi*2, ops+2)
ilog_search (lo, hi, ops)
| hi <= lo = (lo, ops)
| bit mid >= x = ilog_search (lo, mid, ops+3)
| otherwise = ilog_search (mid+1, hi, ops+4)
where mid = (lo + hi) `div` 2

main = do
[n'] <- getArgs
print $ ilog2 $ read n'

--
Ben.

Julio P. Di Egidio

unread,
Feb 12, 2017, 5:39:38 AM2/12/17
to
On 11/02/2017 22:13, Ben Bacarisse wrote:
> "Julio P. Di Egidio" <ju...@diegidio.name> writes:
>> On 11/02/2017 18:22, Ben Bacarisse wrote:
>>> "Julio P. Di Egidio" <ju...@diegidio.name> writes:
>> <snip>
>>>> What about JavaScript? I am not inclined in trying Haskell because
>>>> one thing is reading it, another is writing efficient code.
>>>
>>> JS would be fine though rather messy since JS does not have arbitrary
>>> precision arithmetic built in. Which library should I use? I suggest
>>> node.js's bignumber because it's easy to install, but it is rather slow.
>>
>> You are right, JS won't do: bignums are too slow for anything but toy
>> applications, I have even tried in the past...
>
> The fact that it's slow won't matter much but you don't want to go use
> it, OK.

I haven't tried in node.js, but it's just way too slow in the browser
even for our purposes. That said, I explain and explain, you not only
keep bitching, you also keep acting as if I am the unreasonable one.

> That would work, but I think a common language would be better.

I have explained why the implementation language is irrelevant.

Julio

Julio P. Di Egidio

unread,
Feb 12, 2017, 5:48:07 AM2/12/17
to
On 11/02/2017 19:59, Julio P. Di Egidio wrote:

> 1. Let's write and compare our best algorithms.

I have updated the code on GitHub Gist:
<https://gist.github.com/jp-diegidio/eb2d6c1ae432c0e3c75d9f1e5bba602c>

Two algorithms are implemented:

- ilog2_n/3 : Naive approach.
- ilog2_i/3 : Faster approach that uses integer arithmetic only.

Test results:
-------------

ilog2_i(1267650600228229401496826662165) = 100 (C=32)
ilog2_i(2535301200456458802993529867541) = 101 (C=32)
ilog2_i(5070602400912917605986936278293) = 102 (C=34)
ilog2_i(10141204801825835211973749099797) = 103 (C=30)
ilog2_i(20282409603651670423947374742805) = 104 (C=32)
ilog2_i(40564819207303340847894626028821) = 105 (C=32)
ilog2_i(81129638414606681695789128600853) = 106 (C=34)
ilog2_i(162259276829213363391578133744917) = 107 (C=32)
ilog2_i(324518553658426726783156144033045) = 108 (C=34)
ilog2_i(649037107316853453566312164609301) = 109 (C=34)
ilog2_i(1298074214633706907132624205761813) = 110 (C=36)

% test__case(out x) =>
% for j in [0,100]
% x0 := 2^j
% for i in [0,100]
% x := x0 + i
ilog2_test(ilog2_n) C=1053613
ilog2_test(ilog2_i) C=264235

Julio

Julio P. Di Egidio

unread,
Feb 12, 2017, 6:14:30 AM2/12/17
to
On 11/02/2017 22:13, Ben Bacarisse wrote:

> here is the fastest Haskell version I have
> with counted operations.

OK. Can you now run your code and see what you get for the test cases I
have posted in the meantime?

*Note*: I am only counting the arithmetic operations that create new
numbers, namely I am neither counting comparisons nor I am counting
assignments of return values (every bignums implementation would pass
numbers around by reference). You might have to adjust your code!

Julio

Ben Bacarisse

unread,
Feb 12, 2017, 6:34:40 PM2/12/17
to
"Julio P. Di Egidio" <ju...@diegidio.name> writes:

> On 11/02/2017 22:13, Ben Bacarisse wrote:
>
>> here is the fastest Haskell version I have
>> with counted operations.
>
> OK. Can you now run your code and see what you get for the test cases
> I have posted in the meantime?

I'll post in reply to your numbers.

<snip>
--
Ben.

Ben Bacarisse

unread,
Feb 12, 2017, 6:50:18 PM2/12/17
to
"Julio P. Di Egidio" <ju...@diegidio.name> writes:

> On 11/02/2017 19:59, Julio P. Di Egidio wrote:
>
>> 1. Let's write and compare our best algorithms.
>
> I have updated the code on GitHub Gist:
> <https://gist.github.com/jp-diegidio/eb2d6c1ae432c0e3c75d9f1e5bba602c>
>
> Two algorithms are implemented:
>
> - ilog2_n/3 : Naive approach.
> - ilog2_i/3 : Faster approach that uses integer arithmetic only.

ilog2_i looks like a new algorithm based on the way it performs compared
to the previous code. I am guessing it's O(log log n) but I can't
unravel the Prolog to get to pseudo-code.

I noticed that I was calculating a different log to you (ceiling rather
than floor) so I've had to change my code to:

import Data.Bits

ilog2 :: Integer -> (Int, Int)
ilog2 x = ilog_search (bounds (0, 1, 0))
where
bounds (lo, hi, ops)
| bit hi >= x = (lo, hi, ops+1)
| otherwise = bounds (hi, hi*2, ops+2)
ilog_search (lo, hi, ops)
| hi <= lo = if bit lo > x then (lo-1, ops+2) else (lo, ops+1)
| bit_mid == x = (mid, ops+3)
| bit_mid > x = ilog_search (lo, mid-1, ops+4)
| otherwise = ilog_search (mid+1, hi, ops+4)
where mid = (lo + hi) `div` 2
bit_mid = bit mid


> Test results:
> -------------
>
> ilog2_i(1267650600228229401496826662165) = 100 (C=32)
> ilog2_i(2535301200456458802993529867541) = 101 (C=32)
> ilog2_i(5070602400912917605986936278293) = 102 (C=34)
> ilog2_i(10141204801825835211973749099797) = 103 (C=30)
> ilog2_i(20282409603651670423947374742805) = 104 (C=32)
> ilog2_i(40564819207303340847894626028821) = 105 (C=32)
> ilog2_i(81129638414606681695789128600853) = 106 (C=34)
> ilog2_i(162259276829213363391578133744917) = 107 (C=32)
> ilog2_i(324518553658426726783156144033045) = 108 (C=34)
> ilog2_i(649037107316853453566312164609301) = 109 (C=34)
> ilog2_i(1298074214633706907132624205761813) = 110 (C=36)

I get:

ilog2(1267650600228229401496826662165) = 100 (c=37)
ilog2(2535301200456458802993529867541) = 101 (c=36)
ilog2(5070602400912917605986936278293) = 102 (c=37)
ilog2(10141204801825835211973749099797) = 103 (c=36)
ilog2(20282409603651670423947374742805) = 104 (c=37)
ilog2(40564819207303340847894626028821) = 105 (c=36)
ilog2(81129638414606681695789128600853) = 106 (c=37)
ilog2(162259276829213363391578133744917) = 107 (c=36)
ilog2(324518553658426726783156144033045) = 108 (c=37)
ilog2(649037107316853453566312164609301) = 109 (c=36)
ilog2(1298074214633706907132624205761813) = 110 (c=37)

> % test__case(out x) =>
> % for j in [0,100]
> % x0 := 2^j
> % for i in [0,100]
> % x := x0 + i
> ilog2_test(ilog2_n) C=1053613
> ilog2_test(ilog2_i) C=264235

ilog2_test C=309934

This test does not show much because it just multiplies any absolute
difference in the counts.

The counts can really only help to show the order of growth, and I think
both implementations are O(log log n). The previous code you posted
seemed to be O(log n) as you had stated.

I'd like to know what the algorithm is (and if it's different, what the
old one was!). Using a shared language would help with that.

--
Ben.

Ben Bacarisse

unread,
Feb 12, 2017, 6:55:43 PM2/12/17
to
For you, perhaps, but not for me! For one thing, I'd like to get some
knowledge out of this exchange, and that would mean understanding your
algorithm(s).

If JS is a non-starter, there must surely be some other common language,
no?

--
Ben.

Julio P. Di Egidio

unread,
Feb 12, 2017, 8:36:13 PM2/12/17
to
On 13/02/2017 00:55, Ben Bacarisse wrote:
> "Julio P. Di Egidio" <ju...@diegidio.name> writes:
>> On 11/02/2017 22:13, Ben Bacarisse wrote:
>>> "Julio P. Di Egidio" <ju...@diegidio.name> writes:
>>>> On 11/02/2017 18:22, Ben Bacarisse wrote:
>>>>> "Julio P. Di Egidio" <ju...@diegidio.name> writes:
>>>> <snip>
>>>>>> What about JavaScript? I am not inclined in trying Haskell because
>>>>>> one thing is reading it, another is writing efficient code.
>>>>>
>>>>> JS would be fine though rather messy since JS does not have arbitrary
>>>>> precision arithmetic built in. Which library should I use? I suggest
>>>>> node.js's bignumber because it's easy to install, but it is rather slow.
>>>>
>>>> You are right, JS won't do: bignums are too slow for anything but toy
>>>> applications, I have even tried in the past...
>>>
>>> The fact that it's slow won't matter much but you don't want to go use
>>> it, OK.
>>
>> I haven't tried in node.js, but it's just way too slow in the browser
>> even for our purposes. That said, I explain and explain, you not only
>> keep bitching, you also keep acting as if I am the unreasonable one.
>>
>>> That would work, but I think a common language would be better.
>>
>> I have explained why the implementation language is irrelevant.
>
> For you, perhaps, but not for me! For one thing, I'd like to get some
> knowledge out of this exchange, and that would mean understanding your
> algorithm(s).

Then open my link and read: besides having explained it thoroughly
up-thread, I have now written full pseudo-code along side the Prolog
code, and I have done it JUST FOR YOU: if you cannot even read that, I
wouldn't know what else to try other than telepathy.

> If JS is a non-starter, there must surely be some other common language,
> no?

I am having doubts that the plain English is working here: which is
another way to state that your request (not to mention, your blind
insistence) is ultimately senseless.

That said, I will review your code later: I think you must be
miscounting operations, I would have expected your code to be way slower
than that...

Julio

Ben Bacarisse

unread,
Feb 13, 2017, 7:33:57 AM2/13/17
to
I can read it (and indeed I did read it) but the names (like do and x0,
x1, y and so on) didn't help to make it clear. However, I decided to
persevere and was able to translate it into Haskell:

ilog2_i :: Integer -> (Int, Int)
ilog2_i x = down $ up (x, 1, 0)
where up (x, m, ops) =
let x' = shiftR x m
in if x' /= 0
then up (x', 2*m, ops+2)
else (m, x, ops+1)

down (1, _, ops) = (0, ops)
down (m, x, ops) = log_do (m `div` 2, x, m - 1, ops+2)

log_do (m, x, y, ops) =
let x' = shiftR x m
in if x' /= 0
then log_do(m, x', y + m, ops+2)
else if m /= 1
then log_do(m `div` 2, x, y, ops+2)
else (y, ops+1)

Here I've used x and x-prime (x') where you would have X0 and X1 because
that's more idiomatic in Haskell. I'm pretty sure it's right because I
get exactly the same counts as you so and it runs just about as fast as
I'd expect (i.e. pretty much the same as my version).

(A detail: you have a typo in the first line of the description of _dn.
y0 should be m0 I think.)

I still don't understand the algorithm, because I can't yet see the
overall plan. My code keeps squaring an upper bound for the logarithm
until it is known to be >= the target. Then a conventional binary
search is used to find the log between the upper bound and its square
root (i.e. between the last estimates for the upper bound).

Do you have a similar description for your algorithm? Maybe it is not
conducive to that sort of high-level description.

>> If JS is a non-starter, there must surely be some other common language,
>> no?
>
> I am having doubts that the plain English is working here: which is
> another way to state that your request (not to mention, your blind
> insistence) is ultimately senseless.

What am I blindly insisting on? I seem to be able to annoy you with
everything is do, but I'm just trying to understand and offer up some
ideas.

> That said, I will review your code later: I think you must be
> miscounting operations, I would have expected your code to be way
> slower than that...

It's possible. For the earlier version I wrote all the arithmetic as
separate functions ("addOne", "timesTwo" and so on) and then used
Haskell's profiler to verify that the counts I report are those actually
done, but it messes up the code and I've not done that since my initial
tests.

Both your ilog2_i and my ilog2 seem to be O(log log n).

One thing the counts might obscure is the number of "short" and "long"
arithmetic operations. I think both your code and mine are built from
arbitrary precision binary shifts and some "housekeeping" arithmetic on
small counters. Doing 6 operations where one is a million digit shift
and five are native x+1 operations is going to cheaper than doing one
native x+1 and fine million digits shifts.

--
Ben.

Julio P. Di Egidio

unread,
Feb 17, 2017, 11:05:25 AM2/17/17
to
On 13/02/2017 13:33, Ben Bacarisse wrote:
<snipped>

> I can read it (and indeed I did read it) but the names (like do and x0,
> x1, y and so on) didn't help to make it clear.

x is the input, y is the result, m is the size of the mask (I do not
really mask anything, I use m as the amount by which to shift the
input). Same names for the running values, not just the top-level.

> (A detail: you have a typo in the first line of the description of _dn.
> y0 should be m0 I think.)

Thanks, I have corrected few of those in the meantime.

> I still don't understand the algorithm, because I can't yet see the
> overall plan. My code keeps squaring an upper bound for the logarithm
> until it is known to be >= the target. Then a conventional binary
> search is used to find the log between the upper bound and its square
> root (i.e. between the last estimates for the upper bound).
>
> Do you have a similar description for your algorithm? Maybe it is not
> conducive to that sort of high-level description.

Of course it is, every procedure is amenable to an informal description.

Mine is just a bit longer than that, and it is in two phases, too, which
I have called up and down (dn), but the analogy with your approach I
suppose ends there:

- In the up phase, I start from a mask of size 1 (1 bit) and, at each
iteration, before increasing the size of the mask (but I only double it,
which is one shift, BTW), I also shift to the right the input by the
current size of the mask: this for every iteration, until the input so
reduced is not zero. This phase ends giving back the last reduced input
(which is the remainder after all possible shifts), and the size that
the mask has reached (which at exit is strictly bigger than the size of
the reduced input).

- Then starts the down phase, first with the predicate _dn which only
handles the special case that the size of the mask (as resulting from
the up phase) is 1, in which case the input had to be 1 and the result 0
can be returned; otherwise, _dn first decreases (namely, halves) the
mask (so that it is again lte the size of the input), then calls _do to
perform the actual down loop, passing to it also an accumulator that is
initially set to the (decreased) size of the mask: this will accumulate
the result.

- In the _do (down) loop, the (residual) input is repeatedly shifted
to the right by the current size of the mask, and at every shit the
accumulator itself is increased by the size of the mask, until the
reduced input is not zero; otherwise, if the size of the mask is not yet
equal to 1, decreases (halves) the size of the mask and restarts the
loop; otherwise, the result is set equal to the accumulator.

So, that is the gist of it, though I am pretty sure the structure could
be simplified a little bit, namely, I should rethink the control flow as
I might have messed up things a little bit there. But, more than that,
I should investigate different working parameters, such as squaring
instead of doubling the mask during the up phase, and other possible
variations that certainly would have an impact on the efficiency. -- I
say I should because this is all of course academic, nobody in the real
world would want an integer binary logarithm implemented in user code...

> Both your ilog2_i and my ilog2 seem to be O(log log n).

That is plausible, it is in both cases a double binary search. But I
guess once we have stable code, we can get to the bottom of a precise
analysis, at least in terms of big-O, can't we?

> One thing the counts might obscure is the number of "short" and "long"
> arithmetic operations. I think both your code and mine are built from
> arbitrary precision binary shifts and some "housekeeping" arithmetic on
> small counters. Doing 6 operations where one is a million digit shift
> and five are native x+1 operations is going to cheaper than doing one
> native x+1 and fine million digits shifts.

Yes, that's very right, but... what to do?

On one hand, we could stop counting the "small" operations as well, to
be left with operations on big numbers only, and that should make the
counting at least fairer. Still not completely fair, as for example
multiplication costs way more than addition that costs more than a
shift, etc. We could restrict ourselves to shifts, addition and
subtraction in the present case, but such restrictions would not work in
general (keeping in mind the academic interest, which is indeed about
how to actually analyse algorithms, in general).

Unless we do come to the bottom of what it even means to do a precise
analysis/counting (i.e. where to cut the line between the code and the
underlying platform, and how to cut that line), the only option I see
left is that we just implement all variants in our own respective
systems (you have done it already, I shall implement your idea in
Prolog), then we can do the usual direct comparisons.

BTW, sorry for the delay, I have had a busy week. I should be able to
put my hands on this again across the week-end.

Julio

Julio P. Di Egidio

unread,
Feb 17, 2017, 11:34:14 AM2/17/17
to
On 17/02/2017 17:05, Julio P. Di Egidio wrote:

> - In the _do (down) loop, the (residual) input is repeatedly shifted
> to the right by the current size of the mask, and at every shit the

LOL, the spell checker didn't catch that... My apologies.

Julio

Ben Bacarisse

unread,
Feb 17, 2017, 11:43:25 AM2/17/17
to
"Julio P. Di Egidio" <ju...@diegidio.name> writes:

That appears to a description of the code in English. I was able to
write it in Haskell so I know what it does operationally, but I don't
understand it in any useful way.

> So, that is the gist of it, though I am pretty sure the structure
> could be simplified a little bit, namely, I should rethink the control
> flow as I might have messed up things a little bit there. But, more
> than that, I should investigate different working parameters, such as
> squaring instead of doubling the mask during the up phase, and other
> possible variations that certainly would have an impact on the
> efficiency. -- I say I should because this is all of course
> academic, nobody in the real world would want an integer binary
> logarithm implemented in user code...
>
>> Both your ilog2_i and my ilog2 seem to be O(log log n).
>
> That is plausible, it is in both cases a double binary search.

I know mine does a binary search (I would not call it a double binary
search) but you obviously understand your code better than I do because
I was not be sure what it's big-O was.

> But I
> guess once we have stable code, we can get to the bottom of a precise
> analysis, at least in terms of big-O, can't we?

I'm pretty confident of the O(log log n) result for the simple binary
search method because it does O(log n) steps in a range whose size is
O(log n). The analysis is going to be more fiddly without a high-level
view of the algorithm.

>> One thing the counts might obscure is the number of "short" and "long"
>> arithmetic operations. I think both your code and mine are built from
>> arbitrary precision binary shifts and some "housekeeping" arithmetic on
>> small counters. Doing 6 operations where one is a million digit shift
>> and five are native x+1 operations is going to cheaper than doing one
>> native x+1 and fine million digits shifts.
>
> Yes, that's very right, but... what to do?

I'm OK with this now because I have comparable code in one language.
Both methods perform very much the same as far as I can tell.

> On one hand, we could stop counting the "small" operations as well, to
> be left with operations on big numbers only, and that should make the
> counting at least fairer. Still not completely fair, as for example
> multiplication costs way more than addition that costs more than a
> shift, etc.

Neither code has any multiplications or divisions in that the 2*x and
x//2 can be optimised to shifts (or an addition, in the case of the
doubling). What will matter is big- vs small-number operations. In my
version, the only big-num operation (other than compares) is to set a
bit. The Haskell "bit n" means the same as "1 << n" in other languages,
but it might be more efficient. As it happens, I can't tell the
difference when I replace it with "shiftL 1 n" so I suspect it's not
optimised in any way.

So one could refine the counting by simply splitting out small number
operations from arbitrary precision number operations (in Haskell, the
things that are of type Int vs the things of type Integer).

> We could restrict ourselves to shifts, addition and
> subtraction in the present case, but such restrictions would not work
> in general (keeping in mind the academic interest, which is indeed
> about how to actually analyse algorithms, in general).
>
> Unless we do come to the bottom of what it even means to do a precise
> analysis/counting (i.e. where to cut the line between the code and the
> underlying platform, and how to cut that line), the only option I see
> left is that we just implement all variants in our own respective
> systems (you have done it already, I shall implement your idea in
> Prolog), then we can do the usual direct comparisons.

Sure, though it's not clear that there would be much to gain from doing
another implementation as it's very likely that the relative performance
will be similar. It's possible the Haskell or Prolog penalises (or has
some super fast version of) some operation used by one algorithm and not
the other but since the two algorithms are so similar, that seems
unlikely.

> BTW, sorry for the delay, I have had a busy week. I should be able to
> put my hands on this again across the week-end.

That's OK. I was able to get a common implementation and do some
timings.

--
Ben.

Julio P. Di Egidio

unread,
Feb 17, 2017, 12:22:55 PM2/17/17
to
On 17/02/2017 17:43, Ben Bacarisse wrote:

> That appears to a description of the code in English. I was able to
> write it in Haskell so I know what it does operationally, but I don't
> understand it in any useful way.

It's really your problem, as the whole thing is elementary: while, on
the issue of what constitutes a proper algorithmic analysis, you are
simply dodging... Now, about the informal description I have given you,
maybe I have put too much detail (I just keep thinking the same code in
every language), but, really, how would you rewrite or give a simple
name to "increase the size of the mask while also reducing the input",
which is essentially my "up" phase? I just wouldn't know how to
simplify that further.

> That's OK. I was able to get a common implementation and do some
> timings.

I think I'll do it anyway, I do not trust your counts... ;)

Julio

Ben Bacarisse

unread,
Feb 17, 2017, 4:29:26 PM2/17/17
to
"Julio P. Di Egidio" <ju...@diegidio.name> writes:

> On 17/02/2017 17:43, Ben Bacarisse wrote:
>
>> That appears to a description of the code in English. I was able to
>> write it in Haskell so I know what it does operationally, but I don't
>> understand it in any useful way.
>
> It's really your problem,

Yes of course it's my problem. But if you want people like me to
understand your algorithm, then it becomes yours too. I'm not going to
be upset if you say you only care if people cleverer than I understand
it, or people prepared to put in more work than I, or whatever it is you
think is preventing my grasping it.

> as the whole thing is elementary: while, on
> the issue of what constitutes a proper algorithmic analysis, you are
> simply dodging...

Absolutely agreed -- I am dodging it. You, presumably, are not dodging
it, so what is your proper algorithmic analysis of the code?

<snip>
--
Ben.

Julio P. Di Egidio

unread,
Feb 19, 2017, 5:53:04 AM2/19/17
to
On 13/02/2017 00:50, Ben Bacarisse wrote:
> "Julio P. Di Egidio" <ju...@diegidio.name> writes:
>
>> On 11/02/2017 19:59, Julio P. Di Egidio wrote:
>>
>>> 1. Let's write and compare our best algorithms.
>>
>> I have updated the code on GitHub Gist:
>> <https://gist.github.com/jp-diegidio/eb2d6c1ae432c0e3c75d9f1e5bba602c>

I have uploaded a new version, same link.

>> Two algorithms are implemented:
>>
>> - ilog2_n/3 : Naive approach.
>> - ilog2_i/3 : Faster approach that uses integer arithmetic only.

Now there is also an ilog2_b, which is my implementation of your algorithm.

>> Test results:
>> -------------
>>
>> ilog2_i(1267650600228229401496826662165) = 100 (C=32)
>> ilog2_i(2535301200456458802993529867541) = 101 (C=32)
>> ilog2_i(5070602400912917605986936278293) = 102 (C=34)
>> ilog2_i(10141204801825835211973749099797) = 103 (C=30)
>> ilog2_i(20282409603651670423947374742805) = 104 (C=32)
>> ilog2_i(40564819207303340847894626028821) = 105 (C=32)
>> ilog2_i(81129638414606681695789128600853) = 106 (C=34)
>> ilog2_i(162259276829213363391578133744917) = 107 (C=32)
>> ilog2_i(324518553658426726783156144033045) = 108 (C=34)
>> ilog2_i(649037107316853453566312164609301) = 109 (C=34)
>> ilog2_i(1298074214633706907132624205761813) = 110 (C=36)
>
> ilog2(1267650600228229401496826662165) = 100 (c=37)
> ilog2(2535301200456458802993529867541) = 101 (c=36)
> ilog2(5070602400912917605986936278293) = 102 (c=37)
> ilog2(10141204801825835211973749099797) = 103 (c=36)
> ilog2(20282409603651670423947374742805) = 104 (c=37)
> ilog2(40564819207303340847894626028821) = 105 (c=36)
> ilog2(81129638414606681695789128600853) = 106 (c=37)
> ilog2(162259276829213363391578133744917) = 107 (c=36)
> ilog2(324518553658426726783156144033045) = 108 (c=37)
> ilog2(649037107316853453566312164609301) = 109 (c=36)
> ilog2(1298074214633706907132624205761813) = 110 (c=37)

I get something different for yours:

ilog2_b(1267650600228229401496826662165) = 100 (c=36)
ilog2_b(2535301200456458802993529867541) = 101 (c=37)
ilog2_b(5070602400912917605986936278293) = 102 (c=37)
ilog2_b(10141204801825835211973749099797) = 103 (c=38)
ilog2_b(20282409603651670423947374742805) = 104 (c=36)
ilog2_b(40564819207303340847894626028821) = 105 (c=37)
ilog2_b(81129638414606681695789128600853) = 106 (c=37)
ilog2_b(162259276829213363391578133744917) = 107 (c=38)
ilog2_b(324518553658426726783156144033045) = 108 (c=37)
ilog2_b(649037107316853453566312164609301) = 109 (c=38)
ilog2_b(1298074214633706907132624205761813) = 110 (c=38)

Yours did look suspicious, your counts are pretty much constant.

> ilog2_test C=309934

Again, I get something different for yours: C = 359907.

Here are few more direct comparisons, mainly look at the CPU time, the
number of inferences is not particularly indicative:

% _case(out x) =>
% for j in [0, 100]
% for i in [0, 100]
% x := 2^j + i
?- time(ilog2_test(ilog2_n,C)).
?- time(ilog2_test(ilog2_b,C)).
?- time(ilog2_test(ilog2_i,C)).
C = 1053613.
C = 359907.
C = 264235.
% 2,199,366 inferences, 0.219 CPU in 0.233 seconds (94% CPU, 10054245 Lips)
% 902,684 inferences, 0.094 CPU in 0.092 seconds (102% CPU, 9628629 Lips)
% 690,887 inferences, 0.063 CPU in 0.065 seconds (96% CPU, 11054192 Lips)

?- time(forall((N is 2^2^4, between(1,N,X)), ilog2_n(X,_,_))).
?- time(forall((N is 2^2^4, between(1,N,X)), ilog2_b(X,_,_))).
?- time(forall((N is 2^2^4, between(1,N,X)), ilog2_i(X,_,_))).
% 4,390,901 inferences, 0.313 CPU in 0.299 seconds (105% CPU, 14050883 Lips)
% 3,838,521 inferences, 0.281 CPU in 0.282 seconds (100% CPU, 13648075 Lips)
% 3,132,075 inferences, 0.234 CPU in 0.240 seconds (98% CPU, 13363520 Lips)

?- time(forall((M is 2^2^6,N is M+2^2^4, between(M,N,X)), ilog2_n(X,_,_))).
?- time(forall((M is 2^2^6,N is M+2^2^4, between(M,N,X)), ilog2_b(X,_,_))).
?- time(forall((M is 2^2^6,N is M+2^2^4, between(M,N,X)), ilog2_i(X,_,_))).
% 17,498,377 inferences, 2.125 CPU in 2.125 seconds (100% CPU, 8234530 Lips)
% 7,274,583 inferences, 0.750 CPU in 0.745 seconds (101% CPU, 9699444 Lips)
% 4,522,056 inferences, 0.391 CPU in 0.390 seconds (100% CPU, 11576463 Lips)

As you can see by looking at the CPU time, ilog2_i is significantly
faster and scales better. And I have not implemented it again (it was
in my original version), but an ilog2_f would be even faster than that.

> The counts can really only help to show the order of growth, and I think
> both implementations are O(log log n). The previous code you posted
> seemed to be O(log n) as you had stated.

Nope, my original one was of the same order, in fact even faster as I
was tuning my "initial k" by hand.

Julio

Julio P. Di Egidio

unread,
Feb 19, 2017, 6:06:20 AM2/19/17
to
On 17/02/2017 22:29, Ben Bacarisse wrote:
> "Julio P. Di Egidio" <ju...@diegidio.name> writes:
<snip>
>> It's really your problem,
>
> Yes of course it's my problem. But if you want people like me to
> understand your algorithm, then it becomes yours too.

People like you should be re-educated: you are a systematic liar, that
is the immediate problem you have, the rest is a consequence.

>> while, on
>> the issue of what constitutes a proper algorithmic analysis, you are
>> simply dodging...
>
> Absolutely agreed -- I am dodging it. You, presumably, are not dodging
> it, so what is your proper algorithmic analysis of the code?

Hmm, I thought you were the expert. I will find out and let you know...

EOD.

Julio

Ben Bacarisse

unread,
Feb 19, 2017, 1:15:35 PM2/19/17
to
I re-checked the Haskell counts by replacing all the operations with
functions and running the profiler. The counts check out.

That they are pretty must constant over that range is not surprising for
an O(log log n) algorithm. log log 1267650600228229401496826662165 is
the same as log log 1298074214633706907132624205761813. I'd say all the
counts are "pretty much constant". The variation of one or two is
just noise.

That your counts are different to mine is not surprising in that the two
algorithms are slightly different and will have slightly different
"noise" patterns. (For example, I test for equality and hence always
either add or subtract one to the middle value.)

<snip>
Thanks for the Prolog version. I'll look in more detail later.

--
Ben.

Ben Bacarisse

unread,
Feb 19, 2017, 1:17:10 PM2/19/17
to
"Julio P. Di Egidio" <ju...@diegidio.name> writes:

> On 17/02/2017 22:29, Ben Bacarisse wrote:
>> "Julio P. Di Egidio" <ju...@diegidio.name> writes:
> <snip>
>>> It's really your problem,
>>
>> Yes of course it's my problem. But if you want people like me to
>> understand your algorithm, then it becomes yours too.
>
> People like you should be re-educated: you are a systematic liar,

I replied to your other post before reading this. My mistake.

<snip>
--
Ben.
0 new messages