/*	pbp.h

	This is the reference C++ library implementation of the
	new Parallel Bit Pattern computation model.

	Copyright (C) 2022 by Henry Dietz, All rights reserved.

	This library is free software available under CC BY 4.0,
	https://creativecommons.org/licenses/by/4.0/

	This version is available at http://aggregate.org/PBP
	along with any relevant documentation, such as the
	reference card and citations to publications.

	This code is distributed in the hope that it will be
	useful, but WITHOUT ANY WARRANTY; without even the
	implied warranty of MERCHANTABILITY or FITNESS FOR A
	PARTICULAR PURPOSE.

	Version history:
	20220704: initial version, use at University of Kentucky

	Known bugs:
	* pint ShR, ShL, Div, and Mod have not been validated
*/

#ifndef	_PBPH__
#define	_PBPH__

#include <iostream>
using namespace std;

#define	VERSION	20220704

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <strings.h>
#include <string.h>
#include <stdint.h>
#include <assert.h>

// Array of Bits constants
#define	AOBREGS		(1024)
#define	AOBWAYS		(8)
#define	AOBBITS		(1<<AOBWAYS)
#define	AOBWORDS	(1<<(AOBWAYS-6))
#define	forAOBWORDS(V)	for (int V=0; V<AOBWORDS; ++V)
#define	ACSIZE		(32 * 1024)	// size of applicative cache

// Regular Expression constants
#define	REREGS		(1024)
#ifndef	REWAYS
#define	REWAYS		(10)
#endif
#define	REBITS		(1<<REWAYS)
#define	REAOBS		(1<<(REWAYS-AOBWAYS))
#define	forREAOBS(V)	for (int V=0; V<REAOBS; ++V)

// Pattern Bit constants
// deliberately skips REREGS, cause that can be answer for First(), Ones()
#define	pbit0		(0)
#define	pbit1		(1)
#define	pbitH(E)	((E)+2)
#define	pbitValid(P)	((P)<REREGS)	// a valid pbit value
#define	pbitVOID	(REREGS+1)
#define	pbitNaN		(REREGS+2)

// Pattern Int constants
#define	PINTBITS	(REWAYS)	// restricted by entanglement support
#define	forprec(V)	for (int V=0; V<prec; ++V)
#define	UPINT		0		// unsigned pint type
#define	SPINT		1		// signed pint type

// Array of Bits structures
typedef	uint16_t AoBreg_t;
typedef	uint64_t AoBword_t;
typedef	struct {
	AoBreg_t op;		// opcode
	AoBreg_t a, b, c;	// c = a op b
}	AC_t;

// Regular Expression structures
#if	(REWAYS > 15)
typedef	uint32_t REreg_t;
#else
typedef	uint16_t REreg_t;
#endif

class AC
{
	// Applicative Cache
	AC_t ac[ACSIZE];
	int spot;
	int hit, miss, rehash;

public:
	AC() {
		bzero(((void *) &ac[0]), sizeof(ac));
		hit = 0;
		miss = 0;
		rehash = 0;
	}

	int hash(AoBreg_t op, AoBreg_t a, AoBreg_t b) {
		// should really pick this by benchmarking...
		register int h = (op * (ACSIZE/256));
		h ^= a;
		h ^= b * (ACSIZE/AOBREGS);
		return(h & (ACSIZE-1));
	}

	int avail(AoBreg_t op, AoBreg_t a, AoBreg_t b) {
		register int h = hash(op, a, b);

		for (;;) {
			// empty entry
			if (!ac[h].op) {
				++miss;
				ac[h].op = op;
				ac[h].a = a;
				ac[h].b = b;
				spot = h;
				return(-1);
			}

			// found it
			if ((ac[h].op == op) &&
			    (ac[h].a == a) &&
			    (ac[h].b == b)) {
				++hit;
				return(ac[h].c);
			}

			++rehash;
			h = ((h + 13) & (ACSIZE-1));
		}
	}

	int assoc(AoBreg_t op, AoBreg_t a, AoBreg_t b) {
		// assoiciative operator, so normalize operands
		return((b < a) ? avail(op, a, b) : avail(op, b, a));
	}

	void Stats() {
		fprintf(stderr, "AC: %d hits, %d misses, %d rehashes\n", hit, miss, rehash);
	}

	AoBreg_t cache(AoBreg_t c) {
		return(ac[spot].c = c);
	}
};


class AoB
{
	AC		ac; // applicative cache for AoB values
	AoBword_t	reg[AOBREGS][AOBWORDS];
	AoBreg_t	index[AOBREGS]; // 0 means empty, >0 means reg[index-1]
	AoBreg_t	regs;
	int		hit, rehash;

public:
	// Helpers
	void copy(const AoBreg_t d, const AoBreg_t s) {
		memcpy(&(reg[d][0]), &(reg[s][0]), (AOBWORDS * sizeof(AoBword_t)));
	}

	int match(const AoBreg_t a) {
		forAOBWORDS(i) {
			if (reg[a][i] != reg[regs][i]) return(0);
		}
		return(1);
	}

	void Stats() {
		ac.Stats();
		fprintf(stderr, "AoB: %d regs, %d hits, %d rehashes\n", regs, hit, rehash);
	}

	void Show(const AoBreg_t a) {
		forAOBWORDS(i) {
			for (int b=0; b<64; ++b) {
				putc(((reg[a][i] & (1LL << b)) ? '1' : '0'), stderr);
			}
			putc('\n', stderr);
		}
	}

	// Constructor
	AoB() {
		// initialization patterns for words
		AoBword_t	inits[] = { 0x0000000000000000LL,
				    0xffffffffffffffffLL,
				    0xaaaaaaaaaaaaaaaaLL,
				    0xccccccccccccccccLL,
				    0xf0f0f0f0f0f0f0f0LL,
				    0xff00ff00ff00ff00LL,
				    0xffff0000ffff0000LL,
				    0xffffffff00000000LL };

		// set-up Stats
		hit = 0;
		rehash = 0;

		// set-up 0, 1, H(0), H(1), H(2), ...
		bzero(&(index[0]), (AOBREGS * sizeof(index[0])));
		regs = 0;
		while (regs < (AOBWAYS+2)) {
			forAOBWORDS(i) {
				reg[regs][i] = inits[(regs>7) ? ((i & (1<<(regs-8))) != 0) : regs];
			}

			uniq(); // making unique entry also increments regs
		}
	}

