// ************************************************************************** // // // // eses eses // // eses eses // // eses eseses esesese eses Embedded Systems Group // // ese ese ese ese ese // // ese eseseses eseseses ese Department of Computer Science // // eses eses ese eses // // eses eseses eseseses eses University of Kaiserslautern // // eses eses // // // // ************************************************************************** // // The division of two n-bit radix-2 numbers [N]bool x,y can be reduced to // // the problem to compute the reciprocal z of y with n-bits, i.e., precision // // 2^{-n} which is a number z with 0 <= 1/y - z <= 2^{-n}. Given z, we can // // compute the quotient either as q=floor(x*z) or q=floor(x*z)+1, depending on// // whether the remainder x-q*y satisfies x-q*y<y. // // According to BeCH84, the reciprocal z can be computed as a geometric sum // // 2^j * sum(i=0..n-1) u^i, where u = 1-y*2^{-j} and j>=2 is that uniquely // // defined integer 2^{j-1}<=y<2^j (i.e. ceil(log(y))). The number u must be // // represented by n-bits, and its range for y>1 is 0 < u <= 0.5. Thus, the // // powers are in the range 0< u^i <= 0.5^i. To avoid loss of precision, BeCH84// // suggested to represent u^i with the entire i*n bits that come out of the // // integer multiplication. The arithmetic operations below are based on a // // fixed point representation having 0 bits to the left and N^2 bits to the // // right of the decimal point. // // Assuming that multiplication of M-bit numbers is done in O(log(M)) time, // // the algorithm below computes the reciprocal in time O(log(N)^2), but it is // // only of theoretical interest due to the large bit numbers: it requires the // // multiplication of N^2-bit numbers. // // In the algorithm below, we use fixed point arithmetic where one bit is on // // the left and N*N-1 bits are on the right of the decimal point. The output // // z contains the N bits of the reciprocal approximation right to the decimal // // point with an error up to 2^{-N}. Note that it may still be the case that // // none of the bits is correct, since 0.1 may be approximated as 0.01111111. // // ************************************************************************** // macro N = 8; macro bval(x) = (x ? 1 : 0); macro nval(x) = sum(i=0..N-1) (bval(x{i}) * exp(2,i)); macro fval(x) = nat2real(nval(x)) * exp(2.0,-N); macro fmul(x,y) = (nat2bv(bv2nat(x) * bv2nat(y),2*N*N){-2:-(N*N+1)}); macro fadd(x,y) = (nat2bv(bv2nat(x) + bv2nat(y),N*N){-1:-(N*N)}); macro cmpl(x) = nat2bv(bv2nat(not x)+1,N); macro rsh(x) = false@(x{-1:1}); module ReciprocalBeCH(bv{N} ?y,bv{N} z,event rdy) { bv{N} u; nat{N+1} j; [N]bv{N*N} p; // determine u as maximal position where y{-1:j} are all zero u = y; j = N; while(!u{-1}) { next(u) = u{-2:0}@false; next(j) = j-1; wu: pause; } next(u) = cmpl(u); nc: pause; // powering, i.e. compute the first N powers of u p[0] = true@{false::N*N-1}; for(i=1..N-1) p[i] = false@u@{false::N*(N-1)-1}; for(i=0..log(N)-1) { for(j=exp(2,i)..N-1) { next(p[j]) = fmul(p[j],p[j-exp(2,i)]); } wp: pause; } // summing up the powers of u for(i=0..log(N)-1) { for(j=exp(2,i)..N-1) { next(p[j]) = fadd(p[j],p[j-exp(2,i)]); } ws: pause; } // Extract the desired n-bits from p[N-1]; here we have the problem // that e.g. 0.5 is 0b010..., but we obtain 0b0011.. which is the same. // To avoid difficult rounding, we allow the error to become 2^{-N} // which happens for powers of 2. z = p[N-1]{-1:-N}; while(j>1) { next(z) = rsh(z); next(j) = j-1; wr : pause; } // signal termination of the computation emit(rdy); // check the correctness {real realz,recy,err; realz = fval(z); recy = 1.0/nat2real(nval(y)); err = recy - realz; assert(realz <= recy); assert(err <= exp(2.0,-N)); } } drivenby { y = nat2bv(10,N); await(rdy); }