unum in raku

unum.rakumod

unit module Unum;
=begin LICENSE
Copyright © 2017 John L. Gustafson

Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sub-license, and/or sell copies
of the Software, and to permit persons to whom the Software is furnished to do
so, subject to the following conditions:

This copyright and permission notice shall be included in all copies or
substantial portions of the software.

THE SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES, OR OTHER
LIABILITY, WHETHER IN AN ACTION OR CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
=end LICENSE

constant blue    = "\e[38;2;64;128;255m";
constant brown   = "\e[38;2;153;102;51m";
constant magenta = "\e[38;2;255;0;255m";
constant red     = "\e[38;2;255;32;0m";
constant amber   = "\e[38;2;204;153;51m";
constant reset   = "\e[0m";
constant dim       = "\e[2m";
constant dim-reset = "\e[22m";

multi sqrt(Int $n) {
  (10**($n.chars div 2), { ($_ + $n div $_) div 2 } ... * == *).tail
}

role quire[UInt $] {...}

role posit[UInt $nbits] {
  =for resource
    L<https://posithub.org/docs/Posits4.pdf>

  my $es     := $*posit-exponent-size //
		%*posit<exponent-size> //
		%*ENV<POSIT_EXPONENT_SIZE> //
		2;
  my $npat   := 2**$nbits;
  my $useed  := 2.FatRat**2**$es;
  my $minpos := $useed**(-$nbits + 2);
  my $maxpos := $useed**(+$nbits - 2);
  my $qsize  := 2**(($nbits-2)*2**($es+2) + 5).log2.ceiling;
  my $qextra := $qsize - ($nbits - 2)*2**($es+2);

  sub twoscomp($sign, $p) { ($sign > 0 ?? $p !! $npat - $p) mod $npat }
  sub sign-bit(UInt $p where ^$npat) { +($p ≥ $npat div 2) }
  sub regime-bits(UInt $p where ^$npat) {
    my $q = twoscomp(1 - sign-bit($p), $p);
    my @bits = $q.polymod(2 xx $nbits - 1).reverse;
    my $bit2 = @bits[1];
    my @temp-bits = flat @bits[1..*], 1 - $bit2;
    my $npower = @temp-bits.first(1 - $bit2, :k) - 1;
    @bits[1..($npower+1)];
  }
  sub regime-value(@bits) { @bits.head ?? @bits.elems - 1 !! -@bits; }
  sub exponent-bits(UInt $p where ^$npat) {
    my $q = twoscomp(1 - sign-bit($p), $p);
    my $startbit = regime-bits($q).elems + 2;
    my @bits = $q.polymod(2 xx $nbits - 1).reverse;
    $startbit > @bits.end ?? Empty !!
    @bits[$startbit .. min($startbit + $es - 1, @bits.end)]
  }
  sub fraction-bits(UInt $p where ^$npat) {
    my $q = twoscomp(1 - sign-bit($p), $p);
    my $startbit = regime-bits($q).elems + 2 + $es;
    my @bits = $q.polymod(2 xx $nbits - 1).reverse;
    $startbit > @bits.end ?? Empty !! @bits[$startbit .. *];
  }
  sub p2x(UInt $p where ^$npat) {
    my $s = (-1)**sign-bit($p);
    my $k = regime-value regime-bits $p;
    my @e = exponent-bits $p;
    my @f = fraction-bits $p;
    @e.push: 0 until @e == $es;
    my $e = @e.reduce: 2 * * + *;
    my $f = @f == 0 ?? 1 !! 1.FatRat + @f.reduce(2 * * + *)/2.FatRat**@f;
    given $p {
      when 0          { 0   }
      when $npat div 2 { Inf }
      default         { $s * $useed**$k * 2**$e * $f }
    }
  }
  sub x2p(Real $x) {

    # first, take care of the two exceptional values
    return 0 if $x == 0;
    return $npat div 2 if $x == Inf;

    # working variables
    my UInt ($i, $p);
    my FatRat $y = $x.abs.FatRat;
    my UInt $e = 2**($es - 1);
      
    if $y ≥ 1 { # north-east quadrant
      ($p, $i) = 1, 2;
      # Shift in 1s from the right and scale down.
      ($p, $y, $i) = 2*$p+1, $y/$useed, $i+1 while $y ≥ $useed && $i ≤ $nbits;
      $p *= 2; $i++;
    } else { # south-east quadrant
      ($p, $i) = 0, 1;
      # Shift in 0 s from the right and scale up.
      ($y, $i) = $y*$useed, $i+1 while $y < 1 && $i ≤ $nbits;
      if $i ≥ $nbits {
	$p = 2; $i = $nbits + 1;
      } else { $p = 1; $i++; }
    }
    # Extract exponent bits:
    while $e > 1/2 and $i ≤ $nbits {
      $p *= 2;
      if $y ≥ 2**$e { $y /= 2**$e; $p++ }
      $e div= 2;
      $i++;
    }
    $y--;
    # Fraction bits; substract the hidden bit
    while $y > 0 and $i ≤ $nbits {
      $y *= 2;
      $p = 2*$p + $y.floor;
      $y -= $y.floor;
      $i++
    }
    $p *= 2**($nbits+1-$i);
    $i++;

    # Round to nearest; tie goes to even
    $i = $p +& 1;
    $p = ($p/2).floor;
    $p = $i == 0 ?? $p !!
      $y == 0|1 ?? $p + $p+&1 !!
      $p + 1;

    # Simulate 2's complement
    ($x < 0 ?? $npat - $p !! $p) mod $npat;
  }

  has blob8 $.blob;
  method precision { $nbits }

  multi method UInt { self.blob.read-ubits(0, $nbits) }

  method sign-bit { sign-bit self.UInt }
  method regime-bits { regime-bits self.UInt; }
  method exponent-bits { exponent-bits self.UInt }
  method fraction-bits { fraction-bits self.UInt }

  multi method negate {
    my buf8 $buf .= new;
    $buf.write-ubits:
      0,
      $nbits,
      $npat - (self.blob.read-ubits: 0, $nbits) mod $npat;
    self.new: blob => $buf;
  }

  multi method Real {
    given self.UInt {
      when 0          { 0   }
      when $npat div 2 { Inf }
      default { p2x $_ }
    }
  }
  multi method bits {
    return "zero (posit[$nbits])" if self.UInt == 0;
    return "±∞ (posit[$nbits])"   if self.UInt == $npat div 2;
    my $s = (-1)**self.sign-bit;
    my @r = self.regime-bits;
    my @e = self.exponent-bits;
    my @f = self.fraction-bits;
    [~]
      red,   $s < 0 ?? '-' !! '+',
      amber, @r.join,
      brown, @r ≤ $nbits - 2 ?? 1 - @r.head !! '',
      blue,  @e.join,
      reset,
      dim, @f.join, dim-reset
	;
  }

  multi method gist { "posit[$nbits]:{self.Real}" }

  multi method new(Real $x) {
    my buf8 $buf .= new;
    $buf.write-ubits: 0, $nbits,
      $x == 0 ?? 0 !!
      $x == Inf ?? $npat div 2 !!
      x2p $x;
    return self.new: blob => $buf;
  }

  method sign {
    if    self.Real     == 0 { return 0 }
    elsif self.sign-bit == 1 { return -1 }
    else { return 1 }
  }
  method abs { self.sign < 0 ?? self.negate !! self }
  method sqrt {
    return NaN if self.sign < 0;
    return 0   if self.Real == 0;
    return Inf if self.Real == Inf;
    self.new: sqrt((self.Real*2**(2*$nbits)).Int)/2.FatRat**$nbits;
  }
  method rSqrt {
    return NaN if self.sign < 0;
    return Inf if self.Real == 0;
    return 0   if self.Real == Inf;
    self.new: 1/(sqrt((self.Real*2**(2*$nbits)).Int)/2.FatRat**$nbits);
  }
}

