Note: this is a simple manual HTMLizing of the original post, except that I've removed the ^L's and pagination. See also PDP Systems Architecture., FOLDOC/Jargon's entry for HAKMEM as well as their entries on PDP, PDP-10, and PDP-11.


From: gkn@ucsd.Edu (Gerard K. Newman)
Newsgroups: alt.folklore.computers
Subject: HAKMEM (well, the pieces I have anyway)
Keywords: HAKMEM, PDP-10
Message-ID: <18495@ucsd.Edu>
Date: 6 Sep 90 16:01:18 GMT
Organization: San Diego Supercomputer Center

I got 8 requests for HAKMEM in a period of about 36 hours, so I decided to go ahead and post it. It follows the dotted line below. Enjoy.

Programming hacks for the PDP10/DECSYSTEM-20. The first 36 hacks are taken verbatim from:

HAKMEM - Artificial Intelligence Memo No. 239

Massachusetts Institute of Technology A. I. Laboratory
M. Beeler, R. W. Gosper & R. Schroeppel, February 29, 1972

Items 1 - 36 correspond to items 145 - 180 (pages 72 - 86) in the original document.

WARNING: Numbers in this section are octal (and occasionally binary) unless followed by a decimal point.

105=69.. (And 105.=69 hexadecimal.)

ITEM 1 (Gosper):

Proving that short programs are neither trivial nor exhausted yet, there is the following:
     0/      TLCA 1,1(1)
     1/      see below
     2/      ROT 1,9
     3/      JRST 0

This is a display hack (that is, it makes pretty patterns) with the low 9 bits = Y and the 9 next higher = X; also, it makes interesting, related noises with a stereo amplifier hooked to the X and Y signals. Recommended variations include:

       none            377767,,377767; 757777,,757757; etc.
       TLC 1,2(1)      373777,,0; 300000,,0
       TLC 1,3(1)      -2,,-2; -5,,-1; -6,,-1
       ROT 1,1         7,,7; A0000B,,A0000B
e       ROTC 1,11       ;Can't use TLCA over data.
       AOJA 1,0


Another simple display program: ("munching squares")

It is thought that this was discovered by Jackson Wright on the RLE PDP-1 circa 1962.

     DATAI 2
     ADDB 1,2
     ROTC 2,-22
     XOR 1,2
     JRST .-4

2=X, 3=Y. Try things like 1001002 in data switches. This also does interesting things with operations other than XOR, and rotations other than -22. (Try IOR; AND; TSC; FADR; FDV(!); ROT -14, -9, -20, ...)

ITEM 3 (Schroeppel):

Munching squares is just views of the graph Y = X XOR T for consecutive values for T = time.

ITEM 4 (Cohen, Beeler):

A modification to munching squares which reveals them in frozen states through opening and closing curtains: insert FADR 2,1 before the XOR.

Try data switches =