	// Uniqueness functions
	AoBreg_t hash() {
		// should really pick this by benchmarking...
		AoBreg_t h = 601; // a nice random starting bit pattern
		AoBreg_t *p = ((AoBreg_t *) &(reg[regs][0]));

		for (int i=0; i<(AOBWORDS*(sizeof(AoBword_t)/sizeof(AoBreg_t))); ++i) {
			h = (h << 1) + p[i];
		}

		return((h ^ (h / AOBREGS)) & (AOBREGS-1));
	}

	AoBreg_t uniq() {
		AoBreg_t h = hash();
		for (;;) {
			// empty entry
			if (!index[h]) { index[h] = ++regs; return(regs - 1); }

			// found it
			if (index[h] && match(index[h] - 1)) { ++hit; return(index[h] - 1); }

			++rehash;
			h = ((h + 241) & (AOBREGS-1));
		}
	}

	// Basic operations
	AoBreg_t And(const AoBreg_t a, const AoBreg_t b) {
		// normalize
		register AoBreg_t ta, tb;
		if (a <= b) { ta = a; tb = b; } else { ta = b; tb = a; }

		// check special cases
		if (!ta) return(0);
		if ((ta == 1) || (ta == tb)) return(tb);

		// check associative cache
		int c = ac.avail('&', ta, tb);
		if (c >= 0) return(c);

		// do bitwise operation
		forAOBWORDS(i) reg[regs][i] = reg[ta][i] & reg[tb][i];
		return(ac.cache(uniq()));
	}

	AoBreg_t Or(const AoBreg_t a, const AoBreg_t b) {
		// normalize
		register AoBreg_t ta, tb;
		if (a <= b) { ta = a; tb = b; } else { ta = b; tb = a; }

		// check special cases
		if ((!ta) || (ta == tb)) return(tb);
		if (ta == 1) return(1);

		// check associative cache
		int c = ac.avail('|', ta, tb);
		if (c >= 0) return(c);

		// do bitwise operation
		forAOBWORDS(i) reg[regs][i] = reg[ta][i] | reg[tb][i];
		return(ac.cache(uniq()));
	}

	AoBreg_t Xor(const AoBreg_t a, const AoBreg_t b) {
		// normalize
		register AoBreg_t ta, tb;
		if (a <= b) { ta = a; tb = b; } else { ta = b; tb = a; }

		// check special cases
		if (!ta) return(tb);
		if (ta == tb) return(0);

		// check associative cache
		int c = ac.avail('^', ta, tb);
		if (c >= 0) return(c);

		// do bitwise operation
		forAOBWORDS(i) reg[regs][i] = reg[ta][i] ^ reg[tb][i];
		return(ac.cache(uniq()));
	}

	AoBreg_t Rot(const AoBreg_t a, const AoBreg_t b) {
		// normalize
		register AoBreg_t tb = (b & (AOBBITS-1));

		// check special case
		if (!tb) return(a);

		// check associative cache
		int c = ac.avail('r', a, tb); // doesn't matter that tb isn't a reg
		if (c >= 0) return(c);

		// first, shift by words
		int ws = (tb / (8 * sizeof(AoBword_t)));
		forAOBWORDS(i) {
			reg[regs][i] = reg[a][(i - ws) & ((8 * sizeof(AoBword_t)) - 1)];
		}

		// then by bits if any left
		ws = (tb & ((8 * sizeof(AoBword_t)) - 1));
		if (ws) {
			int rs = ((8 * sizeof(AoBword_t)) - ws);
			AoBword_t lo = (reg[regs][((8 * sizeof(AoBword_t)) - 1)] >> rs);
			forAOBWORDS(i) {
				register AoBword_t hi = (reg[regs][i] >> rs);
				reg[regs][i] = ((reg[regs][i] << ws) | lo);
				lo = hi;
			}
		}

		// cache & return unique
		return(ac.cache(uniq()));
	}

	AoBreg_t Flip(const AoBreg_t a, const AoBreg_t b) {
		// normalize
		register AoBreg_t tb = (b & (AOBBITS-1));

		// check special case
		if (!tb) return(a);

		// check associative cache
		int c = ac.avail('f', a, tb); // doesn't matter that tb isn't a reg
		if (c >= 0) return(c);

		AoBword_t	left[] = { 0x5555555555555555LL,
				   0x3333333333333333LL,
				   0x0f0f0f0f0f0f0f0fLL,
				   0x00ff00ff00ff00ffLL,
				   0x0000ffff0000ffffLL,
				   0x00000000ffffffffLL };
		AoBword_t	right[] = { 0xaaaaaaaaaaaaaaaaLL,
				    0xccccccccccccccccLL,
				    0xf0f0f0f0f0f0f0f0LL,
				    0xff00ff00ff00ff00LL,
				    0xffff0000ffff0000LL,
				    0xffffffff00000000LL };

		// first, copy into work area
		copy(regs, a);

		// do any within-a-word steps
		for (int k=0; k<6; ++k) {
			int bit = (1 << k);
			if (bit & tb) {
				forAOBWORDS(i) {
					reg[regs][i] = (((reg[regs][i] & left[k]) << bit) |
							((reg[regs][i] & right[k]) >> bit));
				}
			}
		}

		// now do across words
		for (int k=0; k<(AOBWAYS-6); ++k) {
			if (((1 << 6) << k) & tb) {
				// copy into regs+1
				copy(regs+1, regs);

				// re-arrange into regs
				forAOBWORDS(i) {
					reg[regs][i] = reg[regs+1][i ^ (1 << k)];
				}
			}
		}

		// cache & return unique
		return(ac.cache(uniq()));
	}

	AoBreg_t Reset(const AoBreg_t a, const AoBreg_t b) {
		// normalize
		register AoBreg_t tb = (b & (AOBBITS-1));

		// check associative cache
		int c = ac.avail('0', a, tb);
		if (c >= 0) return(c);

		// do bitwise operation
		copy(regs, a);
		reg[regs][tb >> 6] &= ~(1LL << (tb & ((1 << 6) - 1)));
		return(ac.cache(uniq()));
	}

	AoBreg_t Set(const AoBreg_t a, const AoBreg_t b) {
		// normalize
		register AoBreg_t tb = (b & (AOBBITS-1));

		// check associative cache
		int c = ac.avail('1', a, tb);
		if (c >= 0) return(c);

		// do bitwise operation
		copy(regs, a);
		reg[regs][tb >> 6] |= (1LL << (tb & ((1 << 6) - 1)));
		return(ac.cache(uniq()));
	}

