Bits of integer

[Albert Chan]

This was a PM that I thought wothy of starting a thread: Number of Bits in a Integer

Albert Chan Wrote:Unless you have correct implementation of log2 (or binary logb), do this instead:

Quote:20 DEF FNB(N)=IP(LN(N+.5)/L2)+1 ! BITS OF INTEGER N

It seems it is harder than it look, even with "perfect" log2

>20 def FNB(N) = IP(LOG2(N))+1
>X = 2^35
>X, FNB(X-1), FNB(X), FNB(X+1)
 34359738368      36      36      36

For binary machine, code is trivial, bit_length(x) = logb(x) + 1

For decimal machine, we may be hit with rounding errors. Any ideas ?

[Werner]

Hi, Albert!
Yes, the method I used in the VA12c challenge will return the correct answer in linear time for all arguments up to 2^35-1 (on a 42S) or to 2^WSIZE-1 on Free42 (WSIZE<45):

Code:

 LBL D
 BINM @ set Binary mode, activates the BASE menu
 CLA
 ARCL ST X @ Alpha reg now contains X in binary
 CLX @ To replace X with the new value - my routine depended on that
 ALENG @ nr of bits
 EXITALL @ exit the BASE menu
 RTN

Sadly, that cannot be used on the 41. Don't know about the 71. On the 48, something similar may be done:

Code:

\<< BIN R\->B \->STR SIZE 3. - \>>

Cheers, Werner

[J-F Garnier]

(12-05-2022 12:12 PM)Albert Chan Wrote:  >20 def FNB(N) = IP(LOG2(N))+1
>X = 2^35
>X, FNB(X-1), FNB(X), FNB(X+1)
 34359738368      36      36      36

For decimal machine, we may be hit with rounding errors. Any ideas ?

The same occurs with LOG10 when trying to get the exponent of a number:
X=10^11
LOG10(X);LOG10(X-1)
11 11
Fortunately, in this case we have the EXPONENT() function.

(12-05-2022 12:41 PM)Werner Wrote:  Yes, the method I used in the VA12c challenge will return the correct answer in linear time for all arguments up to 2^35-1 (on a 42S) or to 2^WSIZE-1 on Free42 (WSIZE<45):

Code:

 LBL D
 BINM @ set Binary mode, activates the BASE menu
 CLA
 ARCL ST X @ Alpha reg now contains X in binary
 CLX @ To replace X with the new value - my routine depended on that
 ALENG @ nr of bits
 EXITALL @ exit the BASE menu
 RTN

Sadly, that cannot be used on the 41. Don't know about the 71. On the 48, something similar may be done:

Code:

\<< BIN R\->B \->STR SIZE 3. - \>>

On the 71B or 75C w/ Math ROM, you can do:
DEF FNB(X)=LEN(BSTR$(X,2))
but I find it quite inelegant to resort to string conversion.

Edit: removed my other proposal. made a mistake in my tests :-(

J-F

[Albert Chan]

(12-05-2022 02:55 PM)J-F Garnier Wrote:  On the 71B or 75C w/ Math ROM, you can do:
DEF FNB(X)=LEN(BSTR$(X,2))
but I find it quite inelegant to resort to string conversion.

>10 DEF FNB(X)=LEN(BSTR$(X,2))
>X = 2^39
>X
549755813888
>FNB(X-1), FNB(X)
 39      40

>Y$ = BSTR$(X,2)
ERR:String Ovfl

FNB(x) work, but intermediate values cannot be stored.
Will the "overflowed" intermediate wreck havoc ?

[Valentin Albillo]

(12-05-2022 04:13 PM)Albert Chan Wrote:  >10 DEF FNB(X)=LEN(BSTR$(X,2))
>X = 2^39
>X
549755813888
>FNB(X-1), FNB(X)
 39      40

>Y$ = BSTR$(X,2)
ERR:String Ovfl

FNB(x) work, but intermediate values cannot be stored.
Will the "overflowed" intermediate wreck havoc ?

You got Overflow because by default undimensioned string variables are limited to 32 characters.

Simply execute DIM Y$[40] and the assignment will work. BSTR$ will never produce a string longer than 40 chars because its argument must be < 999,999,999,999.5

V.

[robve]

(12-05-2022 12:12 PM)Albert Chan Wrote:  This was a PM that I thought wothy of starting a thread: Number of Bits in a Integer

Albert Chan Wrote:Unless you have correct implementation of log2 (or binary logb), do this instead:

It seems it is harder than it look, even with "perfect" log2

>20 def FNB(N) = IP(LOG2(N))+1
>X = 2^35
>X, FNB(X-1), FNB(X), FNB(X+1)
 34359738368      36      36      36

For binary machine, code is trivial, bit_length(x) = logb(x) + 1

For decimal machine, we may be hit with rounding errors. Any ideas ?

For C (and lua etc) I came up with a simple binary search to find the position of the leading bit in a 32 bit integer in 5 search steps:

Code:

  unsigned int n, k, v = <number>;
  k = -!!(v >> 16) & 16;
  n = k;
  v >>= k;
  k = -!!(v >> 8) & 8;
  n += k;
  v >>= k;
  k = -!!(v >> 4) & 4;
  n += k;
  v >>= k;
  k = -!!(v >> 2) & 2;
  n += k;
  v >>= k;
  k = -!!(v >> 1) & 1;
  n += k;
  v >>= k;

It's fast in C, since the compiled code operates with a few registers using bit-wise operations. But on the HP-71B we will need to shift with integer division and integer multiplication. Forth has shift words.

Another trick is to normalize the integer to a clean form before taking the log2, similar to this:

Code:

x |= (x >> 1);
x |= (x >> 2);
x |= (x >> 4);
x |= (x >> 8);
x |= (x >> 16);
n = log2(x & ~(x >> 1)));

