add log2(x [,eps]) builtin function

Added log2(x [,eps]) builtin function.  When x is an integer
power of 2, log2(x) will return an integer, otherwise it will
return the equivalent of ln(x)/ln(2).

CHANGES

@@ -80,6 +80,10 @@ The following are the changes from calc version 2.14.3.5 to date:
     comments in "README.RELEASE", which serves as the contents
     of the calc command "help release".
 
+    Added log2(x [,eps]) builtin function.  When x is an integer
+    power of 2, log2(x) will return an integer, otherwise it will
+    return the equivalent of ln(x)/ln(2).
+
 
 The following are the changes from calc version 2.14.3.4 to 2.14.3.5:
 

LIBRARY

@@ -1,4 +1,4 @@
-	USING THE ARBITRARY PRECISION ROUTINES IN A C PROGRAM
+USING THE ARBITRARY PRECISION ROUTINES IN A C PROGRAM
 
 Part of the calc release consists of an arbitrary precision math link library.
 This link library is used by the calc program to perform its own calculations.
@@ -393,6 +393,12 @@ ZVALUE yourself.  Examples of such checks are:
 	zge32b(z)	(number is >= 2^32)
 	zge64b(z)	(number is >= 2^64)
 
+Return true if z an integer power of 2: set zlog2 to log base 2 of z.
+Return false if a is NOT integer power of 2: set zlog2 to 0.
+The zlog2 arg must be a non-NULL pointer to a ZVALUE.
+
+	zispowerof2(z, &zlog2)
+
 Typically the largest unsigned long is typedefed to FULL.  The following
 macros are useful in dealing with this data type:
 
@@ -524,7 +530,7 @@ macros are:
 	qistwo(q)		(number is 2)
 	qiseven(q)		(number is an even integer)
 	qisodd(q)		(number is an odd integer)
-	qistwopower(q)	(number is a power of 2 >= 1)
+	qisreciprocal(q)	(number is 1 / an integer and q != 0)
 
 The comparisons for NUMBERs are similar to the ones for ZVALUEs.  You use the
 qcmp and qrel functions.

cal/regress.cal