	AoBreg_t Tog(const AoBreg_t a, const AoBreg_t b) {
		// normalize
		register AoBreg_t tb = (b & (AOBBITS-1));

		// check associative cache
		int c = ac.avail('t', a, tb);
		if (c >= 0) return(c);

		// do bitwise operation
		copy(regs, a);
		reg[regs][tb >> 6] ^= (1LL << (tb & ((1 << 6) - 1)));
		return(ac.cache(uniq()));
	}

	AoBreg_t Dom(const AoBreg_t a, const AoBreg_t b) {
		// normalize
		register AoBreg_t tb = (b & (AOBBITS-1));

		// check associative cache
		int c = ac.avail('t', a, tb);
		if (c >= 0) return(c);

		// do bitwise operation
		copy(regs, a);
		for (int i=0; i<(tb >> 6); ++i) reg[regs][i] = ~reg[regs][i];
		reg[regs][tb >> 6] ^= ((2LL << (tb & ((1 << 6) - 1))) - 1LL);
		return(ac.cache(uniq()));
	}

	AoBreg_t Meas(const AoBreg_t a, const AoBreg_t b) {
		// normalize
		register AoBreg_t tb = (b & (AOBBITS-1));
		// don't bother caching, result is 0 or 1
		return((reg[a][tb >> 6] & (1LL << (tb & ((1 << 6) - 1)))) != 0LL);
	}

	AoBreg_t Any(const AoBreg_t a) {
		// could check for Ones(a) == 0, but this is faster
		forAOBWORDS(i) {
			if (reg[regs][i]) return(1);
		}
		return(0);
	}

	AoBreg_t All(const AoBreg_t a) {
		// could check for Ones(a) == AOBBITS, but this is faster
		forAOBWORDS(i) {
			if (reg[regs][i] != 1) return(0);
		}
		return(1);
	}

	AoBreg_t First(const AoBreg_t a) {
		// check special cases... initialized constants
		if (a < (AOBWAYS+2)) {
			if (a == 0) return(1 << AOBWAYS);
			if (a == 1) return(0);
			return(1 << (a - 2));
		}

		// don't cache, result is an integer
		forAOBWORDS(i) {
			register AoBword_t t = reg[a][i];
			if (t) {
				// count trailing zeros
				AoBreg_t pos = (i * (8 * sizeof(AoBword_t)));
				if (!(t & 0x00000000ffffffffLL)) { pos += 32; t >>= 32; }
				if (!(t & 0x000000000000ffffLL)) { pos += 16; t >>= 16; }
				if (!(t & 0x00000000000000ffLL)) { pos += 8; t >>= 8; }
				if (!(t & 0x000000000000000fLL)) { pos += 4; t >>= 4; }
				if (!(t & 0x0000000000000003LL)) { pos += 2; t >>= 2; }
				return((t & 0x0000000000000001LL) ? pos : (pos + 1));
			}
		}
		return(1 << AOBWAYS);
	}

	AoBreg_t Ones(const AoBreg_t a) {
		// check special cases... initialized constants
		if (a < (AOBWAYS+2)) {
			if (a == 0) return(0);
			if (a == 1) return(1 << AOBWAYS);
			return(1 << (AOBWAYS - 1));
		}
		
		// don't cache, result is an integer
		AoBreg_t pop = 0;
		forAOBWORDS(i) {
			register AoBword_t t = reg[a][i];
			if (t) {
				if (~t) {
					// count trailing zeros
					t = (t & 0x5555555555555555LL) + ((t >> 1) & 0x5555555555555555LL);
					t = (t & 0x3333333333333333LL) + ((t >> 2) & 0x3333333333333333LL);
					t = (t & 0x0f0f0f0f0f0f0f0fLL) + ((t >> 4) & 0x0f0f0f0f0f0f0f0fLL);
					t = (t & 0x00ff00ff00ff00ffLL) + ((t >> 8) & 0x00ff00ff00ff00ffLL);
					t = (t & 0x0000ffff0000ffffLL) + ((t >> 16) & 0x0000ffff0000ffffLL);
					t = (t & 0x00000000ffffffffLL) + ((t >> 32) & 0x00000000ffffffffLL);
					pop += t;
				} else	pop += 64;
			}
		}
		return(pop);
	}
};

class RE
{
	AoB		aob; // AoB and applicative cache
	AoBreg_t	reg[REREGS][REAOBS];
	REreg_t		index[AOBREGS]; // 0 means empty, >0 means reg[index-1]
	REreg_t		regs;
	int		hit, rehash;

public:
	// Helpers
	void copy(const REreg_t d, const REreg_t s) {
		memcpy(&(reg[d][0]), &(reg[s][0]), (REAOBS * sizeof(AoBreg_t)));
	}

	int match(const REreg_t a) {
		forREAOBS(i) {
			if (reg[a][i] != reg[regs][i]) return(0);
		}
		return(1);
	}

	void Stats() {
		aob.Stats();
		fprintf(stderr, "RE: %d regs, %d hits, %d rehashes\n", regs, hit, rehash);
	}

	void Show(const REreg_t a) {
		forREAOBS(i) {
			aob.Show(reg[a][i]);
		}
	}

	// Constructor
	RE() {
		// initialization patterns for words
		AoBword_t	inits[] = { 0x0000000000000000LL,
				    0xffffffffffffffffLL,
				    0xaaaaaaaaaaaaaaaaLL,
				    0xccccccccccccccccLL,
				    0xf0f0f0f0f0f0f0f0LL,
				    0xff00ff00ff00ff00LL,
				    0xffff0000ffff0000LL,
				    0xffffffff00000000LL };

		// set-up Stats
		hit = 0;
		rehash = 0;

		// set-up 0, 1, H(0), H(1), H(2), ...
		bzero(&(index[0]), (REREGS * sizeof(index[0])));
		regs = 0;
		while (regs < (REWAYS+2)) {
			forREAOBS(i) {
				reg[regs][i] = ((regs < (AOBWAYS+2)) ? regs : ((i & (1<<(regs-(AOBWAYS+2)))) != 0));
			}

			uniq(); // making unique entry also increments regs
		}
	}

	// Uniqueness functions
	REreg_t hash() {
		// should really pick this by benchmarking...
		REreg_t h = 42; // a nice random starting bit pattern

		for (int i=0; i<REAOBS; ++i) {
			h = (h << 3) | (h & 7);
			h ^= reg[regs][i];
		}

		return((h ^ (h / REREGS)) & (REREGS-1));
	}

