home *** CD-ROM | disk | FTP | other *** search
/ Frozen Fish 1: Amiga / FrozenFish-Apr94.iso / bbs / gnu / libm-src.lha / src / amiga / libm / national / sqrt.s < prev    next >
Text File  |  1993-09-23  |  7KB  |  233 lines
  1. ; Copyright (c) 1985 Regents of the University of California.
  2. ; All rights reserved.
  3. ;
  4. ; Redistribution and use in source and binary forms, with or without
  5. ; modification, are permitted provided that the following conditions
  6. ; are met:
  7. ; 1. Redistributions of source code must retain the above copyright
  8. ;    notice, this list of conditions and the following disclaimer.
  9. ; 2. Redistributions in binary form must reproduce the above copyright
  10. ;    notice, this list of conditions and the following disclaimer in the
  11. ;    documentation and/or other materials provided with the distribution.
  12. ; 3. All advertising materials mentioning features or use of this software
  13. ;    must display the following acknowledgement:
  14. ;    This product includes software developed by the University of
  15. ;    California, Berkeley and its contributors.
  16. ; 4. Neither the name of the University nor the names of its contributors
  17. ;    may be used to endorse or promote products derived from this software
  18. ;    without specific prior written permission.
  19. ;
  20. ; THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  21. ; ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  22. ; IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  23. ; ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  24. ; FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  25. ; DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  26. ; OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  27. ; HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  28. ; LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  29. ; OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  30. ; SUCH DAMAGE.
  31. ;
  32. ;    @(#)sqrt.s    5.4 (Berkeley) 10/9/90
  33. ;
  34.  
  35. ; double sqrt(x)
  36. ; double x;
  37. ; IEEE double precision sqrt
  38. ; code in NSC assembly by K.C. Ng
  39. ; 12/13/85
  40. ;
  41. ; Method:
  42. ;     Use Kahan's trick to get 8 bits initial approximation
  43. ;    by integer shift and add/subtract. Then three Newton
  44. ;    iterations to get down error to within one ulp. Finally
  45. ;    twiddle the last bit to make to correctly rounded 
  46. ;    according to the rounding mode.
  47. ;
  48.     .vers    2
  49.     .text
  50.     .align    2
  51.     .globl    _sqrt
  52. _sqrt:
  53.     enter    [r3,r4,r5,r6,r7],44
  54.     movl    f4,tos
  55.     movl    f6,tos
  56.     movd    2146435072,r2    ; r2 = 0x7ff00000
  57.     movl    8(fp),f0    ; f2 = x 
  58.     movd    12(fp),r3    ; r3 = high part of x
  59.     movd    r3,r4        ; make a copy of high part of x in r4
  60.     andd    r2,r3        ; r3 become the bias exponent of x
  61.     cmpd    r2,r3        ; if r3 = 0x7ff00000 then x is INF or NAN
  62.     bne    L22
  63.                 ; to see if x is INF
  64.     movd    8(fp),r0    ; r0 = low part of x
  65.     movd    r4,r1        ; r1 is high part of x again
  66.     andd    0xfff00000,r1    ; mask off the sign and exponent of x
  67.     ord    r0,r1        ; or with low part, if 0 then x is INF
  68.     cmpqd    0,r1        ; 
  69.     bne    L1        ; not 0; therefore x is NaN; return x.
  70.     cmpqd    0,r4        ; now x is Inf, is it +inf?
  71.     blt    L1        ; +INF, return x 
  72.                 ; -INF, return NaN by doing 0/0
  73. nan:    movl    0f0.0,f0    ; 
  74.     divl    f0,f0
  75.     br    L1
  76. L22:                ; now x is finite
  77.     cmpl    0f0.0,f0    ; x = 0 ?
  78.     beq    L1        ; return x if x is +0 or -0
  79.     cmpqd    0,r4        ; Is x < 0 ?
  80.     bgt    nan        ; if x < 0 return NaN
  81.     movqd    0,r5        ; r5 == scalx initialize to zero
  82.     cmpqd    0,r3        ; is x is subnormal ?  (r3 is the exponent)
  83.     bne    L21        ; if x is normal, goto L21
  84.     movl    L30,f2        ; f2 = 2**54
  85.     mull    f2,f0        ; scale up x by 2**54
  86.     subd    0x1b00000,r5    ; off set the scale factor by -27 in exponent
  87. L21:
  88.                 ; now x is normal
  89.                 ; notations:
  90.                 ;    r1 == copy of fsr
  91.                 ;    r2 == mask of e inexact enable flag
  92.                 ;    r3 == mask of i inexact flag
  93.                 ;    r4 == mask of r rounding mode
  94.                 ;    r5 == x's scale factor (already defined)
  95.  
  96.     movd    0x20,r2
  97.     movd    0x40,r3
  98.     movd    0x180,r4
  99.     sfsr    r0        ; store fsr to r0
  100.     movd    r0,r1        ; make a copy of fsr to r1
  101.     bicd    [5,6,7,8],r0    ; clear e,i, and set r to round to nearest
  102.     lfsr    r0
  103.  
  104.                 ; begin to compute sqrt(x)
  105.     movl    f0,-8(fp)
  106.     movd    -4(fp),r0    ; r0 the high part of modified x
  107.     lshd    -1,r0        ; r0 >> 1
  108.     addd    0x1ff80000,r0    ; add correction to r0  ...got 5 bits approx.
  109.     movd    r0,r6
  110.     lshd    -13,r6        ; r6 = r0>>-15
  111.     andd    0x7c,r6        ; obtain 4*leading 5 bits of r0
  112.     addrd    L29,r7        ; r7 = address of L29 = table[0]
  113.     addd    r6,r7        ; r6 = address of L29[r6] = table[r6]
  114.     subd    0(r7),r0    ; r0 = r0 - table[r6]
  115.     movd    r0,-4(fp)
  116.     movl    -8(fp),f2    ; now f2 = y approximate sqrt(f0) to 8 bits
  117.  
  118.     movl    0f0.5,f6    ; f6 = 0.5
  119.     movl    f0,f4
  120.     divl    f2,f4        ; t = x/y
  121.     addl    f4,f2        ; y = y + x/y
  122.     mull    f6,f2        ; y = 0.5(y+x/y) got 17 bits approx.
  123.     movl    f0,f4
  124.     divl    f2,f4        ; t = x/y
  125.     addl    f4,f2        ; y = y + x/y
  126.     mull    f6,f2        ; y = 0.5(y+x/y) got 35 bits approx.
  127.     movl    f0,f4
  128.     divl    f2,f4        ; t = x/y
  129.     subl    f2,f4        ; t = x/y - y
  130.     mull    f6,f4        ; t = 0.5(x/y-y)
  131.     addl    f4,f2        ; y = y + 0.5(x/y -y) 
  132.                 ; now y approx. sqrt(x) to within 1 ulp
  133.  
  134.                 ; twiddle last bit to force y correctly rounded
  135.     movd    r1,r0        ; restore the old fsr
  136.     bicd    [6,7,8],r0    ; clear inexact bit but retain inexact enable
  137.     ord    0x80,r0        ; set rounding mode to round to zero
  138.     lfsr    r0
  139.  
  140.     movl    f0,f4
  141.     divl    f2,f4        ; f4 = x/y
  142.     sfsr    r0
  143.     andd    r3,r0        ; get the inexact flag
  144.     cmpqd    0,r0
  145.     bne    L18
  146.                 ; if x/y exact, then ...
  147.     cmpl    f2,f4        ; if y == x/y 
  148.     beq    L2
  149.     movl    f4,-8(fp)
  150.     subd    1,-8(fp)
  151.     subcd    0,-4(fp)    
  152.     movl    -8(fp),f4    ; f4 = f4 - ulp
  153. L18:
  154.     bicd    [6],r1
  155.     ord    r3,r1        ; set inexact flag in r1
  156.  
  157.     andd    r1,r4        ; r4 = the old rounding mode
  158.     cmpqd    0,r4        ; round to nearest?
  159.     bne    L17
  160.     movl    f4,-8(fp)
  161.     addd    1,-8(fp)
  162.     addcd    0,-4(fp)
  163.     movl    -8(fp),f4    ; f4 = f4 + ulp
  164.     br    L16
  165. L17:
  166.     cmpd    0x100,r4    ; round to positive inf ?
  167.     bne    L16
  168.     movl    f4,-8(fp)    
  169.     addd    1,-8(fp)
  170.     addcd    0,-4(fp)
  171.     movl    -8(fp),f4    ; f4 = f4 + ulp
  172.  
  173.     movl    f2,-8(fp)    
  174.     addd    1,-8(fp)
  175.     addcd    0,-4(fp)
  176.     movl    -8(fp),f2    ; f2 = f2 + ulp
  177. L16:
  178.     addl    f4,f2        ; y  = y + t
  179.     subd    0x100000,r5    ; scalx = scalx - 1
  180. L2:
  181.     movl    f2,-8(fp)    
  182.     addd    r5,-4(fp)
  183.     movl    -8(fp),f0
  184.     lfsr    r1
  185. L1:
  186.     movl    tos,f6
  187.     movl    tos,f4
  188.     exit    [r3,r4,r5,r6,r7]
  189.     ret    0
  190.     .data
  191. L28:    .byte    64,40,35,41,115,113,114,116,46,99
  192.     .byte    9,49,46,49,32,40,117,99,98,46
  193.     .byte    101,108,101,102,117,110,116,41,32,57
  194.     .byte    47,49,57,47,56,53,0
  195. L29:    .blkb    4
  196.     .double    1204
  197.     .double    3062
  198.     .double    5746
  199.     .double    9193
  200.     .double    13348
  201.     .double    18162
  202.     .double    23592
  203.     .double    29598
  204.     .double    36145
  205.     .double    43202
  206.     .double    50740
  207.     .double    58733
  208.     .double    67158
  209.     .double    75992
  210.     .double    85215
  211.     .double    83599
  212.     .double    71378
  213.     .double    60428
  214.     .double    50647
  215.     .double    41945
  216.     .double    34246
  217.     .double    27478
  218.     .double    21581
  219.     .double    16499
  220.     .double    12183
  221.     .double    8588
  222.     .double    5674
  223.     .double    3403
  224.     .double    1742
  225.     .double    661
  226.     .double    130
  227. L30:    .blkb    4
  228.     .double    1129316352     ;L30:    .double 0,0x43500000
  229. L31:    .blkb    4
  230.     .double 0x1ff00000
  231. L32:    .blkb    4
  232.     .double 0x5ff00000
  233.