role quire[UInt $n] {

  has blob8 $.blob;
  method precision { 16*$n }

  method Real {
    if self.sign-bit == 1
      and self.blob.read-ubits(1, 16*$n - 1) == 0 {
	return NaN
    } else { return self.blob.read-bits(0, 16*$n)*2.FatRat**(16-8*$n); }
  }
  method sign-bit        { self.blob.read-ubits:         0,         1 }
  method carry-guard     { self.blob.read-ubits:         1,        31 }
  method integral-part   { self.blob.read-ubits:        32, 8*$n - 16 }
  method fractional-part { self.blob.read-ubits: 16 + 8*$n, 8*$n - 16 }

  multi method new(Real $x) {
    my buf8 $buf .= new;
    $buf.write-bits: 0, 16*$n, ($x*2.FatRat**(8*$n - 16)).round;
    self.new: blob => $buf;
  }
  method toPosit { posit[$n].new: self.Real }
  method negate {
    my buf8 $buf .= new;
    $buf.write-bits: 0, 16*$n, -self.blob.read-bits: 0, 16*$n;
    self.new: blob => $buf;
  }
  method abs { self.sign-bit ?? self.negate !! self }
}

sub postfix:<p>(Real $r --> posit) is export {
  posit[
    $*posit-precision //
    %*posit<precision> //
    %*ENV<POSIT_PRECISION> //
    8
  ].new: $r
}