As a side note, the LZCNT CPU instruction counts leading zeros and can be used in e.g. C/C++ with __builtin_clz(). The POPCNT CPU instruction counts the number of nonzero bits. That neither work with a HP-71B of course.

- Rob

[Albert Chan]

Hi, robve

We can do binary searches without the costly integer division.

lua> bisect = require'primepi'.bisect
lua> pow2 = {}
lua> for i=1,52 do pow2[i] = 2^i end
lua> function bits(n) return bisect(pow2, n) end

lua> for n = 1, 9 do print(n, bits(n)) end
1      1
2      2
3      2
4      3
5      3
6      3
7      3
8      4
9      4

It would be nice if HP71B has built-in binary search code.
Getting it right is hard.

Are you one of the 10% of programmers who can write a binary search?
Binary Search and Exclusive Upper Bounds

[Thomas Klemm]

(12-06-2022 12:15 AM)Albert Chan Wrote:  Getting it right is hard.

Are you one of the 10% of programmers who can write a binary search?
Binary Search and Exclusive Upper Bounds

Even Java had its problem with it:

• Extra, Extra - Read All About It: Nearly All Binary Searches and Mergesorts are Broken
• JDK-6412541 : Arrays.binarySearch does not work for arrays larger than 1<<30
  

[J-F Garnier]

