TweetFollow Us on Twitter

Nov 94 Challenge
Volume Number:10
Issue Number:11
Column Tag:Programmer’s Challenge

Programmer’s Challenge

By Mike Scanlin, Mountain View, CA

Note: Source code files accompanying article are located on MacTech CD-ROM or source code disks.

The Rules

Here’s how it works: Each month we present a different programming challenge here. First, you write some code that solves the challenge. Second, optimize your code (a lot). Then, submit your solution to MacTech Magazine. We choose a winner based on code correctness, speed, size and elegance (in that order of importance) as well as the postmark of the answer. In the event of multiple equally-desirable solutions, we’ll choose one winner at random (with honorable mention, but no prize, given to the runners up). The prize for each month’s best solution is $50 and a limited-edition “The Winner! MacTech Magazine Programming Challenge” T-shirt (not available in stores).

To help us make fair comparisons, all solutions must be in ANSI compatible C (e.g. don’t use Think’s Object extensions). Use only pure C code. We disqualify any entries with any assembly in them (except for challenges specifically stated to be in assembly). You may call any routine in the Macintosh toolbox (e.g., it doesn’t matter if you use NewPtr instead of malloc). We test entries with the FPU and 68020 flags turned off in THINK C. We time routines with the latest THINK C (with “ANSI Settings”, “Honor ‘register’ first”, and “Use Global Optimizer” turned on), so beware if you optimize for a different C compiler. Limit your code to 60 characters wide. This helps with e-mail gateways and page layout.

We publish the solution and winners for this month’s Programmers’ Challenge two months later. All submissions must be received by the 10th day of the month printed on the front of this issue.

Mark solutions “Attn: Programmers’ Challenge Solution” and send them via e-mail - Internet progchallenge@xplain.com, AppleLink MT.PROGCHAL, CompuServe 71552,174 and America Online MT PRGCHAL. Include the solution, all related files, and your contact info. If you send via snail mail, send a disk with those items on it; see “How to Contact Us” on p. 2.

MacTech Magazine reserves the right to publish any solution entered in the Programming Challenge of the Month. Authors grant MacTech Magazine the non-exclusive right to publish entries without limitation upon submission of each entry. Authors retain copyrights for the code.

Huffman Decoding

Being able to decode a compressed bit stream quickly is important in many applications. This month you’ll get a chance to decode one of the most commonly used compression formats around. Huffman codes are variable length bit strings that represent some other bit string. I’m not going to explain the algorithm here (see any decent book on algorithms for that) but I will explain the format of the symbol table you’ll be given and show you how to decode using it (which is all you need to know to do this Challenge).

The symbol table you are passed consists of an array of elements that look like this:


/* 1 */
typedef struct SymElem {
 unsigned short  symLength;
 unsigned short  sym;
 unsigned short  value;
} SymElem, *SymElemPtr;

where sym is the compressed bit pattern, symLength is the number of bits in sym (from 1 to 16, starting from the least significant bit of sym) and value is the uncompressed output value (16 bits). The symbol table will be sorted smallest to largest first by length and then, within each length, by sym. For example, if you had a table with two SymElems like this:

sym = 3; symLength = 2; value = 0xAAAA
sym = 1; symLength = 3; value = 0xBBBB

and then you were given this compressed bit stream to decode 1100111001001, then the output would be 0xAAAA 0xBBBB 0xAAAA 0xBBBB 0xBBBB because the first two 1 bits are the 2-bit symbol ‘11’ (i.e. 3), and the next 3 bits are the 3-bit symbol ‘001’, and so on.

You will have a chance to create a lookup table in an un-timed init routine but, the amount of memory you can use is variable, from 8K to 256K. Your init routine cannot allocate more than maxMemoryUsage bytes or it will be disqualified (this includes ‘static’ and global data):

void *HuffmanDecodeInit(theSymTable, numSymElems, maxMemoryUsage);
SymElemPtr  theSymTable;
unsigned short   numSymElems;
unsigned long    maxMemoryUsage;

The return value from this init routine will be passed to the actual decode routine (as the privateHuffDataPtr parameter). The decode routine (which is timed) will be called with different sets of compressed data that use the same symbol table:

unsigned long HuffmanDecode(theSymTable, numSymElems, bitsPtr,
 numBits, outputPtr, privateHuffDataPtr)
SymElemPtr  theSymTable;
unsigned short   numSymElems;
char    *bitsPtr;
unsigned long    numBits;
unsigned short   *outputPtr;
void    *privateHuffDataPtr;

You can assume that outputPtr points to a buffer large enough to hold all of the uncompressed data. The return value from HuffmanDecode is the actual number of bytes that were stored in that buffer. The input bits are pointed to by bitsPtr and there are numBits of them. The first bit to decode is the most significant bit of the byte pointed to by bitsPtr. TheSymTable and numSyms are the same parameters that were passed to HuffmanDecodeInit.

Two Months Ago Winner

Wow! The competition was really tight for the Erase Scribble Challenge. So tight, in fact, that the 5th place winner was only about 2% slower than the first place winner. But Challenge champion Bob Boonstra (Westford, MA) was able to implement code that was just a tiny bit more efficient than many other highly efficient entries. And, as a bonus, his entry was smaller than all but one of the other entries. Despite his post-publication disqualification from the Factoring Challenge (see below) Bob remains our champion with four 1st place showings (including this one). Congratulations!

Here are the times and code sizes for each entry. Numbers in parens after a person’s name indicate how many times that person has finished in the top 5 places of all previous Programmer Challenges, not including this one:

Name time code+data

Bob Boonstra (11) 1363 598

Ernst Munter (3) 1378 1434

John Schlack 1391 1482

Tom Elwertowski (1) 1394 910

Mark Chavira 1395 1538

Jim Sokoloff 1481 1094

Allen Stenger (7) 1911 988

Marcel Rivard 2233 2048

Joshua Glazer 168100 466

At least one contestant pointed out that this Challenge was not entirely realistic because: (1) in a real eraser situation the hit-test routine would return the point that was hit to the caller and, (2) the caller would be removing points from the scribble as segments were erased, thus removing the ‘completely static scribble’ characteristic of this Challenge. I agree, it would have been more realistic to have a dynamic scribble but I was trying to limit the complexity of the routine (and I was also trying to give clever people a chance to exploit the static nature of the data by using the init routine).

Bob’s routine is well commented so I won’t discuss it here. He chose an almost identical algorithm to everyone else but he implemented it just a touch better than everyone else.

New Factoring Winner

It seems that Bob Boonstra and I both made mistakes during the Factoring Challenge (June 1994 MacTech): Bob made the mistake of incorrectly handling some input values and I made the mistake of not finding them. Many thanks to Jim Lloyd (Mountain View, CA) for finding this bug and narrowing down the set of inputs that make it happen.

The bug only happens if the number to be factored was created from composite primes where one of them has the high bit set and the other one doesn’t. If they’re both set or if they’re both clear then the bug doesn’t happen. In case you’re using Bob’s code, there is a simple fix (thanks, Bob). Change this line:

*prime2Ptr = (x+y)>>1;

to this:

*prime2Ptr = (x>>1) + (y>>1) + 1;

Because of this bug I’m going to have to retro-actively disqualify Bob’s entry from the Challenge and declare a new winner. However, the new winner is not simply the previous 2nd place winner. The new winner is from a guy whose code was originally sent in a day late (I had to disqualify it) but whose performance is so much better than anyone else who entered (including Bob) that I’m going to allow it to win in the interest of having the best possible factoring code published. It’s about two orders of magnitude faster than the other entries.

So, our new winner is Nick Burgoyne (Berkeley, CA). It turns out that this Challenge was right up Nick’s alley, considering that he has taught factoring math classes in the past. His entry was the only one to use the quadratic sieve algorithm. If you’re interested in learning more about this algorithm then Nick recommends a book by David Bressoud, Factorization and Primality Testing, published by Springer-Verlag in 1989. It covers the underlying mathematics and also gives further references to work on the quadratic sieve. It does not assume an advanced background in math. Nick is willing to discuss factoring with anyone who is interested via e-mail. His internet address is: sbrb@cats.ucsc.edu. Congrats, Nick!

Here’s Bob’s winning solution to the Erase Scribble Challenge, followed by Nick’s winning solution to the Factoring Challenge:

scribble.c