sub postfix:<q>(Real $r --> quire) is export {
  quire[
    $*posit-precision //
    %*posit<precision> //
    %*ENV<POSIT_PRECISION> //
    8
  ].new: $r;
}

# Posit Arithmetic
# ================
multi sub prefix:<->(posit $p --> posit) is export { $p.negate }
multi sub prefix:<+>(posit $p --> posit) is export { $p }

multi sub infix:<+>(posit $p, Real $x) is export { $p.new: $p.Real + $x }
multi sub infix:<+>(Real $x, posit $p) is export { $p.new: $p.Real + $x }
multi sub infix:<*>(Real $x, posit $p) is export { $p.new: $p.Real * $x }
multi sub infix:<*>(posit $p, Real $x) is export { $p.new: $p.Real * $x }
multi sub infix:</>(posit $p, Real $x) is export { $p.new: $p.Real / $x }
multi sub infix:</>(Real $x, posit $p) is export { $p.new: $x / $p.Real }
multi sub infix:<->(posit $p, Real $x) is export { $p.new: $p.Real - $x }

multi sub infix:<+>(posit $p, posit $q where $p.precision == $q.precision) is export { $p.new: $p.Real + $q.Real }
multi sub infix:<->(posit $p, posit $q where $p.precision == $q.precision) is export { $p.new: $p.Real - $q.Real }
multi sub infix:<*>(posit $p, posit $q where $p.precision == $q.precision) is export { $p.new: $p.Real * $q.Real }
multi sub infix:</>(posit $p, posit $q where $p.precision == $q.precision) is export { $p.new: $p.Real / $q.Real }

multi sub infix:<**>(posit $p, 0) is export { $p.Real == 0 ?? NaN !! $p.new: 1 }
multi sub infix:<**>(posit $p, 1) is export { $p }
multi sub infix:<**>(posit $p, 2) is export { $p*$p }
multi sub infix:<**>(posit $p, UInt $n) is export { ($p **($n div 2))**2 * $p**($n mod 2) }

# Quire Arithmetic
# ================
multi sub infix:<*>(Real $x, quire $q) is export { $q.new: $x * $q.Real }
multi sub infix:<*>(quire $q, Real $x) is export { $x*$q }
multi sub infix:</>(quire $q, Real $x) is export { $q*(1/$x) }
multi sub infix:</>(Real $x, quire $q) is export { $q.new: $x / $q.Real }
multi sub infix:<+>(quire $p, quire $q where $p.precision == $q.precision) is export { $p.new: $p.Real + $q.Real }
multi sub infix:<*>(quire $p, quire $q where $p.precision == $q.precision) is export { $p.new: $p.Real * $q.Real }

multi sub infix:<**>(quire $p, 0) is export { $p.Real == 0 ?? NaN !! $p.new: 1 }
multi sub infix:<**>(quire $p, 1) is export { $p }
multi sub infix:<**>(quire $p, 2) is export { $p*$p }
multi sub infix:<**>(quire $p, UInt $n) is export { ($p **($n div 2))**2 * $p**($n mod 2) }