(12-05-2022 02:55 PM)J-F Garnier Wrote:  Edit: removed my other proposal. made a mistake in my tests :-(

Since we are going towards complicate solutions, here is (finally) my proposal for the 71B:

DEF FNB(X)
IF X<1048576 THEN FNB=IP(LOG2(X))+1 ELSE FNB=IP(LOG2(X DIV 1048576))+21
END DEF

J-F

[Albert Chan]

All codes (sorted from slowest to fastest) produced a sum of 388874.
All codes were run combined with below test codes.

Test Code Wrote:10 DESTROY ALL @ RANDOMIZE 1 @ T=0 @ SETTIME 0
20 FOR I=1 TO 10000 @ T=T+FNB(IP(RND*999999999999)+1) @ NEXT I
30 DISP T,TIME

General binary search function turns out much slower, perhaps due to cost of array lookup.
But, binary search code might be useful for other applications, so I post it here.

Below, all timings done on Emu71/DOS WinXP (~ 200X speed)

Note: H=N+1 is an exclusive upper bound

Time = 18.6 seconds Wrote:15 N=39 @ DIM A(N) @ X=1 @ FOR I=1 TO N @ X=X+X @ A(I)=X @ NEXT I

100 DEF FNB(X) @ L=1 @ H=N+1
110 WHILE L<H @ M=IP((L+H)*.5) @ IF X<A(M) THEN H=M ELSE L=M+1
120 END WHILE @ FNB=L @ END DEF

With inclusive upper bound, binary search code look like this:

Code:

100 DEF FNB(X) @ L=1 @ H=N
110 WHILE L<=H @ M=IP((L+H)*.5) @ IF X<A(M) THEN H=M-1 ELSE L=M+1
120 END WHILE @ FNB=L @ END DEF

Avoiding array lookup, even with slow DIV, binary search is faster.

Time = 9.97 seconds Wrote:100 DEF FNB(X) @ B=1
110 IF X>=4294967296 THEN X=X DIV 4294967296 @ B=B+32
120 IF X>=65536 THEN X=X DIV 65536 @ B=B+16
130 IF X>=256 THEN X=X DIV 256 @ B=B+8
140 IF X>=16 THEN X=X DIV 16 @ B=B+4
150 FNB=B+(X>=2)+(X>=4)+(X>=8)
160 END DEF

Hard-coded binary search without DIV is possible, but messy, with many branches.
To give an idea of timing, I did a linear search for hex digit, then zero-in to bits.

Time = 8.25 seconds Wrote:100 DEF FNB(X) @ SELECT X
110 CASE >=68719476736 @ FNB=37+(X>=137438953472)+(X>=274877906944)+(X>=549755813888)
120 CASE >=4294967296 @ FNB=33+(X>=8589934592)+(X>=17179869184)+(X>=34359738368)
130 CASE >=268435456 @ FNB=29+(X>=536870912)+(X>=1073741824)+(X>=2147483648)
140 CASE >=16777216 @ FNB=25+(X>=33554432)+(X>=67108864)+(X>=134217728)
150 CASE >=1048576 @ FNB=21+(X>=2097152)+(X>=4194304)+(X>=8388608)
160 CASE >=65536 @ FNB=17+(X>=131072)+(X>=262144)+(X>=524288)
170 CASE >=4096 @ FNB=13+(X>=8192)+(X>=16384)+(X>=32768)
180 CASE >=256 @ FNB=9+(X>=512)+(X>=1024)+(X>=2048)
190 CASE >=16 @ FNB=5+(X>=32)+(X>=64)+(X>=128)
200 CASE ELSE @ FNB=1+(X>=2)+(X>=4)+(X>=8)
210 END SELECT @ END DEF

J-F Garnier code, limit size of LOG2 argument.

Time = 7.72 seconds Wrote:100 DEF FNB(X)
110 IF X<1048576 THEN FNB=IP(LOG2(X))+1 ELSE FNB=IP(LOG2(X DIV 1048576))+21
120 END DEF

J-F Garnier code, binary string length, fastest of them all! Thanks!

Time = 5.81 seconds Wrote:100 DEF FNB(N)=LEN(BSTR$(N,2))

[John Keith]

(12-05-2022 10:06 PM)robve Wrote:  For C (and lua etc) I came up with a simple binary search to find the position of the leading bit in a 32 bit integer in 5 search steps:

Code:

  unsigned int n, k, v = <number>;
  k = -!!(v >> 16) & 16;
  n = k;
  v >>= k;
  k = -!!(v >> 8) & 8;
  n += k;
  v >>= k;
  k = -!!(v >> 4) & 4;
  n += k;
  v >>= k;
  k = -!!(v >> 2) & 2;
  n += k;
  v >>= k;
  k = -!!(v >> 1) & 1;
  n += k;
  v >>= k;

This is similar to Joe Horn's program from 1-Minute Marvels:

Code:


\<< 0 1 ROT
  FOR c 1 + c
  STEP
\>>

which is especially nice for the 49 and 50 since it is accurate for integers of any size.

[robve]

(12-06-2022 06:49 PM)Albert Chan Wrote:  All codes (sorted from slowest to fastest) produced a sum of 388874.
All codes were run combined with below test codes.

Very nice. But you can avoid divisions (right shifts) solely by using multiplications (left shifts). Do you see what I mean and how to do that?

Also, in general, to find the leading bit in a n-bit binary number is asymptotically faster with binary search taking O(log(n)) asymptotic execution time as opposed to producing a string of 1 to n bits, which takes O(n) time.

The (im)practical aspects of the HP-71B makes binary search slow. Its BASIC is not suitable for this approach, especially when using costly divisions. Using hard coded branches and complicated expressions isn't helping performance either.

(12-06-2022 12:15 AM)Albert Chan Wrote:  Hi, robve
[...]
It would be nice if HP71B has built-in binary search code.
Getting it right is hard.
Are you one of the 10% of programmers who can write a binary search?
Binary Search and Exclusive Upper Bounds

The bit-shift binary search has no issues with edge cases, since we're dealing with integers of 32 bit (or 8, 16, 32, 64 bit for that matter, with a tweak to the code):
1: 0
2: 1
3: 1
4: 2
5: 2
6: 2
7: 2
8: 3
9: 3
10: 3
11: 3
12: 3
13: 3
14: 3
15: 3
16: 4
17: 4
...
511: 8
512: 9
...
4294967292: 31
4294967293: 31
4294967294: 31
4294967295: 31

Sure works like a charm. I like it because the underlying machine code is simple and efficient.

- Rob

[robve]

Here is a cryptic "one liner" in C for binary search with left shifts and a right shift by one bit (for division by two):

unsigned int k, n, u, v, w = <number>; for (k = 16, n = 0, v = 1, u = v << k; k; k >>= 1, u = v << k) if (w >= u) { n += k; v = u; }

or how about the same with multiplications and square roots, in case shifts aren't permitted:

unsigned int k, n, s, u, v, w = <number>; for (k = 16, n = 0, v = 1, u = s = 65536; k; k /= 2, s = sqrt(s), u = v * s) if (w >= u) { n += k; v = u; }

The square roots are simply 65536, 256, 16, 4, 2, so can be precomputed or put in an unrolled loop:

n = 0;
v = 1;
u = 65536;
if (w >= u) { n = 16; v = u; }
u = v * 256;
if (w >= u) { n += 8; v = u; }
u = v * 16;
if (w >= u) { n += 4; v = u; }
u = v * 4;
if (w >= u) { n += 2; v = u; }
u = v * 2;
if (w >= u) { n += 1; }

Add a step to extend to 64 bit numbers or a few more steps to go very far beyond if your machine can handle it.

- Rob

[Albert Chan]

(12-06-2022 07:26 PM)robve Wrote:  Very nice. But you can avoid divisions (right shifts) solely by using multiplications (left shifts).

Yes, it can be done this way too.

The problem is, for HP71B, *both* multiply and divide are costly.
Divide a bit worse, but it is in the same ballpark.

Programming in interpreted BASIC is very different than compiled languages.

Also, pow-of-2 and pow-of-10 numbers never match (except for 2^0 = 10^0 = 1)
This required an extra step to correct for the boundary. (bolded line)

Below could optimize a bit, but likely still slower than binary search with DIV (9.97 seconds)

Time = 11.55 seconds Wrote:10 DESTROY ALL @ RANDOMIZE 1 @ T=0 @ SETTIME 0
20 FOR I=1 TO 10000 @ T=T+FNB(IP(RND*999999999999)+1) @ NEXT I
30 DISP T,TIME

100 DEF FNB(X) @ B=40 ! FNB(x >= 2^39) assumed 40
110 Y=X*4294967296 @ IF Y<1E12 THEN X=Y @ B=B-32
120 Y=X*65536 @ IF Y<1E12 THEN X=Y @ B=B-16
130 Y=X*256 @ IF Y<1E12 THEN X=Y @ B=B-8
140 Y=X*16 @ IF Y<1E12 THEN X=Y @ B=B-4
150 Y=X*4 @ IF Y<1E12 THEN X=Y @ B=B-2
160 Y=X*2 @ IF Y<1E12 THEN X=Y @ B=B-1
170 FNB=B-(X<549755813888)
180 END DEF

[FLISZT]

Is there not a (short) way to do this with an HP-16C, or with any other calculator with a similar or "better" (DM / WP) instruction set?

("better" → for binary computations)

[Paul Dale]

Log₂ is the way here.

The 34S included this function.
The DM43/WP43/C43 will, likely, also include it.


Pauli

[Joe Horn]

Just for fun, here's my RPL solution which appeared first as a "How Does This Work" challenge many years ago. It also was included in Richard Nelson's "One Minute Marvels" collection.

Input: N as an Integer-type object, any size.
Output: Length of binary version of N.

<< 0 1 ROT FOR SWAP 1 + SWAP STEP >>

Not fast, but accurate, and different.

[Werner]

Of course, calling the local variable 'SWAP' is meant to confuse the uninitiated!
But a marvel it is, indeed.
Cheers, Werner

[J-F Garnier]

(12-07-2022 06:33 AM)Paul Dale Wrote:  Log₂ is the way here.

The 34S included this function.
The DM43/WP43/C43 will, likely, also include it.

As noted in the OP, the log2 method can fail when the argument is close to the maximum representable integer.
As does log10 fail to extract the number's exponent on various machines in the same condition:
HP41C, HP15C:
1E9 1 - LOG INT --> 9
Free42:
1E33 1 - LOG IP --> 33

It's just a rounding effect (correct results, no error).

But I agree that, in most cases when the argument range is known and limited, the log2 method is the simplest (but not necessary the fastest) way.

J-F

[Albert Chan]

(12-07-2022 08:41 AM)J-F Garnier Wrote:  But I agree that, in most cases when the argument range is known and limited, the log2 method is the simplest
(but not necessary the fastest) way.

log2(2^n ± x) = n + log2(1 ± x/2^n) ≈ n ± x/2^n/ln(2)

For 12 decimal digits, n <= 39
If last term is going to be round off, last term less than 0.5E-10

x/2^n/ln(2) < 0.5E-10
2^n/x > 0.839759 * 2^35

If x=1, we have n >= 35
If n=39, we have x < 2^4/0.839759 <= 19

>C=2^35
>LOG2(C-2),LOG2(C-1), LOG2(C+1), LOG2(C+2)
34.9999999999      35      35      35.0000000001
>C=2^39
>LOG2(C-20), LOG2(C-19), LOG2(C+19), LOG2(C+20)
38.9999999999      39      39      39.0000000001

The 2 numbers in the middle get rounded to integer, as expected.

--> bits(x) = IP(LOG2(x))+1 only work upto 2^35 - 2 = 34,359,738,366