/* EraseScribble Copyright (c) 1994 J Robert Boonstra

Problem statement: Determine whether a square eraser of diameter eraserSize centered at the thePoint intersects any of the points in the data structure theScribble. Note that while the problem statement refers to line segments, the definition of a "hit" means that only the endpoints matter.

Solution strategy: The ideal approach would be to create a simple bitMap during Initialization indicating whether an eraser at a given location intersected theScribble. The bitMap would be created by stamping a cursor image at the location of each point in theScribble. However, a bitMap covering the required maximum bounding box of 1024 x 1024 would require 2^17 bytes, or four times as much storage as the 32K we are allowed to use.

Therefore, this solution has three cases:

1) If the actual bounding box for theScribble passed to the init function fits in the 32K available, we create a bitMap as above and use it directly.

2) Otherwise, we attempt to create a half-scale bitMap, where each bit represents 4 pixels in the image, 2 in .h and 2 in .v. In the PtInScribble function, we ca then quickly determine in most cases when the eraser does not intersect theScribble, and we have to walk the points in theScribble if the bitMap indicates a possible hit.

3) In the event there is not enough storage for a half-scale bitmap, we create a quarter-scale bitmap, where each bit represents 16 pixels in the image, 4 in .h and 4 in .v. Then we proceed as in case 2.

To optimize examination of the points in theScribble when the bitmap is not full scale, we sort the points in the initialization function and store them in the privateScribbleDataPtr. Although this reduces the amount of storage available for the bitmap by ~2K, it improves worst case performance significantly.

Although the init function is not timed for score, we have written it in assembler, in the spirit of the September Challenge.


/* 2 */
 */

#pragma options(assign_registers,honor_register,mc68020)

Typedefs, defines, and prototoypes
#define ushort unsigned short
#define ulong  unsigned long
#define kTotalStorage  0x8000
/*
 * Layout of privateScribbleData storage:
 *  OFFSET                                 CONTENT
 *  ------                                  -------
 *    0:                                    bitMap
 *    kTotalStorage-kGlobals-4*gNumPoints:  sorted points
 *    kTotalStorage-kGlobals:               global data
 *
 * Globals stored in PrivateScribbleDataPtr
 *   gNumPoints:number of points in the scribble
 *   gHOrigin:  min scribble h - eraserSize/2
 *   gVOrigin:  min scribble v - eraserSize/2
 *   gBMHeight: max scribble h - min scribble h + eraserSize
 *   gBMWidth:  max scribble v - min scribble v + eraserSize
 *   gRowBytes: bitMap rowBytes
 *   gMode:     flag indicating bitmap scale
 */
#define gNumPoints    kTotalStorage-2
#define gHOrigin      kTotalStorage-4
#define gVOrigin      kTotalStorage-6
#define gBMHeight     kTotalStorage-8
#define gBMWidth      kTotalStorage-10
#define gRowBytes     kTotalStorage-12
#define gMode         kTotalStorage-14
#define kGlobals                    14
#define kSortedPoints kTotalStorage-kGlobals
#define kBitMap           0
#define kHalfscaleBitMap  1
#define kQrtrscaleBitMap -1

typedef struct Scribble {
  Point startingPoint;
  Point deltaPoints[1];
} Scribble, *ScribblePtr, **ScribbleHndl;

void *EraseScribbleInit(ScribbleHndl theScribble,
      unsigned short eraserSize);

Boolean PtInScribble(Point thePoint,
      ScribbleHndl theScribble,
      unsigned short eraserSize,
      void *privateScribbleDataPtr);
PtInScribble
#define eraserH       D0
#define eraserV       D1
#define eraserSz      D2
#define pointCt       D6
#define scribblePt    D7

Boolean PtInScribble(Point thePoint,
      ScribbleHndl theScribble,
      unsigned short eraserSize,
      void *privateScribbleDataPtr)
{
  asm 68020 {
    MOVEM.L    D6-D7,-(A7)
    MOVE.L    thePoint,D1
    MOVEA.L   privateScribbleDataPtr,A0
    MOVEQ     #0,D0           ; clear high bits for BFTST
    MOVE.W    D1,D0
    SWAP      D1
; Return noHit if eraser is outside bounding box defined by gVOrigin, 
gHOrigin,
; gVOrigin+gBMHeight, gHOrigin+gBMWIdth.  Note that the bounding box 
has already 
; been expanded by eraserSize/2 in each direction.
    SUB.W     gVOrigin(A0),D1
    BLT       @noHit          ; v < boundingBox.top
    CMP.W     gBMHeight(A0),D1
    BGT       @noHit          ; v > boundingBox.bottom
    SUB.W     gHOrigin(A0),D0
    BLT       @noHit          ; h < boundingBox.left
    CMP.W     gBMWidth(A0),D0
    BGT       @noHit          ; h > boundingBox.right
; Check eraser against bitMap; return noHit if not set
    MOVE.W    gRowBytes(A0),D2
    TST.W     gMode(A0)
    BNE.S     @testQrtrScaleBitMap
; Full-scale bitMap case; can return Hit if bit is set
    MULU.W    D1,D2       ; multiple row by rowBytes
    ADDA      D2,A0       ; A0 points to correct bitMap row
    BFTST     (A0){D0:1}  ; D0 contains bit offset
    BNE       @hit
noHit:
    MOVEQ     #0,D0
    MOVEM.L   (A7)+,D6-D7
    UNLK      A6     ; save one branch by returning directly
    RTS
testQrtrScaleBitMap:
    BGT       @testHalfScaleBitMap
; Qrtr-scale bitMap case; cannot always return Hit if set
    LSR.W     #2,D1      ; * adjust for qrtr-scale bitmap *
    LSR.W     #2,D0      ; * adjust for qrtr-scale bitmap *
    MULU.W    D1,D2      ; multiple row by rowBytes
    ADDA      D2,A0      ; A0 points to correct bitMap row
    BFTST     (A0){D0:1}
    BEQ       @noHit
    BRA.S     @TestPoints
testHalfScaleBitMap:
; Half-scale bitMap case; cannot always return Hit if set
    LSR.W     #1,D1      ; * adjust for half-scale bitmap *
    LSR.W     #1,D0      ; * adjust for half-scale bitmap *
    MULU.W    D1,D2      ; multiple row by rowBytes
    ADDA.L    D2,A0      ; A0 points to correct bitMap row
    BFTST     (A0){D0:1}
    BEQ       @noHit
TestPoints:
; bitMap indicates there might be a hit, need to check
    MOVEA.L   privateScribbleDataPtr,A0
    MOVE.W    eraserSize,eraserSz
    MOVE.L    thePoint,eraserH
    MOVE.L    eraserH,eraserV
    SWAP      eraserV
    LEA       kSortedPoints(A0),A1
; Scan sets of 64 points to find a close match
    MOVE.W    gNumPoints(A0),D7
    LSR.W     #6,D7       ; numPoints/64
    MOVE.W    D7,pointCt
    LSL.W     #8,D7       ; times 4 bytes/pt * 64 pts
    SUBA.L    d7,A1
    SUBQ      #1,pointCt
eightLoop:
    MOVE.L    -(A1),scribblePt
    CMP.W     eraserH,scribblePt
    BLT.S     @hLoop
    ADDI      #65*4,A1
    DBRA      pointCt,@eightLoop;
; Note that the sorted scribblePts have been stored with
; eraserSize/2 already added in h and v
hLoop:
    MOVE.L    -(A1),scribblePt
    CMP.W     eraserH,scribblePt
    BLT.S     @hLoop
; All points from here on have a true .h component >= eraserH-eraserSize/2. 
 Now 
; need to look at those where  .h <= eraserH+eraserSize/2.  Because we 
already added 
; eraserSize/2 to the sorted point values, we compare against eraserH+eraserSz
    ADD.W     eraserSz,eraserH
    ADD.W     eraserV,eraserSz
vLoop:
    CMP.W     eraserH,scribblePt
    BGT.S     @noHit
; scribblePt.h is in range, now check .v
    SWAP      scribblePt
    CMP.W     eraserV,scribblePt
    BLT.S     @vLoopCheck
    SUB.W     eraserSz,scribblePt
    BLE.S     @hit
vLoopCheck:
    MOVE.L    -(A1),scribblePt
    BRA       @vLoop;

hit:
    MOVEQ     #1,D0
    MOVEM.L   (A7)+,D6-D7
  }
}

