about working with fixedpoint calculations

what is it

each number is represented by a 32 bit integer. the fixedpoint position is not explicitly stored. ( as with floating point numbers )

r(A,nA) = A*2^-nA here 'A' is the integer representation of real number r(A,nA) with 'nA' fractional bits.

in the next sections I assume both numbers operated upon have different fractional parts.

• note: I should add notes about loss of precision
• note: I should add notes about rounding

examples in perl

fixedpoint numbers are represented as a hash, with 2 elements,

• 'a' = the integer value.
• 'n' = the number of bits in the fraction

64bit integers are represented as arrays of 2 numbers, first the high word, then the low word.

basic integer operations

sub shl {
    my ($a, $n)= @_;
    return ($a<<$n) if ($n>0);
    return ($a>>-$n) if ($n<0);
    return $a;
}
sub shl64 {
    my ($a, $n)= @_;
    return [ ($a->[0]<<$n)|($a->[1]>>(32-$n)), ($a->[1]<<$n) ] if ($n>0);
    return [ ($a->[0]>>-$n), ($a->[0]<<32+$n)|($a->[1]>>-$n) ] if ($n<0);
    return [ $a->[0], $a->[1] ];
}

sub intadd {
    my ($a, $b)= @_;
    return int($a+$b);
}
sub intsub {
    my ($a, $b)= @_;
    return int($a-$b);
}
sub intmul64 {
    my ($a, $b)= @_;

    my $prod= $a*$b;
    return [ int($prod/2**32), int($prod) ];
}
sub intdiv64 {
    my ($a64, $b)= @_;
    return int(($a64->[0]*2**32 + $a64->[1])/$b);
}

sub fxpshl {
    my ($a, $shift, $n)= @_;
    return { a=>shl($a->{a}, $shift+$n-$a->{n}), n=>$n };
}
# returns  msb such that  2^msb <= $nr < 2*2^msb
sub GetMsb {
    my ($a)= @_;
    my $val= $a->{a};
    my $shifts= 0;
    while($val>1) {
        $val >>=1;
        $shifts++;
    }
    return $shifts;
}

converting between fixedpoint numbers and real numbers

sub fxpdup {
    my ($a, $n)= @_;
    return {a=>shl($a->{a}, $n-$a->{n}), n=>$n};
}

sub tofxp {
    my ($a, $n)= @_;

    return {a=>rint($a*(1<<$n)), n=>$n};
}
sub toval {
    my ($a)= @_;

    return $a->{a}/(1<<$a->{n});
}
sub tostr {
    my ($a)= @_;

    return sprintf("%8.6f[%08lx:%2d]", toval($a), $a->{a}, $a->{n});
}
sub rint {
    my ($a)= @_;
    if ($a>0) {
        return int($a+0.5);
    }
    elsif ($a<0) {
        return int($a-0.5);
    }
    else {
        return 0;
    }

}

arithmetic operations

addition and subtraction

r(C,nC)= r(A,nA) + r(B,nB)

C = ( A+B*2^(nA-nB) ) * 2^(nC-nA)     if nA > nB

C = ( A*2^(nB-nA)+B ) * 2^(nC-nB)     if nA < nB

I think this can be improved upon depending on nC

perl code

sub fxpadd {
    my ($a, $b, $n) = @_;
    if ($b->{n} > $a->{n}) {
        return {a=>shl( intadd(shl($a->{a}, $b->{n}-$a->{n}) , $b->{a}), $n-$b->{n}), n=>$n};
    }
    else {
        return {a=>shl( intadd($a->{a} , shl($b->{a}, $a->{n}-$b->{n})), $n-$a->{n}), n=>$n};
    }
}
sub fxpsub {
    my ($a, $b, $n) = @_;
    if ($b->{n} > $a->{n}) {
        return {a=>shl( intsub(shl($a->{a}, $b->{n}-$a->{n}) , $b->{a}), $n-$b->{n}), n=>$n};
    }
    else {
        return {a=>shl(intsub($a->{a} , shl($b->{a}, $a->{n}-$b->{n})), $n-$a->{n}), n=>$n};
    }
}

multiplication

r(C,nC)= r(A,nA) * r(B,nB)

C = A*B * 2^(nC-nA-nB)

perl code

sub fxpmul {
    my ($a, $b, $n) = @_;

    return {a=>shl64(intmul64($a->{a} , $b->{a}), $n-$a->{n}-$b->{n})->[1], n=>$n};
}

division

r(C,nC)= r(A,nA) / r(B,nB)

C = A / B * 2^(nB+nC-nA)

perl code

sub fxpdiv {
    my ($a, $b, $n)= @_;

    return {a=>intdiv64(shl64([0,$a->{a}], $n-$a->{n}+$b->{n}), $b->{a}), n=>$n};
}

division by multiplication with reciprocal

the reciprocal can be approximated quickly using newton rapson to to solve 'x-1/y=0'

 f(y)= x-1/y   f'(y)= 1/y^2

 y-f(y)/f'(y) = y-xy^2+y = y(2-x*y)

 y[k+1] = y[k] - f(y[k]) / f'(y[k])

  y1= y  *** this initial estimate is done by table lookup.
  y2= y1*(2-x*y1)