	REreg_t uniq() {
		REreg_t h = hash();
		for (;;) {
			// empty entry
			if (!index[h]) { index[h] = ++regs; return(regs - 1); }

			// found it
			if (index[h] && match(index[h] - 1)) { ++hit; return(index[h] - 1); }

			++rehash;
			h = ((h + 347) & (REREGS-1));
		}
	}

	// Basic operations
	REreg_t And(const REreg_t a, const REreg_t b) {
		// normalize
		register REreg_t ta, tb;
		if (a <= b) { ta = a; tb = b; } else { ta = b; tb = a; }

		// check special cases
		if (!ta) return(0);
		if ((ta == 1) || (ta == tb)) return(tb);

		// do AoB operation
		forREAOBS(i) reg[regs][i] = aob.And(reg[ta][i], reg[tb][i]);
		return(uniq());
	}

	REreg_t Or(const REreg_t a, const REreg_t b) {
		// normalize
		register REreg_t ta, tb;
		if (a <= b) { ta = a; tb = b; } else { ta = b; tb = a; }

		// check special cases
		if ((!ta) || (ta == tb)) return(tb);
		if (ta == 1) return(1);

		// do AoB operation
		forREAOBS(i) reg[regs][i] = aob.Or(reg[ta][i],  reg[tb][i]);
		return(uniq());
	}

	REreg_t Xor(const REreg_t a, const REreg_t b) {
		// normalize
		register REreg_t ta, tb;
		if (a <= b) { ta = a; tb = b; } else { ta = b; tb = a; }

		// check special cases
		if (!ta) return(tb);
		if (ta == tb) return(0);

		// do AoB operation
		forREAOBS(i) reg[regs][i] = aob.Xor(reg[ta][i], reg[tb][i]);
		return(uniq());
	}

	REreg_t Rot(const REreg_t a, const REreg_t b) {
		// normalize
		register REreg_t tb = (b & (REBITS-1));

		// check special case
		if (!tb) return(a);

		// first, shift by AoBs
		int ws = (tb / AOBBITS);
		forAOBWORDS(i) {
			reg[regs][i] = reg[a][(i - ws) & (REAOBS-1)];
		}

		// then by bits if any left
		ws = (tb & (AOBBITS - 1));
		if (ws) {
			AoBreg_t lomask = aob.Dom(0, ws-1);
			AoBreg_t himask = aob.Xor(1, ws-1);
			forREAOBS(i) {
				AoBreg_t hi = aob.And(himask, aob.Rot(reg[regs][i], ws));
				AoBreg_t lo = aob.And(lomask, aob.Rot(reg[regs][(i - 1) & (AOBBITS-1)], (AOBBITS-ws)));
				reg[regs+1][i] = aob.Or(hi, lo);
			}

			copy(regs, regs+1);
		}

		// return unique
		return(uniq());
	}

	REreg_t Flip(const REreg_t a, const REreg_t b) {
		// normalize
		register REreg_t tb = (b & (REBITS-1));

		// check special case
		if (!tb) return(a);

		// first, copy into work area
		copy(regs, a);

		// do any within-an-AoB steps
		forREAOBS(i) {
			reg[regs][i] = aob.Flip(reg[a][i], (tb & (AOBBITS-1)));
		}

		// now do across AoBs
		for (int k=0; k<(REWAYS-AOBWAYS); ++k) {
			if (((1 << AOBWAYS) << k) & tb) {
				// copy into regs+1
				copy(regs+1, regs);

				// re-arrange into regs
				forAOBWORDS(i) {
					reg[regs][i] = reg[regs+1][i ^ (1 << k)];
				}
			}
		}

		// return unique
		return(uniq());
	}

	REreg_t Reset(const REreg_t a, const REreg_t b) {
		// normalize
		register REreg_t tb = (b & (REBITS-1));

		// do bitwise operation
		copy(regs, a);
		reg[regs][tb >> AOBWAYS] = aob.Reset(reg[regs][tb >> AOBWAYS], (tb & (AOBBITS-1)));
		return(uniq());
	}

	REreg_t Set(const REreg_t a, const REreg_t b) {
		// normalize
		register REreg_t tb = (b & (REBITS-1));

		// do bitwise operation
		copy(regs, a);
		reg[regs][tb >> AOBWAYS] = aob.Set(reg[regs][tb >> AOBWAYS], (tb & (AOBBITS-1)));
		return(uniq());
	}

	REreg_t Tog(const REreg_t a, const REreg_t b) {
		// normalize
		register REreg_t tb = (b & (REBITS-1));

		// do bitwise operation
		copy(regs, a);
		reg[regs][tb >> AOBWAYS] = aob.Tog(reg[regs][tb >> AOBWAYS], (tb & (AOBBITS-1)));
		return(uniq());
	}

	REreg_t Dom(const REreg_t a, const REreg_t b) {
		// normalize
		register REreg_t tb = (b & (REBITS-1));

		// do bitwise operation
		copy(regs, a);
		for (int i=0; i<(tb >> AOBWAYS); ++i) reg[regs][i] = aob.Xor(1, reg[regs][i]);
		reg[regs][tb >> AOBWAYS] = aob.Dom(reg[regs][tb >> AOBWAYS], (tb & (AOBBITS-1)));
		return(uniq());
	}

	REreg_t Meas(const REreg_t a, const REreg_t b) {
		// normalize
		register REreg_t tb = (b & (REBITS-1));
		// result is 0 or 1; no need to do uniq()
		return(aob.Meas(reg[a][tb >> AOBWAYS], (tb & (AOBBITS-1))));
	}

	REreg_t Any(const REreg_t a) {
		// result is 0 or 1; no need to do uniq()
		forREAOBS(i) {
			if (aob.Any(reg[a][i])) return(1);
		}
		return(0);
	}

	REreg_t All(const REreg_t a) {
		// result is 0 or 1; no need to do uniq()
		forREAOBS(i) {
			if (!aob.All(reg[a][i])) return(0);
		}
		return(1);
	}

	REreg_t First(const REreg_t a) {
		// check special cases... initialized constants
		if (a < (REWAYS+2)) {
			if (a == 0) return(1 << REWAYS);
			if (a == 1) return(0);
			return(1 << (a - 2));
		}

		// result is an integer
		register REreg_t r = 0;
		forREAOBS(i) {
			register REreg_t t = aob.First(reg[a][i]);
			r += t;
			if (t < AOBBITS) {
				return(r);
			}
		}
		return(r);
	}