EraseScribbleInit
void *EraseScribbleInit(ScribbleHndl theScribble,
      unsigned short eraserSize)
{
ulong bitMapStorageAvail;
  asm 68020 {
    MOVEM.L   D3-D7/A2,-(A7)
; Allocate storage - use full 32K.
; Note that it is the responsibility of the caller to release this storage.
    MOVE.L    #kTotalStorage,D0
    NewPtr    CLEAR
    MOVE.L    A0,A2    ; A0 is privateDataPtr
; Scan thru the scribble to find min and max in h and v
    MOVEA.L   theScribble,A1
    MOVEA.L   (A1),A2  ; A2 = *theScribble
; Initialize mins and maxes to be the starting point
    MOVE.L    (A2)+,D7 ; current scribble point
    MOVEQ     #0,D5    ; clear high bits for ADDA.L later
    MOVE.W    D7,D5    ; D5 = max h
    MOVE.W    D7,D4    ; D4 = min h
    MOVE.W    D7,D1    ; D1 = current h
    SWAP      D7       ; D7 = max v
    MOVE.W    D7,D6    ; D6 = min
    MOVE.W    D7,D2    ; D2 = current v
    LEA       kSortedPoints(A0),A1 ; sorted point storage
; Set up eraserSize/2
    MOVE.W    eraserSize,D3
    LSR.W     #1,D3    ; D3 = eraserSize/2
minMaxLoop:
; Store point for subsequent sorting
    MOVE.W    D1,D0
    ADD.W     D3,D0
    MOVE.W    D0,-(A1) ; store current h + eraserSize/2
    MOVE.W    D2,D0
    ADD.W     D3,D0
    MOVE.W    D0,-(A1) ; store current v + eraserSize/2
; Fetch next point
    MOVE.L    (A2)+,D0 ; fetch deltaPoint
    BEQ       @minMaxDone
    ADD.W     D0,D1    ; update current h
    CMP.W     D1,D5
    BGE       @noNewHMax
    MOVE.W    D1,D5    ; store new max h
    BRA.S     @noNewHMin
noNewHMax:
    CMP.W     D1,D4
    BLE       @noNewHMin
    MOVE.W    D1,D4    ; store new min h
noNewHMin:
    SWAP      D0
    ADD.W     D0,D2    ; update current v
    CMP.W     D2,D7
    BGE       @noNewVMax
    MOVE.W    D2,D7    ; store new max v
    BRA.S     @noNewVMin
noNewVMax:
    CMP.W     D2,D6
    BLE       @noNewVMin
    MOVE.W    D2,D6    ; store new min v
noNewVMin:
    BRA.S     @minMaxLoop
minMaxDone:

; Calculate number of points
    LEA       kSortedPoints(A0),A2
    MOVE.L    A2,D0
    SUB.L     A1,D0
    LSR.W     #2,D0
    MOVE.W    D0,gNumPoints(A0)
; Calculate bitmap storage available
    SUBQ.L    #8,A1    ; Reserve room for sentinals
    MOVE.L    A1,D0
    SUB.L     A0,D0
    MOVE.L    D0,bitMapStorageAvail
; Calculate origin = min h/v minus eraserSize/2
    MOVE.W    eraserSize,D0
; Calculate number of columns
    SUB.W     D3,D4    ; adjust h origin for eraserSize/2
    MOVE.W    D4,gHOrigin(A0) ; D4 = H Origin
    ADD.W     D0,D5    ; add eraserSize to width
    SUB.W     D4,D5
    MOVE.W    D5,gBMWidth(A0)
    ADDQ      #7,D5    ; round up to byte level
    LSR.W     #3,D5
    MOVE.W    D5,gRowBytes(A0) ; D5 = rowBytes
; Calculate number of rows
    SUB.W     D3,D6    ; adjust v origin for eraserSize/2
    MOVE.W    D6,gVOrigin(A0) ; D6 = V Origin
    ADD.W     D0,D7    ; add eraserSize to height

    SUB.W     D6,D7
    MOVE.W    D7,gBMHeight(A0) ; D7 = number of rows
    ADDQ      #1,D7
; Calculate number of bytes needed to store bitmap.
    MULU.W    D5,D7
    CMP.L     bitMapStorageAvail,D7
    BLE.S     @haveEnoughStorage

; Not enough storage, so we try a half-scale bitMap
    MOVEQ     #kHalfscaleBitMap,D1
    MOVE.W    D1,gMode(A0)
; Adjust gRowBytes for half-scale bitMap
    MOVE.W    gBMWidth(A0),D5
    ADDI      #15,D5   ; round up to byte level
    LSR.W     #4,D5
    MOVE.W    D5,gRowBytes(A0) ; D5 = rowBytes
    LSR.W     #1,D0    ; adjust eraser size for half-scale bitmap
    MOVE.W    gBMHeight(A0),D7
    ADDQ      #1,D7
    LSR.W     #1,D7
    ADDQ      #1,D7
    MULU.W    D5,D7
    CMP.L     bitMapStorageAvail,D7
    BLE.S     @haveEnoughStorage

; Still not enough storage, so we set up a qrtr-scale bitMap
    MOVEQ     #kQrtrscaleBitMap,D1
    MOVE.W    D1,gMode(A0)
; Adjust gRowBytes for qrtr-scale bitMap
    MOVE.W    gBMWidth(A0),D5
    ADDI      #31,D5    ; round up to byte level
    LSR.W     #5,D5
    MOVE.W    D5,gRowBytes(A0) ; D5=rowBytes
    LSR.W #1,D0 ; adjust eraser size for qrtr-scale bitmap

haveEnoughStorage:
; Create bitMap
    MOVEA.L   theScribble,A1
    MOVEA.L   (A1),A2

; Initial location to stamp eraser image in bitmap is the
; startingPoint, minus the bitmap origin, minus eraserSize/2;
    MOVE.L    (A2)+,D2 ; Fetch startingPoint
    MOVE.W    D2,D1    ; D1 = current scribble point h
    SWAP      D2       ; D2 = current scribble point v
    SUB.W     D4,D1    ; D1 = scribble point h - H origin
    SUB.W     D3,D1    ; offset h by -eraserSize/2
    SUB.W     D6,D2    ; D2 = scribble point v - V origin
    SUB.W     D3,D2    ; offset v by -eraserSize/2
    MOVEQ     #0,D4    ; clear high bits for BFINS later

    MOVE.W    D0,D7    ; save eraser width

    MOVEQ     #-1,D3   ; set bits for insertion using BFINS;
; Loop through all points in the scribble
stampPointLoop:
    MOVEA.L   A0,A1
    MOVE      D7,D6    ; D6 = row counter for DBRA smear
    MOVE.W    D2,D0    ; D0 = v - origin

    TST.W     gMode(A0)
    BEQ       @L1
    BGT       @L0
    LSR.W     #1,D0    ; adjust for qrtr-scale bitmap
L0: LSR.W     #1,D0    ; adjust for half-scale bitmap
L1:
    MULU.W    D5,D0    ; D0 = (v - origin) * rowBytes
    ADDA.L    D0,A1    ; A1 = privateStorage - rowOffset
    MOVE.W    D7,D0
    ADDQ      #1,D0    ; D0=eraserSize/2+1, the number of
                       ;   bits to set in bitMap;

;   Loop through eraserSize+1 rows in bitMap for this point
stampRowLoop:
    MOVE.W    D1,D4    ; D4 = scribble h - origin

    TST.W     gMode(A0)
    BEQ       @L3
    BGT       @L2
    LSR.W     #1,D4    ; adjust for qrtr-scale bitmap
L2: LSR.W     #1,D4    ; adjust for half-scale bitmap
L3:
    BFINS     D3,(A1){D4:D0};insert mask D3 of width D0
                            ;  into bitmap A1 at offset D4
    ADDA.L    D5,A1         ; increment by rowBytes
    DBRA      D6,@stampRowLoop

;   Fetch next deltaPoint, quit if done
    MOVE.L    (A2)+,D0
    BEQ       @bitMapDone
    ADD.W     D0,D1    ; update current scribble h
    SWAP      D0
    ADD.W     D0,D2    ; update current scribble v
    BRA.S     @stampPointLoop

bitMapDone:
; Sort scribble points for fast lookup.  Simple exchange sort will do.
outerSortLoop:
    MOVEQ     #0,D6       ; exchange flag = no exchanges
    LEA       kSortedPoints(A0),A2 ; sorted pt storage
    MOVE.W    gNumPoints(A0),D7  ;outer loop counter
    SUBQ     #2,D7      ; DBRA adjustment
    MOVE.L    -(A2),D0  ; D0 = compare value - h
    MOVE.L    D0,D1
    SWAP      D1        ; D1 = compare value - v
innerSortLoop:
    MOVE.L    -(A2),D2
    MOVE.L    D2,D3
    SWAP      D3
    CMP.W     D0,D2     ; primary sort by h
    BLT.S     @doSwap
    BGT.S     @noSwap
    CMP.W     D1,D3     ; secondary sort by v
    BGE.S     @noSwap
doSwap:
    MOVE.L    D2,4(A2)
    MOVE.L    D0,(A2)
    MOVEQ     #1,D6
    BRA.S     @testInnerLoop
noSwap:
    MOVE.L    D2,D0
    MOVE.W    D3,D1
testInnerLoop:
    DBRA      D7,@innerSortLoop
    TST.W     D6
    BNE.S     @outerSortLoop;
    MOVE.L    #0x80008000,-(a2) ; Sentinal                     
       ; to end point loop
    MOVE.L    #0x7FFF7FFF,-(a2) ; Sentinal
 ; to end point loop

    MOVE.L    A0,D0 ; return ptr to private data
    MOVEM.L   (A7)+,D3-D7/A2
  }
}

