home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 3 / TheARMClub_PDCD3.iso / hensa / documentation / documents / a252sqrt < prev    next >
Internet Message Format  |  1999-04-27  |  20KB
  1. From: chughes@maths.tcd.ie (Conrad Hughes)
  2. Subject: Re: 32/64 bit square root.
  3. Date: 25 Sep 91 16:30:11 GMT
  4.  
  5. Here's the best code I've managed yet; it manages to fit into four
  6. instructions per bit. If anyone has managed better, please let me
  7. know.. If you do use the code, it'd be nice to receive a mention
  8. :-)
  9. DEFFNsqrt(N,R,T,S):LOCALA%:IFS<>N [OPTpass:MOV S,N:]
  10. [OPTpass:MOV R,#0:]:FORA%=7TO0STEP-1:[OPTpass:ADD T,R,#1<<A%+15
  11. CMP S,T,LSL#A%+1:SUBGE S,S,T,LSL#A%+1:ADDGE R,R,#1<<A%+16:]NEXT:[OPTpass
  12. ADD T,R,#1<<14:CMP S,T:SUBGE S,S,T:ADDGE R,R,#1<<15:]FORA%=2TO15:[OPTpass
  13. ADD T,R,#1<<15-A%:RSBS S,T,S,LSL#1:ADDLT S,S,T:ADDGE R,R,#1<<16-A%:]NEXT
  14. [OPTpass:ADD T,R,R:ADD T,T,#1:RSBS S,T,S,LSL#2:ADDLT S,S,T:ADDGE R,R,#1
  15. CMP S,R:ADDGE R,R,#1:]=0
  16. ..defines a function which calculates the square root of a number
  17. stored in one register, with the fixed point between bits 15 and 16;
  18. I didn't write down behaviour for negative numbers, but you should
  19. be able to test this out for yourself, and it shouldn't be too
  20. difficult to extend to different formats. To convert a number into
  21. this format, you just multiply it by 65536, and to convert back, you
  22. divide by 65536 (i.e., 1.000 would be represented as 65536 in a
  23. register).
  24.  
  25. So, for example, you want R1=SQRT(R0):
  26. FNsqrt(0,1,2,3) will preserve the contents of R0, and use R2 and R3
  27. for data - both will be corrupted.
  28. FNsqrt(0,1,2,0) will corrupt R0, but doesn't need R3.
  29. There are no other circumstances under which two of the registers
  30. the routine uses can safely be the same.
  31.  
  32. The loops are unwrapped for speed (actually using loops would take
  33. up much less space, but would probably use up another register, and
  34. make the routine twice as slow). The routine takes up 400 or 404
  35. bytes, depending on whether N=S or not.
  36.  
  37. One thing - this is slow :-}, so wherever possible you should avoid
  38. using square roots.. Similarly division.
  39.  
  40. Conrad
  41.  
  42.  
  43.  
  44. From: gcwillia@daisy.waterloo.edu (Graeme Williams)
  45. Subject: Re: The (non-charity) ray tracer
  46. Date: 20 Nov 91 19:21:07 GMT
  47.  
  48. The algorithm I'd use would be a Newton-Raphson iteration scheme designed
  49. to solve x^2=z. (Whether this is the fastest method I don't know)
  50.  
  51. The way this works is you make a guess for sqrt(z), say x_0. (The guess
  52. can be very rough, I'd do something like find out how many bits in z
  53. (say there were 23) half that (11 or 12) and then guess x_0 (2^11 or 2^12.)
  54.  
  55. One then makes successive approximations via 
  56.  
  57.  x_{n+1} = x_{n} - (x_{n}^2 - z)/2/x_{n}
  58.  
  59. This scheme will converge very rapidly to sqrt(z) if your initial guess
  60. is at all close to z (give or take a factor of 2 or 4 is close enough).
  61. You might want to write (x_{n}^2-z) as (x_{n}+z)*(x_{n}-z) it might be
  62. quicker (I've no idea).
  63.  
  64. As an example: if z=529020918
  65.  
  66. and we guess     x0=16384
  67. then we find     x1=24336.43719
  68.                  x2=23037.12504
  69.                  x3=23000.48392
  70.                  x4=23000.45473
  71.                  x5=23000.45473
  72.                  x6=23000.45473
  73.  
  74. so you can see 4 or maybe 5 iterations at most is enough. Each iteration
  75. (if you arrange things *carefully*) will need only a 16x16 bit multiply
  76. and a 32x16 bit divide.
  77.  
  78. Graeme Williams
  79.  
  80.  
  81.  
  82. From: torbenm@diku.dk (Torben AEgidius Mogensen)
  83. Date: 21 Nov 91 15:59:37 GMT
  84. Sender: torbenm@freke.diku.dk
  85.  
  86. gcwillia@daisy.waterloo.edu (Graeme Williams) writes:
  87. >The algorithm I'd use would be a Newton-Raphson iteration scheme designed
  88. >to solve x^2=z. (Whether this is the fastest method I don't know)
  89. >...
  90. >so you can see 4 or maybe 5 iterations at most is enough. Each iteration
  91. >(if you arrange things *carefully*) will need only a 16x16 bit multiply
  92. >and a 32x16 bit divide.
  93.  
  94. There is a *much* faster way. If your number is of the form a*2^n,
  95. where 0.5 <= a < 1.0 then do the following:
  96.  
  97.     sqrt(a*2^(2n)) = sqrt(a)*2^n
  98.     sqrt(a*2^(2n+1)) = sqrt(a/2)*2^(n+1)
  99.  
  100. This reduces the problem to finding a square root of a number a,
  101. 0.25 <= a < 1.0. This can be done by this method, which should have
  102. a cost comparable to 2 or 3 divisions.
  103.  
  104.     { 0.25 <= a < 1.0 }
  105.     x := 4*(a-0.25);
  106.     r := 1;
  107.     n := 1;
  108.     while (n<number-of-bits) and (a<>0) do begin
  109.       { invariant: r^2+x = 4^n*a, x<2*r+1 }
  110.       n := n+1;
  111.       t := 4*r+1;
  112.       x := 4*x;
  113.       r := 2*r;
  114.       if x >= t then begin
  115.         x := x-t;
  116.         r := r+1
  117.       end
  118.     end;
  119.     r := r/2^n;
  120.     { sqrt(a)-r < 2^-n }
  121.  
  122. Note that only addition, subtractions and multiplications by 2 or 4
  123. are used in each part of the loop. At the end a single n-bit shift is
  124. done. Note that r and t are integers, while x is a real.
  125.  
  126. The (a<>0) test allows you to exit the loop if an exact root is found,
  127. but it can be omitted without causing error (so the loop can be
  128. unwound).
  129.  
  130. We can reformulate the algorithm to keep x < 1.0 throughout the
  131. algorithm:
  132.  
  133.     { 0.25 <= a < 1.0 }
  134.     x := (a-0.25)/2;
  135.     r := 1;
  136.     for n := 2 to number-of-bits do begin
  137.       { invariant: r^2+2^n*x= 4^n*a, x<(2*r+1)*2^-n }
  138.       t := 4*r+1;
  139.       x := 2*x;
  140.       r := 2*r;
  141.       if x >= t/2^(n+1) then begin
  142.         x := x-t/2^(n+1);
  143.         r := r+1
  144.       end
  145.     end;
  146.     r := r/2^number-of-bits;
  147.     { sqrt(a)-r < 2^-number-of-bits }
  148.  
  149. Assuming we have a 32 bit mantissa which is stored with the most
  150. significant bit as 0.5, we would get the following ARM code, which
  151. uses 213 cycles to get a 32 bit result. A test of x=0 could be added
  152. at every iteration. This would make minimum time shorter, but assuming
  153. there is no large number of squares in the input, this would make the
  154. average time longer. 
  155.  
  156.     \ on entry Ra holds a
  157.     \ on exit Rr holds sqrt(a)
  158.     \ uses Rx and Rt (either can be equal to Ra)
  159.     SUB Rx,Ra,#2^30
  160.     MOV Rx,Rx LSR #1    \ x := (a-0.25)/2;
  161.     MOV Rr,#1        \ r := 1;
  162.  
  163.     MOV Rt,Rr ASL #2
  164.     ADD Rt,Rt,#1        \ t := 4*r+1;
  165.     MOV Rx,Rx ASL #1    \ x := 2*x;
  166.     MOV Rr,Rr, ASL #1    \ r := 2*r;
  167.     SUBS Rt,Rx,Rt ASL #29    \ if x >= t/2^(n+1) then begin
  168.     MOVPL Rx,Rt        \   x := x-t/2^(n+1);
  169.     ADDPL Rr,Rr,#1        \   r := r+1
  170.  
  171.     ...
  172.  
  173.     MOV Rt,Rr ASL #2
  174.     ADD Rt,Rt,#1        \ t := 4*r+1;
  175.     MOV Rx,Rx ASL #1    \ x := 2*x;
  176.     MOV Rr,Rr, ASL #1    \ r := 2*r;
  177.     SUBS Rt,Rx,Rt ASL #0    \ if x >= t/2^(n+1) then begin
  178.     MOVPL Rx,Rt        \   x := x-t/2^(n+1);
  179.     ADDPL Rr,Rr,#1        \   r := r+1
  180.  
  181.     \ at exit: 0.5 <= Rr < 1.0 when considered as a mantissa
  182.  
  183. We can actually get this better by avoiding the doubling of r in
  184. each iteration. This way r and t are no longer integers, but the code
  185. is better anyway: 
  186.  
  187.     { 0.25 <= a < 1.0 }
  188.     x := a-0.25;
  189.     r := 0.5;
  190.     for n := 2 to number-of-bits do begin
  191.       { invariant: r^2+x/2^(n-1) = a, x<r+2^-(n+1) }
  192.       t := r+2^-(n+2);
  193.       x := 2*x
  194.       if x >= t then begin
  195.         x := x-t;
  196.         r := r+2^-(n+1)
  197.       end
  198.     end;
  199.     { sqrt(a)-r < 2^-number-of-bits }
  200.  
  201. In ARM this uses 149 cycles to get a 32 bit result. Try to better that!
  202.  
  203.     \ on entry Ra holds a
  204.     \ on exit Rr holds sqrt(a)
  205.     \ uses Rx and Rt (either can be equal to Ra)
  206.     SUB Rx,Ra,#2^30        \ x := a-0.25;
  207.     MOV Rr,#2^31        \ r := 0.5;
  208.  
  209.     ADD Rt,Rr,#2^28        \ t := r+2^-(n+2);
  210.     MOV Rx,Rx ASL #1    \ x := 2*x;
  211.     SUBS Rt,Rx,Rt        \ if x >= t then begin
  212.     MOVPL Rx,Rt        \   x := x-t;
  213.     ADDPL Rr,Rr,#2^29    \   r := r+2^-(n+1)
  214.  
  215.     ...
  216.  
  217.     ADD Rt,Rr,#2^0        \ t := r+2^-(n+2);
  218.     MOV Rx,Rx ASL #1    \ x := 2*x;
  219.     SUBS Rt,Rx,Rt        \ if x >= t then begin
  220.     MOVPL Rx,Rt        \   x := x-t;
  221.     ADDPL Rr,Rr,#2^1    \   r := r+2^-(n+1)
  222.  
  223.                 \ t := r+2^-(n+2); 
  224.                 \ { 2^-(n+2) is less than precision }
  225.                 \ x := 2*x;
  226.     RSBS Rt,Rr,Rx ASL #1    \ if x >= t then begin
  227.                 \   x := x-t; { x not needed }
  228.     ADDPL Rr,Rr,#2^0    \   r := r+2^-(n+1)
  229.  
  230.     \ at exit: 0.5 <= Rr < 1.0 when considered as a mantissa
  231.  
  232.  
  233. Note: I haven't tried this code on my Archimedes, as it is at home and
  234. I am at work. I have, however, tried the pseudo-Pascal code.
  235.  
  236.  
  237. From: torbenm@diku.dk (Torben AEgidius Mogensen)
  238. Date: 22 Nov 91 10:46:26 GMT
  239.  
  240. I made a small error in this. The correct code is: 
  241.  
  242.     { 0.25 <= a < 1.0 }
  243.     x := a-0.25;
  244.     r := 0.5;
  245.     for n := 2 to number-of-bits do begin
  246.       { invariant: r^2+x/2^(n-2) = a, x<r+2^-(n+1) }
  247.       t := r+2^-(n+1);
  248.       x := 2*x
  249.       if x >= t then begin
  250.         x := x-t;
  251.         r := r+2^-n
  252.       end
  253.     end;
  254.     { sqrt(a)-r < 2^-number-of-bits }
  255.  
  256. There are, however, still a small problem: x might be greater than 1
  257. after the doubling. We can solve this by delaying the doubling:
  258.  
  259.     { 0.25 <= a < 1.0 }
  260.     x := a-0.25;
  261.     r := 0.5;
  262.     for n := 2 to number-of-bits do begin
  263.       { invariant: r^2+x/2^(n-2) = a, x<r+2^-(n+1) }
  264.       t := r+2^-(n+1)
  265.       if x >= t/2 then begin
  266.         x := x-t/2;
  267.         r := r+2^-n
  268.       end;
  269.       x := 2*x
  270.     end;
  271.     { sqrt(a)-r < 2^-number-of-bits }
  272.  
  273. Also, my ARM code used unsigned numbers but compared them as if they
  274. were signed (using PL). Comparison of unsigned numbers requires HS.
  275. The code is still untested, so there might be further bugs. Note also
  276. that you can't write the constant 2^31 in BASIC assembler. Use 1<<31
  277. instead. The number of cycles have increased to 154, as an extra pass
  278. through the loop is required to get 32 bit precision:
  279.  
  280.     \ on entry Ra holds a
  281.     \ on exit Rr holds sqrt(a)
  282.     \ uses Rx and Rt (either can be equal to Ra)
  283.     SUB Rx,Ra,#2^30        \ x := a-0.25;
  284.     MOV Rr,#2^31        \ r := 0.5;
  285.  
  286.     ADD Rt,Rr,#2^29        \ t := r+2^-(n+1);
  287.     SUBS Rt,Rx,Rt,LSR #1    \ if x >= t/2 then begin
  288.     MOVHS Rx,Rt        \   x := x-t/2;
  289.     ADDHS Rr,Rr,#2^30    \   r := r+2^-n
  290.     MOV Rx,Rx,ASL #1    \ x := 2*x;
  291.  
  292.     ...
  293.  
  294.     ADD Rt,Rr,#2^0        \ t := r+2^-(n+1);
  295.     SUBS Rt,Rx,Rt,LSR #1    \ if x >= t/2 then begin
  296.     MOVHS Rx,Rt        \   x := x-t/2;
  297.     ADDHS Rr,Rr,#2^1    \   r := r+2^-n
  298.     MOV Rx,Rx,ASL #1    \ x := 2*x;
  299.  
  300.                 \ t := r+2^-(n+1); 
  301.                 \ { 2^-(n+1) is less than precision }
  302.     SUBS Rt,Rx,Rr,LSR #1    \ if x >= t/2 then begin
  303.                 \   x := x-t/2; { x not needed }
  304.     ADDHS Rr,Rr,#2^0    \   r := r+2^-n
  305.                 \ x := 2*x; { x not needed }
  306.     \ at exit: 0.5 <= Rr < 1.0 when considered as a mantissa
  307.  
  308.  
  309. Torben Mogensen (torbenm@diku.dk)
  310.  
  311.  
  312.  
  313. From: dseal@armltd.co.uk (David Seal)
  314. Subject: Square roots (was Re: The (non-charity) ray tracer)
  315. Date: 21 Nov 91 15:52:14 GMT
  316.  
  317. The Newton-Raphson technique is a good one when you've got a fast division -
  318. e.g. if you've got division in hardware (and haven't got square root in
  319. hardware :-). However, if you're going to do division by a software routine
  320. (as you must on an ARM), there is a "long square root" algorithm which is
  321. very similar to the "long division" algorithms used to do division in
  322. software. It takes only a bit more time than one division *in total*,
  323. compared with a division, an addition and a shift *for each of several
  324. iterations* in Newton-Raphson. I would therefore recommend it for ARM square
  325. root routines when speed is important.
  326.  
  327. I won't supply code, since (a) I don't know what type of numbers you want to
  328. take square roots of, how accurate a result you want, etc. - these details
  329. will obviously affect the code; (b) I don't have time to write code, whereas
  330. I have a written description of the algorithm on file which I can simply
  331. include here!
  332.  
  333. I'll illustrate the "long square root" algorithm in base 10 - suppose we
  334. want to take the square roots of 2 and 20:
  335.  
  336. (A) First write down the number you want to take the square root of, split
  337.     into groups of two digits, with one split at the decimal point. In these
  338.     particular cases, this gives us:
  339.  
  340.        -----------------                     -----------------
  341.       V 02. 00 00 00 ...                    V 20. 00 00 00 ...
  342.  
  343.     (Note that the two numbers chosen have the same pattern of digits, but
  344.     different splits because of the different positions of the decimal
  345.     point. This is the only reason they produce different answers.)
  346.  
  347. (B) Initialise by finding the largest square which is smaller than the first
  348.     pair of digits: in the first case the first pair of digits is 02 = 2, so
  349.     this square is 1^2 = 1, and in the second the first pair of digits is
  350.     20, so this square is 4^2 = 16. Subtract the square from the first pair
  351.     of digits, and put its square root in as the first digit of the result:
  352.  
  353.          1.                                    4.
  354.        -----------------                     -----------------
  355.       V 02. 00 00 00 ...                    V 20. 00 00 00 ...
  356.        - 1                                   -16
  357.         --                                    --
  358.          1                                     4
  359.  
  360.     Then bring the next pair of digits down:
  361.  
  362.          1.                                    4.
  363.        -----------------                     -----------------
  364.       V 02. 00 00 00 ...                    V 20. 00 00 00 ...
  365.        - 1                                   -16
  366.         ------                                ------
  367.            100                                   400
  368.  
  369. (C) To find the next digit of the result: take twice the answer so far. Look
  370.     at numbers of the form N * (twice answer so far concatenated with N) for
  371.     N a digit. The largest one of these that is still less than or equal to
  372.     the remainder we've got gives the next digit of the result, and the
  373.     value of N * (twice answer so far concatenated with N) is what we should
  374.     subtract from the remainder before bringing down the next pair of
  375.     digits:
  376.  
  377.       Twice answer so far is 2*1=2.         Twice answer so far is 2*4=8.
  378.       0*20 = 0, less than 100               0*80 = 0, less than 400
  379.       1*21 = 21, less than 100              1*81 = 81, less than 400
  380.       2*22 = 44, less than 100              2*82 = 164, less than 400
  381.       3*23 = 69, less than 100              3*83 = 249, less than 400
  382.       4*24 = 96, less than 100              4*84 = 336, less than 400
  383.       5*25 = 125, more than 100             5*85 = 425, more than 400
  384.       So the next digit of the answer       So the next digit of the answer
  385.         is 4, and we must subtract 96         is 4, and we must subtract 336
  386.  
  387.          1.  4                                 4.  4
  388.        -----------------                     -----------------
  389.       V 02. 00 00 00 ...                    V 20. 00 00 00 ...
  390.        - 1                                   -16
  391.         ------                                ------
  392.            100                                   400
  393.           - 96                                  -336
  394.            ------                                ------
  395.               400                                  6400
  396.  
  397. (D) Repeat step (C) as many times as you need result digits - I'll do it
  398.     once more for the examples:
  399.  
  400.       Twice answer so far is 2*14=28.       Twice answer so far is 2*44=88.
  401.       0*280 = 0, less than 400              0*880 = 0, less than 6400
  402.       1*281 = 281, less than 400            1*881 = 881, less than 6400
  403.       2*282 = 564, more than 400            2*882 = 1764, less than 6400
  404.       So the next digit of the answer       :
  405.         is 1, and we must subtract 281      :
  406.                                             7*887 = 6209, less than 6400
  407.                                             8*888 = 7104, more than 6400
  408.                                             So the next digit of the answer
  409.                                               is 7, and we must subtract 6209
  410.  
  411.          1.  4  1                              4.  4  7
  412.        -----------------                     -----------------
  413.       V 02. 00 00 00 ...                    V 20. 00 00 00 ...
  414.        - 1                                   -16
  415.         ------                                ------
  416.            100                                   400
  417.           - 96                                  -336
  418.            ------                                ------
  419.               400                                  6400
  420.             - 281                                 -6209
  421.              -------                               -------
  422.                11900                                 19100
  423.  
  424. As to why the method works, I'll give some hints and leave you to work out
  425. the details yourself. To give the rough reason for it, note that the
  426. "remainder" in the above calculations is actually the difference between the
  427. number you are trying to take the square root of and the square of the
  428. answer so far:
  429.  
  430.       2                                     20
  431.       = 1^2 + 1                             = 4^2 + 4
  432.       = 1.4^2 + 0.04                        = 4.4^2 + 0.64
  433.       = 1.41^2 + 0.0119                     = 4.47^2 + 0.0191
  434.       = ...                                 = ...
  435.  
  436. [ Aside: this fact yields another useful property of the algorithm, which is
  437.   that (as in long division) your answer is exactly right if the remainder
  438.   becomes 0 and there are no further non-zeros to bring down. This can be
  439.   very helpful e.g. for calculating square roots in floating point systems
  440.   that conform to IEEE standard 754. ]
  441.  
  442. Then prove this - it depends on the identity:
  443.  
  444.      (10x+y)^2 - (10x)^2  =  20xy + y^2
  445.                           =  (20x + y) * y
  446.  
  447. The quantities (20x + y) * y are what the digit calculation in part (C) is
  448. actually determining.
  449.  
  450. One important observation: this algorithm simplifies immensely in binary!
  451. I'll run through it, illustrating it for the case of the square root of two
  452. again:
  453.  
  454. (A) Write down the argument, split into pairs aligned with the binary point:
  455.  
  456.        --------------------
  457.       V 10. 00 00 00 00 ...
  458.  
  459. (B) If the first pair is 00, the first digit of the answer is 0; otherwise
  460.     it is 1. Subtract the square of this first digit from the first pair and
  461.     bring down the next pair:
  462.  
  463.          1.
  464.        --------------------
  465.       V 10. 00 00 00 00 ...
  466.        - 1
  467.         ------
  468.            100
  469.  
  470. (C) There are only two possible next digits, namely 0 and 1. But 0 * (twice
  471.     the answer so far concatenated with 0) is always 0, which must be less
  472.     than or equal to the remainder. So all that matters is whether 1 *
  473.     (twice the answer so far concatenated with 1) is less than or equal to
  474.     the remainder: if it is, the next result digit is 1; if it is greater
  475.     than the remainder, the next result digit is 0. One more thing: twice
  476.     the answer so far concatenated with 1 is the answer so far concatenated
  477.     with 01.
  478.  
  479.     So the result digit selection procedure becomes: concatenate the answer
  480.     so far with 01, and compare with the remainder. If greater, leave the
  481.     remainder unchanged and the result digit is 0; if less than or equal,
  482.     subtract it from the remainder and the result digit is 1.
  483.  
  484.     To continue the example:
  485.  
  486.       Result so far concatenated with 01 is 101, which is greater than 100,
  487.       so result digit is 0 and subtract 0:
  488.  
  489.          1.  0
  490.        --------------------
  491.       V 10. 00 00 00 00 ...
  492.        - 1
  493.         ------
  494.            100
  495.           -  0
  496.            ------
  497.             10000
  498.  
  499.       Result so far concatenated with 01 is 1001, which is less than 10000,
  500.       so result digit is 1 and subtract 1001:
  501.  
  502.          1.  0  1
  503.        --------------------
  504.       V 10. 00 00 00 00 ...
  505.       -  1
  506.        -------
  507.            100
  508.          -   0
  509.           -------
  510.             10000
  511.            - 1001
  512.             --------
  513.                11100
  514.  
  515.       Result so far concatenated with 01 is 10101, which is less than 11100,
  516.       so result digit is 1 and subtract 10101:
  517.  
  518.          1.  0  1  1
  519.        --------------------
  520.       V 10. 00 00 00 00 ...
  521.       -  1
  522.        -------
  523.            100
  524.          -   0
  525.           -------
  526.             10000
  527.            - 1001
  528.             --------
  529.                11100
  530.              - 10101
  531.               ---------
  532.                   11100
  533.  
  534.       Result so far concatenated with 01 is 101101, which is greater than
  535.       11100, so result digit is 0 and subtract 0:
  536.  
  537.          1.  0  1  1  0
  538.        --------------------
  539.       V 10. 00 00 00 00 ...
  540.       -  1
  541.        -------
  542.            100
  543.          -   0
  544.           -------
  545.             10000
  546.            - 1001
  547.             --------
  548.                11100
  549.              - 10101
  550.               ---------
  551.                   11100
  552.                -      0
  553.                 ----------
  554.                    1110000
  555.  
  556. One further point: Step (B) is just a special case of step (C) if we
  557.   initialise the "square root so far" to 0 - i.e. you don't need special
  558.   case hardware for the first iteration.
  559.  
  560. Hope this helps.
  561.  
  562. David Seal
  563. dseal@armltd.co.uk
  564.