@@ -1,7 +1,7 @@
 /*
  * regress - calc regression and correctness test suite
  *
- * Copyright (C) 1999-2017,2021  David I. Bell and Landon Curt Noll
+ * Copyright (C) 1999-2017,2021,2023  David I. Bell and Landon Curt Noll
  *
  * Calc is open software; you can redistribute it and/or modify it under
  * the terms of the version 2.1 of the GNU Lesser General Public License
@@ -2997,6 +2997,45 @@ define test_2600()
 	    strcat(str(tnum++),
 		   ': round(log(1.2+1.2i,1e-10),10) == ',
 		   '0.2296962439+0.3410940885i'));
+	vrfy(log2(2) == 1,
+			    strcat(str(tnum++), ': log2(2) == 1'));
+	vrfy(log2(4) == 2,
+			    strcat(str(tnum++), ': log2(4) == 2'));
+	vrfy(log2(1024) == 10,
+			    strcat(str(tnum++), ': log2(1024) == 10'));
+	vrfy(log2(2^500) == 500,
+			    strcat(str(tnum++), ': log2(2^500) == 500'));
+	vrfy(log2(1/2^23209) == -23209,
+			    strcat(str(tnum++), ': log2(1/2^23209) == -23209'));
+	vrfy(isint(log2(1/2^23209)),
+			    strcat(str(tnum++), ': isint(log2(1/2^23209))'));
+	vrfy(round(log2(127),10) == 6.9886846868,
+	    strcat(str(tnum++),
+		   ': round(log2(127),10) == 6.9886846868'));
+	vrfy(round(log2(23209),10) == 14.5023967426,
+	    strcat(str(tnum++),
+		   ': round(log2(23209),10) == 14.5023967426'));
+	vrfy(round(log2(2^17-19),10) == 16.9997908539,
+	    strcat(str(tnum++),
+		   ': round(log2(2^17-19),10) == 16.9997908539'));
+	vrfy(round(log2(-127),10) == 6.9886846868+4.5323601418i,
+	    strcat(str(tnum++),
+		   ': round(log2(-127),10) == 6.9886846868+4.5323601418i'));
+	vrfy(round(log2(2+3i),10) == 1.8502198591+1.4178716307i,
+	    strcat(str(tnum++),
+		   ': round(log2(2+3i),10) == 1.8502198591+1.4178716307i'));
+	vrfy(round(log2(-2+3i),10) == 1.8502198591+3.1144885111i,
+	    strcat(str(tnum++),
+		   ': round(log2(-2+3i),10) == 1.8502198591+3.1144885111i'));
+	vrfy(round(log2(2-3i),10) == 1.8502198591-1.4178716307i,
+	    strcat(str(tnum++),
+		   ': round(log2(2-3i),10) == 1.8502198591-1.4178716307i'));
+	vrfy(round(log2(-2-3i),10) == 1.8502198591-3.1144885111i,
+	    strcat(str(tnum++),
+		   ': round(log2(-2-3i),10) == 1.8502198591-3.1144885111i'));
+	vrfy(round(log2(17+0.3i),10) == 4.0876874474+0.0254566819i,
+	    strcat(str(tnum++),
+		   ': round(log2(17+0.3i),10) == 4.0876874474+0.0254566819i'));
         epsilon(i),;
 	print tnum++: ': epsilon(i),;';
 

calcerr.tbl

@@ -1,7 +1,7 @@
 #
 # calcerr - error codes and messages
 #
-# Copyright (C) 1999-2006,2021  Ernest Bowen
+# Copyright (C) 1999-2006,2021,2023  Ernest Bowen and Landon Curt Noll
 #
 # Calc is open software; you can redistribute it and/or modify it under
 # the terms of the version 2.1 of the GNU Lesser General Public License
@@ -542,3 +542,6 @@ E_HMS2H1	Non-real-number arguments 1, 2 or 3 for hms2h
 E_HMS2H2	Invalid rounding arg 4 for hms2h
 E_HM2H1		Non-real-number arguments 1 or 2 for hm2h
 E_HM2H2		Invalid rounding arg 4 for hm2h
+E_LOG2_1	Bad epsilon argument for log2
+E_LOG2_2	Non-numeric first argument for log2
+E_LOG2_3	Cannot calculate log2 for this value

cmath.h

@@ -1,7 +1,7 @@
 /*
  * cmath - data structures for extended precision complex arithmetic
  *
- * Copyright (C) 1999-2007,2014  David I. Bell
+ * Copyright (C) 1999-2007,2014,2023  David I. Bell and Landon Curt Noll
  *
  * Calc is open software; you can redistribute it and/or modify it under
  * the terms of the version 2.1 of the GNU Lesser General Public License
@@ -97,6 +97,7 @@ E_FUNC COMPLEX *c_root(COMPLEX *c, NUMBER *q, NUMBER *epsilon);
 E_FUNC COMPLEX *c_exp(COMPLEX *c, NUMBER *epsilon);
 E_FUNC COMPLEX *c_ln(COMPLEX *c, NUMBER *epsilon);
 E_FUNC COMPLEX *c_log(COMPLEX *c, NUMBER *epsilon);
+E_FUNC COMPLEX *c_log2(COMPLEX *c, NUMBER *epsilon);
 E_FUNC COMPLEX *c_cos(COMPLEX *c, NUMBER *epsilon);
 E_FUNC COMPLEX *c_sin(COMPLEX *c, NUMBER *epsilon);
 E_FUNC COMPLEX *c_cosh(COMPLEX *c, NUMBER *epsilon);

comfunc.c

@@ -1,7 +1,7 @@
 /*
  * comfunc - extended precision complex arithmetic non-primitive routines
  *
- * Copyright (C) 1999-2007,2021-2023  David I. Bell and Ernest Bowen
+ * Copyright (C) 1999-2007,2021-2023  David I. Bell, Landon Curt Noll and Ernest Bowen
  *
  * Primary author:  David I. Bell
  *
@@ -39,9 +39,13 @@
  */
 STATIC COMPLEX *cln_10 = NULL;
 STATIC NUMBER *cln_10_epsilon = NULL;
