Israel Koren

UNIVERSITY OF MASSACHUSETTS Dept. of Electrical & Computer Engineering Digital Computer Arithmetic ECE 666 Part 9a Evaluation of Elementary Functions - I

• Israel Koren

Efficient Algorithms for Evaluating Elementary Functions

• Hardware algorithms which are accurate and fast compared to software algorithms

• Useful for scientific hand-held calculators and numerical processors

• Examples: e , ln x, sin x, cos x

• Straightforward method: look-up table

• x,y of length n bits

• To evaluate y=F(x) - ROM of size 2  n

• Prohibitive for n  20 - size  16 M

• Second method: Taylor series expansion

Taylor Series Expansion

• Example:

• Converges rapidly when x is a fraction

• Converges slowly for x close to 1

• Hardware implementation complex - separate logical networks for different functions

• Algorithms requiring a smaller number of computational steps for same precision exist

• Most commonly used methods: based on either polynomial or rational approximations

Example - Exponential Function

• Exponent divided into integer part I and fractional part f : x log 2 e = I+f

• Calculation of 2 in either fixed-point or floating-point - straightforward

• Calculation of 2 - rational approximation - ratio between 2 polynomials

• Example: Two degree-5 polynomials

• ai,bi (i=0,1,...,5) known constants - 10 multiplications, 10 additions, 1 division

• Approximations with smaller execution time exist

Alternative Methods for Evaluating Elementary Functions

• Similar to division-by-convergence

• 2 (or more) recursive formulas - one is forced to a constant while other yields result

• Example: Evaluating y=e for some fraction x0

• bi's selection - x0,x1,…,xm approach 0

• m - number of iterations needed to ensure xm=0

• Linear convergence - m linear function of number of bits n

Exponential Function

• In particular:

• xm=0 

• Similar to division-by-convergence - product yi e remains constant for all bi's

• Need a simple way to select bi's to ensure convergence of xi+1 to 0

• Also, multiplication by bi should be simple and fast

• Use the form bi =(1+si 2 ), where si  {-1,0,1}
• Multiplication reduced to shift and add

• The term ln bi=ln (1+si2 ) can be positive or negative

• Must ensure convergence of positive or negative xi+1 to 0

Precalculated ln(12 )

• Online calculation too slow

• Stored in a look-up table (ROM of size 2nn)

• (i) xi+1=xi -ln (1+si 2 )

• (ii) yi+1=yi (1+si 2 )

• To calculate e - find the vector

• s={s0,s1,…,sm-1} so that xm=0 or

• If sl restricted to {-1,0,1} - smallest/largest x0 correspond to all sl =-1 / sl =1

• For every x0 in interval there exists vector s guaranteeing convergence to 0

Domain of x0

• For large m (20) bounds converge: -1.24  x0  1.56

• If x0 positive fraction - simple scheme for selecting si - one-sided selection rule: si  {0,1}

• Similar to quotient-bit selection in restoring division

• In step i+1 form difference D=xi-ln (1+2 )

• If D positive or zero - set si=1 and xi+1=D
• If D negative - set si=0 and xi+1=xi

• Analyze subtract in xi+1 - examine Taylor series:

• At step i can cancel bit of weight 2 by subtracting ln (1+si 2 ) from xi
• Up to n steps needed to ensure convergence of n-bit fraction to 0 - convergence linear and m=n

Modify Selection Rule

• Improves performance of algorithm

• If bit in ith position =1 - select si=1 and calculate xi+1=xi-ln (1+si 2 )

• If 0 - select si=0, go to next step without subtract - capability of skipping over 0’s

• More complex selection rules can be used as well

Example: e in 10-bit precision

• One-sided selection rule

• x0=0.2510= 0.01000 0000002

• x0

• Select s2=1 and obtain x3=x2-0.00111 00100 = 0.00000 11100

• Set s3=s4=0  s5=0 , s6=s7=s8=1, and s9=s10=0

• Yielding x11=0

Example - Cont.

• In parallel :y11=y0(1+2 )(1+2 )(1+2 )(1+2 ) yielding y11=1.01001 000112 = 1.2841810

• Exact value (to 5 five decimal digits)=1.2840310

• Approximation error =0.158·2

• If si=-1 allowed - selected values (s0,…,s10)=0,1,-1,1,0,0,0,1,1,1,1

