MiniSoftFP32 - A simple small software 32-bit single precision floating point implementation

MiniSoftFP32.pas

unit MiniSoftFP32; // Copyright (C) 2018, Benjamin "BeRo" Rosseaux (benjamin@rosseaux.de) - License: CC0
// Declaimer / Notice of caution:
// Attention, this code implements only the basic functions, but for example not the correct handling of
// Infinity, NaN, division-by-zero special cases and so on.
// In short, this code is only intended for demystifying the base floating point arithmetics (using 32-bit
// single precision floating point values in this implementation).
{$ifdef fpc}
 {$mode delphi}
 {$if defined(cpu386) or defined(cpuamd64)}
  {$asmmode intel}
 {$ifend}
{$endif}

interface

type PPmsfp32Int8=^Pmsfp32Int8;
     Pmsfp32Int8=^Tmsfp32Int8;
     Tmsfp32Int8={$ifdef fpc}Int8{$else}shortint{$endif};

     PPmsfp32UInt8=^Pmsfp32UInt8;
     Pmsfp32UInt8=^Tmsfp32UInt8;
     Tmsfp32UInt8={$ifdef fpc}UInt8{$else}byte{$endif};

     PPmsfp32Int16=^Pmsfp32Int16;
     Pmsfp32Int16=^Tmsfp32Int16;
     Tmsfp32Int16={$ifdef fpc}Int16{$else}smallint{$endif};

     PPmsfp32UInt16=^Pmsfp32UInt16;
     Pmsfp32UInt16=^Tmsfp32UInt16;
     Tmsfp32UInt16={$ifdef fpc}UInt16{$else}word{$endif};

     PPmsfp32Int32=^Pmsfp32Int32;
     Pmsfp32Int32=^Tmsfp32Int32;
     Tmsfp32Int32={$ifdef fpc}Int32{$else}longint{$endif};

     PPmsfp32UInt32=^Pmsfp32UInt32;
     Pmsfp32UInt32=^Tmsfp32UInt32;
     Tmsfp32UInt32={$ifdef fpc}UInt32{$else}longword{$endif};

     PPmsfp32Int64=^Pmsfp32Int64;
     Pmsfp32Int64=^Tmsfp32Int64;
     Tmsfp32Int64=Int64;

     PPmsfp32UInt64=^Pmsfp32UInt64;
     Pmsfp32UInt64=^Tmsfp32UInt64;
     Tmsfp32UInt64=UInt64;

     PPmsfp32Void=^Pmsfp32Void;
     Pmsfp32Void=Pointer;

     PPmsfp32Float=^Pmsfp32Float;
     Pmsfp32Float=^Tmsfp32Float;
     Tmsfp32Float=Single;

function MiniSoftFP32Add(const aLeft,aRight:Tmsfp32Float):Tmsfp32Float;
function MiniSoftFP32Sub(const aLeft,aRight:Tmsfp32Float):Tmsfp32Float;
function MiniSoftFP32Mul(const aLeft,aRight:Tmsfp32Float):Tmsfp32Float;
function MiniSoftFP32Div(const aLeft,aRight:Tmsfp32Float):Tmsfp32Float;

implementation

{$ifndef fpc}
function SARLongint(Value,Shift:Tmsfp32Int32):Tmsfp32Int32;
{$ifdef cpu386}
{$ifdef fpc} assembler; inline;
asm
 mov ecx,edx
 sar eax,cl
end ['eax','edx','ecx'];
{$else} assembler; register;
asm
 mov ecx,edx
 sar eax,cl
end;
{$endif}
{$else}
{$ifdef cpuarm} assembler; inline;
asm
 mov r0,r0,asr R1