+STATIC COMPLEX *cln_2 = NULL;
+STATIC NUMBER *cln_2_epsilon = NULL;
 STATIC NUMBER _q10_ = { { _tenval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
+STATIC NUMBER _q2_ = { { _twoval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
 STATIC NUMBER _q0_ = { { _zeroval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
 COMPLEX _cten_ = { &_q10_, &_q0_, 1 };
+COMPLEX _ctwo_ = { &_q2_, &_q0_, 1 };
 
 
 /*
@@ -530,6 +534,7 @@ c_ln(COMPLEX *c, NUMBER *epsilon)
 	return r;
 }
 
+
 /*
  * Calculate base 10 logarithm by:
  *
@@ -546,7 +551,7 @@ c_log(COMPLEX *c, NUMBER *epsilon)
 	 * compute ln(c) first
 	 */
 	ln_c = c_ln(c, epsilon);
-	/* log(1) == 0 */
+	/* quick return for log(1) == 0 */
 	if (ciszero(ln_c)) {
 		return ln_c;
 	}
@@ -585,6 +590,61 @@ c_log(COMPLEX *c, NUMBER *epsilon)
 }
 
 
+/*
+ * Calculate base 2 logarithm by:
+ *
+ *	log(c) = ln(c) / ln(2)
+ */
+COMPLEX *
+c_log2(COMPLEX *c, NUMBER *epsilon)
+{
+	int need_new_cln_2 = true;	/* false => use cached cln_2 value */
+	COMPLEX *ln_c;			/* ln(x) */
+	COMPLEX *ret;			/* base 2 of x */
+
+	/*
+	 * compute ln(c) first
+	 */
+	ln_c = c_ln(c, epsilon);
+	/* quick return for log(1) == 0 */
+	if (ciszero(ln_c)) {
+		return ln_c;
+	}
+
+	/*
+	 * save epsilon for ln(2) if needed
+	 */
+	if (cln_2_epsilon == NULL) {
+		/* first time call */
+		cln_2_epsilon = qcopy(epsilon);
+	} else if (qcmp(cln_2_epsilon, epsilon) == true) {
+		/* replaced cached value with epsilon arg */
+		qfree(cln_2_epsilon);
+		cln_2_epsilon = qcopy(epsilon);
+	} else if (cln_2 != NULL) {
+		/* the previously computed ln(2) is OK to use */
+		need_new_cln_2 = false;
+	}
+
+	/*
+	 * compute ln(2) if needed
+	 */
+	if (need_new_cln_2 == true) {
+		if (cln_2 != NULL) {
+			comfree(cln_2);
+		}
+		cln_2 = c_ln(&_ctwo_, cln_2_epsilon);
+	}
+
+	/*
+	 * return ln(c) / ln(2)
+	 */
+	ret = c_div(ln_c, cln_2);
+	comfree(ln_c);
+	return ret;
+}
+
+
 /*
  * Calculate the complex cosine within the specified accuracy.
  * This uses the formula:

func.c

@@ -2192,6 +2192,55 @@ f_log(int count, VALUE **vals)
 }
 
 
+S_FUNC VALUE
+f_log2(int count, VALUE **vals)
+{
+	VALUE result;
+	COMPLEX ctmp, *c;
+	NUMBER *err;
+
+	/* initialize VALUE */
+	result.v_subtype = V_NOSUBTYPE;
+
+	err = conf->epsilon;
+	if (count == 2) {
+		if (vals[1]->v_type != V_NUM)
+			return error_value(E_LOG2_1);
+		err = vals[1]->v_num;
+	}
+	switch (vals[0]->v_type) {
+		case V_NUM:
+			if (!qisneg(vals[0]->v_num) &&
+			    !qiszero(vals[0]->v_num)) {
+				result.v_num = qlog2(vals[0]->v_num, err);
+				result.v_type = V_NUM;
+				return result;
+			}
+			ctmp.real = vals[0]->v_num;
+			ctmp.imag = qlink(&_qzero_);
+			ctmp.links = 1;
+			c = c_log2(&ctmp, err);
+			break;
+		case V_COM:
+			c = c_log2(vals[0]->v_com, err);
+			break;
+		default:
+			return error_value(E_LOG2_2);
+	}
+	if (c == NULL) {
+		return error_value(E_LOG2_3);
+	}
+	result.v_type = V_COM;
+	result.v_com = c;
+	if (cisreal(c)) {
+		result.v_num = qlink(c->real);
+		result.v_type = V_NUM;
+		comfree(c);
+	}
+	return result;
+}
+
+
 S_FUNC VALUE
 f_cos(int count, VALUE **vals)
 {
@@ -10165,10 +10214,8 @@ STATIC CONST struct builtin builtins[] = {
 	 "hypotenuse of right triangle within accuracy c"},
 	{"ilog", 2, 2, 0, OP_NOP, {.null = NULL}, {.valfunc_2 = f_ilog},
 	 "integral log of a to integral base b"},
-#if 0 /* XXX - to be added later */
 	{"ilogn", 2, 2, 0, OP_NOP, {.null = NULL}, {.valfunc_2 = f_ilog},
 	 "same is ilog"},
-#endif /* XXX */
 	{"ilog10", 1, 1, 0, OP_NOP, {.null = NULL}, {.valfunc_1 = f_ilog10},
 	 "integral log of a number base 10"},
 	{"ilog2", 1, 1, 0, OP_NOP, {.null = NULL}, {.valfunc_1 = f_ilog2},
@@ -10268,6 +10315,8 @@ STATIC CONST struct builtin builtins[] = {
 	 "natural logarithm of value a within accuracy b"},
 	{"log", 1, 2, 0, OP_NOP, {.null = NULL}, {.valfunc_cnt = f_log},
 	 "base 10 logarithm of value a within accuracy b"},
+	{"log2", 1, 2, 0, OP_NOP, {.null = NULL}, {.valfunc_cnt = f_log2},
+	 "base 2 logarithm of value a within accuracy b"},
 	{"lowbit", 1, 1, 0, OP_LOWBIT, {.null = NULL}, {.null = NULL},
 	 "low bit number in base 2 representation"},
 	{"ltol", 1, 2, FE, OP_NOP, {.numfunc_2 = f_legtoleg}, {.null = NULL},

help/log2

@@ -14,21 +14,30 @@ DESCRIPTION
     Approximate the base 2 logarithm function of x by a multiple of
     epsilon, the error having absolute value less than 0.75 * eps.
 
+    When x is an integer power of 2, log2(x) will return an integer
+    regardless of the value of eps or epsilon().
+
     If y is a positive integer, log(x, 2^-y) will usually be correct
     to the y-th decimal place.
 
-    When x if a power of 2, log2(x) will return an integer regardless
-    of the value of eps or epsilon().
-
 EXAMPLE
-    ; print log2(2), log2(4), log2(1024), log2(2^500)
+    ; print log2(2), log2(4), log2(1024), log2(2^500), log2(1/2^23209)
-    1 2 10 500
+    1 2 10 500 -23209
 
-    ; print log(127), log(23209), log(2^17-19)
+    ; print log2(127), log2(23209), log2(2^17-19)
-    ~6.98868468677216585326 ~14.50239674255407864468 ~16.99979085393743521984
+    6.98868468677216585326 14.50239674255407864468 16.99979085393743521984
 
-    ; print log(2+3i, 1e-5)
+    ; print log2(-127)
-    ~1.85020558320709803073+~1.41786049195700786266i
+    6.98868468677216585326+4.53236014182719380961i
+
+    ; print log2(2+3i), log2(-2+3i)
+    1.85021985907054608020+1.41787163074572199658i 1.85021985907054608020+3.11448851108147181304i
+
+    ; print log2(2-3i), log2(-2-3i)
+    1.85021985907054608020-1.41787163074572199658i 1.85021985907054608020-3.11448851108147181304i
+
+    ; print log2(17+0.3i, 1e-75), log2(-17-0.3i, 1e-75)
+    4.08768744737521449955+0.02545668190826163773i 4.08768744737521449955-4.50690345991893217190i
 
 LIMITS
     x != 0

qmath.c

@@ -1,7 +1,7 @@
 /*
  * qmath - extended precision rational arithmetic primitive routines
  *
- * Copyright (C) 1999-2007,2014,2021-2023  David I. Bell and Ernest Bowen
+ * Copyright (C) 1999-2007,2014,2021-2023  David I. Bell, Landon Curt Noll and Ernest Bowen
  *
  * Primary author:  David I. Bell
  *

qmath.h

@@ -1,7 +1,7 @@
 /*
  * qmath - declarations for extended precision rational arithmetic
  *
- * Copyright (C) 1999-2007,2014,2021  David I. Bell
+ * Copyright (C) 1999-2007,2014,2021,2023  David I. Bell and Landon Curt Noll
  *
  * Calc is open software; you can redistribute it and/or modify it under
  * the terms of the version 2.1 of the GNU Lesser General Public License
@@ -185,6 +185,7 @@ E_FUNC NUMBER *qsin(NUMBER *q, NUMBER *epsilon);
 E_FUNC NUMBER *qexp(NUMBER *q, NUMBER *epsilon);
 E_FUNC NUMBER *qln(NUMBER *q, NUMBER *epsilon);
 E_FUNC NUMBER *qlog(NUMBER *q, NUMBER *epsilon);
+E_FUNC NUMBER *qlog2(NUMBER *q, NUMBER *epsilon);
 E_FUNC NUMBER *qtan(NUMBER *q, NUMBER *epsilon);
 E_FUNC NUMBER *qsec(NUMBER *q, NUMBER *epsilon);
 E_FUNC NUMBER *qcot(NUMBER *q, NUMBER *epsilon);
@@ -248,12 +249,12 @@ E_FUNC NUMBER *swap_HALF_in_NUMBER(NUMBER *dest, NUMBER *src, bool all);
 #define qistwo(q)	(zistwo((q)->num) && zisunit((q)->den))
 #define qiseven(q)	(zisunit((q)->den) && ziseven((q)->num))
 #define qisodd(q)	(zisunit((q)->den) && zisodd((q)->num))
-#define qistwopower(q)	(zisunit((q)->den) && zistwopower((q)->num))
 #define qistiny(q)	(zistiny((q)->num))
 
 #define qhighbit(q)	(zhighbit((q)->num))
 #define qlowbit(q)	(zlowbit((q)->num))
 #define qdivcount(q1, q2)	(zdivcount((q1)->num, (q2)->num))
+#define qisreciprocal(q)	(zisunit((q)->num) && !ziszero((q)->den))
 /* operation on #q may be undefined, so replace with an inline-function */
 /* was: #define qlink(q)	((q)->links++, (q)) */
 static inline NUMBER* qlink(NUMBER* q) { if(q) { (q)->links++; } return q; }

qtrans.c

@@ -1,7 +1,7 @@
 /*
  * qtrans - transcendental functions for real numbers
  *
- * Copyright (C) 1999-2007,2021-2023  David I. Bell and Ernest Bowen
+ * Copyright (C) 1999-2007,2021-2023  David I. Bell, Landon Curt Noll and Ernest Bowen
  *
  * Primary author:  David I. Bell
  *
@@ -43,10 +43,12 @@ HALF _qlgeden_[] = { 25469 };
 NUMBER _qlge_ = { { _qlgenum_, 1, 0 }, { _qlgeden_, 1, 0 }, 1, NULL };
 
 /*
- * cache the natural logarithm of 10
+ * cache the natural logarithm of 10 and 2
  */
 STATIC NUMBER *ln_10 = NULL;
 STATIC NUMBER *ln_10_epsilon = NULL;
+STATIC NUMBER *ln_2 = NULL;
+STATIC NUMBER *ln_2_epsilon = NULL;
 
 /*
  * cache pi
@@ -1184,7 +1186,7 @@ qlog(NUMBER *q, NUMBER *epsilon)
 	 * compute ln(c) first
 	 */
 	ln_q = qln(q, epsilon);
-	/* log(1) == 0 */
+	/* quick return for log(1) == 0 */
 	if (qiszero(ln_q)) {
 		return ln_q;
 	}
@@ -1223,6 +1225,142 @@ qlog(NUMBER *q, NUMBER *epsilon)
 }
 
 
+/*
+ * Calculate the base 2 logarithm
+ *
+ *	log(q) = ln(q) / ln(2)
+ */
+NUMBER *
+qlog2(NUMBER *q, NUMBER *epsilon)
+{
+	int need_new_ln_2 = true;	/* false => use cached ln_2 value */
+	NUMBER *ln_q;			/* ln(x) */
+	NUMBER *ret;			/* base 2 logarithm of x */
+
+	/* firewall */
+	if (qiszero(q) || qiszero(epsilon)) {
+		math_error("logarithm of 0");
+		not_reached();
+	}
+
+	/*
+	 * special case: q is integer power of 2
+	 *
+	 * When q is integer power of 2, the base 2 logarithm is an integer.
+	 * We return a base 2 logarithm that is in integer in this case.
+	 *
+	 * From above we know that q != 0, so we need to check that q>0.
+	 *
+	 * We have two cases for a power of two: a positive power of 2, and a
+	 * negative power of 2 (i.e., 1 over a power of 2).
+	 */
+	if (qispos(q)) {
+		ZVALUE zlog2;		/* base 2 logarithm when q is power of 2 */
+
+		/*
+		 * case: q>0 is an integer
+		 */
+		if (qisint(q)) {
+
+			/*
+			 * check if q is an integer power of 2
+			 */
+			if (zispowerof2(q->num, &zlog2)) {
+
+				/*
+				 * case: q>0 is an integer power of 2
+				 *
+				 * Return zlog2, which is an integer power of 2 as a NUMBER.
+				 */
+				ret = qalloc();
+				zcopy(zlog2, &ret->num);
+				zfree(zlog2);
+				return ret;
+
+			/*
+			 * case: integer q is not an integer power of 2
+			 */
+			} else {
+
+				/* free the zispowerof2() 2nd arg and proceed to calculate ln(q)/ln(2) */
+				zfree(zlog2);
+			}
+
+		/*
+		 * case: q>0 is 1 over an integer
+		 */
+		} else if (qisreciprocal(q)) {
+
+			/*
+			 * check if q is 1 over an integer power of 2
+			 */
+			if (zispowerof2(q->den, &zlog2)) {
+
+				/*
+				 * case: q>0 is an integer power of 2
+				 *
+				 * Return zlog2, which is an negative integer power of 2 as a NUMBER.
+				 */
+				ret = qalloc();
+				zlog2.sign = !zlog2.sign;	/* zlog2 = -zlog2 */
+				zcopy(zlog2, &ret->num);
+				zfree(zlog2);
+				return ret;
+
+			/*
+			 * case: reciprocal q is not an integer power of 2
+			 */
+			} else {
+
+				/* free the zispowerof2() 2nd arg and proceed to calculate ln(q)/ln(2) */
+				zfree(zlog2);
+			}
+		}
+	}
+
+	/*
+	 * compute ln(c) first
+	 */
+	ln_q = qln(q, epsilon);
+	/* quick return for log(1) == 0 */
+	if (qiszero(ln_q)) {
+		return ln_q;
+	}
+
+	/*
+	 * save epsilon for ln(2) if needed
+	 */
+	if (ln_2_epsilon == NULL) {
+		/* first time call */
+		ln_2_epsilon = qcopy(epsilon);
+	} else if (qcmp(ln_2_epsilon, epsilon) == true) {
+		/* replaced cached value with epsilon arg */
+		qfree(ln_2_epsilon);
+		ln_2_epsilon = qcopy(epsilon);
+	} else if (ln_2 != NULL) {
+		/* the previously computed ln(2) is OK to use */
+		need_new_ln_2 = false;
+	}
+
+	/*
+	 * compute ln(2) if needed
+	 */
+	if (need_new_ln_2 == true) {
+		if (ln_2 != NULL) {
+			qfree(ln_2);
+		}
+		ln_2 = qln(&_qtwo_, ln_2_epsilon);
+	}
+
+	/*
+	 * return ln(q) / ln(2)
+	 */
+	ret = qqdiv(ln_q, ln_2);
+	qfree(ln_q);
+	return ret;
+}
+
+
 /*
  * Calculate the result of raising one number to the power of another.
  * The result is calculated to the nearest or next to nearest multiple of

version.c

@@ -3,7 +3,7 @@
  *
  * Copyright (C) 1999-2023  David I. Bell and Landon Curt Noll
  *
- * Primary author:  David I. Bell
+ * See "version.h" for the actual calc version constants.
  *
  * Calc is open software; you can redistribute it and/or modify it under
  * the terms of the version 2.1 of the GNU Lesser General Public License

zfunc.c

@@ -2236,3 +2236,81 @@ zissquare(ZVALUE z)
 	zfree(tmp);
 	return (n ? true : false);
 }
+
+
+/*
+ * test if a number is an integer power of 2
+ *
+ * given:
+ *	z	value to check if it is a power of 2
+ *	zlog2	>= 0 ==> *zlog2 power of 2 of z (return true)
+ *		< 0 z is not a power of 2 (return false)
+ *
+ * returns:
+ *	true	z is a power of 2
+ *	false	z is not a power of 2
+ */
+bool
+zispowerof2(ZVALUE z, ZVALUE *zlog2)
+{
+	FULL ilogz=0;		/* potential log base 2 return value or -1 */
+	HALF tophalf;		/* most significant HALF in z */
+	LEN len;		/* length of z in HALFs */
+	int i;
+
+	/* firewall */
+	if (zlog2 == NULL) {
+		math_error("%s: zlog2 NULL", __func__);
+		not_reached();
+	}
+
+	/* zero and negative values are never integer powers of 2 */
+	if (ziszero(z) || zisneg(z)) {
+		*zlog2 = _neg_one_;
+		return false;
+	}
+
+	/*
+	 * trim z just in case
+	 *
+	 * An untrimmed z will give incorrect results.
+	 */
+	ztrim(&z);
+
+	/*
+	 * all HALFs below the top HALF must be zero
+	 */
+	len = z.len;
+	for (i=0, ilogz=0; i < len-1; ++i, ilogz+=BASEB) {
+		if (z.v[i] != 0) {
+			*zlog2 = _neg_one_;
+			return false;
+		}
+	}
+
+	/*
+	 * top HALF must be a power of 2
+	 *
+	 * For non-zero values of tophalf,
+	 * (tophalf & (tophalf-1)) == 0 ==> tophalf is a power of 2,
+	 * (tophalf & (tophalf-1)) != 0 ==> tophalf is NOT a power of 2.
+	 */
+	tophalf = z.v[len-1];
+	if ((tophalf == 0) || ((tophalf & (tophalf-1)) != 0)) {
+		*zlog2 = _neg_one_;
+		return false;
+	}
+
+	/*
+	 * count the bits in the top HALF which is a power of 2
+	 */
+	for (i=0; i < BASEB && tophalf != ((HALF)1 << i); ++i) {
+		++ilogz;
+	}
+
+	/*
+	 * return power of 2
+	 */
+	utoz(ilogz, zlog2);
+	return true;
+}

zmath.c

@@ -1,7 +1,7 @@
 /*
  * zmath - extended precision integral arithmetic primitives
  *
- * Copyright (C) 1999-2007,2021-2023  David I. Bell and Ernest Bowen
+ * Copyright (C) 1999-2007,2021-2023  David I. Bell, Landon Curt Noll and Ernest Bowen
  *
  * Primary author:  David I. Bell
  *
@@ -2163,34 +2163,3 @@ zpopcnt(ZVALUE z, int bitval)
 	 */
 	return cnt;
 }
-
-
-#if 0 /* XXX - to be added later */
-/*
- * test if a number is a power of 2
- *
- * given:
- *	z	value to check if it is a power of 2
- *	zlog2	set to log base 2 of z if z is a power of 2, 0 otherwise
- *
- * returns:
- *	true	z is a power of 2
- *	false	z is not a power of 2
- */
-bool
-zispowerof2(ZVALUE z, ZVALUE *zlog2)
-{
-	/* firewall */
-	if (zlog2 == NULL) {
-		math_error("%s: zlog2 NULL", __func__);
-		not_reached();
-	}
-
-	/* zero and negative values are never powers of 2 */
-	if (ziszero(z) || zisneg(z)) {
-		return false;
-	}
-
-	/* XXX - add code here - XXX */
-}
-#endif /* XXX */

zmath.h

@@ -1,7 +1,7 @@
 /*
  * zmath - declarations for extended precision integer arithmetic
  *
- * Copyright (C) 1999-2007,2014,2021,2023  David I. Bell
+ * Copyright (C) 1999-2007,2014,2021,2023  David I. Bell and Landon Curt Noll
  *
  * Calc is open software; you can redistribute it and/or modify it under
  * the terms of the version 2.1 of the GNU Lesser General Public License
@@ -415,6 +415,7 @@ E_FUNC FLAG zsqrt(ZVALUE z1, ZVALUE *dest, long R);
 E_FUNC void zroot(ZVALUE z1, ZVALUE z2, ZVALUE *dest);
 E_FUNC bool zissquare(ZVALUE z);
 E_FUNC void zhnrmod(ZVALUE v, ZVALUE h, ZVALUE zn, ZVALUE zr, ZVALUE *res);
+E_FUNC bool zispowerof2(ZVALUE z, ZVALUE *zlog2);
 
 
 /*