• x11=0 ; y11=1.01001 000102

• Approximation error=0.842·2

Extending Argument of Exponential Function

• x = x log 2 e ln 2

• x log 2 e = I + f (I - integer ; 0  f < 1)

• Calculate y=e =e =e =2 e

• Set x0=f ln 2  0  x0 < ln 2 =0.693

• 2 is easily dealt with

• Incorporation into exponent part of floating-point number
• Through a shift operation for fixed-point arithmetic

The Logarithm Function

• Calculating e - continued summation of terms ln (1+si 2 ) to force convergence of xi to 0 - additive normalization

• Multiplicative normalization - xi is forced to 1 (or some other nonzero constant) by continued multiplication with precalculated factors

• bi selected so that

• After m steps (multiplicative inverse)

Domain of Algorithm

• Multiplication factor bi: 1+si 2 for si{-1,0,1}

• Domain:

• For large m:

• Domain for x0: 0.21  x0  3.45

• Positive normalized fractions are in domain

• If argument not normalized floating-point number - rewrite as x=x0·2 ; 0.5  x0<1

• ln x=ln x0+Ex ln 2

One Sided Selection Rule si{0,1}

• has i leading 1's  xi+1=xibi=xi+xisi2 (with si=1) has (i+1) leading 1's

• In yi+1 formula:

• g(bl)=ln(bl) 

• xi+1 approaches 1  yi+1 converges to

Calculating Natural Logarithm

• Same table of ln (1+si 2 ) needed here

• Base e can be replaced by any other base

• Base 10 - used in hand-held calculators

• Finally:

• yi+1 + ln xi+1 = (yi-ln bi)+(ln xi+ln bi)= yi+ln xi is a constant

• xi+1 approaches 1  yi+1 approaches y0+ln x0

Example - ln (0.50) in 10-bit precision

• y11=-0.10110 001012=-0.6923810

• Exact result (in 5 decimal digits) = -0.69315

• Approximation error = 0.783·2

Operation in (ii) - Multiplication

• g(bi)=bi (k - any real number):

• Example:

• k=1 : yi+1  y0 / x0 - divide operation
• k=1/2 : yi+1  y0 /  x0 - reciprocal of square root
• y0=x0 : yi+1   x0 - square root calculated

• g(bi)= bi ; for simplicity g(bi)=(1+si 2 )

• Multiplication by bi more complex - not efficient for square root extraction

Trigonometric Functions

• E = cos x + j sin x where j= -1

• For exponential function:

• For trigonometric functions select bi=(1+jsi 2 )

• This complex number can be written as

• where

• Polar form of complex numbers

• where

Domain for Trigonometric Functions

• si's selected so that x0,x1,…,xm  0

• 

• Denoting:

• Then:

• Convergence domain for m  20 : (including 0  x0 /2 = 1.57)

Equations for x and y

• xi's real numbers - yi's complex numbers

• yi=Zi + jWi (Zi - real part ; Wi - imaginary part)

• Initial value y0 can be complex - Z0+jW0

• Setting W0=0, Z0=y0 - desired values obtained:

• Volder developed differently using rotations in polar system - CORDIC (COordinated Rotation DIgital Computer)

Selection Rule for Si

• tan (si 2 )= si tan (2 )

• n constants stored in ROM

• For si {-1,0,1}

• K not constant - depends on s which depends on x0

• Also, full-length comparison between xi & tan (2 )

• If si restricted to {-1,1} -

• Constant that can be precalculated

• For m > 16, K=1.6468 - set y0=1/K Zm =cos x0 ; Wm = sin x0

Selection Rule - Modified

• Selection rule becomes

• Similar to rule for nonrestoring division

• Examine convergence rate of xi+1 to 0 - Taylor series expansion of tan (si2 ) (in radians):

• Linear convergence

• For i>n/3, all terms except first negligible - reduce size of ROM

Example: sin(/4) in 10-bit precision

• /4=0.78539810=0.11001001002

• Initially Z0=1/K=0.10011011102 (=0.607410)

• Set si=sign(xi)

• Final results: Z11=0.10110100112 and W11=0.10110101002

• “exact” results in 10-bit precision: both 0.10110101002= 0.707110

• Approximation error not larger than 2