Newsgroups: comp.sys.mac.programmer

I wrote a routine that seems almost correct, and that allows you to
quickly calculate a set of shifts and add/subs that approximate the
division by any divisor. I know there's some stuff wrong with it, but
thought it's worth sending out now, while the iron is hot.

It's in HP RMB (BASIC), which is great for prototyping in, though
a little non-standard, and lacks a 32-bit integer.

Herewith.

ScottB

10    ! RE-SAVE "SHIFT-ADD-DIVIDE"
20    !
30    ! Purpose
40    ! -------
50    ! Find a continued fraction whose denominators are
60    ! all powers of 2, and that approximates the
70    ! desired divisor using a given number of terms, or
80    ! including all terms whose denominator < 2^32
90    !
100   ! Can be used to evaluate which shifts and adds
110   ! will compute a division arbitrarily close to
120   ! a given divisor.
130   !
140   ! Results returned in A(*). Values in A(*) are exponents
150   ! of 2 (right shifts); the signs of the values indicate
160   ! whether to add or subtract the shifted value to the
170   ! continued fraction sum.
180   !
190   ! E.g. asking for a divisor of 1000 to 4 terms gives the
200   ! values: 10, 15, -17, 21.
210   ! This implies the algorithm for X / 1000 approximation
220   ! with 4 terms is:
230   ! X/1000 = X>>10 + X>>15 - X>>17 + X>>21
240   ! The actual division this performs is very close to:
250   ! X/1000.0724845 (as shown by "approximation" printout.)
260   !
270   ! I haven't done the numeric analysis to find out exactly
280   ! how close the approximation is, but it should be near
290   ! the approximation value I print.
300   ! Give me a break; I did it between phone calls in tech
310   ! support
320   !
330   ! I chose HP'S "Rocky Mountain BASIC" because it's quick
340   ! and easy to prototype in. I hope most of the stuff
350   ! is obvious. PROUND(N,M) means round the (decimal) number
360   ! N to the (decimal) place M. So PROUND(X,-1) means
370   ! round X to the .1 position, and PROUND(X,0) means ROUND(X).
380   !
390   ! INTEGERs are only 2 bytes (2's complement). REALS are
400   ! 8-byte IEEE floating.
410   !
420   ! Scott Bayes
430   ! bayes@hpisla.lvld.hp.com
440   !
450   !
460   INTEGER A(99),I      ! A(*) STORES TERMS OF THE FRACTION
470   REAL T,D,S,E
480   INPUT "Divisor",D
490   INPUT "# OF TERMS (0 TO GET ALL POSSIBLE TERMS)",N
500   E=D                  ! E IS THE "ERROR" RESIDUAL
510   I=0                  ! I COUNTS TERMS GENERATED
520   LOOP
530     A(I)=FNClosest2(E) ! FIND N SUCH THAT 2^N CLOSEST TO E
540     ! THE ABOVE COULD PROBABLY BE IMPROVED SINCE IT FINDS THE
550     ! ARITHMETICALLY CLOSEST POWER OF 2. i THINK (WITHOUT ANY
560     ! ANALYSIS TO BACK IT UP) THAT THE GEOMETRICALLY CLOSEST
570     ! WOULD BE OPTIMAL.
580     !
590   EXIT IF ABS(A(I))>31 OR ((I>=N-1) AND (N<>0)) ! EXITS THE LOOP
600     S=FNSum(A(*),I)    ! EVALUATOR FUNCTION
610     I=I+1
620   EXIT IF S=D          ! IN CASE D=2^N FOR SOME N
630     E=D*S/(S-D)        ! NEW RESIDUAL TERM
640   END LOOP
650   PRINT I+1;"TERMS:";
660   FOR J=0 TO I
670     PRINT A(J);
680   NEXT J
690   PRINT 
700   T=FNSum(A(*),I)      ! T IS THE EFFECTIVE DENOMINATOR
710   PRINT "APPROXIMATION=";T
720   END
730   !
740   DEF FNClosest2(Y)    ! FIND N SUCH THAT 2^N IS CLOSEST TO ABS(Y)
750     X=ABS(Y)           ! SIGN OF THE RESULTS FLAGS THAT WE SHOULD
760                        ! SUBTRACT THE TERM
770     RETURN SGN(Y)*PROUND(LOG(X)/LOG(2),0)
780   FNEND
790   !
800   DEF FNSum(INTEGER X(*),N) ! EVALUATE THE CONTINUED FRACTION
810     REAL S
820     S=0                     ! SUMS SHIFTED TERMS 1/2^A(I)
830     FOR I=0 TO N
840       S=S+2^(-ABS(X(I)))*SGN(X(I))
850     NEXT I
860     RETURN 1/S
870   FNEND