home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1993 #3 / NN_1993_3.iso / spool / comp / sys / mac / programm / 22280 < prev    next >
Encoding:
Internet Message Format  |  1993-01-27  |  4.0 KB
  1. From: bayes@hplvec.LVLD.HP.COM (Scott Bayes)
  2. Date: Mon, 25 Jan 1993 23:50:39 GMT
  3. Subject: Re: faster division routine wanted.
  4. Message-ID: <830027@hplvec.LVLD.HP.COM>
  5. Organization: Hewlett-Packard Co., Loveland, CO
  6. Path: sparky!uunet!pageworks.com!world!eff!sol.ctr.columbia.edu!spool.mu.edu!sdd.hp.com!hpscit.sc.hp.com!hplextra!hpfcso!hplvec!bayes
  7. Newsgroups: comp.sys.mac.programmer
  8. References: <5662@daily-planet.concordia.ca>
  9. Lines: 99
  10.  
  11. I wrote a routine that seems almost correct, and that allows you to
  12. quickly calculate a set of shifts and add/subs that approximate the
  13. division by any divisor. I know there's some stuff wrong with it, but
  14. thought it's worth sending out now, while the iron is hot.
  15.  
  16. It's in HP RMB (BASIC), which is great for prototyping in, though
  17. a little non-standard, and lacks a 32-bit integer.
  18.  
  19. Herewith.
  20.  
  21. ScottB
  22.  
  23. 10    ! RE-SAVE "SHIFT-ADD-DIVIDE"
  24. 20    !
  25. 30    ! Purpose
  26. 40    ! -------
  27. 50    ! Find a continued fraction whose denominators are
  28. 60    ! all powers of 2, and that approximates the
  29. 70    ! desired divisor using a given number of terms, or
  30. 80    ! including all terms whose denominator < 2^32
  31. 90    !
  32. 100   ! Can be used to evaluate which shifts and adds
  33. 110   ! will compute a division arbitrarily close to
  34. 120   ! a given divisor.
  35. 130   !
  36. 140   ! Results returned in A(*). Values in A(*) are exponents
  37. 150   ! of 2 (right shifts); the signs of the values indicate
  38. 160   ! whether to add or subtract the shifted value to the
  39. 170   ! continued fraction sum.
  40. 180   !
  41. 190   ! E.g. asking for a divisor of 1000 to 4 terms gives the
  42. 200   ! values: 10, 15, -17, 21.
  43. 210   ! This implies the algorithm for X / 1000 approximation
  44. 220   ! with 4 terms is:
  45. 230   ! X/1000 = X>>10 + X>>15 - X>>17 + X>>21
  46. 240   ! The actual division this performs is very close to:
  47. 250   ! X/1000.0724845 (as shown by "approximation" printout.)
  48. 260   !
  49. 270   ! I haven't done the numeric analysis to find out exactly
  50. 280   ! how close the approximation is, but it should be near
  51. 290   ! the approximation value I print.
  52. 300   ! Give me a break; I did it between phone calls in tech
  53. 310   ! support
  54. 320   !
  55. 330   ! I chose HP'S "Rocky Mountain BASIC" because it's quick
  56. 340   ! and easy to prototype in. I hope most of the stuff
  57. 350   ! is obvious. PROUND(N,M) means round the (decimal) number
  58. 360   ! N to the (decimal) place M. So PROUND(X,-1) means
  59. 370   ! round X to the .1 position, and PROUND(X,0) means ROUND(X).
  60. 380   !
  61. 390   ! INTEGERs are only 2 bytes (2's complement). REALS are
  62. 400   ! 8-byte IEEE floating.
  63. 410   !
  64. 420   ! Scott Bayes
  65. 430   ! bayes@hpisla.lvld.hp.com
  66. 440   !
  67. 450   !
  68. 460   INTEGER A(99),I      ! A(*) STORES TERMS OF THE FRACTION
  69. 470   REAL T,D,S,E
  70. 480   INPUT "Divisor",D
  71. 490   INPUT "# OF TERMS (0 TO GET ALL POSSIBLE TERMS)",N
  72. 500   E=D                  ! E IS THE "ERROR" RESIDUAL
  73. 510   I=0                  ! I COUNTS TERMS GENERATED
  74. 520   LOOP
  75. 530     A(I)=FNClosest2(E) ! FIND N SUCH THAT 2^N CLOSEST TO E
  76. 540     ! THE ABOVE COULD PROBABLY BE IMPROVED SINCE IT FINDS THE
  77. 550     ! ARITHMETICALLY CLOSEST POWER OF 2. i THINK (WITHOUT ANY
  78. 560     ! ANALYSIS TO BACK IT UP) THAT THE GEOMETRICALLY CLOSEST
  79. 570     ! WOULD BE OPTIMAL.
  80. 580     !
  81. 590   EXIT IF ABS(A(I))>31 OR ((I>=N-1) AND (N<>0)) ! EXITS THE LOOP
  82. 600     S=FNSum(A(*),I)    ! EVALUATOR FUNCTION
  83. 610     I=I+1
  84. 620   EXIT IF S=D          ! IN CASE D=2^N FOR SOME N
  85. 630     E=D*S/(S-D)        ! NEW RESIDUAL TERM
  86. 640   END LOOP
  87. 650   PRINT I+1;"TERMS:";
  88. 660   FOR J=0 TO I
  89. 670     PRINT A(J);
  90. 680   NEXT J
  91. 690   PRINT 
  92. 700   T=FNSum(A(*),I)      ! T IS THE EFFECTIVE DENOMINATOR
  93. 710   PRINT "APPROXIMATION=";T
  94. 720   END
  95. 730   !
  96. 740   DEF FNClosest2(Y)    ! FIND N SUCH THAT 2^N IS CLOSEST TO ABS(Y)
  97. 750     X=ABS(Y)           ! SIGN OF THE RESULTS FLAGS THAT WE SHOULD
  98. 760                        ! SUBTRACT THE TERM
  99. 770     RETURN SGN(Y)*PROUND(LOG(X)/LOG(2),0)
  100. 780   FNEND
  101. 790   !
  102. 800   DEF FNSum(INTEGER X(*),N) ! EVALUATE THE CONTINUED FRACTION
  103. 810     REAL S
  104. 820     S=0                     ! SUMS SHIFTED TERMS 1/2^A(I)
  105. 830     FOR I=0 TO N
  106. 840       S=S+2^(-ABS(X(I)))*SGN(X(I))
  107. 850     NEXT I
  108. 860     RETURN 1/S
  109. 870   FNEND
  110.