end ['r0','R1'];
{$else}{$ifdef fpc}inline;{$endif}
begin
{$ifdef HasSAR}
 result:=SARLongint(Value,Shift);
{$else}
 Shift:=Shift and 31;
 result:=(Tmsfp32UInt32(Value) shr Shift) or (Tmsfp32UInt32(Tmsfp32Int32(Tmsfp32UInt32(0-Tmsfp32UInt32(Tmsfp32UInt32(Value) shr 31)) and Tmsfp32UInt32(0-Tmsfp32UInt32(ord(Shift<>0))))) shl (32-Shift));
{$endif}
end;
{$endif}
{$endif}
{$endif}

function MiniSoftFP32Add(const aLeft,aRight:Tmsfp32Float):Tmsfp32Float;
var Left,Right,
    SignLeft,SignRight,
    AbsoluteLeft,AbsoluteRight,
    BaseLeft,BaseRight,
    ExponentLeft,ExponentRight,
    MantissaLeft,MantissaRight,
    Exponent,
    Mantissa,
    Shift,SignMagnitude,NormalizedMantissa,
    Value:Tmsfp32Int32;
begin

 Left:=Tmsfp32Int32(pointer(@aLeft)^);
 Right:=Tmsfp32Int32(pointer(@aRight)^);

 SignLeft:=Left shr 31;
 SignRight:=Right shr 31;

 ExponentLeft:=(Left shr 23) and $ff;
 ExponentRight:=(Right shr 23) and $ff;

 MantissaLeft:=((Left and $7fffff) shl 1) or $1000000;
 MantissaRight:=((Right and $7fffff) shl 1) or $1000000;

 if SignLeft<>0 then begin
  AbsoluteLeft:=-MantissaLeft;
 end else begin
  AbsoluteLeft:=MantissaLeft;
 end;

 if SignRight<>0 then begin
  AbsoluteRight:=-MantissaRight;
 end else begin
  AbsoluteRight:=MantissaRight;
 end;

 if ExponentRight>ExponentLeft then begin
  Shift:=ExponentRight-ExponentLeft;
  if Shift>31 then begin
   Shift:=31;
  end;
  Exponent:=ExponentRight;
  BaseLeft:=SARLongint(AbsoluteLeft,Shift);
  BaseRight:=AbsoluteRight;
 end else begin
  Shift:=ExponentLeft-ExponentRight;
  if Shift>31 then begin
   Shift:=31;
  end;
  Exponent:=ExponentLeft;
  BaseLeft:=AbsoluteLeft;
  BaseRight:=SARLongint(AbsoluteRight,Shift);
 end;
 
 Mantissa:=((SignLeft shl 26) or (SignLeft shl 25) or (BaseLeft and $1ffffff))+
           ((SignRight shl 26) or (SignRight shl 25) or (BaseRight and $1ffffff));

 if (Mantissa and (1 shl 26))<>0 then begin
  SignMagnitude:=(1-Mantissa) and $7ffffff;
 end else begin
  SignMagnitude:=(Mantissa+1) and $7ffffff;
 end;

 inc(Exponent);
 NormalizedMantissa:=SignMagnitude shr 1;

 if (SignMagnitude and $3fffffc)<>0 then begin
  while (NormalizedMantissa and (1 shl 24))=0 do begin
   NormalizedMantissa:=NormalizedMantissa shl 1;
   dec(Exponent);
  end;
 end else begin
  NormalizedMantissa:=NormalizedMantissa shl 24;
  dec(Exponent,24);
 end;

 if (Left and $7fffffff)=0 then begin
  if (Right and $7fffffff)=0 then begin
   Value:=0;
  end else begin
   Value:=Right;
  end;
 end else if (Right and $7fffffff)=0 then begin
  Value:=Left;
 end else if ((NormalizedMantissa and $1ffffff)=0) or ((Exponent and $100)<>0) then begin
  Value:=0;
 end else begin
  Value:=((Mantissa and $4000000) shl 5) or (Exponent shl 23) or ((NormalizedMantissa shr 1) and $7fffff);
 end;

 result:=Tmsfp32Float(pointer(@Value)^);

end;

