Interesting article on new number type that promises better precision
than float
lazarus / 64-FPC development

----------------------------------------------------------------------
edvard · 2016-05-26 01:23

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:

QuoteBetter 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...

----------------------------------------------------------------------
skalogryz · 2016-05-26 02:51

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: 10..$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

----------------------------------------------------------------------
Eugene Loza · 2016-05-26 09:25

Quote from: edvard on May 26, 2016, 03:23:21 amMaybe we don't need it,
maybe we doWell... 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...

----------------------------------------------------------------------
MathMan · 2016-05-31 14:26

Quote from: Eugene Loza on May 26, 2016, 11:25:49 amQuote from: edvard
on May 26, 2016, 03:23:21 amMaybe we don't need it, maybe we doWell...
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