Recent

Author Topic: Interesting article on new number type that promises better precision than float  (Read 3012 times)

edvard

  • Full Member
  • ***
  • Posts: 165
I read through this article thinking "Hey, we could use a Pascal port of this", and looking at the Python port, it doesn't look simple:

Quote
Better Than Floating - New Number Format Avoids Imprecision
...
A new format called unum saves bits in all senses of the word and promises better precision.
http://www.i-programmer.info/news/202-number-crunching/9707-better-than-floating-new-number-format-avoids-imprecision.html

Maybe we don't need it, maybe we do, and the specifics are over my head. It seems to be laid out in the article pretty clearly though, and gee whiz, if it can be done in Python...  8-)
All children left unattended will be given a mocha and a puppy.

Debian 'Testing' 64-bit | Xfce 4.12 | FreePascal 3.0 + Lazarus 1.6.2, GTK2+ and Qt

skalogryz

  • Global Moderator
  • Hero Member
  • *****
  • Posts: 2345
    • havefunsoft.com
it seems like unum has been around for awhile.

Here's the code to start with :)
Code: [Select]
unit unumtype;

{$mode objfpc}{$H+}

interface

type
  unumbits = bitpacked record
  {$IFDEF ENDIAN_LITTLE}
    fracsize: 0..$0F;  // 4 bit
    expsize : 0..$07;  // 3 bit
    ubit    : 0..$01;  // 1 bit
    frac    : 0..$0FFF;//12 bit
    exp     : 0..$FF;  // 8 bit
    sign    : 0..$01;  // 1 bit
  {$ENDIF}
  end;

  unum = record
  case byte of
    0: (b: unumbits);
  end;

procedure i2u(i: integer; out u: unumbits); overload;
procedure i2u(i: integer; out u: unum); inline; overload;
procedure u2f(const u: unum; out s: single); inline; overload;
function u2f(const u: unum): single; inline; overload;
function isuexact(const u: unum): Boolean; inline;
procedure f2u(f: single; out u: unum);

implementation

const
  exp_ofs   = 127;
  frac_max  = 12;

type
  ieeesingle = bitpacked record
    {$IFDEF ENDIAN_LITTLE}
    frac : 0..$7FFFFF;
    exp  : 0..$FF;
    sign : 0..$1;
    {$ENDIF}
  end;

procedure _i2u(lu: LongWord; out u: unumbits; asign: byte);
var
  l  : Longword;
  ex : byte;
  f  : longword;
  ubit: byte;
  m  : LongWord;
  ln : integer;
begin
  if lu=0 then begin
    ex:=0;
    l:=0;
    f:=0;
    u.exp:=0;
  end else begin
    l:=1 shl 31;
    ex:=31;
    while (l and lu) =0 do begin
      l:=l shr 1;
      dec(ex);
    end;
    f:=lu and (l-1);
    u.exp:=ex+exp_ofs;
  end;
  u.sign:=asign;

  ubit:=0;
  if ex<frac_max then begin
    f:=f shl (frac_max-ex);
  end else if ex>frac_max then begin
    ln:=ex-frac_max;
    m:=(1 shl ln) - 1;
    if f and m = 0
      then ubit:=0
      else ubit:=1;
    f:=f shr ln;
  end;
  u.frac:=f;
  u.ubit:=ubit;
  u.expsize:=7; // hard-coding... whatever...
  u.fracsize:=12;
end;

procedure i2u(i: integer; out u: unumbits);
var
  s : byte;
begin
  if i<0 then begin
    s:=1;
    i:=abs(i);
  end else
    s:=0;
  _i2u(longWord(i), u, s);
end;

procedure i2u(i: integer; out u: unum);
begin
  i2u(i, u.b);
end;

procedure u2f(const u: unum; out s: single); inline;
var
  ie : ieeesingle absolute s;
begin
  ie.exp:=u.b.exp;
  ie.frac:=(u.b.frac shl 11);
  ie.sign:=u.b.sign;
end;

function u2f(const u: unum): single; inline;
begin
  u2f(u,Result);
end;

function isuexact(const u: unum): Boolean;
begin
  Result:=u.b.ubit=0;
end;

procedure f2u(f: single; out u: unum);
var
  ie : ieeesingle absolute f;
begin
  u.b.sign:=ie.sign;
  u.b.exp:=ie.exp;
  u.b.frac:=ie.frac shr 11;
  if ie.frac and ((1 shl 11) - 1)=0
    then u.b.ubit:=0
    else u.b.ubit:=1;
  u.b.expsize:=7;
  u.b.fracsize:=12;
end;

end.

sample:
Code: [Select]
uses unumtype;
var
  u  : unum;
  f  : single;
begin
  i2u(2335, u);
  f:=u2f(u);
  writeln(f:0:2,' ',u.b.ubit);

  i2u(23554, u);
  f:=u2f(u);
  writeln(f:0:2,' ',u.b.ubit);

  i2u(65536, u);
  f:=u2f(u);
  writeln(f:0:2,' ',u.b.ubit);
end.
no arithmetic operations are implemented, obviously
« Last Edit: May 26, 2016, 06:13:22 am by skalogryz »
Patron Cocoa Widgetset development https://www.patreon.com/skalogryz

Eugene Loza

  • Hero Member
  • *****
  • Posts: 563
    • My "almost daily" development blog
Maybe we don't need it, maybe we do
Well... it doesn't seem very promising to me for complex numerical calculations:
1. At the moment I understand that normal numerical co-processors don't seem to support it. So complex calculations will be slower.
2. Non-standard length of 29 bits... looks weird to me. I am really looking forward to 32 or 64 bit floats with full co-processor acceleration.
3. It fails just the same way as double does. E.g. I have an infinite integral in my dissertation consisting of 2 almost equal "halves" with opposite signs and it's the difference I need which occasionally might be 10 to 30 digits below least significant digit of extended type. unum doesn't seem to be able warn me when I'm getting an absolutely wrong result due to accuracy lack. Variable length float would have been very nice, but unum is fixed :)
So... I see it simply as something a bit more accurate than extended type but significantly slower and memory-inefficient.
Maybe I'm missing something...
Lazarus 1.9 + FPC 3.1.1 Debian Jessie 64 bit.

My Free and Open Source games in Lazarus/FreePascal/CastleGameEngine:
https://decoherence.itch.io/
(and some ancient games in Turbo Pascal too)
Sources are here: https://github.com/eugeneloza?tab=repositories

MathMan

  • Full Member
  • ***
  • Posts: 166
Maybe we don't need it, maybe we do
Well... it doesn't seem very promising to me for complex numerical calculations:
...
So... I see it simply as something a bit more accurate than extended type but significantly slower and memory-inefficient.
Maybe I'm missing something...

Looks the same to me. Regarding your other point maybe you should take a look at MPFR - you'd need a Pascal wrapper though as this is a C library. There are also some libs avail for variable length floating point with precision checking - but I am missing the names of the top of my head.

MathMan