home *** CD-ROM | disk | FTP | other *** search
/ Amiga ISO Collection / AmigaUtilCD2.iso / Misc / The_SOURCE#7.LHA / Source / NumericalMethods.lha / NumericalMethods / sqrt.s < prev   
Encoding:
Text File  |  1994-05-27  |  4.6 KB  |  137 lines
  1. Newsgroups: comp.graphics.algorithms
  2. Path: usenet.ee.pdx.edu!cs.uoregon.edu!sgiblab!swrinde!cs.utexas.edu!howland.reston.ans.net!pipex!sunic!umdac!cs.umu.se!christer
  3. From: christer@cs.umu.se (Christer Ericson)
  4. Subject: Re: Fastest long integer square root algorithm?
  5. Message-ID: <Cnn27G.DEv@cs.umu.se>
  6. Sender: news@cs.umu.se (News Administrator)
  7. Organization: Dept. of Computing Science, Umea Univ., 901 87 Umea, Sweden
  8. References: <2neb6t$mcq@geraldo.cc.utexas.edu>
  9. Date: Sat, 2 Apr 1994 15:40:25 GMT
  10. Lines: 124
  11.  
  12. In <2neb6t$mcq@geraldo.cc.utexas.edu> Steve Mariotti <stevem@zaphod.ee.pitt.edu> writes:
  13. >
  14. >My goal in life right now is to find the fastest integer square root
  15. >algorithm that is humanly possible.  I know, I should probably set my
  16. >sights a little higher.
  17. >[...]
  18.  
  19. This piece of code is written to work on all 680x0-processors (which
  20. seemed to be your target architechture). You can't go faster than this
  21. on a 68000 without resorting to some sort of table lookup.
  22.  
  23. It is the basic Newton/Babylonian method, but with some interesting
  24. twists.
  25.  
  26.  
  27. --- cut here ---
  28.  
  29. unsigned short ce_quick_sqrt(n)
  30. register unsigned long n;
  31. {
  32.     asm {
  33. ;-------------------------
  34. ;
  35. ; Routine       : ISQRT; integer square root
  36. ; I/O parameters: d0.w = sqrt(d0.l)
  37. ; Registers used: d0-d2/a0
  38. ;
  39. ; This is a highly optimized integer square root routine, based
  40. ; on the iterative Newton/Babylonian method
  41. ;
  42. ;    r(n + 1) = (r(n) + A / R(n)) / 2
  43. ;
  44. ; Verified over the full interval [0x0L,0xFFFFFFFFL]
  45. ;
  46. ; Some minor compromises have been made to make it perform well
  47. ; on the 68000 as well as 68020/030/040. This routine outperforms
  48. ; the best of all other algorithms I've seen (except for a table
  49. ; lookup).
  50. ;
  51. ; Copyright (c) Christer Ericson, 1993. All rights reserved.
  52. ;
  53. ; Christer Ericson, Dept. of Computer Science, University of Umea,
  54. ; S-90187 UMEA, SWEDEN. Internet: christer@cs.umu.se
  55. ;
  56.  
  57. ;-------------------------
  58. ;
  59. ; Macintosh preamble since THINK C passes first register param
  60. ; in d7. Remove this section on any other machine
  61. ;
  62.     move.l    d7,d0
  63. ;-------------------------
  64. ;
  65. ; Actual integer square root routine starts here
  66. ;
  67.     move.l    d0,a0        ; save original input value for the iteration
  68.     beq.s    @exit        ; unfortunately we must special case zero
  69.     moveq    #2,d1        ; init square root guess
  70. ;-------------------------
  71. ;
  72. ; Speed up the process of finding a power of two approximation
  73. ; to the actual square root by shifting an appropriate amount
  74. ; when the input value is large enough
  75. ;
  76. ; If input values are word only, this section can be removed
  77. ;
  78.     move.l    d0,d2
  79.     swap    d2        ; do [and.l 0xFFFF0000,d2] this way to...
  80.     tst.w    d2        ; go faster on 68000 and to avoid having to...
  81.     beq.s    @skip8        ; reload d2 for the next test below
  82.     move.w    #0x200,d1    ; faster than lsl.w #8,d1 (68000)
  83.     lsr.l    #8,d0
  84. @skip8    and.w    #0xFE00,d2    ; this value and shift by 5 are magic
  85.     beq.s    @skip4
  86.     lsl.w    #5,d1
  87.     lsr.l    #5,d0
  88. @skip4
  89.  
  90. ;-------------------------
  91. ;
  92. ; This finds the power of two approximation to the actual square root
  93. ; by doing the integer equivalent to sqrt(x) = 2 ^ (log2(x) / 2). This
  94. ; is done by shifting the input value down while shifting the root
  95. ; approximation value up until they meet in the middle. Better precision
  96. ; (in the step described below) is gained by starting the approximation
  97. ; at 2 instead of 1 (this means that the approximation will be a power
  98. ; of two too large so remember to shift it down).
  99. ;
  100. ; Finally (and here's the clever thing) since, by previously shifting the
  101. ; input value down, it has actually been divided by the root approximation
  102. ; value already so the first iteration is "for free". Not bad, eh?
  103. ;
  104. @loop    add.l    d1,d1
  105.     lsr.l    #1,d0
  106.     cmp.l    d0,d1
  107.     bcs.s    @loop
  108. @skip    lsr.l    #1,d1        ; adjust the approximation
  109.     add.l    d0,d1        ; here we just add and shift to...
  110.     lsr.l    #1,d1        ; get the first iteration "for free"!
  111. ;-------------------------
  112. ;
  113. ; The iterative root value is guaranteed to be larger than (or equal to)
  114. ; the actual square root, except for the initial guess. But since the first
  115. ; iteration was done above, the loop test can be simplified below
  116. ; (Oh, without the bvs.s the routine will fail on most large numbers, like
  117. ; for instance, 0xFFFF0000. Could there be a better way of handling these,
  118. ; so the bvs.s can be removed? Nah...)
  119. ;
  120. @loop2    move.l    a0,d2        ; get original input value
  121.     move.w    d1,d0        ; save current guess
  122.     divu.w    d1,d2        ; do the Newton method thing
  123.     bvs.s    @exit        ; if div overflows, exit with current guess
  124.     add.w    d2,d1
  125.     roxr.w    #1,d1        ; roxr ensures shifting back carry overflow
  126.     cmp.w    d0,d1
  127.     bcs.s    @loop2        ; exit with result in d0.w
  128. @exit
  129.     }
  130. }
  131.  
  132. --- and here ---
  133.  
  134. Christer Ericson --- Internet: christer@cs.umu.se --- tel: +46-90-166794
  135. Department of Computer Science, University of Umea, S-90187 UMEA, SWEDEN
  136.  
  137.