home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #27 / NN_1992_27.iso / spool / comp / lang / rexx / 1326 < prev    next >
Encoding:
Internet Message Format  |  1992-11-19  |  2.1 KB
  1. Path: sparky!uunet!mcsun!uknet!comlab.ox.ac.uk!comlab.ox.ac.uk!imc
  2. From: imc@comlab.ox.ac.uk (Ian Collier)
  3. Newsgroups: comp.lang.rexx
  4. Subject: Square roots and things
  5. Message-ID: <2741.imc@uk.ac.ox.prg>
  6. Date: 19 Nov 92 18:27:23 GMT
  7. Sender: imc@comlab.ox.ac.uk (Ian Collier)
  8. Organization: Oxford University Computing Laboratory, UK
  9. Lines: 52
  10. X-Local-Date: Thursday, 19th November 1992 at 6:27pm GMT
  11. Originator: imc@msc5.comlab
  12.  
  13. This is a follow-up to a long article I posted previously about calculating
  14. a square root.  Sorry if this bores you...
  15.  
  16. Anyway, my digit-by-digit method has been beaten fair-and-square by Newton's
  17. method, because I wrote a sqrt function which
  18.  
  19.  * makes a good estimate based on the mantissa and exponent of the argument
  20.  * starts off at low precision and, after the low precision answer has been
  21.    found, doubles the precision at each step.
  22.  
  23. These together make a sqrt function which has the same complexity as a
  24. divide instruction.  So, in cases where decimal divide is a native
  25. instruction (including any code written in REXX) this method will almost
  26. always be faster than my digit-by-digit method.
  27.  
  28. Tests gave (in elapsed seconds):
  29.  
  30. digits   digit_method   Newton's_method
  31.   10       0.59           0.13
  32.   50       3.68           0.72
  33.  100       8.63           2.19
  34.  200      26.29           7.87
  35.  400      89.26          30.31
  36.  800     304.14         125.23
  37.  
  38. Here is the new sqrt() function:
  39.  
  40. sqrt:procedure
  41. parse numeric d . /* or d=digits() */
  42. numeric form scientific
  43. parse arg x
  44. if x=0 then return 0
  45. parse value format(x,1,1,,0) with n +1 'E' exp /* error if x<0 */
  46. if exp='' then exp=0
  47. r=n/2||'E'||(exp+1)%2 /* divide both mantissa and exponent by 2 for initial
  48.                          estimate */
  49. d1=5
  50. numeric digits d1+2
  51. numeric fuzz 1
  52. do until r = r1       /* Get estimate accurate up to d1 digits */
  53.    parse value r (x/r+r)/2  with r1 r
  54. end
  55. do while d1<d
  56.    d1=min(d1+d1,d)    /* Double up until required accuracy is reached */
  57.    numeric digits d1+2
  58.    parse value r (x/r+r)/2  with r1 r
  59. end
  60. numeric digits d
  61. return format(r)
  62.  
  63. Ian Collier
  64. Ian.Collier@prg.ox.ac.uk | imc@ecs.ox.ac.uk
  65.