home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / magazine / drdobbs / 1986 / 08 / letter.aug next >
Text File  |  1986-08-31  |  24KB  |  733 lines
  1. Listing One:
  2. ;DeSmet function isqrt(source);
  3. ;                    ...source is a long integer (32 bit)...
  4. ;Returns square root of source in ax, using Newton's method; see Scanlon's
  5. ;8086 book for similar function (Scanlon's is not sufficiently general, and has
  6. ;an error)...
  7. ;
  8. ;REGISTER USAGE       cx:bx   ...stores copy of 32-bit source throughout...
  9. ;                     di      ...``last''estimate of isqrt(source)...
  10. ;                     si      ...current estimate of isqrt(source)...
  11. ;                     dx:ax   ...used to divide source by di...
  12.  
  13. cseg
  14. public isqrt_
  15. isqrt_:          push bp
  16.           mov bp,sp
  17.                 mov bx, [bp+4]       ;Store a copy of source in cx:bx;
  18.                       mov cx, [bp+6]       ;cx:bx preserved till almostdone...
  19. ;-----Start block to determine initial estimate; base estimate on most---------
  20. ;-----significant non-zero byte of cx:bx---------------------------------------
  21.           cmp ch,0               ;Note 46341 covers largest positive
  22.                 je test_cl           ;long that most compilers will pass.
  23.            mov di,46341         ;If you use 32-bit unsigned, replace
  24.              jmp load_si          ;with 65535, and if cx>=FFFEH, jump
  25.           ;------------        ;out and set result= 65535.
  26.   test_cl:      cmp cl,0
  27.           je text_bh
  28.           mov di,2896
  29.           jmp load_si
  30.           ;------------
  31.   test_bh:      cmp bh,0
  32.           je its_bl
  33.           mov di,181
  34.           jmp load_si
  35.           ;------------
  36.   its_bl:       mov di,8
  37.   load_si:      mov si,di
  38. ;-----End block to determine initial estimate----------------------------------
  39.  
  40. ;-----Begin loop to refine the estimate----------------------------------------
  41.   refine:       mov dx,cx            ;Load dx:ax pair with source in prep
  42.           mov ax,bx       ;for divide by di=estimate...
  43.           div di
  44.           ;------Block to average quotient and last estimate-------
  45.           shr ax,1       ;We can't just add di,ax then
  46.           adc di,0       ;shr di,1 because sum of di and ax
  47.           shr di,1       ;may exceed 65535...
  48.           ;-----------
  49.           sub si,di       ;Obtain difference betw. old (si)
  50.           jz almostdone     ;and new estimates; if 0, we're
  51.           cmp si,1       ;almost done...Else if diff. is
  52.           je almostdone     ;1 or -1 we're almost done...
  53.           cmp si,-1
  54.           je almostdone
  55.           mov si,di       ;Store current value di in si as
  56.                   ;``old'' estimate for next iteration.
  57.           jmp refine
  58. ;-----End loop to refine the estimate------------------------------------------
  59.   almostdone:   mov ax,di       ;Check to see if estimate*estimate
  60.           mul di       ;is less than cx:bx; this step is
  61.           sub bx,ax       ;for fussbudgets who demand final
  62.           sbb cx,dx       ;integer sqrt be < real sqrt; ditch
  63.           jns done       ;it and save approx. 60 clocks...
  64.           dec di       ;If product >cx:bx, subtract 1...
  65.   done:        mov ax,di       ;DeSmet looks for result in ax...
  66.           mov sp,bp
  67.           pop bp
  68.           ret
  69.  
  70. /* Test driver for isqrt()... */
  71.  
  72. main()
  73. {
  74. long start,i,NumSqrts;
  75. printf(``ENTER start: '');
  76. scanf(``%D'',&start);
  77. printf(``ENTER NumSqrts: '');
  78. scanf(``%D'',&NumSqrts);
  79. printf(``isqrt(%D)= %u\n'',start,isqrt(start));
  80. putchar('\007');
  81. for(i=start;i<start+NumSqrts;++i)isqrt(i); /* Remove this isqrt() to find */
  82. putchar('\007');          /* time taken by loop itself...*/
  83. }
  84.  
  85.  
  86. /* Another short driver; this one better for verifying algorithm */
  87.  
  88. main()
  89. {
  90. long source;
  91. unsigned result;
  92. printf(``ENTER # for sqrt; negative exits\n'');
  93. while(printf(``ENTER #: ''),scanf(``%D'',&source),source>=OL){
  94.   result=isqrt(source);
  95.   printf(``result= SQRT(%D)= %u\n'',source,result);
  96.   printf(``result*result= %D\n'',((long)result)*((long)result));
  97.   printf(``(result+1)*(result+1)= %D\n\n'',(result+1L)*(result+1L));
  98.   }
  99. }
  100. -----------------------------------------------------------------------------
  101. Listing Two: Bit-Shifting Method (Slower)
  102.  
  103. ;DeSmet function 1sqrt(); takes a long (32-bit) integer as an argument, returns
  104. ;a short (16-bit) integer square root.  Function result returned in ax.
  105. ;Modified after 68000 code published in DDJ #109, Nov. 1985, p. 90. Comments
  106. ;give roughly analogous 68000 instructions; correspondence with 68000 registers
  107. ;is:
  108. ;   D0 = sp:si  (initially holds argument)
  109. ;  D1 = dx:di  (Error term)
  110. ;  D2 = ax    (Running estimate)
  111. ;  D3 = bx    (High bracket; may exceed 16 bits on last iteration)
  112. ;  D4 = cx    (Loop counter)
  113. ;
  114. ;Note sp is used as a general register, so can't push or pop between
  115. ;``mov bp,sp'' and ``mov sp,bp''...also, we must disable DOS's timer
  116. ;interrupt, because it manipulates the stack every 55 milliseconds...
  117.  
  118. cseg
  119. public  lsqrt_
  120.  
  121.   lsqrt_:        push bp
  122.           mov bp,sp
  123.  
  124. ;------Now can address stack; 4-byte argument starts at bp+4------------------
  125.            cli       ;Clear interrupts (lock out)...
  126.           mov si,word [bp+4]   ;Place argument in sp:si= ``DO''...
  127.           mov sp,word [bp+6] 
  128.   
  129.           xor di,di       ;Zero ``D1'' and ``D2''...
  130.           mov dx,di
  131.           mov ax,di
  132.           mov cx,16       ;Loop counter cx= ``D4''...
  133.  
  134.   sqrt1:        shl si,1
  135.           rcl sp,1       ;``asl.l  #1,D0''...
  136.           rcl di,1       ;``roxl.l #1,D1''...
  137.           rcl dx,1
  138.  
  139.           shl si,1
  140.           rcl sp,1       ;Repeat shift and rotate...
  141.           rcl di,1
  142.           rcl dx,1
  143.  
  144.           shl ax,1       ;``asl.l #1,D2''...
  145.  
  146.           mov bx,ax       ;``move.1 D2,D3''...
  147.  
  148.           shl bx,1       ;``asl.l #1,D3''...
  149.           jc is_carry     ;Jump out of loop if new ``D3''exceeds
  150.              ;16 bits (may happen on last
  151.              ;iteration)...
  152.           cmp dx,0
  153.           jb sqrt2       ;``cmp.l D3,D1''...
  154.           ja past1         ;``bls  sqrt2''...
  155.           cmp di,bx
  156.           jbe sqrt2
  157.  
  158.   past1:        inc ax       ;``addq.1 #1,D2''...
  159.           inc bx       ;``addq.1 #1,D3''...
  160.  
  161.           sub di,bx       ;``sub.1 D3,D1''...
  162.           sbb dx,0
  163.  
  164.   sqrt2:        loop sqrt1     ;``dbra D4,sqrt1''...
  165.  
  166.           jmp past3       ;Skip 'is_carry' block if we finished
  167.              ;loop through 16 iterations (no jc,
  168.              ;``D3'' stayed <17 bits)...
  169. ;------------------------------------------------------------------------------
  170.   is_carry;     cmp dx,0001H         ;If we got here, there was a carry
  171.           ja past2             ;for shl bx,1 on last iteration of
  172.                 jb past3             ;loop; compare upper word ``D1''
  173.           cmp di,bx            ;against upper word ``D3'', which is
  174.           jbe past3       ;now 0001H; then compare lower words
  175.              ;of ``D1'' and ``D3''...
  176.   past2:        inc ax               ;``addq.1 #1,D2''...
  177. ;------------------------------------------------------------------------------
  178.   past3:        sti         ;Restore interrupts...
  179.           mov sp,bp       ;Restore frame, function result
  180.           pop bp       ;is already in ax....
  181.           ret
  182.  
  183. *************************************************************************
  184. Listing Three:
  185.  
  186.  
  187. *                  *
  188. *  Integer Square Root (32 to 16 bit).        *
  189. *                  *
  190. *  (Exact method, not approximate).        *
  191. *                  *
  192. *  Call with:              *
  193. *    DO.L = Unsigned number.          *
  194. *                  *
  195. *  Returns:              *
  196. *    DO.L = SQRT(DO.L)          *
  197. *                  *
  198. *  Notes:  Result fits in DO.W, but is valid in longword.    *
  199. *    Takes from 122 to 1272 cycles (including rts).    *
  200. *    Averages 610 cycles measured over first 65535 roots.  *
  201. *    Averages 1104 cycles measured over first 500000 roots.  *
  202. *                  *
  203. *************************************************************************
  204.  
  205.   .globl lsqrt
  206. *      Cycles
  207. lsqrt  tst.l d0  (4)  ; Skip doing zero.
  208.   beg.s done   (10/8)
  209.   cmp.l #$10000,d0 (14)  ; If is a longword, use the long routine.
  210.   bhs.s glsqrt  (10/8)
  211.   cmp.w #625,d0  (8)  ; Would the short word routine be quicker?
  212.   bhi.s gsqrt  (10/8)  ; No, use general purpose word routine.
  213. *        ; Otherwise fall into special routine.
  214. *
  215. * For speed, we use three exit points.
  216. * This is cheesy, but this is a speed-optimized subroutine!
  217.  
  218. *************************************************************************
  219. *                  *
  220. *  Faster Integer Square Root (16 to 8 bit).  For small arguments. *
  221. *                  *
  222. *  (Exact method, not approximate).        *
  223. *                  *
  224. *  Call with:              *
  225. *    DO.W = Unsigned number.          *
  226. *                  *
  227. *  Returns:              *
  228. *    DO.W = SQRT(DO.W)          *
  229. *                  *
  230. *  Notes:  Result fits in DO.B, but is valid in word.    *
  231. *    Takes from 72 (d0=1) to 504 (d0=625) cycles    *
  232. *    (including rts).          *
  233. *                  *
  234. *  Algorithm supplied by Motorola.          *
  235. *                  *
  236. *************************************************************************
  237.  
  238. * Use the theorem that a perfect square is the sum of the first 
  239. * sqrt(arg) number of odd integers.
  240.  
  241. *       Cycles
  242.   move.w d1,-(sp)  (8)
  243.   move.w #-1,d1  (8)
  244. qsqrt1  addq.w #2,d1  (4)
  245.   sub.w d1,d0  (4)
  246.   bpl qsqrt1  (10/8)
  247.   asr.w #1,d1  (8)
  248.   move.w d1,d0  (4)
  249.   move.w (sp)+,d1 (12)
  250. done  rts    (16)
  251.  
  252. *************************************************************************
  253. *                  *
  254. *  Integer Square Root (16 to 8 bit).        *
  255. *                  *
  256. *  (Exact method, not approximate).        *
  257. *                  *
  258. *  Call with:              *
  259. *    DO.W = Unsigned number.          *
  260. *                  *
  261. *  Returns:              *
  262. *    DO.L = SQRT(DO.W)          *
  263. *                  *
  264. *  Uses:  D1-D4 as temporaries --          *
  265. *    D1 = Error term;          *
  266. *    D2 = Running estimate;          *
  267. *    D3 = High bracket;          *
  268. *    D4 = Loop counter          *
  269. *                  *
  270. *  Notes: Result fits in DO.B, but is valid in word.    *
  271. *                  *
  272. *    Takes from 512 to 592 cycles (including rts).    *
  273. *                  *
  274. *    Instruction times for branch-type instructions     *
  275. *    listed as (X/Y) are for (taken/not taken).    *
  276. *                  *
  277. *************************************************************************
  278.  
  279. *      Cycles
  280. gsqrt  movem.w d1-d4,-(sp) (24)
  281.   move.w #7,d4   (8)  ; Loop count (bits-1 of result).
  282.   clr.w d1  (4)  ; Error term in D1.
  283.   clr.w d2  (4)
  284. sqrt1  add.w d0,d0  (4)  ; Get 2 leading bits a time and add
  285.   addx.w d1,d1  (4)  ; into Error term for interpolation.
  286.   add.w d0,d0  (4)  ; (Classical method, easy in binary).
  287.   addx.w d1,d1   (4)
  288.   add.w d2,d2   (4)  ; Running estimate * 2.
  289.   move.w d2,d3  (4)
  290.   add.w d3,d3   (4)
  291.   cmp.w d3,d1  (4)
  292.   bls.s sqrt2  (10/8)  ; New Error term > 2* Running estimate?
  293.   addq.w #1,d2  (4)  ; Yes, we want a `1' bit then.
  294.   addq.w #1,d3  (4)  ; Fix up new Error term.
  295.   sub.w d3,d1   (4)
  296. sqrt2  dbra d4,sqrt1  (10/14)  ; Do all 8 bit-pairs.
  297.   move.w d2,d0  (4)
  298.   movem.w (sp)+,d1-d4 (28)
  299.   rts    (16)
  300.  
  301. *************************************************************************
  302. *                  *
  303. *    Integer Square Root (32 to 16 bit).        *
  304. *                  *
  305. *  (Exact Method, not approximate).        *
  306. *                  *
  307. *  Call with:              *
  308. *    DO.L = Unsigned number.          *
  309. *                  *
  310. *  Returns:              *
  311. *    DO.L = SQRT(DO.L)          *
  312. *                  *
  313. *  Uses:  D1-D4 as temporaries --          *
  314. *    D1 = Error term;          *
  315. *    D2 = Running estimate;          *
  316. *    D3 = High bracket;          *
  317. *    D4 = Loop counter.          *
  318. *                  *
  319. *  Notes:  Result fits in DO.W, but is valid in longword.    *
  320. *                  *
  321. *    Takes from 1080 to 1236 cycles (including rts.)    *
  322. *                  *
  323. *    Two of the 16 passes are unrolled from the loop so that *
  324. *    quicker instructions may be used where there is no   *
  325. *    danger of overflow (in the early passes).    *
  326. *                  *
  327. *    Instruction times for branch-type instructions    *
  328. *    listed as (X/Y) are for (taken/not taken).    *
  329. *                  *
  330. *************************************************************************
  331.  
  332. *      Cycles
  333. glsqrt  movem.l  d1-d4,-(sp) (40)
  334.   moveq #13,d4  (4)  ; Loop count (bits-1 of result).
  335.   moveq #0,d1  (4)  ; Error term in D1.
  336.   moveq #0,d2  (4)
  337. lsqrt1  add.l d0,d0  (8)  ; Get 2 leading bits a time and add
  338.   addx.w d1,d1  (4)  ; into Error term for interpolation.
  339.   add.l d0,d0   (8)  ; (Classical method, easy in binary).
  340.   addx.w d1,d1   (4)  
  341.   add.w d2,d2   (4)  ; Running estimate * 2.
  342.   move.w d2,d3  (4)
  343.   add.w d3,d3   (4)
  344.   cmp.w d3,d1  (4)
  345.   bls.s lsqrt2  (10/8)   ; New Error term > 2* Running estimate?
  346.   addq.w #1,d2   (4)  ; Yes, we want a `1' bit then.
  347.   addq.w #1,d3  (4)  ; Fix up new Error term.
  348.   sub.w d3,d1  (4)
  349. lsqrt2  dbra d4,lsqrt1  (10/14) ; Do first 14 bit-pairs.
  350.  
  351.   add.l d0,d0  (8)  ; Do 15-th bit-pair.
  352.   addx.w d1,d1  (4)
  353.   add.l d0,d0  (8)
  354.   addx.l d1,d1  (8)
  355.   add.w d2,d2   (4)
  356.   move.l d2,d3   (4)
  357.   add.w d3,d3   (4)
  358.   cmp.l d3,d1   (6)
  359.   bls.s lsqrt3  (10/8)
  360.   addq.w #1,d2  (4)
  361.   addq.w #1,d3   (4)
  362.   sub.l d3,d1  (8)
  363.  
  364. lsqrt3  add.l d0,d0  (8)  ; Do 16-th bit-pair.
  365.   addx.l d1,d1   (8)
  366.   add.l d0,d0  (8)
  367.   addx.l d1,d1   (8)
  368.   add.w d2,d2  (4)
  369.   move.l d2,d3   (4)
  370.   add.l d3,d3   (8)
  371.   cmp.l d3,d1   (6)
  372.   bls.s lsqrt4  (10/8)
  373.   addq.w #1,d2  (4)
  374. lsqrt4  move.w d2,d0  (4)
  375.   movem.l (sp)+,d1-d4 (44)
  376.   rts    (16)
  377.  
  378.   end
  379.  
  380.  
  381.  
  382. *************************************************************************
  383. Listing Four:
  384.  
  385. *                  *
  386. *  Integer Square Root (32 to 16 bit).        *
  387. *                  *
  388. *  (Newton-Raphson method).          *
  389. *                  *
  390. *  Call with:              *
  391. *    DO.L = Unsigned number.          *
  392. *                  *
  393. *  Returns:              *
  394. *    DO.L = SQRT(DO.L)          *
  395. *                  *
  396. *  Notes:  Result fits in DO.W, but is valid in longword.    *
  397. *    Takes from 338 cycles (1 shift, 1 division) to    *
  398. *    1580 cycles (16 shifts, 4 divisions)  (including rts).  *
  399. *    Averages 854 cycles measured over first 65535 roots.  *
  400. *    Averages 992 cycles measured over first 500000 roots.  *
  401. *                  *
  402. *************************************************************************
  403.  
  404.   .globl lsqrt
  405. *      Cycles
  406. lsqrt  movem.l d1-d2,-(sp) (24)
  407.   move.l d0,d1  (4)  ; Set up for guessing algorithm.
  408.   beq.s return  (10/8)  ; Don't process zero.
  409.   moveq #1,d2  (4)
  410.  
  411. guess  cmp.l d2,d1  (6)  ; Get a guess that is guaranteed to be
  412.   bls.s newton   (10/8)  ; too high, but not by much, by dividing the
  413.   add.l d2,d2   (8)  ; argument by two and multiplying a 1 by 2
  414.   lsr,l #1,d1  (10)  ; until the power of two passes the modified
  415.   bra.s guess   (10)  ; argument, then average these two numbers.
  416.  
  417. newton  add.l d1,d2   (8)  ; Average the two guesses.
  418.   lsr.l #1,d2  (10)
  419.   move.l d0,d1  (4)  ; Generate the next approximation(s)
  420.   divu d2,d1  (140)  ; via the Newton-Raphson method.
  421.   bvs.s done  (10/8)  ; Handle out-of-range input (cheats!)
  422.   cmp.w d1,d2  (4)  ; Have we converged?
  423.   bls.s done  (10/8)  
  424.   swap d1    (4)  ; No, kill the remainder so the
  425.   clr.w d1  (4)  ; next average comes out right.
  426.   swap d1    (4)  
  427.   bra.s newton  (10)
  428.  
  429. done  clr.w d0  (4)  ; Return a word answer in longword.
  430.   swap d0    (4)
  431.   move.w d2,d0  (4)
  432. return  movem.l (sp)+,d1-d2 (28)
  433.   rts    (16)
  434.  
  435.   end
  436. **************************************************************************
  437. Listing Five:
  438.  
  439. Program Squareroot;
  440. {
  441.   Squareroot algorithm & testprogram; DDJ March 1986, p.122
  442.  
  443.   Features:  - sqrt routine in 68000 machine language;
  444.        - long integer loopcount;
  445. }
  446.  
  447. const  { Iteration count for test loop }
  448.   NNR = 6E4;  { real, for printing of statistics }
  449.   NNS = '60000';  { string, for assignment to long integer }
  450.  
  451. type  long = record
  452.       low,high : integer;
  453.     end;
  454.  
  455. var  finished  : boolean;  { flag for loop }
  456.   number, limit   : long;    { loop count, loop limit } 
  457.   sqrrt,        { square root result }
  458.   sqrrto,        { last displayed square root }
  459.   t1,t2    : integer;  { parameters for system time }
  460.   time1, time2  : real;   { start, end time }
  461.   
  462. { --- ROUTINES FOR LONG (32-BIT) INTEGER SUPPORT --- }
  463.  
  464. procedure lg_clr(var l:long);   external;
  465. { Clears long integer l }
  466.  
  467. procedure lg_asn_DU(s:string; var l:long);   external;
  468. { Assigns the unsigned decimal string s to the long integer l }
  469.  
  470. procedure lg_inc1(var l:long);   external;
  471. { Increments long integer l by 1 }
  472.  
  473. procedure lg_grt(a,b:long; var flag:boolean);   external;
  474. { Compares long integers a and b and assign result to flag }
  475.  
  476. { --- SQUAREROOT ROUTINE --- }
  477.  
  478. procedure sqrt(number:long; var result:integer);  external;
  479. { Calculates square root of 'number' and returns it in 'result' }
  480.  
  481. begin { main }
  482.  
  483.   sqrrto := 100;
  484.  
  485.   lg_clr(number);  { number := 0 }
  486.   lg_asn_DU(NNS,limit); {limit := NNS }
  487.   finished := false;
  488.  
  489.   write('Start...');
  490.   time(t1,t2);  time1 := t2;  { get start time }
  491.  
  492.   while not finished do { calculate sqrt(number) }
  493.     begin
  494.       sqrt(number,sqrrt);
  495.     {
  496.       if sqrrt <> sqrrto
  497.    then begin
  498.      write ('Number = ',number.high:6,' | ',number.low:6);
  499.      writeln('  ---  Sqrt = ',sqrrt:4);
  500.      sqrrto := sqrrt;
  501.    end;
  502.     }
  503.       lg_inc1(number);      { number := number + 1 }
  504.       lg_grt(number,limit,finished);  { finished := (number > limit) }
  505.     end;
  506.    
  507.     time(t1,t2);   time2 := t2;  { get end time }
  508.  
  509.     writeln('finished !'); writeln;
  510.     writeln('Time: ',(time2-time1)/60,' seconds');
  511.  
  512. end.
  513. *******************************************************************************
  514. Listing Six:
  515.  
  516. 1 *
  517. 2 * Squareroot algorithm;  DDJ March 1986, p.122
  518. 3 *
  519. 4 * 68000 assembly language version
  520. 5 *
  521. 6 * Features: - equivalent to compiler-generated code;
  522. 7 *
  523. 8 *
  524. 9 * procedure sqrt(number:long; var result:integer);
  525. 10 *
  526. 11 * Calculates integer square root of 'number' and returns it in 'result';
  527. 12 * 
  528. 13 *
  529. 14 * Register usage:
  530. 15 * --------------
  531. 16 *
  532. 17 * D0 : word register    A0 : parameter stack pointer
  533. 18 * D1 : number              A1 : scratch register
  534. 19 * D2 : guess1
  535. 20 * D3 : guess2
  536. 21 * D4 : error
  537. 22 *
  538. 23   proc  sqrt,2    ;2 words of parameters
  539. 24 *
  540. 25 * Get parameters from stack
  541. 26 *
  542. 27  move.l   (sp)+,a0  ;get return address  
  543. 28  move.w  2(sp),a1  ;get ^number
  544. 29  move.l  (a1),d1    ;get number
  545. 30 *
  546. 31 * bra Q15 ;--- for timing only
  547. 32 *
  548. 33  beq  Q8    ;if number=0, skip
  549. 34 *
  550. 35 * Set initial values
  551. 36 *
  552. 37 Q7  moveq  #1,d2    ;guess1 := 1
  553. 38  move.1  d1,d3    ;guess2 := number
  554. 39 *
  555. 40 * Do shifts until guess1 ~~~  guess2
  556. 41 *
  557. 42 Q9  move.l  d2,d0    ;guess1 to work register
  558. 43  asl.l  #1,d0    ;guess1 * 2
  559. 44  cmp.l  d3,d0    ;compare with guess2
  560. 45  bge  Q11    ;branch if guess1 * 2 >= guess2
  561. 46  asl.l   #1,d2    ;guess1 := guess1 * 2
  562. 47  asr.l  #1,d3    ;guess2 := guess2 / 2
  563. 48  bra  Q9
  564. 49 *
  565. 50 * Now do divisions
  566. 51 *
  567. 52 Q11
  568. 53 Q13  add.l  d3,d2    ;guess := guess1 + guess2
  569. 54  asr.l  #1,d2    ;guess1 := (guess1+guess2)/2
  570. 55  move.l  d1,d0    ;number to work register
  571. 56  divs  d2,d0    ;number / guess1
  572. 57  move.w  d0,d3    ;guess2 := number/guess1
  573. 58  ext.l  d3    ;extend to 32 bits
  574. 59  move.l  d2,d0    ;guess1 to work register
  575. 60  sub.l  d3,d0    ;guess1-guess2
  576. 61  move.l  d0,d4    ;error := guess1-guess2
  577. 62  ble  Q14    ;if error <= 0
  578. 63  bra  Q13    ;loop back if error > 0
  579. 64 Q14  move.l  d2,d0    ;result := guess1
  580. 65  bra  Q15
  581. 66 Q8  moveq   #0,d0    ;result := 0
  582. 67 *
  583. 68 * Set result & return to caller
  584. 69 *
  585. 70 Q15  movea.w (sp)+,a1  ;get ^result
  586. 71  move.w  d0,(a1)    ;store result
  587. 72 *
  588. 73  addq.l  #2,sp    ;drop ^number
  589. 74  jmp  (a0)    ;return to caller
  590. 75 *
  591. 76  .nolist
  592.  
  593.  
  594. ******************************************************************************
  595. Listing Seven:
  596.  
  597. 1 *
  598. 2 * Squareroot algorithm; DDJ March 1986, p.122
  599. 3 *
  600. 4 * 68000 assembly language version
  601. 5 *
  602. 6 * Features:  - hand-optimized machine code;
  603. 7 *
  604. 8 * 
  605. 9 * procedure sqrt(number:integer; var result:integer);
  606. 10 *
  607. 11 * Calculates square root of 'number' and returns it in 'result';
  608. 12 *
  609. 13 *
  610. 14 * Register usage:
  611. 15 * --------------
  612. 16 *
  613. 17 * D0 : work register, error    A0 : ^result
  614. 18 * D1 : number      A1 : ^number
  615. 19 * D2 : guess1,result
  616. 20 * D3 : guess2
  617. 21 * D4 : temporary storage for return address
  618. 22 *
  619. 23 *
  620. 24   .proc  sqrt,2    ;2 words of parameters
  621. 25 *
  622. 26 * Get parameters from stack
  623. 27 * 
  624. 28  moveq #0,d2    ;result := 0 --- remove for timing
  625. 29 *
  626. 30   move.l   (sp)+,d4  ;get return address
  627. 31   move.w  (sp)+,a0  ;get ^result
  628. 32  move.w (sp)+,a1    ;get ^number
  629. 33  move.l  (a1),d1    ;get number
  630. 34 *
  631. 35 * bra Q15 ;--- for timing only
  632. 36 *
  633. 37  beq   Q15    ;if number=0, then result=0
  634. 38 *
  635. 39 * Set initial values
  636. 40 *
  637. 41 Q7  moveq  #1,d2    ;guess1 := 1
  638. 42  move.l  d1,d3    ;guess2 := number
  639. 43 *
  640. 44 * Do shifts until guess1 ~  guess2
  641. 45 *
  642. 46 Q9  asl.l  #1,d2    ;guess1 := guess1 * 2
  643. 47  cmp.l  d3,d2    ;compare with guess2
  644. 48  bge  Q11    ;branch if guess1 * 2 >= guess2
  645. 49  asr.l  #1,d3    ;guess2 := guess2 / 2
  646. 50  bra  Q9
  647. 51 *
  648. 52 Q11  asr.l  #1,d2    ;adjust guess1
  649. 53 *
  650. 54 * Now do divisions
  651. 55 *
  652. 56 Q13  add.l  d3,d2    ;guess1 := guess1 + guess2
  653. 57   asr.l  #1,d2    ;guess1 := (guess1+guess2)/2
  654. 58  move.l  d1,d3    ;guess2 := number
  655. 59  divs  d2,d3    ;guess2 := number / guess1
  656. 60  ext.l  d3    ;extend guess2 to 32 bits
  657. 61  move.l  d2,d0    ;guess1 to work register
  658. 62  sub.l  d3,d0    ;error := guess1 - guess2
  659. 63  bgt  Q13    ;loop back if error > 0
  660. 64 *
  661. 65 * Store result and return to caller
  662. 66 *
  663. 67 Q15  move.w   d2,(a0)    ;store result   
  664. 68 *
  665. 69   movea.l  d4,a0    ;move return address to adr-reg
  666. 70  jmp  (a0)    ;return to caller
  667. 71 *
  668. 72  .nolist
  669.  
  670.  
  671. ***************************************************************************
  672.  
  673. Listing Eight:
  674.  
  675. 1 *
  676. 2 * Squareroot algorithm; DDJ November 1985, p.88
  677. 3 *
  678. 4 * 68000 assembly language
  679. 5 *
  680. 6 * procedure sqrt(number:long; var result:integer);
  681. 7 *
  682. 8 * Calculates the integer squareroot of 'number' and returns it in 'result'
  683. 9 *
  684. 10 *
  685. 11 * Register usage
  686. 12 * --------------
  687. 13 *
  688. 14 * D0 : number    A0 : scratch, for pointers
  689. 15 * D1 : error term    A1 : scratch, for pointers
  690. 16 * D2 : estimate
  691. 17 * D3 : corrector term
  692. 18 * D4 : loop counter
  693. 19 *
  694. 20 *
  695. 21  .proc  sqrt,2    ;2 words of parameters
  696. 22 *
  697. 23   move.w   6(sp),a0  ;get ^number
  698. 24  move.l  (a0),d0    ;get number into d0
  699. 25 *
  700. 26 * bra exit ;--- for timing only
  701. 27 *
  702. 28  moveq  #16-1,d4  ;set loopcount,16*2 = 32 bits
  703. 29  moveq  #0,d1    ;clear error term
  704. 30   moveq  #0,d2    ;clear estimate
  705. 31 *
  706. 32 * Calculate squareroot
  707. 33 *
  708. 34 sqrtl asl.l #1,d0    ;get 2 leading bits,
  709. 35  roxl.l  #1,d1    ;one at a time, into
  710. 36  asl.l  #1,d0    ;the error term
  711. 37  roxl.l  #1,d1
  712. 38  asl.l   #1,d2    ;estimate * 2
  713. 39  move.l  d2,d3    
  714. 40  asl.l  #1,d3    ;corrector = 2 * new estimate
  715. 41  cmp.l  d3,d1
  716. 42  bls  sqrt2    ;branch if error term <= corrector
  717. 43  addq.l  #1,d2    ;otherwise, add low bit to estimate
  718. 44  addq.l  #1,d3
  719. 45  sub.l  d3,d1    ;and calculate new error term
  720. 46 *
  721. 47 sqrt2 dbra  d4,sqrtl  ;do all 16 bit-pairs
  722. 48 *
  723. 49 * Set result & return to caller
  724. 50 *
  725. 51 exit  move.l   (sp)+,a0  ;get return address
  726. 52  move.w  (sp)+,a1  ;get ^result
  727. 53  move.w  d2,(a1)    ;store integer result
  728. 54  addq.l  #2,sp    ;drop ^number
  729. 55  jmp  (a0)    ;return to caller
  730. 56 *
  731. 57   .nolist
  732.  
  733.