	REreg_t Ones(const REreg_t a) {
		// check special cases... initialized constants
		if (a < (REWAYS+2)) {
			if (a == 0) return(0);
			if (a == 1) return(1 << REWAYS);
			return(1 << (REWAYS - 1));
		}
		
		// result is an integer
		REreg_t pop = 0;
		forREAOBS(i) {
			pop += aob.Ones(reg[a][i]);
		}
		return(pop);
	}
};

RE re;	// instantiate everything

class pbit
{
public:
	REreg_t v;  // value

	// Constructor
	pbit(int value=pbitVOID) {
		v = (((value < 0) || (value > pbitVOID)) ? pbitNaN : value);
	}

	// Helpers
	void Show() const {
		if (v < pbitVOID) re.Show(v);
		else if (v == pbitVOID) fprintf(stderr, "pbitVOID\n");
		else fprintf(stderr, "pbitNaN\n");
	}

	int Valid() const {
		return((v >= 0) && (v < pbitVOID));
	}

	// Operators
	void operator=(const pbit& a) {
		v = (a.Valid() ? a.v : pbitNaN);
	}

	pbit And(const pbit& a) const {
		return((a.Valid() && Valid()) ? pbit(re.And(a.v, v)) : pbit(pbitNaN));
	}
	pbit operator&(const pbit& a) const { return(And(a)); }
	pbit operator*(const pbit& a) const { return(And(a)); }

	pbit Or(const pbit& a) const {
		return((a.Valid() && Valid()) ? pbit(re.Or(a.v, v)) : pbit(pbitNaN));
	}
	pbit operator|(const pbit& a) const { return(Or(a)); }

	pbit Xor(const pbit& a) const {
		return((a.Valid() && Valid()) ? pbit(re.Xor(a.v, v)) : pbit(pbitNaN));
	}
	pbit operator^(const pbit& a) const { return(Xor(a)); }
	pbit operator+(const pbit& a) const { return(Xor(a)); }

	pbit Not() const { return(Valid() ? pbit(re.Xor(1, v)) : pbit(pbitNaN)); }
	pbit operator~() const { return(Not()); }

	pbit LNot() const { return(Valid() ? pbit(v == 0) : pbit(pbitNaN)); }
	pbit operator!() const { return(LNot()); }

	pbit Sub(const pbit& a) const {
		// this - a means this & ~a
		return((a.Valid() && Valid()) ? pbit(re.And(v, re.Xor(1, a.v))) : pbit(pbitNaN));
	}
	pbit operator-(const pbit& a) const { return(Sub(a)); }

	pbit LT(const pbit& a) const {
		// this < a means a & ~this
		return((a.Valid() && Valid()) ? pbit(re.And(a.v, re.Xor(1, v))) : pbit(pbitNaN));
	}
	pbit operator<(const pbit& a) const { return(LT(a)); }

	pbit GT(const pbit& a) const { return(a < pbit(v)); }
	pbit operator>(const pbit& a) const { return(GT(a)); }

	pbit LE(const pbit& a) const {
		// this <= a means ~(this > a)
		return((a.Valid() && Valid()) ? pbit(re.Xor(1, GT(a).v)) : pbit(pbitNaN));
	}
	pbit operator<=(const pbit& a) const { return(LE(a)); }

	pbit GE(const pbit& a) const { return(a <= pbit(v)); }
	pbit operator>=(const pbit& a) const { return(GE(a)); }

	pbit NE(const pbit& a) const { return(Xor(a)); }
	pbit operator!=(const pbit& a) const { return(NE(a)); }

	pbit EQ(const pbit& a) const {
		// this == a means ~(this != a)
		return((a.Valid() && Valid()) ? pbit(re.Xor(1, NE(a).v)) : pbit(pbitNaN));
	}
	pbit operator==(const pbit& a) const { return(EQ(a)); }

	pbit Rot(const int a) const {
		return(Valid() ? pbit(re.Rot(v, a)) : pbit(pbitNaN));
	}
	pbit operator<<(const int a) const { return(Rot(a)); }

	pbit Flip(const int a) const {
		return(Valid() ? pbit(re.Flip(v, a)) : pbit(pbitNaN));
	}

	pbit Reset(const int a) const { return(Valid() ? pbit(re.Reset(v, a)) : pbit(pbitNaN)); }

	pbit Set(const int a) const { return(Valid() ? pbit(re.Set(v, a)) : pbit(pbitNaN)); }

	pbit Dom(const int a) const { return(Valid() ? pbit(re.Dom(v, a)) : pbit(pbitNaN)); }

	pbit Meas(const int a) const { return(Valid() ? pbit(re.Meas(v, a)) : pbit(pbitNaN)); }

	pbit Any() const { return(Valid() ? pbit(re.Any(v)) : pbit(pbitNaN)); }

	pbit All() const { return(Valid() ? pbit(re.All(v)) : pbit(pbitNaN)); }

	REreg_t First() const { return(Valid() ? re.First(v) : pbitNaN); }

	REreg_t Ones() const { return(Valid() ? re.Ones(v) : pbitNaN); }

	pbit Mux(const pbit& a, const pbit& b) const {
		// a few simplifications
		if (v == 0) return(b);
		if (v == 1) return(a);
		if (a.v == b.v) return(a);
		if ((!Valid()) || (!a.Valid()) || (!b.Valid())) return(pbitNaN);

		// do the bit ops
		return((b & *this) | (a & ~*this));
	}

};

// quantum-style pbit operators
void NOT(pbit& q) { q = pbit(1) ^ q; }
void CNOT(pbit& c, pbit& t) { t = c ^ t; }
void SWAP(pbit& i0, pbit& i1) { pbit t = i0; i0 = i1; i1 = t; }
void CCNOT(pbit& a, pbit& b, pbit& c) { c = c ^ (a & b); }
void CSWAP(pbit& c, pbit& i0, pbit& i1) { pbit t = c.Mux(i0,i1); i1 = c.Mux(i1,i0); i0 = t; }
void H(pbit& q, int c) { q = q ^ pbit(c + 2); }
int MEASat = 0;
void SETMEAS(int m) { MEASat = (m & (REBITS - 1)); }
void SETMEAS() { SETMEAS(rand()); }
int MEAS(pbit& q) { q = q.Meas(MEASat); return(q.v); }

	

class pint
{
	uint8_t	type;		// 0 unsigned, 1 signed
	uint8_t	prec;		// precision in pbits
	pbit	bit[PINTBITS];	// actual pbits

public:
	// constructors
	pint() {
		// default to 1 unsigned pbit that's NaN
		type = UPINT;
		prec = 1;
		bit[0] = pbit(pbitNaN);
	}