function MiniSoftFP32Sub(const aLeft,aRight:Tmsfp32Float):Tmsfp32Float;
var Right:Tmsfp32UInt32;
begin
 Right:=Tmsfp32UInt32(pointer(@aRight)^) xor $80000000;
 result:=MiniSoftFP32Add(aLeft,Tmsfp32Float(pointer(@Right)^));
end;

function MiniSoftFP32Mul(const aLeft,aRight:Tmsfp32Float):Tmsfp32Float;
const Bit47:Tmsfp32Int64=Tmsfp32Int64(1) shl 47;
var Left,Right,
    Sign,
    ExponentLeft,ExponentRight,
    MantissaLeft,MantissaRight,
    Exponent,Mantissa,
    Value:Tmsfp32Int32;
    Mantissa64:Tmsfp32Int64;
begin

 Left:=Tmsfp32Int32(pointer(@aLeft)^);
 Right:=Tmsfp32Int32(pointer(@aRight)^);

 Sign:=(Left xor Right) and $80000000;

 ExponentLeft:=(Left shr 23) and $ff;
 ExponentRight:=(Right shr 23) and $ff;

 MantissaLeft:=(Left and $7fffff) or $800000;
 MantissaRight:=(Right and $7fffff) or $800000;
 
 Mantissa64:=Tmsfp32Int64(MantissaLeft)*MantissaRight;

 Exponent:=(ExponentLeft+ExponentRight)-127;
 if (Mantissa64 and Bit47)<>0 then begin
  inc(Exponent);
  Mantissa:=(Tmsfp32Int32(Mantissa64 shr 23)+1) and $ffffff;
 end else begin
  Mantissa:=(Tmsfp32Int32(Mantissa64 shr 22)+1) and $ffffff;
 end;

 if (ExponentLeft=0) or (ExponentRight=0) then begin
  Value:=0;
 end else if (Exponent and $100)=0 then begin
  Value:=Sign or ((Exponent and $ff) shl 23) or (Mantissa shr 1);
 end else if (Exponent and $80)=0 then begin
  Value:=Sign or ($ff shl 23) or (Mantissa shr 1);
 end else begin
  Value:=0;
 end;

 result:=Tmsfp32Float(pointer(@Value)^);

end;

function MiniSoftFP32Div(const aLeft,aRight:Tmsfp32Float):Tmsfp32Float;
const Bit25:Tmsfp32Int64=Tmsfp32Int64(1) shl 25;
var Left,Right,
    Sign,
    ExponentLeft,ExponentRight,
    MantissaLeft,MantissaRight,
    Mantissa,Exponent,
    Value:Tmsfp32Int32;
begin

 Left:=Tmsfp32Int32(pointer(@aLeft)^);
 Right:=Tmsfp32Int32(pointer(@aRight)^);

 Sign:=(Left xor Right) and $80000000;

 ExponentLeft:=(Left shr 23) and $ff;
 ExponentRight:=(Right shr 23) and $ff;

 MantissaLeft:=(Left and $7fffff) or $800000;
 MantissaRight:=(Right and $7fffff) or $800000;

 Mantissa:=(MantissaLeft*Bit25) div MantissaRight;

 Exponent:=(ExponentLeft-ExponentRight)+126;

 if (Mantissa and Bit25)<>0 then begin
  inc(Exponent);
  Mantissa:=(Mantissa shr 1) and $ffffff;
 end else begin
  Mantissa:=Mantissa and $ffffff;
 end;

 if ExponentLeft=0 then begin
  Value:=0;
 end else if ExponentRight=0 then begin
  Value:=Sign or ($ff shl 23);
 end else if (Exponent and $100)=0 then begin
  Value:=Sign or ((Exponent and $ff) shl 23) or ((Mantissa+1) shr 1);
 end else if (Exponent and $80)=0 then begin
  Value:=Sign or ($ff shl 23) or (Mantissa shr 1);
 end else begin
  Value:=0;
 end;

 result:=Tmsfp32Float(pointer(@Value)^);

end;

end.