Nick’s Factoring Solution

Factor64.c


/* 3 */
// Factor the product N of two primes each ¾ 32 bits
#include "Factor64.h"

// Factor64() is a simplified version of the quadratic 
//  sieve using multiple polynomials
Factor64
void Factor64(ulong nL,ulong nH,
   ulong *p1,ulong *p2) {

 register ushort ls,k;
 register ulong  i,j;
 
 ushort p,s,qi,hi,r,c,sP,sN,ts
 ushort b,m,bm,br,m2,S[5];
 ulong  d,e,sX,sQ,sq;

 Int    B,C,Q,R,T,X,Y;
 ushort  *Ptr,*pf,*sf,*lg,*r1,*r2;
 ushort *lv,**hm,*hp,**gm,*gp,*gi
 ushort *pv,**Xv,*Yv;

// Check whether N is square
 N[1] = nL & 0xffff; N[2] = nL >> 16;
 N[3] = nH & 0xffff; N[4] = nH >> 16;
 k = 4; while (N[k] == 0) k--; N[0] = k;
 
 sq = floor(sqrt(nL + 4294967296.0*nH));
 
 Set(S,sq); Mul(S,S,S);

 if (Comp(S,N) == 0) {
 *p1 = sq;
 *p2 = sq;
 return;
 } 
// Check N for small factors (up to 0x800)   
 for (k = 0; k < 616; k += 2) {
 p = Prm[k];
 if (Modq(N,p) == 0) { 
 Divq(N,p);
 *p1 = p;
 *p2 = Unset(N);
 return;
 }
 } 
//  Allocate memory
 Ptr = malloc(0x20000);
 if (Ptr == 0) {
 printf(" malloc failed \n");
 exit(0);
 }
 X     = (ushort *) Ptr;  Ptr += 20;
 Y     = (ushort *) Ptr;  Ptr += 20;
 B     = (ushort *) Ptr;  Ptr += 20;
 C     = (ushort *) Ptr;  Ptr += 20;
 Q     = (ushort *) Ptr;  Ptr += 20;
 R     = (ushort *) Ptr;  Ptr += 20;
 T     = (ushort *) Ptr;  Ptr += 20;
 pf   = (ushort *) Ptr;   Ptr += 120;
 sf   = (ushort *) Ptr;   Ptr += 120;
 lg   = (ushort *) Ptr;     Ptr += 120;
 r1   = (ushort *) Ptr;   Ptr +=   120;
 r2   = (ushort *) Ptr;   Ptr +=   120;
 lv   = (ushort *) Ptr;   Ptr += 12000;
 hm    = (ushort **) Ptr; Ptr +=   240;
 hm[0] = (ushort *)  Ptr; Ptr += 14400;
 gm    = (ushort **) Ptr; Ptr +=   240;
 gm[0] = (ushort *)  Ptr; Ptr += 24000;
 pv    = (ushort *)  Ptr; Ptr +=   120;
 Yv    = (ushort *)  Ptr; Ptr +=   120;
 Xv    = (ushort **) Ptr; Ptr +=   240;
 Xv[0] = (ushort *)  Ptr; 
 
 for (k = 1; k < 120; k++) {
 hm[k] = hm[k-1] + 120;
   gm[k] = gm[k-1] + 120 + k;
 Xv[k] = Xv[k-1] + 6;
 }

/* sieve parameters: ts and qs are cutoffs, b is size of prime base and 
m2 is size of sieve*/
 k = Bitsize(N);
 b = 2 + (5*k)/4;
 bm = b + 3;
 
 if (k > 40) m =  50*k +  400;
 else   m = 100*k - 1600;
 m2 = m + m; 
 
 if (k > 40) ts = 196 + k/2;
 else   ts = 11*k - 24 - (k*k)/8;
 sq = ceil(sqrt(sq/m));
 
// Construct the prime base
L0:k = 2;
 i = 0;
 while (k < b) {
 p = Prm[i];
 i++;
 s = Modq(N,p);
 if (qrs(s,p) == 1) {
 pf[k] = p;
 sf[k] = mrt(s,p);
 lg[k] = Prm[i];
 k++;
 }
 i++;
 }
 pf[1] = 2;

// Construct quadratic polynomials as needed             
 while (Prm[i] < sq) i += 2;
 hi = i;
 qi = 0;

L1: do {if (hi < 616) sP = Prm[hi];
      else        sP = npr(sP);
      while ((sP&3) == 1) {
   hi += 2;
 if (hi < 616) sP = Prm[hi];
        else        sP = npr(sP);
 }
 sN = Modq(N,sP);
 hi += 2;
 } while (qrs(sN,sP) == -1);

 Set(T,sN);
 Dif(T,T,N); 
 Divq(T,sP); 
 d = Modq(T,sP); d *= sP; d += sN;
 sX = hrt(d,sP);
 Set(X,sX);

 e = (ulong) sP*sP;
 Set(B,e);
 Mulq(B,m);  Dif(B,B,X);
 Mul(C,B,B); Dif(C,C,N);
 Divq(C,sP); Divq(C,sP);
// Do the sieve for current polynoial
   if ((B[1] & 1) == 0) {i = 1; j = 0;}
   else        {i = 0; j = 1;}

   if      ((N[1] & 7) == 0)  ls = 21;
   else if ((N[1] & 3) == 0)  ls = 14;
   else                   ls =  7;
 
 while (i < m2) {
      lv[i] = ls; i += 2;
 lv[j] =  0; j += 2;
 }

 for (k = 2; k < b; k++)  {
      p = pf[k]; 
      s = sf[k]; 
      e = sX % p; 
      d = sP % p; d *= d; d %= p;  d = inv(d,p);

 if (e < s) e += p;
 i = e-s; i *= d; i %= p; i = m-i; i %= p;
 r1[k] = i;
  
 j = e+s; j *= d; j %= p; j = m-j; j %= p;
 r2[k] = j;
      
      if (i > j) {
        e = i; 
        i = j; 
        j = e;
      }
       ls = lg[k];
        while (j < m2) {
        lv[i] += ls; i += p;
        lv[j] += ls; j += p;
        }
        if (i < m2) lv[i] += ls;
    }
   
// Factor polynomial to find rows of matrix hm[qi,k]
    for (i = 0; i < m2; i++) 
 if (lv[i] > ts) {
      hp = hm[qi];
      d = (ulong) i*sP;
      Set(T,d); Mul(Q,T,T); Add(Q,Q,C);
      R[0] = B[0]; R[1] = B[1];
      R[2] = B[2]; R[3] = B[3];
      s = i+i; Mulq(R,s);
      if (Comp(Q,R) == 1) hp[0] = 0;
        else     hp[0] = 1;
      Dif(Q,Q,R);

      hp[1] = 0;
     while ((Q[1] & 1) == 0) {
         hp[1] += 1;
        Shiftr(Q);
      }

 k = b;
 while (Q[0] > 2) {
 k--;
 if (k == 1) goto L2;
        p = pf[k]; 
        j = i % p;
        if (j == r1[k] || j == r2[k]) {
        hp[k] = 1;
        Divq(Q,p);
              while (Modq(Q,p) == 0) {
               hp[k] += 1; 
               Divq(Q,p);
              }
          }
            else hp[k] = 0;
      }
        sQ = Unset(Q);
        while (sQ > 1) {
        k--;
        if (k == 1) goto L2;
        p = pf[k]; 
        j = i % p;
        if (j == r1[k] || j == r2[k]) {
        hp[k] = 1;
        sQ /= p;
               while (sQ%p == 0) {
                 hp[k] += 1;
                 sQ /= p;
               }
          }
 else hp[k] = 0;
 }

   while (k > 2) {
   k--; 
   hp[k] = 0;
   }
        Mulq(T,sP); 
        Dif(Xv[qi],T,B);
        Yv[qi] = sP;
        qi += 1;
        if (qi == bm) goto L3;
L2:;
   }
 goto L1;

// Row reduce gm = hm % 2 and find X^2 = Y^2 
// (mod N)
L3:k = N[0];
 g = (double) N[k];
 if (d > 1) g += N[k-1]/65536.0;
 if (d > 2) g += N[k-2]/4294967296.0;
 if (d > 3) g += 1/4294967296.0;
 
 for (r = 0; r < bm; r++) { 
      hp = hm[r]; gp = gm[r];
        for (c = 0; c < b; c++) 
            gp[c]  = hp[c] & 1;
        br = b + r;
 for (c = b; c < br;c++) gp[c] = 0;
 gp[br] = 1;
    }
   for (r = 0; r < bm; r++) {
        br = b + r;
        gp = gm[r];
        c = b - 1;  
        while (gp[c] == 0 && c>0) c--;

        if (c > 0 || gp[0] == 1) { 
            for (i = r+1; i<bm; i++) {
            gi = gm[i];
                if (gi[c] == 1) {
                 gi[c] = 0;
                 for (k = 0; k<c; k++)
 gi[k] ^= gp[k]; 
                  for (k=b; k<=br; k++)
 gi[k] ^= gp[k];
        }
        }
        }
        else {
            for (i = 1; i < b; i++)
 pv[i] = 0;
            Set(X,1); Set(Y,1);
            for (k = b; k <= br; k++) 
                if (gp[k] == 1) {
                 j = k - b;
                 Mul(X,X,Xv[j]); 
                 if (X[0]>12) Mod(X);
                 Mulq(Y,Yv[j]);  
                 if (Y[0]>12) Mod(Y);
                    for (i=1;i<b; i++)
               pv[i] += (long) hm[j][i];
                }
   Mod(X);
   
   k = pv[1] >> 1;
   while (k > 15) {
   for (i = Y[0]; i > 0; i--)
   Y[i+1] = Y[i];
   Y[1] = 0; 
   Y[0] += 1;
   k -= 16;
   if (Y[0] > 12) Mod(Y);
   }
   Mulq(Y,1 << k);
 
            for (i = 2; i < b; i++) {
            j = pv[i];
            if (j == 0) continue;
            if (j == 2) Mulq(Y,pf[i]);
                else {
                    T[1] = pf[i]; 
                    T[0] = 1;
                 while (j > 2) {
                 j -= 2; 
                 Mulq(T,pf[i]);
                 }
                 Mul(Y,Y,T);
                }
                if (Y[0] > 12) Mod(Y);
            }
            Mod(Y);
            
            Dif(T,X,Y); 
            Add(R,X,Y);
            if (Comp(N,R)<=0) Dif(R,R,N);
            if (T[0] == 0 || R[0] == 0) 
 continue;
 *p1 = Gcd(T); 
 *p2 = Gcd(R);
        return;  
   }
 } 
    bm += 5; 
    ts -= 2;
    goto L0;     
}
End of the function Factor64().

