//Headers
#define _USE_MATH_DEFINES

#include <iostream>
#include <set>
#include <map>
#include <algorithm>
#include <cmath>
#include <deque>
#include <list>
#include <vector>
#include <cstdlib>
#include <cstdio>
#include <fstream>
#include <climits>
#include <string>
#include <cstring>
#include <cfloat>
#include <sstream>
#include <utility>
#include <stack>

using namespace std;

//cmath constants
M_E 		e 		2.71828182845904523536
M_LOG2E		log2(e)		1.44269504088896340736
M_LOG10E	log10(e)	0.434294481903251827651
M_LN2		ln(2)		0.693147180559945309417
M_LN10		ln(10)		2.30258509299404568402
M_PI		pi		3.14159265358979323846
M_PI_2		pi/2		1.57079632679489661923
M_PI_4		pi/4		0.785398163397448309616
M_1_PI		1/pi		0.318309886183790671538
M_2_PI		2/pi		0.636619772367581343076
M_2_SQRTPI	2/sqrt(pi)	1.12837916709551257390
M_SQRT2		sqrt(2)		1.41421356237309504880
M_SQRT1_2	1/sqrt(2)	0.707106781186547524401

//This check if x is a perfect square
bool isPerfectSquare(int x)
{
	int s = sqrt(x);
	return s*s == x;
}
 
//This check if n is a Fibonacci Number
bool isFibonacci(int n)
{
	return isPerfectSquare(5*n*n + 4) || isPerfectSquare(5*n*n - 4);
}

//v will contains all primes between 2 and size (inclusive)
void sieve(vector<size_t>& v, size_t size)
{
	if(size < 1)
		return;

	v.resize(size, 1);
	for(int i = 2; i * i < size; ++i)
		if(v[i])
			for(int j = i + i; j < size; j += i)
				v[j] = 0;
}
//How to get the primes
/*for(int i = 0; i < v.size(); ++i)
	if(v[i] == 1)
		printf("%zu\n", i);*/

//Compute n!
size_t factorial(size_t n)
{
	if(n < 2)
		return 1;

	size_t prod = n;
	while(n > 1)
		prod *= --n;

	return prod;
}

//(n!)/(n-r)!
//Where n is the number of things to choose from, and we choose r of them (No repetition, order matters)
size_t permutation(size_t total, size_t amount)
{
	return factorial(total)/factorial(total - amount);
}

//(n!)/(r!(n-r)!)
//Where n is the number of things to choose from, and we choose r of them (No repetition, order doesn't matter)
size_t combinaison(size_t total, size_t amount)
{
		size_t top = factorial(total);
		size_t bottom = factorial(amount) * factorial(total - amount);
		return top/bottom;
}

//Combinations with repetition ::: 1 , 1 , 1 ::: 1 , 1 , 2 ::: ... ::: 3 , 3 , 3
//(n+r-1)!/(r!(n-r)!)
//Where n is the number of things to choose from, and we choose r of them (Repetition allowed, order doesn't matter) 
size_t combinaison_rep(size_t total, size_t amount)
{
	size_t top = factorial(total + amount - 1);
	size_t bottom = factorial(amount) * factorial(total - amount);
	return top/bottom;
}

//Log base n
float logbn(int val, int base)
{
		return log(val)/log(base);
}

//Greatest common divider
size_t gcd(size_t a, size_t b)
{
	while(b)
	{
		size_t x = a % b;
		a = b;
		b = x;
	}

	return a;
}

//Distance between two points
double distance(double x1, double y1, double x2, double y2)
{
	return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}

//Check if the point c is on the line segment defined by a and b
bool isonline(pair<int, int> a, pair<int, int> b, pair<int, int> c)
{
	return acos(abs(a.first-b.first)/dist(a, b)) == acos(abs(a.first-c.first)/dist(a, c));
}

//This function can be modified if you already have the length of the three sides
double AreaTriangle(double x1, double y1, double x2, double y2, double x3, double y3)
{
	double dab = distance(x1, y1, x2, y2);
	double dbc = distance(x2, y2, x3, y3);
	double dca = distance(x3, y3, x1, y1);

	double halfP = (dab+dbc+dca)/2.0;
	return sqrt(halfP*(halfP-dab)*(halfP-dbc)*(halfP-dca));
}