	pint(const int konst) {
		// auto sizes to konst, also determines sign from konst
		type = (konst < 0);
		prec = PINTBITS;

		// install the bits
		forprec(i) {
			bit[i] = pbit(((1 << i) & konst) ? 1 : 0);
		}

		// get rid of any unnecessary bits
		*this = Minimize();
	}

	pint(const int konst, const int bits) {
		// make an integer with constant value konst
		// signed if bits is negative
		if (bits >= 0) {
			type = 0; // unsigned
			prec = bits;
		} else {
			type = 1; // signed
			prec = -bits;
		}

		// install bits
		forprec(i) {
			bit[i] = pbit(((1 << i) & konst) ? 1 : 0);
		}
	}

	// Helpers
	void Summary() const {
		// summarize this pint
		fprintf(stderr, "%s%d_t {", (type ? "int" : "uint"), prec);
		forprec (i) {
			fprintf(stderr, "%d%s", bit[i].v, ((i==(prec-1)) ? "}\n" : ","));
		}
	}

	void Show() const {
		// show bits in this pint
		Summary();
		for (int i=0; i<prec; ++i) bit[i].Show();
	}

	void Enumerate() const {
		// enumerate integer values in this pint
		const char *f = ((type == SPINT) ? "%d%s" : "%u%s");
		for (int i=0; i<REBITS; ++i) {
			fprintf(stderr, f, Meas(i), (((i%8)==7) ? "\n" : "\t"));
		}
	}

	int Valid() const {
		// check to see if this pbit is valid
		if (type > 1) return(0);
		if ((prec < 1) || (prec >= PINTBITS)) return(0);
		forprec (i) {
			if (!bit[i].Valid()) return(0);
		}
		return(1);
	}

	pint Minimize() const {
		// make value as small as possible
		pint r = *this;

		if (r.type) {
			// signed
			for (int i=r.prec-1; i>0; --i) {
				if (r.bit[i].v != r.bit[i-1].v) return(r);
				--r.prec;
			}
		} else {
			// unsigned
			for (int i=r.prec-1; i>0; --i) {
				if (r.bit[i].v) return(r);
				--r.prec;
			}
		}
		return(r);
	}

	pint Extend(const int p) const {
		// Extend to precision p, truncate if already bigger
		int pnew = ((p > PINTBITS) ? PINTBITS : p);
		pint r = *this;

		// Signed/unsigned extend
		if (r.prec < pnew) {
			for (int i=r.prec; i<pnew; ++i) {
				r.bit[i].v = (r.type ? r.bit[i-1].v : 0);
			}
		}
		r.prec = pnew;
		return(r);
	}

	pint Promote(const pint& b) const {
		// promote this to covering type for this, b
		register int p = ((prec < b.prec) ? b.prec : prec);
		pint r = *this;

		// if the signs differ and the signed one isn't bigger
		// than the unsigned one, add a bit and make both signed
		if (((type && (!b.type)) && (prec <= b.prec)) ||
		    ((b.type && (!type)) && (b.prec <= prec))) {
			++p;
		}

		// Extend precision
		r = r.Extend(p);

		// Correct sign
		r.type = (type | b.type);
		return(r);
	}

	pint Logic() const {
		// reduce to a logical value, i.e., a single pbit
		pint r = *this;
		for (int i=1; i<prec; ++i) {
			r.bit[0] = r.bit[0] | r.bit[i];
		}
		r.type = 0; // always unsigned
		r.prec = 1; // single bit
		return(r);
	}

	// Operators
	void operator=(const pint& a) {
		type = a.type;
		prec = a.prec;
		forprec (i) bit[i] = a.bit[i];
	}
	void operator=(const int a) {
		*this = pint(a);
	}

	pint And(pint b) const {
		pint r = Promote(b);
		b = b.Promote(r);
		for (int i=0; i<r.prec; ++i) {
			r.bit[i] = r.bit[i] & b.bit[i];
		}
		r = r.Minimize();
		return(r);
	}
	pint operator&(const pint& b) const { return(And(b)); }
	pint operator&=(pint b) { *this = And(b); return(*this); }

	pint LAnd(pint b) const {
		pint r = Logic();
		b = b.Logic();
		r.bit[0] = r.bit[0] & b.bit[0];
		return(r);
	}
	pint operator&&(const pint& b) const { return(LAnd(b)); }

	pint Or(pint b) const {
		pint r = Promote(b);
		b = b.Promote(r);
		for (int i=0; i<r.prec; ++i) {
			r.bit[i] = r.bit[i] | b.bit[i];
		}
		r = r.Minimize();
		return(r);
	}
	pint operator|(const pint& b) const { return(Or(b)); }
	pint operator|=(pint b) { *this = Or(b); return(*this); }

	pint LOr(pint b) const {
		pint r = Logic();
		b = b.Logic();
		r.bit[0] = r.bit[0] | b.bit[0];
		return(r);
	}
	pint operator||(const pint& b) const { return(LOr(b)); }

	pint Xor(pint b) const {
		pint r = Promote(b);
		b = b.Promote(r);
		for (int i=0; i<r.prec; ++i) {
			r.bit[i] = r.bit[i] ^ b.bit[i];
		}
		r = r.Minimize();
		return(r);
	}
	pint operator^(const pint& b) const { return(Xor(b)); }
	pint operator^=(pint b) { *this = Xor(b); return(*this); }

	pint LXor(pint b) const {
		pint r = Logic();
		b = b.Logic();
		r.bit[0] = r.bit[0] ^ b.bit[0];
		return(r);
	}
	// oddly, there is no C++ ^^ operator...

	pint Not() const {
		pint r = *this;
		forprec(i) r.bit[i] = ~bit[i];
		return(r);
	}
	pint operator~() const { return(Not()); }

	pint LNot() const {
		pint r = *this;
		r = r.Logic();
		r.bit[0] = ~r.bit[0];
		return(r);
	}
	pint operator!() const { return(LNot()); }

	pint EQ(pint b) const {
		pint r = Xor(b);
		r = r.Logic();
		r.bit[0] = ~r.bit[0];
		return(r);
	}
	pint operator==(const pint& b) const { return(EQ(b)); }

	pint NE(pint b) const {
		pint r = Xor(b);
		r = r.Logic();
		return(r);
	}
	pint operator!=(const pint& b) const { return(NE(b)); }

	pint Rot(const int b) const {
		pint r = *this;
		forprec(i) r.bit[i] = bit[i].Rot(b);
		return(r);
	}