// Inverse of b modulo m  (Euclid’s algorithm)
inv
ulong inv(ulong b,ushort m) {

 register ulong u,v,t,n;  
 u = 1;
 v = 0;
 t = 1;
 n = m;
 
 for (;;) { 
 while ((b & 1) == 0) {   
 v <<= 1; 
 if (t & 1) t += m;
 t >>= 1; 
 b >>= 1;
 }
 if (b == 1) break;
 
 if (b > n) {
 b -= n;
 u += v;
 }
 else {
 n -= b;
 v += u;
 }
 while ((n & 1) == 0) {
 u <<= 1; 
 if (t & 1) t += m;
 t >>= 1; 
 n >>= 1;
 }
 }
 t *= u; t %= m;
 return t;
}

// Determine if n is square modulo p  (QR algorithm)
qrs
short qrs(ushort n,ushort p) {
 
     register short  j;
     register ushort n2,n4;

     if (n == 1) return 1;

     j = 1;
     for (;;) {
       n2 = n + n;
       n4 = n2 + n2;
         p %= n4;
         while (p > n) {
          if (n & 2)  j = -j;
   if (p > n2) p -= n2;
                 else   p = n2 - p;
         }
         if (p == 1) return j;
         n %= p;
         if (n == 1) return j;
      }
}