• I am not sure how to actually get the initial estimate
• write this in terms of int operations

perl code

sub fxpnewtdiv {
    my ($a, $b, $n)= @_;

    return fxpmul($a, newtreciprocal($b, $n), $n);
}

#  todo: write this in terms of int operations
sub newtreciprocal {
    my ($x, $n)= @_;

    my $y= GetInitialEstimate($x);   ## cheating here

    for (my $i= 0 ; $i < 10 ; $i++)
    {
        $y= fxpmul($y, fxpsub(tofxp(2, $n), fxpmul($x, $y, $n), $n), $n);
    }
    return $y;
}
# initial estimate for reciprocal calc: use highest bit.
#   ( if x=2^n + 2^(n-?) + ... -> return 2^-n
sub GetInitialEstimate {
    my ($x)= @_;

    my $val= toval($x);
    my $l2= int(log($val)/log(2));

    return tofxp(2**-$l2, $x->{n});
}

division by repeated multiplication

• first both numbers are normalized, such that the divisor is between 1 and 2.
• ....... [ finish description ]

perl code

sub fxpmuldiv {
    my ($a, $b, $n)= @_;
    my $msbit=GetMsb($b);
    my $shiftcount= $b->{n}-$msbit;

    my $normb= fxpshl($b, $shiftcount, $b->{n});
    my $norma= fxpshl($a, $shiftcount, $a->{n});

    for (1..6) {
        $b= fxpmul(fxpsub(tofxp(2,$n), $normb, $n), $b, $n);
        $a= fxpmul(fxpsub(tofxp(2,$n), $norma, $n), $a, $n);
    }

    return $a;
}

calculating the natural logarithm

a(i) = 1/(1-2^-i) ,  i=1...

r(L,nL)= ... r(A,nA)

r(A,nA)= 2^scale * prod( a(i) )
     where 'scale' is such that  2<= r(A,nA)*2^-scale < 1

r(L,nL)= log(2^scale * prod( a(i) ) )
       = scale*log(2) + sum( log(a(i)) )

L*2^-nL = log(A*2^-nA) = log(A)-nA*log(2)

perl code

sub fxplog {
    my ($a, $n)= @_;

    if ($a->{a}<=0) {
        return;
    }

    my $base= $n;
    # lg2 contains values [ 2*2^n .. 2^n >
    my @va2= map { rint(2**$base/(1-2**(-$_))) } (1..$base); unshift @va2, 0;
    my @lg2= map { -rint(2**$n*log(1-2**(-$_))) } (1..$base); unshift @lg2, 0;


    my $val= $a->{a};

    # scale $a such that it is just over 2*2^n : bit (n+1) = 1
    my $msb= GetMsb($a);
    $val = shl($val, 1+$base-$msb);
    my $l= -rint((1+$a->{n}-$msb)*log(2)*2**$n); #-$a->{n};

    for (1..$base)
    {
        while ($val>=$va2[$_])
        {
            $val -= ($val >> $_);
            $l += $lg2[$_];
        }
    }

    my $res= {a=>$l, n=>$n};
    return $res;
}

calculating a square root

r(L,nL)= sqrt(r(A,nA))

perl code

sub fxpsqrt {
    my ($a, $n)= @_;

    # a= (a+1)/2  - initial estimate
    my $sqrt= fxpshl(fxpadd($a, tofxp(1,$n), $n), -1, $n);

    for (0..7) {
        $sqrt= fxpshl(fxpadd($sqrt, fxpdiv($a, $sqrt, $n), $n), -1, $n);
    }

    return $sqrt;
}

calculating the exponential

r(L,nL)= exp(r(A,nA))

perl code

sub fxpexp {
    my ($a, $n)= @_;

    my @fact= (1);
    push @fact, $fact[-1]*$_ for (1..8);
    my @fxfact= map { tofxp($_, $n) } @fact;


    my $result= tofxp(1, $n);
    my $xn= fxpdup($a, $n);

    for (1..8) {
        $result= fxpadd($result, fxpdiv($xn, $fxfact[$_], $n), $n);
        $xn= fxpmul($xn, $a, $n);
    }

    return $result;
}

calculating a exponantiation

# pow can be implemented as follows:
#
#
# a^sum(2^b(i)) = prod ( a^(2^b(i)) )
#
# res= 1;
# aa= a;
# for (i=fixedpoint ; i<32 ; i++)
# {
#    if (b&bit(i))
#    {
#       res *= aa;
#    }
#    aa *= aa;
# }
# aa= sqrt(a);
# for (i=fixedpoint-1 ; i>=0 ; i--)
# {
#    if (b&bit(i))
#    {
#       res *= aa;
#    }
#    aa = sqrt(aa);
# }
#
#
# a^b can then be implemented as  exp2(b*log2(a))
#

perl code

- none yet.

calculating trigoniometric functions

perl code

- none yet