//Make sure that the line is not vertical
double DistenceFromLineWithSlope(double m1, double b1, double x, double y)
{
	return fabs(y - m1*x - b1)/sqrt(m1*m1 + 1);	
}

//Backup in case of problem
double DistenceFromLineWithSlope2(double m, double b, double x, double y)
{
	if(m == 0.0)
		return fabs(b - y);

	double A = fabs((m*x + b) - y);

	double B = fabs((y/m - b/m) - x);

	if(A == 0.0 && B == 0.0)
		return 0.0;

	return sqrt(A*A*B*B/(A*A+B*B));
}

//Backup in case of problem
double DistenceFromLineWithSlope3(double m1, double b1, double x, double y)
{
	if(m1 == 0.0)
		return fabs(b1 - y);

	double m2 = -1.0/m1;
	double b2 = y - m2*x;

	double xi = (b2 - b1)/(m1 - m2);
	double yi = (m2*b1 - m1*b2)/(m2 - m1);

	return distance(xi, yi, x, y);
}

//Horizontal distance between the point (x, y) and the line segment defined by (x1, y1), (x2, y2)
//Note: this contition must be true for this to work : (y1 >= y2 && y >= y1 && y <= y2) 
double xdistance(double x, double y, double x1, double y1, double x2, double y2)
{
	if(x1 == x2 && y1 == y2)
		return distance(x, y, x1, y1);

	if(x1 == x2)
		return fabs(x - x1);

	if(y1 == y2)
	{
		double d1 = fabs(x - x1);
		double d2 = fabs(x - x2);

		return (d1 < d2) ? d1 : d2;
	}

	double m = (y2-y1)/(x2-x1);
	double b = y1 - m*x1;

	double xi = (y - b)/m;

	return fabs(x - xi);
}

//Shortest distance between the point (x, y) and the line segment *defined* by (x1, y1), (x2, y2)
//Note: this is for a *finite* line segment
double DistenceFromFinLineSeg(double x, double y, double x1, double y1, double x2, double y2)
{
	if(x1 == x2 && y1 == y2)
		return distance(x, y, x1, y1);

	if(x1 == x2)
	{
		if((y >= y1 && y <= y2) || (y <= y1 && y >= y2))		
			return fabs(x1 - x);
		if((y > y1 && y1 > y2) || (y < y1 && y1 < y2))
			return distance(x1, y1, x, y);
		if((y > y2 && y2 > y1) || (y < y2 && y2 < y1))
			return distance(x2, y2, x, y);
	}

	//Copy the if block above and change all the "y" for "x" and
	//all the "x" for "y" except in the distance function calls
	if(y1 == y2)
	{
		if((x >= x1 && x <= x2) || (x <= x1 && x >= x2))		
			return fabs(y1 - y);
		if((x > x1 && x1 > x2) || (x < x1 && x1 < x2))
			return distance(x1, y1, x, y);
		if((x > x2 && x2 > x1) || (x < x2 && x2 < x1))
			return distance(x2, y2, x, y);
	}

	double m1 = (y2 - y1)/(x2 - x1);
	double b1 = y1-(m1*x1);
	double m2 = -1.0/m1;
	double b2 = y-(m2*x);

	double xi = (b2 - b1)/(m1 - m2);
	double yi = (b1*m2 - b2*m1)/(m2 - m1);

	if(
	((x1 < x2 && y1 < y2) && (xi >= x1 && x1 <= x2) && (yi >= y1 && yi <= y2)) ||
	((x1 > x2 && y1 > y2) && (xi <= x1 && x1 >= x2) && (yi <= y1 && yi >= y2)) ||
	((x1 < x2 && y1 > y2) && (xi >= x1 && x1 <= x2) && (yi <= y1 && yi >= y2)) ||
	((x1 > x2 && y1 < y2) && (xi <= x1 && x1 >= x2) && (yi >= y1 && yi <= y2))
	)
		return DistenceFromLineWithSlope(m1, b1, x, y);
	else
	{
		double d1 = distance(x, y, x1, y1);
		double d2 = distance(x, y, x2, y2);

		return (d1 < d2) ? d1 : d2;
	}
}

