about working with fixedpoint calculations

WARNING code on this page may not work properly, also the algorithms described may be incorrect, this is a work in progress.

what is it

In most embedded or older systems there is no floatingpoint processor present, to be able to still do calculations with 'real' numbers, you can pretend that each integer is divided by a fixed number, shifting the position of 'bit0' to the left.

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.

examples in perl

all operations have perl code as example implementation. here is the 'integer' library that this code is based on.

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

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);
}
here are some generic fixed point operations:
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
Subtraction is exactly the same.

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

in order not to loose precision, It is advisable that the multiplication is done on a double size word ( ie. 64 bits, if all numbers are stored as 32 bits )
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

It may sometimes be faster to calculate the inverse of the divisor, and then multiply the dividend by this inverse, than to do the actual division, with shifts, and substracts.

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)

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

another way to do a division, is to multiplying both the dividend and divisor by the same number, such that the divisor gets closer to '1'. and repeat this until we have the desired precision.

perl code

WARNING:!!! I don't think the code below is correct. !!!!
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

[todo- add more description]
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

a square root is also approximated using newton rapson
r(L,nL)= sqrt(r(A,nA))

[todo- more description]

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))
this code uses a taylor expansion, to calculate the exponential. unfortunately it is not very accurate.

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