// Square root of n modulo p
mrt
ushort mrt(ushort n,ushort p) {

   register ulong q,s;
   register long  x,u,v;
   ulong m;
   long  t,r;

    if (n == 1) return 1;
    q = (p+1) >> 1;

 m = n;
 if ((q & 1) == 0) {
        q >>= 1;
        s   = 1;
        while (q > 0) {
            if (q & 1) {
            s *= m; s %= p;
            }
            m *= m; m %= p;
            q >>= 1;
        }
        if (s+s > p) s = p-s;
      return s;
   }

   s = 0;
   t = n;
   while (qrs(t,p) == 1) {
        s += 1;
        t += 1-s-s;
        if (t%p == 0) {
        if (s+s > p) s = p-s;
        return s;
      }
      if (t < p) t += p;
   }

    q >>= 1;
    n = t;
    r = s;
    u = 1;
    v = 1;
    while (q > 0) {
        x = n*u; x %= p; x *= u; x %= p;
        t = r*r; t %= p;
        if (t >= x) x = t-x;
        else        x = t-x+p;
        u *= r; u %= p; u += u; u %= p;
        r = x;
        if (q & 1) {
          x = n*u; x %= p; x *= v;
            x %= p;
            t = s*r; t %= p;
            if (t >= x) x = t-x;
            else        x = t-x+p;
            v *= r; v %= p; v += s*u; v %= p;
            s = x;
        }
        q >>= 1;
   
 if (s+s > p) s = p-s;
   return s;
}

// Square root of n modulo p^2 for p = 3 (mod 4)     
hrt
ulong hrt(ulong n,ushort p) {

  register ulong  s,t,x,y;
  ulong  m;

 m = n%p;
 s = m;
 y = 1;
 t = p-3; t >>= 2;
    while (t > 0) {
        if (t&1) {
        y *= s; y %= p;
        }
        s *= s; s %= p;
        t >>= 1;
   }
 x = m*y; x %= p;
 s = (ulong) p*p;
 t = x*x;
 if (n < t) n += s;
 n -= t; n /= p; n *= y; n %= p;
 t = p; t += 1; t >>= 1;
 n *= t; n %= p; n *= p; n += x;
 if (n+n > s) n = s-n;
 return n;
}

// Get next prime  
npr
ushort npr(ushort p) {

 register ushort d,s,k;

   do {p += 2;
   s = floor(sqrt(p));
   k = 0;
   d = 3;
 while (d <= s) {
 if (p%d == 0) break;
 k += 2;
 d = Prm[k];
 }
     }
     while (d <= s);
     return p;
}

// Addition  S = A + B
Add
void Add(Int S,Int A,Int B) {

 register ushort *pH, *pL;
 register ulong  t;
 register ushort c,k;
 ushort s,dH,dL;
 
 if (A[0] > B[0] ) {
 pH = A; dH = A[0]; pL = B; dL = B[0];
 }
 else {
 pH = B; dH = B[0]; pL = A; dL = A[0];
 }
 if (dL == 0) {
 if (S != pH) 
 for (k=0;k<=dH; k++) S[k]=pH[k];
 return;
 }
 
 k = 0;
 c = 0;
 while (k < dL) { 
 k++;
 t = (ulong) pH[k] + pL[k] + c;
 if (t >= 0x10000) {
 t -= 0x10000; 
 c = 1;
 }
 else c = 0;
 S[k] = t;
 }
 while (c == 1 && k < dH) {
 k++;
 s = pH[k];
 if (s == 0xFFFF) {
 S[k] = 0; 
 c = 1;
 }
 else {
 S[k] = s + 1;
 c = 0;
 }
 }
 while (k < dH) {
 k++;
 S[k] = pH[k];
 }
 if (c == 1) {
 dH += 1;
 S[dH] = 1;
 }
 S[0] = dH; 
}

// Difference  D = |A - B|
Dif
void Dif(Int D,Int A,Int B) {

 register ushort *pH, *pL;
 register long   t;
 register ushort c,k;
 ushort   s,dH,dL;
 short    e;

 k = A[0];
 if (k > B[0]) e = 1;
 else {
 if (k < B[0]) e = -1;
 else {
 while (A[k]==B[k] && k > 0) k--;
 if (k == 0) {
 D[0] = 0;
 return;
 }
 if (A[k] > B[k]) e = 1;
 else e = -1;
 }
 }
 if (e == 1) {
 pH = A; dH = A[0];
 pL = B; dL = B[0];
 }
 else {
 pH = B; dH = B[0];
 pL = A; dL = A[0];
 }
 if (dL == 0) {
 if (D != pH) 
 for (k=0;k<=dH;k++) D[k]=pH[k];
 return;
 }

 c = 0;
 k = 0;
 while (k < dL) {
 k++;   
 t = (long) pH[k] - pL[k] - c;
 if (t < 0) {
 t += 0x10000; 
 c = 1;
 }
 else c = 0;
 D[k] = t;
 }
 while (c == 1 && k < dH) {
 k++;
 s = pH[k];
 if (s == 0) {
 D[k] = 0xFFFF; 
 c = 1;
 }
 else {
 D[k] = s - 1;
 c = 0;
 }
 }
 while (k < dH) {
 k++;
 D[k] = pH[k];
 }
 k = dH;
 while (D[k] == 0 && k > 0) k--;
 D[0] = k;
}

// Multiply  P = A * B
Mul
void Mul(Int P,Int A,Int B) {

 register ushort *pB;
 register ulong  t;
 register ushort s,n,k;
 ushort d,j;
 
 if (A[0] == 0 || B[0] == 0) {
 P[0] = 0;
 return;
 }
 d = A[0] + B[0];
 for (k = 1; k <= d; k++) bufm[k]= 0;
 
 pB = B;
 j = 0;
 while (j < A[0]) {
 j++;
 s = A[j];
 k = 0;
 n = j;
 while (k < pB[0]) {
 k++;
 t = (ulong) s * pB[k];
 bufm[n] += t & 0xFFFF;
 n++;
 bufm[n] += t >> 16;
 } 
 }

 k = 1;
 while (k < d) {
 t = bufm[k];
 bufm[k] = t & 0xFFFF;
 k++;
 bufm[k] += t >> 16;
 }
 if (bufm[d] == 0) d -= 1;
 
 for (k = 1; k<=d; k++) P[k]=bufm[k];
 P[0] = d;
}

// Quick multiply  A = A * b
Mulq
void Mulq(Int A,ushort b) {

 register ushort *pA;
 register ulong  t,w,c;
 register ushort d,k;
 
 d = A[0];
 if (d == 0) return;
 
 pA = A;
 if (b == 0) {
 pA[0] = 0;
 return;
 }
 
 c = 0;
 k = 0;
 while (k < d) {
 k++;
 t = (ulong) pA[k] * b;
 w = (t & 0xFFFF) + c;
 c = t >> 16;
 if (w >= 0x10000) {
 w -= 0x10000;
 c += 1;
 }
 pA[k] = w;
 }
 if (c > 0) {
 d += 1;
 pA[d] = c;
 }
 A[0] = d;
}

// Modulo  R = R % N 
Mod
void Mod(Int R) {

 register ulong  t;
 register long   w;
 register ushort k,j;
 ushort d,n,tL,tH;
 ulong  c;
 short  e;
 double f;

 d = N[0];
 n = R[0];
 if (d > n) return;

 for (k = 1; k<=n; k++) buf[k]=R[k];

 while (n >= d) {
 
 f = buf[n]*65536.0;
 if (n > 1) f += buf[n-1];
 t = f/g;
 e = n - d;
 if (e > 0) {
 tL = t & 0xFFFF;
 tH = t >> 16;
 }
 else {
 tL = t >> 16;
 if (tL == 0) {
 if (t < 0xFFFF) break;
 else {
 k = d;
 while (buf[k] == N[k] 
 && k > 0) k--;
 if (k > 0 && buf[k] < N[k])
 break;
 tL = 1;
 }
 }
 tH = 0;
 e = 1;
 }

 c = 0;
 j = e;
 for (k = 1; k <= d; k++) {
 t = (ulong) N[k]*tL + c;
 w = (ulong) buf[j]-(t & 0xFFFF);
 t >>= 16;
 if (w < 0) {
 buf[j] = 0x10000 + w;
 t += 1;
 }
 else buf[j] = w;
 j++;
 w = (ulong) buf[j] - t;
 if (w < 0) {
 buf[j] = 0x10000 + w;
 c = 0x10000;
 }
 else {
 buf[j] = w;
 c = 0;
 }
 }
 
 if (tH > 0) {
 c = 0;
 j = e + 1;
 for (k = 1; k <= d; k++) {
 t = (ulong) N[k]*tH + c;
 w = (ulong) buf[j]-(t&0xFFFF);
 t >>= 16;
 if (w < 0) {
 buf[j] = 0x10000 + w;
 t += 1;
 }
 else buf[j] = w;

 if (k == d) break;
 j++;
 w = (ulong) buf[j] - t;
 if (w < 0) {
 buf[j] = 0x10000 + w;
   c = 0x10000;
 }
 else {
 buf[j] = w;
 c = 0;
 }
 }
 }
 while (buf[n] == 0 && n > 0) n--;
 if (n == d && buf[d] < N[d]) break;
 }

 for (k = 1; k <= n; k++) R[k] = buf[k];
 R[0] = n;
}

// Quick modulo  A = A % m
Modq
ushort Modq(Int A,ushort m) {

 register ulong  z,t;
 ushort n,e;

 n = A[0];
 if (n == 0) return 0;

 for (e = 1; e <= n; e++) buf[e] = A[e];

 while (n > 1) {
 e = n - 1;
 t = ((ulong) buf[n] << 16) + buf[e];
 z = t/m; 
 t -= z*m;
 buf[e] = t;
 if (t == 0) do e--;
 while (buf[e] == 0 && e > 0);
 n = e;
 }
 return buf[1] % m;
}

// Quick divide  A = A / m
Divq
void Divq(Int A,ushort m) {

 register ulong  z,t;
 ushort n,e;

 n = A[0];
 if (n == 0) return;

 for (e = 1; e <= n; e++) {
 buf[e] = A[e];
 A[e] = 0;
 }

 while (n > 1) {
 e = n - 1;
 t = ((ulong) buf[n] << 16) + buf[e];
 z = t/m;

 if (z < 0x10000) A[e] = z;
 else {
 A[e] = z & 0xFFFF;
 A[n] = z >> 16;
 }
 t -= z*m;
 buf[e] = t;
 if (t == 0) do e--;
 while (buf[e] == 0 && e > 0);
 n = e;
 }
 n = buf[1];
 if (n >= m) A[1] = n/m;
 if (A[A[0]] == 0) A[0] -= 1;
}

// Shift right  X = X >> 1
Shiftr
void Shiftr(Int X) {

 register ulong  t;
 register ushort i,j,c,d;

 d = X[0];
 if (d == 0) return;

 i = 0;
 j = 1;
 c = X[1] >> 1;
 while (j < d) {
 j++;
 t = (ulong) X[j] << 15;
 t += c;
 i++;
 X[i] = t & 0xFFFF;
 c = t >> 16;
 }
 if (c == 0) d -= 1;
 else X[d] = c;
 X[0] = d;
}

// Compare : {+1  0  -1} as {X > Y   X == Y   X < Y}
Comp
short Comp(Int X,Int Y) {

 register ushort d;

 d = X[0];
 if (d > Y[0]) return  1;
 if (d < Y[0]) return -1;

 while (X[d] == Y[d] && d > 0) d--;

 if (d == 0) return 0;
 if (X[d] > Y[d]) return 1;
 return -1;
}

// Convert from unsigned long to Int 
Set
void Set(Int X, ulong n) {

 if (n == 0) {
 X[0] = 0;
 return;
 }
 if (n < 0x10000) {
 X[0] = 1;
 X[1] = n;
 return;
 }
 X[0] = 2;
 X[1] = n & 0xffff;
 X[2] = n >> 16;
}

// Convert from Int to unsigned long 
Unset
ulong Unset(Int X) {

 register ulong  n;
 register ushort d;

 d = X[0];
 if (d == 0) return 0;
 n = (ulong) X[1];
 if (d == 1) return n;
 return (n + ((ulong) X[2] << 16));
}

// Number of Bits in X
Bitsize
ulong Bitsize(Int X) {

 register ushort d,t;
 register ulong  n;

 d = X[0];
 if (d == 0) return 0;

 n = (ulong) d << 4;
 t = 0x8000;
 while ((t & X[d]) == 0) {
 n  -= 1;
 t >>= 1;
 }
 return n;
}

// The greatest common divisor of A and N                
Gcd
ulong Gcd(Int A) {

 register long k;

 for (k = 0; k < 5; k++) buf[k] = N[k];

 while ((A[1]&1) == 0) Shiftr(A);

 for (;;)
 switch (Comp(A,buf)) {
 case  1 :  Dif(A,A,buf);
 while ((A[1]&1) == 0) Shiftr(A);
 break;
 case -1 :  Dif(buf,buf,A);
 while ((buf[1]&1) == 0) Shiftr(buf);
 break;
 case  0 :  return Unset(A);
 }
}
End of file Factor64.c

Factor64.h


/* 4 */
// Header file for 64 bit factorization program.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define  ulong    unsigned long
#define  ushort  unsigned short
typedef  ushort *  Int;

ulong    inv(ulong,ushort);
short    qrs(ushort,ushort);
ushort   mrt(ushort,ushort);
ulong    hrt(ulong,ushort);
ushort   npr(ushort);

void     Add(Int,Int,Int);
void     Dif(Int,Int,Int);
void     Mul(Int,Int,Int);
void     Mulq(Int,ushort);
void     Mod(Int);
ushort   Modq(Int,ushort);
void     Divq(Int,ushort);
void     Shiftr(Int);
short    Comp(Int,Int);
void     Set(Int,ulong);
ulong    Unset(Int);
ulong    Bitsize(Int);
ulong    Gcd(Int);

ushort   buf[20], N[5];
ulong    bufm[20];
double   g;

// Odd primes below 0x800 (and their 10*logs)
ushort Prm[] = {
   3, 11,   5, 17,   7, 20,  11, 24,  13, 26,
  17, 29,  19, 30,  23, 32,  29, 34,  31, 35,
  37, 37,  41, 38,  43, 38,  47, 39,  53, 40,
  59, 41,  61, 42,  67, 43,  71, 43,  73, 43,
  79, 44,  83, 45,  89, 45,  97, 46, 101, 46,
 103, 46, 107, 46, 109, 46, 113, 47, 127, 48,
 131, 48, 137, 49, 139, 49, 149, 50, 151, 50,
 157, 50, 163, 50, 167, 51, 173, 51, 179, 51,
 181, 51, 191, 52, 193, 52, 197, 52, 199, 52,
 211, 53, 223, 54, 227, 54, 229, 54, 233, 54,
 239, 54, 241, 54, 251, 55, 257, 55, 263, 55,
 269, 55, 271, 56, 277, 56, 281, 56, 283, 56,
 293, 56, 307, 57, 311, 57, 313, 57, 317, 57,
 331, 58, 337, 58, 347, 58, 349, 58, 353, 58,
 359, 58, 367, 59, 373, 59, 379, 59, 383, 59,
 389, 59, 397, 59, 401, 59, 409, 60, 419, 60,
 421, 60, 431, 60, 433, 60, 439, 60, 443, 60,
 449, 61, 457, 61, 461, 61, 463, 61, 467, 61,
 479, 61, 487, 61, 491, 61, 499, 62, 503, 62,
 509, 62, 521, 62, 523, 62, 541, 62, 547, 63,
 557, 63, 563, 63, 569, 63, 571, 63, 577, 63,
 587, 63, 593, 63, 599, 63, 601, 63, 607, 64,
 613, 64, 617, 64, 619, 64, 631, 64, 641, 64,
 643, 64, 647, 64, 653, 64, 659, 64, 661, 64,
 673, 65, 677, 65, 683, 65, 691, 65, 701, 65,
 709, 65, 719, 65, 727, 65, 733, 65, 739, 66,
 743, 66, 751, 66, 757, 66, 761, 66, 769, 66,
 773, 66, 787, 66, 797, 66, 809, 66, 811, 66,
 821, 67, 823, 67, 827, 67, 829, 67, 839, 67,
 853, 67, 857, 67, 859, 67, 863, 67, 877, 67,
 881, 67, 883, 67, 887, 67, 907, 68, 911, 68,
 919, 68, 929, 68, 937, 68, 941, 68, 947, 68,
 953, 68, 967, 68, 971, 68, 977, 68, 983, 68,
 991, 68, 997, 69,1009, 69,1013, 69,1019, 69,
1021, 69,1031, 69,1033, 69,1039, 69,1049, 69,
1051, 69,1061, 69,1063, 69,1069, 69,1087, 69,
1091, 69,1093, 69,1097, 70,1103, 70,1109, 70,
1117, 70,1123, 70,1129, 70,1151, 70,1153, 70,
1163, 70,1171, 70,1181, 70,1187, 70,1193, 70,
1201, 70,1213, 71,1217, 71,1223, 71,1229, 71,
1231, 71,1237, 71,1249, 71,1259, 71,1277, 71,
1279, 71,1283, 71,1289, 71,1291, 71,1297, 71,
1301, 71,1303, 71,1307, 71,1319, 71,1321, 71,
1327, 71,1361, 72,1367, 72,1373, 72,1381, 72,
1399, 72,1409, 72,1423, 72,1427, 72,1429, 72,
1433, 72,1439, 72,1447, 72,1451, 72,1453, 72,
1459, 72,1471, 72,1481, 73,1483, 73,1487, 73,
1489, 73,1493, 73,1499, 73,1511, 73,1523, 73,
1531, 73,1543, 73,1549, 73,1553, 73,1559, 73,
1567, 73,1571, 73,1579, 73,1583, 73,1597, 73,
1601, 73,1607, 73,1609, 73,1613, 73,1619, 73,
1621, 73,1627, 73,1637, 74,1657, 74,1663, 74,
1667, 74,1669, 74,1693, 74,1697, 74,1699, 74,
1709, 74,1721, 74,1723, 74,1733, 74,1741, 74,
1747, 74,1753, 74,1759, 74,1777, 74,1783, 74,
1787, 74,1789, 74,1801, 74,1811, 75,1823, 75,
1831, 75,1847, 75,1861, 75,1867, 75,1871, 75,
1873, 75,1877, 75,1879, 75,1889, 75,1901, 75,
1907, 75,1913, 75,1931, 75,1933, 75,1949, 75,
1951, 75,1973, 75,1979, 75,1987, 75,1993, 75,
1997, 75,1999, 76,2003, 76,2011, 76,2017, 76,
2027, 76,2029, 76,2039, 76,   0,  0 };

End of file Factor64.h







  
 
AAPL
$119.00
Apple Inc.
+1.40
MSFT
$47.75
Microsoft Corpora
+0.28
GOOG
$540.37
Google Inc.
-0.71

MacTech Search:
Community Search:

Software Updates via MacUpdate

Skype 7.2.0.412 - Voice-over-internet ph...
Skype allows you to talk to friends, family and co-workers across the Internet without the inconvenience of long distance telephone charges. Using peer-to-peer data transmission technology, Skype... Read more
HoudahSpot 3.9.6 - Advanced file search...
HoudahSpot is a powerful file search tool built upon MacOS X Spotlight. Spotlight unleashed Create detailed queries to locate the exact file you need Narrow down searches. Zero in on files Save... Read more
RapidWeaver 6.0.3 - Create template-base...
RapidWeaver is a next-generation Web design application to help you easily create professional-looking Web sites in minutes. No knowledge of complex code is required, RapidWeaver will take care of... Read more
iPhoto Library Manager 4.1.10 - Manage m...
iPhoto Library Manager lets you organize your photos into multiple iPhoto libraries. Separate your high school and college photos from your latest summer vacation pictures. Or keep some photo... Read more
iExplorer 3.5.1.9 - View and transfer al...
iExplorer is an iPhone browser for Mac lets you view the files on your iOS device. By using a drag and drop interface, you can quickly copy files and folders between your Mac and your iPhone or... Read more
MacUpdate Desktop 6.0.3 - Discover and i...
MacUpdate Desktop 6 brings seamless 1-click installs and version updates to your Mac. With a free MacUpdate account and MacUpdate Desktop 6, Mac users can now install almost any Mac app on macupdate.... Read more
SteerMouse 4.2.2 - Powerful third-party...
SteerMouse is an advanced driver for USB and Bluetooth mice. It also supports Apple Mighty Mouse very well. SteerMouse can assign various functions to buttons that Apple's software does not allow,... Read more
iMazing 1.1 - Complete iOS device manage...
iMazing (was DiskAid) is the ultimate iOS device manager with capabilities far beyond what iTunes offers. With iMazing and your iOS device (iPhone, iPad, or iPod), you can: Copy music to and from... Read more
PopChar X 7.0 - Floating window shows av...
PopChar X helps you get the most out of your font collection. With its crystal-clear interface, PopChar X provides a frustration-free way to access any font's special characters. Expanded... Read more
OneNote 15.4 - Free digital notebook fro...
OneNote is your very own digital notebook. With OneNote, you can capture that flash of genius, that moment of inspiration, or that list of errands that's too important to forget. Whether you're at... Read more

Latest Forum Discussions

See All

Raby (Games)
Raby 1.0.3 Device: iOS Universal Category: Games Price: $2.99, Version: 1.0.3 (iTunes) Description: ***WARNING - Raby runs on: iPhone 5, iPhone 5C, iPhone 5S, iPhone 6, iPhone 6 Plus, iPad Mini Retina, iPad Mini 3, iPad 4, iPad Air,... | Read more »
Oddworld: Stranger's Wrath (Games)
Oddworld: Stranger's Wrath 1.0 Device: iOS Universal Category: Games Price: $5.99, Version: 1.0 (iTunes) Description: ** PLEASE NOTE: Oddworld Stranger's Wrath requires at least an iPhone 4S, iPad 2, iPad Mini or iPod Touch 5th gen... | Read more »
Bounce On Back (Games)
Bounce On Back 1.0 Device: iOS Universal Category: Games Price: $2.99, Version: 1.0 (iTunes) Description: | Read more »
Dwelp (Games)
Dwelp 1.0 Device: iOS Universal Category: Games Price: $.99, Version: 1.0 (iTunes) Description: === 50% off for a limited time, to celebrate release === Dwelp is an elegant little puzzler with a brand new game mechanic. To complete a... | Read more »
Make Way for Fat Chicken, from the Maker...
Make Way for Fat Chicken, from the Makers of Scrap Squad Posted by Jessica Fisher on November 26th, 2014 [ permalink ] Relevant Games has announced they will be releasing their reverse tower defense game, | Read more »
Tripnary Review
Tripnary Review By Jennifer Allen on November 26th, 2014 Our Rating: :: TRAVEL BUCKET LISTiPhone App - Designed for the iPhone, compatible with the iPad Want to create a travel bucket list? Tripnary is a fun way to do exactly that... | Read more »
Ossian Studios’ RPG, The Shadow Sun, is...
Ossian Studios’ RPG, The Shadow Sun, is Now Available for $4.99 Posted by Jessica Fisher on November 26th, 2014 [ permalink ] Universal App - Designed for iPhone and iPad | Read more »
Mmmm, Tasty – Having the Angry Birds for...
The very first Angry Birds debuted on iOS back in 2009. When you sit back and tally up the number of Angry Birds games out there and the impact they’ve had on pop culture as a whole, you just need to ask yourself: “How would the birds taste... | Read more »
Rescue Quest Review
Rescue Quest Review By Jennifer Allen on November 26th, 2014 Our Rating: :: PATH BASED MATCH-3Universal App - Designed for iPhone and iPad Guide a wizard to safety by matching gems. Rescue Quest might not be an entirely original... | Read more »
You Can Play the Final Chapter of Lone W...
You Can Play the Final Chapter of Lone Wolf: Dawn Over V’taag Right Now Posted by Jessica Fisher on November 26th, 2014 [ permalink ] Universal App - Designed for iPhone and iPad | Read more »

Price Scanner via MacPrices.net

BEVL Releases Dock Tailored for iPhone 6 and...
Seattle based BEVL has released their first product: an iPhone dock that is divergent in build quality, rock-solid function and visual simplicity to complement the iPhone. BEVL is now accepting... Read more
Black Friday: $150 off 13-inch Retina MacBook...
 Best Buy has 13-inch 2.6GHz Retina MacBook Pros on sale for $150 off MSRP on their online store as part of their Black Friday sale. Choose free shipping or free local store pickup (if available).... Read more
Black Friday: $300 off 15-inch Retina MacBook...
 B&H Photo has the new 2014 15″ Retina MacBook Pros on sale for $300 off MSRP as part of their Black Friday sale. Shipping is free, and B&H charges NY sales tax only: - 15″ 2.2GHz Retina... Read more
Black Friday: Up to $140 off MacBook Airs, fr...
 B&H Photo has 2014 MacBook Airs on sale for up to $140 off MSRP as part of their Black Friday sale. Shipping is free, and B&H charges NY sales tax only: - 11″ 128GB MacBook Air: $799 $100... Read more
Black Friday: 13-inch 2.5GHz MacBook Pro on s...
 Best Buy has the 13″ 2.5GHz MacBook Pro on sale for $899.99 on their online store as part of their Black Friday sale. Choose free shipping or free instant local store pickup (if available). Their... Read more
2014 1.4GHz Mac mini on sale for $449, save $...
 B&H Photo has the new 1.4GHz Mac mini on sale for $449.99 including free shipping plus NY tax only. Their price is $50 off MSRP, and it’s the lowest price available for this new model. Adorama... Read more
Early Black Friday pricing on 27-inch 5K iMac...
 B&H Photo continues to offer Black Friday sale prices on the 27″ 3.5GHz 5K iMac, in stock today and on sale for $2299 including free shipping plus NY sales tax only. Their price is $200 off MSRP... Read more
Early Black Friday sale prices on iPad Air 2,...
 MacMall is discounting iPad Air 2s by up to $75 off MSRP as part of their Black Friday sale. Shipping is free: - 16GB iPad Air WiFi: $459 $40 off - 64GB iPad Air WiFi: $559 $40 off - 128GB iPad Air... Read more
Early Black Friday MacBook Air sale prices, $...
 MacMall has posted early Black Friday MacBook Air sale prices. Save $101 on all models for a limited time: - 11″ 1.4GHz/128GB MacBook Air: $798 - 11″ 1.4GHz/256GB MacBook Air: $998 - 13″ 1.4GHz/... Read more
Why iPhone 6 Tablet/Laptop Cannibalization Is...
247wallst.com blogger Douglas A. McIntyre noted last week that according to research posted on the Applovin blog site the iPhone 6 is outselling the iPhone 6 Plus by a wide margin . Hardly a surprise... Read more

Jobs Board

*Apple* Solutions Consultant (ASC) - Apple (...
**Job Summary** The ASC is an Apple employee who serves as an Apple brand ambassador and influencer in a Reseller's store. The ASC's role is to grow Apple Read more
Senior Event Manager, *Apple* Retail Market...
…This senior level position is responsible for leading and imagining the Apple Retail Team's global event strategy. Delivering an overarching brand story; in-store, Read more
*Apple* Retail - Multiple Positions (US) - A...
Sales Specialist - Retail Customer Service and Sales Transform Apple Store visitors into loyal Apple customers. When customers enter the store, you're also the Read more
*Apple* Solutions Consultant (ASC) - Apple (...
**Job Summary** The ASC is an Apple employee who serves as an Apple brand ambassador and influencer in a Reseller's store. The ASC's role is to grow Apple Read more
*Apple* Solutions Consultant (ASC) - Apple (...
**Job Summary** The ASC is an Apple employee who serves as an Apple brand ambassador and influencer in a Reseller's store. The ASC's role is to grow Apple Read more
All contents are Copyright 1984-2011 by Xplain Corporation. All rights reserved. Theme designed by Icreon.