	pint Flip(const int b) const {
		pint r = *this;
		forprec(i) r.bit[i] = bit[i].Flip(b);
		return(r);
	}

	pint Reset(const int b) const {
		pint r = *this;
		forprec(i) r.bit[i] = bit[i].Reset(b);
		r = r.Minimize();
		return(r);
	}

	pint Set(const int b) const {
		pint r = *this;
		forprec(i) r.bit[i] = bit[i].Set(b);
		r = r.Minimize();
		return(r);
	}

	pint Dom(const int b) const {
		pint r = *this;
		forprec(i) r.bit[i] = bit[i].Dom(b);
		return(r);
	}

	int Meas(const int b) const {
		// note that Measurement always produces 0 or 1
		// bit values, never superpositions, which is why
		// this can return an ordinary int value
		register int v = 0;
		forprec(i) if ((bit[i].Meas(b)).v) v |= (1<<i);
		if (type == SPINT) {
			// perform sign extension?
			if (v & (1<<(prec-1))) {
				// yes, negative value
				v |= ~((1<<(prec-1))-1);
			}
		}
		return(v);
	}
	int Meas() const {
		// Quantum-like random measurement
		return(Meas(rand() & (REBITS - 1)));
	}

	REreg_t First() const {
		// first is first 1 in any pbit
		register REreg_t first = pbitNaN;
		forprec(i) {
			register REreg_t t = bit[i].First();
			if (t < first) first = t;
		}
		return(first);
	}

	REreg_t Ones() const {
		// ones is number of non-zero multi-bit values
		pint r = *this;
		r = r.Logic();
		return(bit[0].Ones());
	}

	pint GT(pint b) const {
		pbit gt(0);
		pbit ep(1);
		pint r = *this;
		r = r.Promote(b);
		b = b.Promote(r);

		if (type) {
			// signed > flips bits if they were negative
			gt = (~r.bit[prec-1]) & b.bit[prec-1];
			ep = ~(r.bit[prec-1] ^ b.bit[prec-1]);
			pbit s = r.bit[prec-1] & b.bit[prec-1];
			for (int i=0; i<r.prec; ++i) {
				r.bit[i] = s ^ r.bit[i];
				b.bit[i] = s ^ b.bit[i];
			}
		}

		// unsigned > looks for msb that differs
		for (int i=r.prec-1; i>=0; --i) {
			gt = (gt | (ep & (r.bit[i] & ~b.bit[i])));
			ep = (ep & ~(r.bit[i] ^ b.bit[i]));
		}

		// always returns 1 bit unsigned
		r.bit[0] = gt;
		r.prec = 1;
		r.type = UPINT;
		return(r);
	}
	pint operator>(const pint& b) const { return(GT(b)); }

	pint LT(pint b) const {
		b = (b > *this);
		return(b);
	}
	pint operator<(const pint& b) const { return(LT(b)); }

	pint GE(pint b) const {
		b = ~(*this < b);
		return(b);
	}
	pint operator>=(const pint& b) const { return(GE(b)); }

	pint LE(pint b) const {
		b = ~(*this > b);
		return(b);
	}
	pint operator<=(const pint& b) const { return(LE(b)); }

	pint Min(pint b) const {
		pint r = *this;
		r = r.Promote(b);
		b = b.Promote(r);
		pint gt = (r > b);
		for (int i=0; i<b.prec; ++i) {
			r.bit[i] = gt.bit[0].Mux(r.bit[i], b.bit[i]);
		}
		return(r.Minimize());
	}

	pint Max(pint b) const {
		pint r = *this;
		r = r.Promote(b);
		b = b.Promote(r);
		pint gt = (r > b);
		for (int i=0; i<b.prec; ++i) {
			r.bit[i] = gt.bit[0].Mux(b.bit[i], r.bit[i]);
		}
		return(r.Minimize());
	}

	pint ShR(pint b) const {
		// shift right, treats shift amount as unsigned
		pint r = *this;
		for (int i=(b.prec-1); i>=0; --i) {
			forprec(j) {
				pbit x;

				if ((j + (1 << i)) >= r.prec) {
					x = (type ? r.bit[prec-1] : pbit(0));
				} else {
					x = r.bit[j + (1 << i)];
				}

				r.bit[j] = b.bit[i].Mux(x, r.bit[j]);
			}
		}
		return(r.Minimize());
	}
	pint operator>>(pint b) const { return(ShR(b)); }
	pint operator>>=(pint b) { *this = ShR(b); return(*this); }

	pint ShL(pint b) const {
		// shift left, treats shift amount as unsigned
		pint r = *this;
		r = r.Extend(r.prec + (1 << b.prec) - 1);
		for (int i=0; i<b.prec; ++i) {
			for (int j=prec-1; j>=0; --j) {
				r.bit[j] = b.bit[i].Mux(((j >= (1<<i)) ? bit[j-(1<<i)] : pbit(0)), bit[j]);
			}
		}
		return(r.Minimize());
	}
	pint operator<<(pint b) const { ShL(b); return(*this); }
	pint operator<<=(pint b) { *this = ShL(b); return(*this); }

	pint Add(pint b) const {
		pint r = *this;
		r = r.Promote(b);
		r = r.Extend(r.prec + 1); // k bits + k bits = k+1 bits
		b = b.Promote(r);

		// ripple carry adder
		pbit c = pbit(0);
		for (int i=0; i<b.prec; ++i) {
			pbit x = r.bit[i] ^ b.bit[i];
			pbit n = r.bit[i] & b.bit[i];
			r.bit[i] = x ^ c;
			c = n | (x & c);
		}
		return(r.Minimize());
	}
	pint operator+(pint b) const { return(Add(b)); }
	pint operator+=(pint b) { *this = Add(b); return(*this); }
	// prefix ++
	pint& operator++() { *this = Add(pint(1)); return(*this); }
	// postfix ++ note that int arg is a strange C++ ism
	pint operator++(int) { pint r = *this; *this = Add(pint(1)); return(r); }

	pint Sub(pint b) const {
		pint r = *this;
		r = r.Promote(b);
		r = r.Extend(r.prec + 1); // k bits + k bits = k+1 bits
		b = b.Promote(r);

		// ripple carry subtract
		pbit c = pbit(1);
		for (int i=0; i<b.prec; ++i) {
			pbit n = ~b.bit[i];
			pbit x = r.bit[i] ^ n;
			n = r.bit[i] & n;
			r.bit[i] = x ^ c;
			c = n | (x & c);
		}
		return(r.Minimize());
	}
	pint operator-(pint b) const { return(Sub(b)); }
	pint operator-=(pint b) { *this = Sub(b); return(*this); }
	// prefix --
	pint& operator--() { *this = Sub(pint(1)); return(*this); }
	// postfix -- note that int arg is a strange C++ ism
	pint operator--(int) { pint r = *this; *this = Sub(pint(1)); return(r); }

