How to init an array in Posix bc?

[ewmayer]

Originally Posted by bgbeuning

bc seems to be doing base 10 arithmetic.
gmp and most "advanced" packages do base 2^32 (or close) arithmetic.

Don't know what internals bc uses, no disputing it's shite compared to any halfway-decent compiled library - like I said, it would simply be nice to have some option that satisfies the following criteria:

1. Freeware;
2. No compile/install (with all the attendant versioning issues, nontrivial for us Luddites using older OSes and such) needed;
3. No internet pipe needed;
4. Decent set of number-theoretic tools such as gcd, powmod, modest factoring capability (say numbers somewhere in the 50-100 digit range);
5. User can add/modify functionality as desired.

No disputing that bc has some major drawbacks - lack of bitwise operations and binary shift being an obvious one. That needs workarounds to get acceptably fast functionality, for, say bit-counting - if you try to use the 1-bit-a-time bits() function in my above post on any input beyond 1000 bits or so, things slow to a crawl. So here are 2 alternatives - on a 500000-bit input, on my Core2-based Macbook classic, bits2() needs around 60 seconds to count the bits; bits3() needs under 5, still slow, but a nice speedup, and none of this is intended for serious work in the megabit size range:

/* Power-of-2 powers-of-2 needed by fast bit-counting function: pow2[i] = 2^2^i: */
pow2[0] = 2;
pow2[1] = pow2[0]^2;
pow2[2] = pow2[1]^2;
[etc]
/* Recursive divide-and-conquer version asymptotically faster than 1-bit-at-a-time: */
define bits2(n) {
	auto i,j,tmp;
	tmp = abs(n);
	if(tmp < 2) return tmp;
	i = 0; j = 0;	/* j has current bitcount */
	while(tmp >= pow2[i]) {
		i += 1;
	}
	j += 2^(i-1);
	tmp /= pow2[i-1];
	j += bits2(tmp);
	return(j);
}

/* Multiply much faster than divide, so this version faster than recursive div-based one: */
define bits3(n) {
	auto i,j,tmp,tm2;
	n = abs(n);
	if(n < 2) return n;
	i = 0; j = 1; tmp = 1;	/* j has current bitcount */
	while(n >= pow2[i++]) { }
	while(i--) {
		tm2 = tmp*pow2[i];
		if(tm2 < n) {
			j += 2^i;
			tmp = tm2;
		}
	}
	return(j);
}

[paulunderwood]

I found this page which might be of interest to you: http://www.numbertheory.org/gnubc/bc_programs.html

[paulunderwood]

Another way to find how many bits a number is, is to load mathlib and use logs, and afterwards use "scale=0" to get back to using integers

bc -lq
l(2^50000)/l(2)+1
50001.00000000000000052166

[paulunderwood]

Originally Posted by paulunderwood

Another way to find how many bits a number is, is to load mathlib and use logs, and afterwards use "scale=0" to get back to using integers

bc -lq
l(2^50000)/l(2)+1
50001.00000000000000052166

l(2^44497-1)/l(2)+1
44498.00000000000000046424
l(2^44497)/l(2)+1
44498.00000000000000046424

This method is not exact.

Here is my bit-counter:

define count_bits(n) {
        auto m, i;
        m = 1;
        i = 0;
        while( m <= n ) {
                m *= 2;
                i += 1;
        }
        return(i);
}

[paulunderwood]

length(n)/l(2)*l(10)

gives a ballpark figure that can tweaked for exactness

define fast_count_bits(n) {
        auto r;
        scale=0;
        r = length(n);
        scale=20;
        r = r/l(2)*l(10);
        scale = 0;
        r = r - r%1;
        while (2^r > n) { r -= 1 }
        return(r+1)
}

This gives some runtime warnings and returns a non-integer.

[paulunderwood]

Originally Posted by paulunderwood

define fast_count_bits(n) {
        auto r;
        scale=0;
        r = length(n);
        scale=20;
        r = r/l(2)*l(10);
        scale = 0;
        r = r - r%1;
        while (2^r > n) { r -= 1 }
        return(r+1)
}

This gives some runtime warnings and returns a non-integer.

Use r = (r - r%1)/1; instead!

[ewmayer]

Originally Posted by paulunderwood

I found this page which might be of interest to you: http://www.numbertheory.org/gnubc/bc_programs.html

Thanks - that looks very useful, I would've been surprised to have been the first person to envision such a load-text-code-and-go basic-comp-NT system built on bc.

Re. loading the mathlib, it would definitely be nice to have both pure-int and arbitrary-prec float avaialable and playing nice side-by-side, but doesn't that wreck some key number-theoretic functions? E.g. mod no longer works as expected.

[paulunderwood]

Originally Posted by ewmayer

Thanks - that looks very useful, I would've been surprised to have been the first person to envision such a load-text-code-and-go basic-comp-NT system built on bc.

Re. loading the mathlib, it would definitely be nice to have both pure-int and arbitrary-prec float avaialable and playing nice side-by-side, but doesn't that wreck some key number-theoretic functions? E.g. mod no longer works as expected.

Setting scale=0; ensures mod works.

[ewmayer]

Originally Posted by paulunderwood

Setting scale=0; ensures mod works.

So when mixing int and float operations we just need to make sure the former save and restore scale - nice. On large (by bc standards) inputs your bit-counter is 10-20% faster than my pure-int one; this sort of function may be a worst-case scenario here due to bc's lack of built-in bitwise functionality.

[ewmayer]

Paul, you were just one little step away from 'virtually instantaneous' in your FP-based bit-counting trials ... Since length() works in both int ('bc') and float ('bc -l') modes to count decimal digits and is very fast (validating the surmise that bc's internals are decimal-multiprecision-based), in either mode we can simply supplement that with a fixed-point-mul/div by a scaled version of log2(10), whose rightmost digit is a rounded-up truncation of said logarithm:

define bits(n) {
return length(n)*3321928095/1000000000;
}

[paulunderwood]

Originally Posted by ewmayer

Paul, you were just one little step away from 'virtually instantaneous' in your FP-based bit-counting trials ... Since length() works in both int ('bc') and float ('bc -l') modes to count decimal digits and is very fast (validating the surmise that bc's internals are decimal-multiprecision-based), in either mode we can simply supplement that with a fixed-point-mul/div by a scaled version of log2(10), whose rightmost digit is a rounded-up truncation of said logarithm:

define bits(n) {
return length(n)*3321928095/1000000000;
}

Without -l, just bc -q ...

bits(2^100)
102
bits(2^101)
102
bits(2^102)
102

You still need to home in from a ballpark figure. Right?

The revised function, assuming no mathlib loaded or mathlib loaded and scale=0:

define fast_count_bits(n) {
        auto r;
        r = length(n)*3321928095/1000000000;
        while ( 2^r > n ) { r -= 1; }
        return(r+1)
}