//Shortest distance between the point (x, y) and the line segment *passing* by (x1, y1), (x2, y2)
//Note: use this if you can consider the line segment as infinite
double DistenceFromInfLineSeg(double x, double y, double x1, double y1, double x2, double y2)
{
	//if(x1 == x2 && y1 == y2)
		//Undefined behavior

	if(x1 == x2)
	{
		if(x == x1)
			return 0.0;
		
		return fabs(x- x1);
	}

	//Copy the if block above and change all the "y" for "x" and all the "x" for "y"
	if(y1 == y2)
	{
		if(y == y1)
			return 0.0;
		
		return fabs(y- y1);
	}

	double m = (y2 - y1)/(x2 - x1);
	double b = y1-(m*x1);

	return DistenceFromLineWithSlope(m, b, x, y);

}

//The next three functions are for regular polygons
//By definition, all sides of a regular polygon are equal in length.

//Area of a regular polygon knowing one side length
double AreaRegPolySide(double nbside, double length)
{
	return (length*length*nbside)/(4*tan(M_PI/nbside));
}

//Area of a regular polygon knowing the circumradius
//Note: the circumradius is distance from the center to a vertex
double AreaRegPolyRad(double nbside, double r)
{
	return (r*r*nbside*sin((2*M_PI)/nbside))/2;
}

//Area of a regular polygon knowing the apothem/inradius
//Note: the apothem/inradius is the perpendicular distance from center to a side
double AreaRegPolyApo(double nbside, double a)
{
	return a*a*nbside*tan(M_PI/nbside);
}

//Area of a polygon knowing the coordinates of the vertices
//Note: This works with regular, irregular, convex and 
//concave polygons but not with self-crossing polygons
//Note: The vertices must be in order, either clockwise or 
//counter-clockwise starting at any vertex. If the order is
//clockwise the result is positive and if the order is
//counter-clockwise the result is negative (*-1 to fix)
double AreaPolyCoord(vector<double> &x, vector<double> &y)
{
	double sum = 0.0;

	size_t j = x.size() - 1;
	for(size_t i = 0; i < x.size(); j = i, ++i)
		sum += (x[j] + x[i])*(y[j] - y[i]);

	return sum/2.0;
}

//This function check if a point is in a polygon using the Ray-casting algorithm
//Note: if the point is on the border of the polygon the result is an undefined behavior 
bool pnpoly(vector<double> &xp, vector<double> &yp, double x, double y)
{
	bool c = false;
	size_t j = xp.size() - 1;

	for(size_t i = 0; i < xp.size(); j = i, ++i)
		if(((yp[i] > y) != (yp[j] > y)) && (x < (xp[j] - xp[i]) * (y - yp[i])/(yp[j] - yp[i]) + xp[i]))
			c = !c;

	return c;
}

//This version of the algorithm consider points on the borders as inside the polygon
//Note: you can change all the "return false" for "return true" if you want to exclude the borders 
bool pnpolyb(vector<double> &xp, vector<double> &yp, double x, double y)
{
	bool c = false;
	size_t j = xp.size() - 1;

	for(size_t i = 0; i < xp.size(); j = i, ++i)
	{
		if(xp[i] == x && yp[i] == y)
			return true;
		else if(yp[i] == yp[j])
		{
			if(((x <= xp[j] && x >= xp[i]) || (x <= xp[i] && x >= xp[j])) && y == yp[i])
				return true;
		}
		else if(((yp[i] <= y) && (y < yp[j])) || ((yp[j] <= y) && (y < yp[i])))
		{
			if(x == (xp[j] - xp[i])*(y-yp[i])/(yp[j] - yp[i]) + xp[i])
				return true;
			else if(x < (xp[j] - xp[i])*(y-yp[i])/(yp[j] - yp[i]) + xp[i])
				c = !c;
		}
	}

	return c;
}

//Used in ConvexHull
typedef pair<long, long> point;

//Helper function for ConvexHull
bool cw(const point &a, const point &b, const point &c)
{
	return (b.first - a.first) * (c.second - a.second) - (b.second - a.second) * (c.first - a.first) < 0;
}

//This function compute the convex hull for a set of points
//Note: Convex hull mean the smallest convex set of points in which a given set of points is contained
vector<point> ConvexHull(vector<point> p) 
{
	if(p.size() <= 1)
		return p;

	sort(p.begin(), p.end());

	size_t cnt = 0;
	vector<point> q(p.size() * 2);

	for(size_t i = 0; i < p.size(); q[cnt++] = p[i++])
		for (; cnt >= 2 && !cw(q[cnt - 2], q[cnt - 1], p[i]); --cnt);

	for(ssize_t i = p.size() - 2, t = cnt; i >= 0; q[cnt++] = p[i--])
		for(; cnt > t && !cw(q[cnt - 2], q[cnt - 1], p[i]); --cnt);

	q.resize(cnt - 1 - (q[0] == q[1]));

	return q;
}

