Beeler, M., Gosper, R.W., and Schroeppel, R.

The original examples were in PDP-10 assembler language. The PDP-10 had 36-bit data registers; I have tried to generalise the tricks so as to work on any word length (subject to restrictions as noted). Moreover I have chosen octal, decimal and hexadecimal notation to ease reading the magic numbers involved. I have also collected additional (non-HAKMEM) tricks.

w = w * 0xffff0001;

u ^= v; v ^= u; u ^= v;

/* evenp7() adds even parity to its 7 bit character argument... */ unsigned evenp7(unsigned36 a) { return ((a * 002010040201) /* 5 adjacent copies */ & 021042104377) /* every 4th bit of left 4 copies + right copy */ % (15<<7); /* casting out 15.'s in hexadecimal shifted 7 */ } /* odd parity on 7 bits (Schroeppel) */ unsigned oddp7(unsigned36 a) { return ((a * 000010040201) /* 4 adjacent copies */ | 007555555400) /* leaves every 3rd bit+offset+right copy */ % (9<<7); /* powers of 2^3 are +-1 mod 9 */ /* changing 7555555400 to 27555555400 gives even parity */ } /* count number of 1's in 9-bit argument (Schroeppel) */ unsigned count_ones(unsigned36 a) { return ((a * 01001001001) /* 4 adjacent copies */ & 042104210421) /* every 4th bit */ % 15; /* casting out 15.'s in hexadecimal */ } /* if argument is 6 bit quantity, return 6 bits reversed (Schroeppel) */ unsigned reverse_6bits(unsigned36 a) { return ((a * 02020202) /* 4 adjacent copies shifted */ & 0104422010) /* where bits coincide with reverse repeated base 2^8 */ % 255; /* casting out 2^8 - 1's */ } /* reverse 7 bits (Schroeppel) */ unsigned reverse_7bits(unsigned36 a) { return ((a * 010004002001) /* 4 copies sep by 000's base 2 */ & 0210210210010) /* where bits coincide with reverse repeated base 2^8 */ % 255; /* casting out 2^8 - 1's */ /* reverse 8 bits (Schroeppel) */ unsigned reverse_8bits(unsigned41 a) { return ((a * 0x000202020202) /* 5 copies in 40 bits */ & 0x010884422010) /* where bits coincide with reverse repeated base 2^10 */ /* PDP-10: 041(6 bits):020420420020(35 bits) */ % 1023; /* casting out 2^10 - 1's */ }

unsigned count_ones(unsigned36 a) { unsigned36 b; b = (a>>>1) & 0333333333333; a -= b; b = (b>>>1) & 0333333333333; a -= b; /* each octal digit is replaced by number of 1's in it */ a = ((a>>3) + a) & 0070707070707; return a % 63; /* casting out 63's */ /* This code (10 PDP-10 instructions), with constants extended, would work on word lengths up to 62; eleven suffice up to 254.

unsigned nexthi_same_count_ones(unsigned a) { /* works for any word length */ unsigned c = (a & -a); unsigned r = a+c; return (((r ^ a) >> 2) / c) | r); }

union { int i; float f; } u; if ((float)u.i == u.f) /* then we have u.i = 0x00000000 */

NEW X = OLD X - *epsilon* * OLD Y

NEW Y = OLD Y + *epsilon* * NEW(!) X

This makes a very round ellipse centered at the origin with its size
determined by the initial point. *epsilon* determines the angular
velocity of the circulating point, and slightly affects the
eccentricity. If *epsilon* is a power of 2, then we don't even need
multiplication, let alone square roots, sines, and cosines! The
"circle" will be perfectly stable because the points soon become
periodic.

The circle algorithm was invented by mistake when I tried to save one register in a display hack! Ben Gurley had an amazing display hack using only about six or seven instructions, and it was a great wonder. But it was basically line-oriented. It occurred to me that it would be exciting to have curves, and I was trying to get a curve display hack with minimal instructions.

X(N+1) = (2 - *epsilon*^2) * X(N) - X(N-1)

Y(N+1) = (2 - *epsilon*) * Y(N) - Y(N-1).

These are just the Chebychev recurrence with cos *theta* (the angular
increment) = 1-*epsilon*^2/2. Thus X(N) and Y(N) are expressible in the
form R cos(N *theta* + *phi*). The *phi*'s and R for X(N) and Y(N) can be
found from N=0,1. The *phi*'s will differ by less than *pi*/2 so that the
curve is not really a circle. The algorithm is useful nevertheless,
because it needs no sine or square root function, even to get
started.

X(N) and Y(N) are also expressible in closed form in the algebra of ordered pairs described under linear recurrences, but they lack the remarkable numerical stability of the "simultaneous" form of the recurrence.

X^2 - *epsilon* X Y + Y^2 = constant,

where the constant is determined by the initial point. This ellipse
has its major axis at 45 degrees (if *epsilon* > 0) or 135 degrees
(if *epsilon* < 0) and has eccentricity

sqrt(*epsilon*/(1 + *epsilon*/2)).

- If the result loops with period = 1 with sign +, you are on a sign-magnitude machine.
- If the result loops with period = 1 at -1, you are on a twos-complement machine.
- If the result loops with period > 1, including the beginning, you are on a ones-complement machine.
- If the result loops with period > 1, not including the beginning, your machine isn't binary -- the pattern should tell you the base.
- If you run out of memory, you are on a string or Bignum system.
- If arithmetic overflow is a fatal error, some fascist pig with a read-only mind is trying to enforce machine independence. But the very ability to trap overflow is machine dependent.

let X = the sum of many powers of two = ...111111

now add X to itself; X + X = ...111110

thus, 2X = X - 1 so X = -1

therefore algebra is run on a machine (the universe) which is twos-complement.

(SETQ X (PROG2 0 Y (SETQ Y X)))