4000,,4      1000,,2002      2000,,4      0,,1002
(Notation:  <left half>,,<right half>
Also try the FADR after the XOR, switches = 1001,,1.

ITEM 5 (Minsky):

Here is an elegant way to draw almost circles on a point-plotting display. CIRCLE ALGORITHM:
     NEW X = OLD X - e * OLD Y
     NEW Y = OLD Y + e * NEW(!) X

This makes a very round ellipse centered at the origin with its size determined by the initial point. e determines the angular velocity of the circulating point, and slightly affects the eccentricity. If e 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.

ITEM 6 (Schroeppel):

PROBLEM: Although the reason for the circle algorithm's stability is unclear, what is the number of distinct sets of radii? (Note: algorithm is invertible, so all points have predecessors.)

ITEM 7 (Gosper):

Separating X from Y in the above recurrence,
     X(N+1) = (2-e^2)*X(N) - X(N-1)
     Y(N+1) = (2-e^2)*Y(N) - Y(N-1).
These are just the Chebychev recurrence with cos A (the angular increment) = 1-(e^2)/2. Thus X(N) and Y(N) are expressible in the form R cos(N A + B). The B's and R for X(N) and Y(N) can be found from N=0,1. The B'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.

ITEM 8 (Salamin):

With exact arithmetic, the circle algorithm is stable if |e| < 2. In this case, all points lie on the ellipse
     X^2 - e * X * Y + Y^2 = constant,
where the constant is determined by the initial point. This ellipse has its major axis at 45 degrees (if e > 0) or 135 degrees (if e < 0) and has eccentricity
     (e/(1 + e/2)^2.

ITEM 9 (Minsky):

To portray a 3-dimensional solid on a 2-dimensional display, we can use a single circle algorithm to compute orbits for the corners to follow. The (positive or negative) radius of each orbit is determined by the distance (forward or backward) from some origin to that corner. The solid will appear to wobble rigidly about the origin, instead of simply rotating.

ITEM 10 (Gosper):

The myth that any given programming language is machine independent is easily exploded by computing the sum of powers of 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 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.

By this strategy, consider the universe, or, more precisely algebra:

  • 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.
  • ITEM 11 (Liknaitzky):

    To subtract the right half of an accumulator from the left (as in restarting an AOBJN counter):
         IMUL A,[377777,,1]

    ITEM 12 (Mitchell):

    To make an AOBJN pointer when the origin is fixed and the length is a variable in A:
         HRLOI A,-1(A)
         EQVI A,ORIGIN

    ITEM 13 (Freiberg):

    If instead, A is a pointer to the last word
         HRLOI A,-ORIGIN(A)
         EQVI A,ORIGIN

    Slightly faster: change the HRLOIs to MOVSIs and the EQVI addresses to -ORIGIN-1. These two routines are clearly adjustable for BLKOs and other fenceposts.

    ITEM 14 (Gosper, Salamin, Schroeppel):

    A miniature (recursive) sine and cosine routine follows.
         COS:    FADR A,[1.57079632679] ;Pi/2
         SIN:    MOVM B,A        ;argument in A
                 CAMG B,[.00017] ;< 3^(-3) / 2^(13)
                 POPJ P,         ;sin X = X, within 27. bits
                 FDVRI A,(-3.0)
                 PUSHJ P,SIN     ;sin -X/3
                 FMPR B,B
                 FSC B,2
                 FADRI B,(-3.0)
                 FMPRB A,B       ;sin X = 4(sin -X/3)^3 -3(sin -X/3)
                 POPJ P,         ;sin in A, sin or |sin| in B
         ;|sin| in B occurs when angle is smaller than end test

    Changing both -3.0's to +3.0's gives sinh:

         sinh X = 3 sinh X/3 + 4 (sinh X/3)^3.

    Changing the first -3.0 to a +9.0, then inserting PUSHJ P,.+1 after PUSHJ P,SIN gains about 20% in speed and uses half the pushdown space (< 5 levels in the first 4 quadrants). PUSHJ P,.+1 is a nice way to have something happen twice. Other useful angle multiplying formulas are:

         tanh X = (2 tanh X/2)/(1 + (tanh X/2)^2)
         tan X = (2 tan X/2)/(1 - (tan X/2)^2)
    if infinity is handled correctly. For cos and cosh, one can use:
         cos X = 1 - 2 (sin X/2)^2, cosh X = 1 + 2 (sinh X/2)^2.

    In general, to compute functions like e^X, cos X, elliptic functions, etc. by iterated application of double and triple argument formulas, it is necessary to subtract out the constant in the Taylor series and transform the range reduction formula accordingly. Thus:

         F(X) = cos(X)-1    F(2X) = 2F*(F+2)     F(e) = -e^2/2
         G(X) = e^X-1       G(2X) = G*(G+2)      G(e) = e

    This is to prevent the destruction of the information in the range-reduced argument by the addition of a quantity near 1 upon the success of the e test. The addition of such a quantity in the actual recurrences is OK since the information is restored by the multiply. In fact, a cheap and dirty test for F(e) sufficiently small is to see if the addition step has no effect. People lucky enough to have a square root instruction can get natural log by iterating X <- X/((1+X)^2 + 1) until 1+X = 1.

    Then multiply by 2^(number of iterations). Here, a LSH or FSC would work.

    ITEM 15 (Gosper, Schroeppel):

    (Numbers herein are decimal.)

    The correct epsilon test in such functions as the foregoing SIN are generally the largest argument for which addition of the second term has no effect on the first. In SIN, the first term is x and the second is -(x^3)/6, so the answer is roughly the x which makes the ratio of those terms (1/(2^27); so x = (3^(-3))/(2^(13)). But this is not exact, since the precise cutoff is where the neglected term is the power of 2 whose 1 bit coincides with the first neglected (28th) bit of the fraction. Thus, (x^3)/6 = 1/(2^(27)) * 1/(2^(13)), so x = (3^(-3)) / (2^(13).

    ITEM 16 (Gosper):

    Here is a way to get log base 2. A and B are consecutive. Call by PUSHJ P,LOG2 with a floating point argument in A.
         LOG2:   LSHC A,-33
                 MOVSI C,-201(A)
                 TLC C,211000    ;Speciner's bum
                 MOVEI A,200     ;exponent and sign sentinel
         LOGL:   LSH B,-9
         REPEAT 7, FMPR B,B      ;moby flunderflo
                 LSH B,2
                 LSHC A,7
                 SOJG A,LOGL     ;fails on 4th try
                 LSH A,-1
                 FADR A,C
                 POPJ P,         ;answer in A

    Basically, you just square seven times and use the low seven bits of the exponent as the next seven bits of the log.

    ITEM 17 (Gosper):

    To swap the contents of two locations in memory:
         EXCH A,LOC1
         EXCH A,LOC2
         EXCH A,LOC1
    Note: LOC1 must not equal LOC2! If this can happen use MOVE-EXCH-MOVEM, clobbering A.

    ITEM 18 (Gosper):

    To swap two bits in an accumulator:
         TRCE A,BITS
         TRCE A,BITS
         TRCE A,BITS
    Note (Nelson): last TRCE never skips, and used to be TRC, but TRCE is less forgettable. Also, use TLCE or TDCE if the bits are not in the right half.

    ITEM 19 (Sussman):

    To exchange two variables in LISP without using a third variable:
         (SETQ X (PROG2 O Y (SETQ Y X)))

    ITEM 20 (Samson):

    To take MAX in A of two byte pointers (where A and B are consecutive accumulators):
         ROTC A,6
         CAMG A,B
         EXCH A,B
         ROTC A,-6

    ITEM 21 (Freiberg):

    A byte pointer can be converted to a character address < 2^(18) by:
         MULI A,<# bytes/word>
         SUBI B,1-<# b/w>(A)

    To get full word character address, use SUB into magic table.

    ITEM 22 (Gosper, Liknaitzky):

    To rotate three consecutive accumulators N < 37. places:
         ROTC A,N
         ROT B,-N
         ROTC B,N

    Thus M AC's can be ROTC'ed in 2M-3 instructions. (Stallman): For 73. > N > 35.:

         ROTC A,N-36
         EXCH A,C
         ROT B,36.-N
         ROTC A,N-72.

    ITEM 23 (Gosper, Freiberg):

    ;B gets 7 bit character in A with even parity
            IMUL A,[2010040201]     ;5 adjacent copies
            AND A,[21042104377]     ;every 4th bit of left 4 copies + right copy
            IDIVI A,17              ;casting out 15.'s in  hexadecimal shifted 7
    ;odd parity on 7 bits (Schroeppel)
            IMUL A,[10040201]       ;4 adjacent copies
            IOR A,[755555540]       ;leaves every 3rd bit+offset+right copy
            IDIVI A,9               ;powers of 2^3 are +- mod 9
    ;changing 7555555400 to 27555555400 gives even parity
    ;if A is a 9 bit quantity, B gets number of 1's (Schroeppel)
            IMUL A,[1001001001]     ;4 copies
            AND A,[42104210421]     ;every 4th bit
            IDIVI A,17              ;casting out 15.'s in hexadecimal
    ;if A is 6 bit quantity, B gets 6 bits reversed (Schroeppel)
            IMUL A,[2020202]        ;4 copies shifted
            AND A,[104422010]       ;where bits coincide with reverse repeated base 2^8
            IDIVI A,377             ;casting out 2^8 -1's
    ;reverse 7 bits (Schroeppel)
            IMUL A,[10004002001]    ;4 copies set by 000's base 2 (may set arith. o'flow)
            AND A,[210210210010]    ;where bits coincide with reverse repeated base 2[8]
            IDIVI A,377             ;casting out 377's
    ;reverse 8 bits (Schroeppel)
            MUL A,[100200401002]    ;5 copies in A and B
            AND B,[20420420020]     ;where bits coincide with reverse repeated base 2[10]
            ANDI A,41               ;"
            DIVI A,1777             ;casting out 2[10]-1's

    ITEM 24 (PDP-1 hackers):

    foo,    lat       /DATAI switches
            adm a     /ADDB
            and (707070
            adm b
            iot 14    /output AC sign bit to a music flip-flop
            jmp foo

    Makes startling chords, arpeggios, and slides, with just the sign of the AC. This translates to the PDP-6 (roughly) as:

    FOO:    DATAI 2
            ADDB 1,2
            AND 2,[707070707070]    ;or 171717171717, 363636363636, 454545454545, ...
            ADDB 2,3
            LDB 0,[360600,,2]
            JRST FOO
    Listen to the square waves from the low bits of 0.

    ITEM 25 (in order of one-ups-manship: Gosper, Mann, Lenard, [Root and Mann]):

    To count the ones in a PDP-6/10 word:
            LDB B,[014300,,A]       ;or MOVE B,A then LSH B,-1
            AND B,[333333,,333333]
            SUB A,B
            LSH B,-1
            AND B,[333333,,333333]
            SUBB A,B   ;each octal digit is replaced by number of 1's in it
            LSH B,-3
            ADD A,B
            AND A,[070707,070707]
            IDIVI A,77              ;casting out 63.'s
    These ten instructions, with constants extended, would work on word lengths up to 62.; eleven suffice up to 254..

    ITEM 26 (Jensen):

    Useful strings of non-digits and zeros can arise when carefully chosen negative numbers are fed to unsuspecting decimal print routines. Different sets arise from different methods of character-to-digit conversion.

    Example (Gosper):

    DPT:    IDIVI F,12
            HRLM G,(P)      ;tuck remainder on pushdown list
            SKIPE F
            PUSHJ P,DPT
            LDB G,[220600,,(P)]   ;retrieve low 6 bits of remainder
            TRCE G,"0       ;convert digit to character
            SETOM CCT       ;that was no digit!
    TYO:    .IOT TYOCHN,G   ;or DATAO or IDPB ...
            AOS G,CCT
            POPJ P,

    This is the standard recursive decimal print of the positive number in F, but with a LDB instead of a HLRZ. It falls into the typeout routine which returns in G the number of characters since the last carriage return. When called with -36., DPT types carriage return, line feed, and resets CCT, the character position counter.

    ITEM 27 (Gosper):

    Since integer division can never produce a larger quotient than dividend, doubling the dividend and divisor beforehand will distinguish division by zero from division by 1 or anything else, in situations where division by zero does nothing.

    ITEM 28 (Gosper):

    The fundamental operation for building list structure, called CONS, is defined to: find a free cell in memory, store the argument in it, remove it from the set of free cells, return a pointer to it, and call the garbage collector when the set is empty. This can be done in two instructions:
    CONS:   EXCH A,[EXCH A,[...[PUSHJ P,GC]...]]
            EXCH A,CONS
    Of course, the address-linked chain of EXCH's indicated by the nested brackets is concocted by the garbage collector. This method has the additional advantage of not constraining an accumulator for the free storage pointer.
            EXCH A,CONS
            EXCH A,@CONS
    Returns cell addressed by A to free storage list; returns former cell contents in A.

    ITEM 29 (Gosper):

    The incantation to fix a floating number is usually
         MULI A,400          ;exponent to A, fraction to A+1
         TSC A,A             ;1's complement magnitude of excess 200 exponent
         ASH A+1,-200-27.-8(A) ;answer in A+1

    If number is known positive, you can omit the TSC. On the PDP-10

         UFA A,[+-233000,,]  ;not in PDP-6 repertoire
         TLC A+1,233000      ;if those bits really bother you
    When you know the sign of A, and |A| < 2^(26), you can
         FAD A,[+-233400,,]  ;or FADR for rounded fix!
         TLC A,233400        ;if those bits are relevant
    where the sign of the constant must match A's. This works on both machines and doesn't involve A+1. On the 10, FADRI saves a cycle and a constant, and rounds.

    ITEM 30 (Gosper, Nelson):

    21963283741. = 243507216435 is a fixed point of the float function on the PDP-6/10, i.e., it is the only positive number whose floating point representation equals its fixed.

    ITEM 31 (Gosper):

    To get the next higher number (in A) with the same number of 1 bits: (A, B, C, D do not have to be consecutive)
         MOVE B,A
         MOVN C,B
         AND C,B
         ADD A,C
         MOVE D,A
         XOR D,B
         LSH D,-2
         IDIVM D,C
         IOR A,C

    ITEM 32 (Gosper):

    The "banana phenomenon" was encountered when processing a character string by taking the last 3 letters typed out, searching for a random occurrence of that sequence in the text, taking the letter following that occurrence, typing it out, and iterating. This ensures that every 4-letter string output occurs in the original. The program typed ANANANANANANANA.... We note an ambiguity in the phrase, "the Nth occurrence of." In one sense, there are five 00's in 0000000000; in another there are nine. The editing program TECO finds five. Thus it finds only the first ANA in BANANA, and is thus obligated to type N next. By Murphy's Law, there is but one NAN, thus forcing A, and thus a loop. An option to find overlapped instances would be useful, although it would require backing up N-1 characters before seeking the next N character string.


    Certain plotters and displays are constrained to approximate curves by a sequence of king-moves between points on a lattice.

    Many curves and contours are definable by F(X,Y) = 0 with F changing sign on opposite sides of the curve. The following algorithm will draw most such curves more accurately than polygonal approximations and more easily than techniques which search for a "next" X and Y just one move away.

    We observe that a good choice of lattice points is just those for which F, when evaluated on one of them, has opposite sign and smaller magnitude than on one or more of its four immediate neighbors.[+] This tends to choose the nearer endpoint of each graph paper line segment which the curve crosses, if near the curve F is monotone with distance from the curve.

    First, divide the curve into arcs within which the curve's tangent lies within one 45 degree semiquadrant. We can show that for reasonable F, only two different increments (say north and northwest) are needed to visit the desired points.

    Thus, we will be changing one coordinate (incrementing Y) every step, and we have only to check whether changing the other (decrementing X) will reduce the magnitude of F. (If F increases with Y, F(X,Y+1) > -F(X-1,Y+1) means decrement X.) F can often be manipulated so that the inequality simplifies and so that F is easily computed incrementally from X and Y.

    As an example, the following computes the first semiquadrant of the circle

                         F = X^2+Y^2-R^2 = 0.
         C0:     F <- 0, Y <- 0, X <- R
         C1:     F <- F+2Y+1, Y <- Y+1
         C2:     if+ F +> X, F <- F-2X+1, X <- X-1
         C3:     if Y < X-1, go to C1
         C4:     (Link to next arc) if Y = X-1, Y <- Y+1, X <- X-1

    This can be bummed by maintaining Z = 2Y+1 instead of Y. Symmetry may be used to compute all eight semiquadrants at once, or the loop may be closed at C2 and C3 with two PUSHJ's to provide the palindrome of decisions for the first quadrant. There is an expression for the number of steps per quadrant, but it has a three-way conditional dependent upon the midpoint geometry. Knowing this value, however, we can replace C3 and C4 with a simple loop count and an odd-even test for C4.

    The loop must be top-tested (C3 before C1) if the "circle" R=1, with four diagonal segments, is possible.

    All this suggests that displays might be designed with an increment mode which accepts bit strings along with declarations of the form: "0 means north, 1 means northwest". 1100 (or 0011) will not occur with a curve of limited curvature; thus, it could be used as an escape code, but this would be an annoying restriction.

    See the following illustration of circles drawn this way.

    [+] In case of a tie, i.e., F has equal magnitudes with opposite signs on adjacent points, do not choose both points but rather have some arbitrary yet consistent preference for, say, the outer one. The problem can't arise for C2 in the example because the inequality F >= X is really F > -(F-2X+1) or F > X-.5.

    ITEM 34 (Schroeppel, Salamin):

    Suppose Y satisfies a differential equation of the form
         P(X)Y(Nth derivative) + .....  + Q(X) =R(X)
    where P, ..... Q, and R are polynomials in X (for example, Bessel's equation, X^2Y''+XY'+(X^2-N^2)Y = 0) and A is an algebraic number. Then Y(A) can be evaluated to N places in time proportional to N(ln N)^3.

    Further, e^X and ln X or any elementary function can be evaluated to N places in n(ln N)^2 for X a real number. If F(X) (by Newton's method), and the first derivative of F(X). Also zeta(3) and gamma can be done in N(ln N)^3.

    ITEM 35 (Gosper):

    A program which searches a character string for a given substring can always be written by iterating the sequence fetch-compare-transfer (ILDB-CAIE-JRST on the PDP6/10) once for each character in the sought string. The destinations of the transfers (address fields of the JRST's) must, however, be computed as functions of the sought string. Let
         0 1 2 3 4
           S A S S Y
           0 1 0 2 2
    stand for the program
         T0:     ILDB C,A       ;C gets next char from pointer in A
         T1:     CAIE C,"S      ;skip if it's an S
                 JRST TO        ;loop back on failure
                 ILDB C,A       ;next
         T2:     CAIE C,"A      ;skip if A
                 JRST T1        ;could be an S
                 ILDB C,A
         T3:     CAIE C,"S
                 JRST TO        ;S, A, non S, so start over
                 ILDB C,A       ;next
         T4:     CAIE C,"S
                 JRST T2        ;could be SAS.ASSY
                 ILDB C,A
                 CAIE C,"Y
                 JRST T2        ;could be SASS.ASSY
         ;found SASSY

    In other words, a number >0 in the top row is a location in the program where the corresponding letter of the middle row is compared with a character of th input string. If it differs, the number in the bottom row indicates the location where comparison is to resume. If it matches, the next character of the middle row is compared with the next character of the input string.

    Let J be a number in the top row and K be the number below J, so that TK is the address of the Jth JRST. For each J = 1, 2, ... we compute K(J) as follows:

         K(1) = 0.  Let P be a counter, initially 0.

    For each succeeding J, increment P. If the Pth letter = the Jth, K(J) = K(P). Otherwise, K(J) = P, and P is reset to 0. (P(J) is the largest number such that the first P characters match the last P characters in the first J characters of the sought string.)

         J=      0 1                        0 1 2 3 4 5
                   M I S S I S S I P P I      I S S I S S I P P I
         K(J)=     0 1 1 1 1 1 1 1 1 1 1      0 1 1 0 1 1 0 5 1 0
                 0 1 2 3                    0 1 2 3
                   C O C A C O L A            S A S S A F R A S
                   0 1 0 2 0 1 3 1            0 1 0 2 1 3 1 1 0

    To generalize this method to search for N strings at once, we produce a program of ILDB-CAIE=JRST's for each of the sought strings, omitting the initial ILDB from all but the first. We must compute the destination of the Jth JRST in the Ith program, TKM(I,J) which is the location of the Kth compare in the Mth program.

    It might be reasonable to compile such an instruction sequence whenever a search is initiated, since alternative schemes usually require saving or backing up the character pointer.

    ITEM 36 (Gosper):

    A problem which may arise in machine processing of visual information is the identification of corners on a noisy boundary of a polygon. Assume you have a broken line. If it is a closed loop, find the vertex furthest from the centroid (or any place). Open the loop by making this place both endpoints and calling it a corner. We define the corner of a broken line segment to be the point the sum of whose distances from the endpoints is maximal. This will divide the segment in two, allowing us to proceed recursively, until our corner isn't much cornerier than the others along the line.

    The perpendicular distance which the vector C lies from the line connecting vectors A and B is just

          (C - A) x (B - A)
         ------------------- ,
              2 |A - B|
    but maximizing this can lose on very pointy V's. The distance sum hack can lose on very squashed Z's.

    ITEM 37 (@Hurley: origin unknown)

    Random number generation:
            ;Initialize random-number generation
            GTAD                    ;Get current time and date
            MOVEM T1,SEED           ;Make it the seed
            CALL RANDOM             ;Go get a random number
            ;Return with psuedo random number in T2
    RANDOM: MOVE T1,SEED            ;Get the seed
            MUL T1,RAN              ;Make next random number
            MOVEM T2,SEED           ;Save low order as new seed
            IDIVI T1,10             ;Get random number in range 0-7
    RAN:    343277,,244615          ;Low order portion of 5^15
    SEED:   0

    gkn Gerard K. Newman        gkn@sds.sdsc.edu    619.534.5076
    San Diego Supercomputer Center  gkn@sdsc.bitnet         619.534.5152 FAX
    PO Box 85608                    sdsc::gkn (27.1/span)
       San Diego, CA  92186-9784   ucsd!gkn