//The folowing functions compute the approximate area under the curve (integral)
//Note: Rectangle is less precise than trapezoid, trapezoid is less precise than simpson13 and so on
//Note: Smaller the STEP is, slower the execution will be, in the same way, more precision mean slower execution
#define STEP 0.000001

template<typename T>
T rectangle(T (*f)(T), T a, T b)
{
	T sum = (T)0.0;
	T tmp = (T)STEP/(T)2.0;

	b += numeric_limits<T>::epsilon();
	for(T i = a; i <= b; i += (T)STEP)
		sum += f(i + tmp);

	sum *= (T)STEP;

	return sum;
}

template<typename T>
T trapezoid(T (*f)(T), T a, T b)
{
	T sum = (T)0.0;

	b += numeric_limits<T>::epsilon();
	for(T i = a; i <= b;)
	{
		sum += f(i);
		i += (T)STEP;
		sum += f(i);
	}

	sum *= (T)STEP/(T)2.0;

	return sum;
}

template<typename T>
T simpson13(T (*f)(T), T a, T b)
{
	T sum = (T)0.0;
	T inStep = (T)STEP/(T)2.0;

	b += numeric_limits<T>::epsilon();
	for(T i = a; i <= b;)
	{
		sum += f(i) 
			+ (T)4.0*f(i + inStep);

		i += (T)STEP;
		sum += f(i);
	}

	sum *= inStep/(T)3.0;

	return sum;
}

template<typename T>
T simpson38(T (*f)(T), T a, T b)
{
	T sum = (T)0.0;
	T inStep = (T)STEP/(T)3.0;

	b += numeric_limits<T>::epsilon();
	for(T i = a; i <= b; i += (T)STEP)
		sum += (f(i) + f(i + (T)STEP))
			+ (T)3.0*(f(i + inStep) + f(i + (T)2.0*inStep));

	sum *= (T)3.0*inStep/(T)8.0;

	return sum;
}

template<typename T>
T boole(T (*f)(T), T a, T b)
{
	T sum = (T)0.0;
	T inStep = (T)STEP/(T)4.0;

	b += numeric_limits<T>::epsilon();
	for(T i = a; i <= b; i += (T)STEP)
		sum += (T)7.0*(f(i) + f(i + (T)STEP))
			+ (T)32.0*(f(i + inStep) + f(i + (T)3.0*inStep))
			+ (T)12.0*f(i + (T)2.0*inStep);

	sum *= (T)2.0*inStep/(T)45.0;

	return sum;
}

template<typename T>
T closedNC5(T (*f)(T), T a, T b)
{
	T sum = (T)0.0;
	T inStep = (T)STEP/(T)5.0;

	b += numeric_limits<T>::epsilon();
	for(T i = a; i <= b; i += (T)STEP)
		sum += (T)19.0*(f(i) + f(i + (T)STEP))
			+ (T)75.0*(f(i + inStep) + f(i + (T)4.0*inStep))
			+ (T)50.0*(f(i + (T)2.0*inStep) + f(i + (T)3.0*inStep));

	sum *= (T)5.0*inStep/(T)288.0;

	return sum;
}


template<typename T>
T closedNC6(T (*f)(T), T a, T b)
{
	T sum = (T)0.0;
	T inStep = (T)STEP/(T)6.0;

	b += numeric_limits<T>::epsilon();
	for(T i = a; i <= b; i += (T)STEP)
		sum += (T)41.0*(f(i) + f(i + (T)STEP))
			+ (T)216.0*(f(i + inStep) + f(i + (T)5.0*inStep)) 
			+ (T)27.0*(f(i + (T)2.0*inStep) + f(i + (T)4.0*inStep))
			+ (T)272.0*f(i + (T)3.0*inStep);

	sum *= inStep/(T)140.0;

	return sum;
}

template<typename T>
T closedNC7(T (*f)(T), T a, T b)
{
	T sum = (T)0.0;
	T inStep = (T)STEP/(T)7.0;

	b += numeric_limits<T>::epsilon();
	for(T i = a; i <= b; i += (T)STEP)
		sum += (T)751.0*(f(i) + f(i + (T)STEP))
			+ (T)3577.0*(f(i + inStep) + f(i + (T)6.0*inStep))
			+ (T)1323.0*(f(i + (T)2.0*inStep) + f(i + (T)5.0*inStep))
			+ (T)2989.0*(f(i + (T)3.0*inStep) + f(i + (T)4.0*inStep));

	sum *= (T)7.0*inStep/(T)17280.0;

	return sum;
}