	pint Neg() const {
		pint r = pint(0).Sub(*this);
		return(r.Minimize());
	}
	pint operator-() const { return(Neg()); }

	pint Mul(pint b) const {
		// multiply
		pint v(0);
		pint one(1);

		// faster if first arg has fewer bits
		forprec(i) {
			pint c(0,PINTBITS);
			for (int j=0; ((j<b.prec)&&((i+j)<PINTBITS)); ++j) {
				c.bit[i+j] = b.bit[j] & bit[i];
			}
			c = c.Minimize();
			v = v + c;
		}

		return(v.Minimize());
	}
	pint operator*(pint b) const { return(Mul(b)); }
	pint operator*=(pint b) { *this = Mul(b); return(*this); }

	pint Mul(pint b, const int bits) const {
		// multiply, but only get least significant bits
		pint v(0);
		pint one(1);
		if (b.prec > bits) b.prec = bits;
		int rounds = ((prec > bits) ? bits : prec);

		// faster if first arg has fewer bits
		for (int i=0; i<rounds; ++i) {
			pint c = b;

			for (int j=0; j<c.prec; ++j) {
				c.bit[j] = c.bit[j] & bit[i];
			}
			v = v + c;
			b = b << one;

			// discard any extra bits
			if (v.prec > bits) v.prec = bits;
			if (b.prec > bits) b.prec = bits;
		}

		return(v.Minimize());
	}

	pint Abs() const {
		// absolute value (trivial if unsigned)
		if (type == UPINT) {
			return(*this);
		}

		pint s(0);
		s.bit[0] = bit[prec-1];
		s = (s & -*this) | ((~s) & *this);
		s.type = UPINT; // result of Abs() is unsigned
		return(s.Minimize());
	}

	pint Signed() const {
		// force signed; preserves bits, not value
		pint r = *this;
		r.type = SPINT;
		return(r);
	}

	pint UnSigned() const {
		// force unsigned
		pint r = *this;
		r.type = UPINT;
		return(r);
	}

	pint Div(pint b) const {
		// divide
		pint s(0);
		pint a = *this;

		// signed?
		if (a.type || b.type) {
			// save the result sign
			a = a.Promote(b);
			b = b.Promote(a);
			b.bit[0] = a.bit[prec-1] ^ b.bit[prec-1];
			a = a.Abs();
			b = b.Abs();
		}

		// do the unsigned divide
		a = a.Promote(b);
		b = b.Promote(a);
		pint rem(a.prec, 0);
		pint quo = rem;

		// classic shift & subtract loop
		for (int i=a.prec-1; i>=0; --i) {
			// shift left 1 bit
			for (int j=a.prec-1; j>0; --j) {
				rem.bit[j] = rem.bit[j-1];
			}
			rem.bit[0] = a.bit[i];

			pint ge = (rem >= b);
			rem = (ge & (rem - b)) | ((!ge) & rem);
			rem = rem.Extend(a.prec);
			quo.bit[i] = ge.bit[0];
		}

		// fix sign if possibly negative
		if (s.bit[0].v != 0) {
			// make it signed again
			quo = quo.Promote(a);
			quo = (s & -quo) | ((~s) & quo);
		}

		return(quo.Minimize());
	}
	pint operator/(pint b) const { return(Div(b)); }
	pint operator/=(pint b) { *this = Div(b); return(*this); }

	pint Rem(pint b) const {
		// C's % operator is called modulus,
		// but modulus should result in sign of divisor;
		// C's % is really remainder, with sign of dividend
		pint s(0);
		pint a = *this;

		// signed?
		if (a.type || b.type) {
			// save the result sign
			a = a.Promote(b);
			b = b.Promote(a);
			b.bit[0] = a.bit[prec-1] ^ b.bit[prec-1];
			a = a.Abs();
			b = b.Abs();
		}

		// do the unsigned divide
		a = a.Promote(b);
		b = b.Promote(a);
		pint rem(a.prec, 0);

		// classic shift & subtract loop
		for (int i=a.prec-1; i>=0; --i) {
			// shift left 1 bit
			for (int j=a.prec-1; j>0; --j) {
				rem.bit[j] = rem.bit[j-1];
			}
			rem.bit[0] = a.bit[i];

			pint ge = (rem >= b);
			rem = (ge & (rem - b)) | ((!ge) & rem);
			rem = rem.Extend(a.prec);
		}

		// fix sign if possibly negative
		if (s.bit[0].v != 0) {
			// make it signed again
			rem = rem.Promote(a);
			rem = (s & -rem) | ((~s) & rem);
		}

		return(rem.Minimize());
	}
	pint operator%(pint b) const { return(Rem(b)); }
	pint operator%=(pint b) { *this = Rem(b); return(*this); }

	friend pint H(const int bits, const int mask);

	int Any() const {
		// this is the classic SIMD test
		// true if any entanglement channel value is non-0
		pint r = *this;
		r = r.Logic();
		return(r.bit[0].v ? 1 : 0);
	}

	int All() const {
		// this is the classic SIMD test
		// true is all entanglement channels values are non-0
		pint r = *this;
		r = r.Logic();
		return((r.bit[0].v == 1) ? 1 : 0);
	}

};

pint
H(const int bits, const int mask) {
	// make a bits-pbit Hadamard value using entanglement channels per mask
	// this is sort-of a Hadamard pint constructor...
	// but declaring it as a friend seems the only way to name it differently
	register int m = (mask ? mask : (REBITS - 1));
	register int h = 0;
	register int b = ((bits > 0) ? bits : 1); // force at least 1 bit
	pint r;

	// set precision and signedness from bits
	r.prec = ((r.type = (b < 0)) ? -b : b);
	for (int i=0; i<r.prec; ++i) {
		while (!(m & 1)) {
			++h;
			if (!(m >>= 1)) {
				// not enough bits in mask
				r.prec = i;
				return(r);
			}
		}
		r.bit[i] = pbit(pbitH(h));
		m >>= 1;
		++h;
	}
	return(r);
}

pint
H(const int bits)
{
	// Hadamard constructor across all entanglement channels
	return(H(bits, ((1 << PINTBITS) - 1)));
}

#endif