template<typename T>
T closedNC8(T (*f)(T), T a, T b)
{
	T sum = (T)0.0;
	T inStep = (T)STEP/(T)8.0;

	b += numeric_limits<T>::epsilon();
	for(T i = a; i <= b; i += (T)STEP)
		sum += (T)989.0*(f(i) + f(i + (T)STEP))
			+ (T)5888.0*(f(i + inStep) + f(i + (T)7.0*inStep))
			- (T)928.0*(f(i + (T)2.0*inStep) + f(i + (T)6.0*inStep))
			+ (T)10496.0*(f(i + (T)3.0*inStep) + f(i + (T)5.0*inStep))
			- (T)4540.0*f(i + (T)4.0*inStep);

	sum *= (T)4.0*inStep/(T)14175.0;

	return sum;
}

template<typename T>
T closedNC9(T (*f)(T), T a, T b)
{
	T sum = (T)0.0;
	T inStep = (T)STEP/(T)9.0;

	b += numeric_limits<T>::epsilon();
	for(T i = a; i <= b; i += (T)STEP)
		sum += (T)2857.0*(f(i) + f(i + (T)STEP))
			+ (T)15741.0*(f(i + inStep) + f(i + (T)8.0*inStep))
			+ (T)1080.0*(f(i + (T)2.0*inStep) + f(i + (T)7.0*inStep))
			+ (T)19344.0*(f(i + (T)3.0*inStep) + f(i + (T)6.0*inStep))
			+ (T)5778.0*(f(i + (T)4.0*inStep) + f(i + (T)5.0*inStep));

	sum *= (T)9.0*inStep/(T)89600.0;

	return sum;
}

template<typename T>
T closedNC10(T (*f)(T), T a, T b)
{
	T sum = (T)0.0;
	T inStep = (T)STEP/(T)10.0;

	b += numeric_limits<T>::epsilon();
	for(T i = a; i <= b; i += (T)STEP)
		sum += (T)16067.0*(f(i) + f(i + (T)STEP))
			+ (T)106300.0*(f(i + inStep) + f(i + (T)9.0*inStep))
			- (T)48525.0*(f(i + (T)2.0*inStep) + f(i + (T)8.0*inStep))
			+ (T)272400.0*(f(i + (T)3.0*inStep) + f(i + (T)7.0*inStep))
			- (T)260550.0*(f(i + (T)4.0*inStep) + f(i + (T)6.0*inStep))
			+ (T)427368.0*f(i + (T)5.0*inStep);

	sum *= (T)5.0*inStep/(T)299376.0;

	return sum;
}

//This will return a string which represents the sum of the 
//two NON-NEGATIVE integers represented by the strings a and b
string Add( const string& a, const string& b )
{
	string ans;
	int carry = 0;
	int indexA = a.size()-1, indexB = b.size()-1;

	while(indexA >= 0 && indexB >= 0)
	{
		int t = int(a[indexA]) - '0' + b[indexB] - '0' + carry;
		carry = t/10;
		ans = char('0' + t%10) + ans;
		indexA--, indexB--;
	}
	while(indexA >= 0)
	{
		int t = int(a[indexA]) - '0' + carry;
		carry = t/10;
		ans = char('0' + t%10) + ans;
		indexA--;
	}
	while(indexB >= 0)
	{
		int t = int(b[indexB]) - '0' + carry;
		carry = t/10;
		ans = char('0' + t%10) + ans;
		indexB--;
	}
	if(carry)
		ans = char('0' + carry) + ans;

	return ans;
}

//This convert an integer to a string
string itoa(int n)
{
	if(n == 0)
		return "0";

	bool negative = (n < 0);
	string ans;

	if(negative)
		n = -n;

	while(n != 0)
	{
		ans = char('0' + n%10) + ans;
		n /= 10;
	}

	if(negative)
		ans = '-' + ans;

	return ans;
}

//Multiply intergers represented with strings
string Multiply(string a, string b)
{
	if(a.size() < 5 && b.size() < 5)
	{
		int intA = atoi(a.c_str());
		int intB = atoi(b.c_str());

		return itoa(intA*intB);
	}

	while(a.size() < b.size())
		a = '0' + a;
	while(a.size() > b.size())
		b = '0' + b;

	int size = a.size()/2;
	int r = a.size()%2;

	string a1 = Multiply(a.substr(0, size), b.substr(0, size)) + string((size + r)*2, '0');
	string a2 = Multiply(a.substr(0, size), b.substr(size)) + string((size + r), '0');
	string a3 = Multiply(a.substr(size), b.substr(0, size)) + string((size + r), '0');
	string a4 = Multiply(a.substr(size), b.substr(size));

	return Add(Add(a1, a2), Add(a3, a4));
}

//The following pieces of code are important bit twiddling hacks

//Detect if two integers have opposite signs
bool f = ((x ^ y) < 0);

//Determining if an integer is a power of 2
f = (v & (v - 1)) == 0;
//Note: 0 is incorrectly considered a power of 2 here. To remedy this, use:
f = v && !(v & (v - 1));

//Counting bits set (naive way)
for(c = 0; v; v >>= 1)
	c += v & 1;
	
//Counting bits set, Brian Kernighan's way
for(c = 0; v; c++)
	v &= v - 1;
	
//The best method for counting bits in a 32-bit integer v is the following:
v = v - ((v >> 1) & 0x55555555);
v = (v & 0x33333333) + ((v >> 2) & 0x33333333);
c = ((v + (v >> 4) & 0xF0F0F0F) * 0x1010101) >> 24;

//A generalization of the best bit counting method to integers 
//of bit-widths upto 128 (parameterized by type T) is this:
v = v - ((v >> 1) & (T)~(T)0/3);
v = (v & (T)~(T)0/15*3) + ((v >> 2) & (T)~(T)0/15*3);
v = (v + (v >> 4)) & (T)~(T)0/255*15;
c = (T)(v * ((T)~(T)0/255)) >> (sizeof(T) - 1) * 8;

//Computing parity the naive way
unsigned int v;       // word value to compute the parity of
bool parity = false;  // parity will be the parity of v

while(v)
{
	parity = !parity;
	v = v & (v - 1);
}

//Reverse bits the obvious way
unsigned int v;     // input bits to be reversed
unsigned int r = v; // r will be reversed bits of v; first get LSB of v
int s = sizeof(v) * 7; // extra shift needed at end

for(v >>= 1; v; v >>= 1)
{   
	r <<= 1;
	r |= v & 1;
	s--;
}
r <<= s; // shift when v's highest bits are zero

//Finding integer log base 2 of an integer (aka the position of the highest bit set)
//Find the log base 2 of an integer with the MSB N set in O(N) operations (the obvious way)
unsigned int v; // 32-bit word to find the log base 2 of
unsigned int r = 0; // r will be lg(v)

while(v >>= 1) // unroll for more speed...
	r++;

//Count the consecutive zero bits (trailing) on the right linearly
unsigned int v;  // input to count trailing zero bits
int c;  // output: c will count v's trailing zero bits,
        // so if v is 1101000 (base 2), then c will be 3
if(v)
{
	v = (v ^ (v - 1)) >> 1;  // Set v's trailing 0s to 1s and zero rest
	for (c = 0; v; c++)
		v >>= 1;
}
else
	c = 8 * sizeof(v);

//Count the consecutive zero bits (trailing) on the right in parallel
unsigned int v;      // 32-bit word input to count zero bits on right
unsigned int c = 32; // c will be the number of zero bits on the right
v &= -signed(v);
if (v) c--;
if (v & 0x0000FFFF) c -= 16;
if (v & 0x00FF00FF) c -= 8;
if (v & 0x0F0F0F0F) c -= 4;
if (v & 0x33333333) c -= 2;
if (v & 0x55555555) c -= 1;

//Count the consecutive zero bits (trailing) on the right by casting to a float
unsigned int v;            // find the number of trailing zeros in v
int r;                     // the result goes here
float f = (float)(v & -v); // cast the least significant bit in v to a float
r = (*(uint32_t *)&f >> 23) - 0x7f;

//Round up to the next highest power of 2
unsigned int v; // compute the next highest power of 2 of 32-bit v

v--;
v |= v >> 1;
v |= v >> 2;
v |= v >> 4;
v |= v >> 8;
v |= v >> 16;
v++;