package FixedPoint;
# simplistic module for adding and halving stringified
# new, add, halve; no negatives
# new() does NOT properly handle a floating point argument that default-stringifies to scientific notation

use warnings;
use strict;


sub new($;$) {
    my $class = shift;
    my $val = '' . (shift || '0.0');
    my $self = bless \$val, $class;
    $$self .= '.0' if $self->dotpos() >= length($$self);
    return $self;
}

sub copy($) {
    my $orig = shift;
    my $val = '' . $$orig;
    my $self = bless \$val, ref($orig);
    $$self .= '.0' if $self->dotpos() >= length($$self);
    return $self;
}

sub dotpos($) {
    my $val = ${$_[0]};
    my $dot = index($val, '.');
    $dot = index($val.='.0','.') if $dot<0;
    return $dot;
}

sub ilen($) {   # the length of the integer portion
    my $v = ${$_[0]}; # just the value; don't want side effects
    $v =~ s/^0+(?!\.)|(?<!\.)0+$//g;      # strip leading zeroes (except 0.) and trailing zeroes (except .0);
    return dotpos(\$v) || 1;        # if dotpos==0, there's an implied 0 before the dot
}

sub flen($) {   # the length of the fractional portion
    my $v = ${$_[0]}; # just the value; do not affect $self
    $v =~ s/^0+(?!\.)|(?<!\.)0+$//g;      # strip leading zeroes (except 0.) and trailing zeroes (except .0);
    my $dot = dotpos(\$v);
    return 0 if $dot>=length($v);   # no fractional length if the dot is imaginary
    return length($v)-$dot-1;
}

sub lengths($) {
    my $self = shift;
    sprintf '<<%d,%d>>', $self->ilen(), $self->flen();
}

sub halve($) {      # in-place division
    my $self = shift;
    my $rem = 0;
    my $pos = 0;
    my $dot = $self->dotpos;
    if($dot<0) {
        $$self .= '.0';
        $dot = $self->dotpos;
    }
    if(0==$dot) {
        substr $$self, 0, 0, '0';   # precede ".xxx" with "0" => "0.xxx"
        ++$dot;
    }
    while ( $pos < length($$self) ) {
        next if $pos eq $dot;
        my $d = substr($$self, $pos, 1);
        $d += 10*$rem;
        my $h = int($d / 2);
        $rem = $d % 2;
        substr $$self, $pos, 1, $h;
        if(0){
            local $, = "\t";
            local $\ = $/;
            print $pos, $d, $h, $rem,,$$self;
        }
        $$self .= '0' if $rem && ($pos+1 == length($$self));
    } continue {
        ++$pos;
    }
    $$self =~ s/^0+(?!\.)|(?<!\.)0+$//g;      # strip leading zeroes (except 0.) and trailing zeroes (except .0);
}

sub add($) {    # in-place addition
    my($self, $addend) = @_;
    my $sdot = $self->dotpos;
    my $adot = $addend->dotpos;
    my $delt = 0;

    if($sdot >= length($$self))   { $$self   .= ".0"; $sdot = $self->ilen; }
    if($adot >= length($$addend)) { $$addend .= ".0"; $adot = $addend->ilen; }

    # line up fixed points
    if($sdot < $adot) {
        $delt = $adot - $sdot;
        substr $$self, 0, 0, '0'x$delt;
        $sdot = $self->dotpos;
    } elsif ($adot < $sdot ) {
        $delt = $sdot - $adot;
        substr $$addend, 0, 0, '0'x$delt;
        $adot = $addend->dotpos;
    }

    # line up end of fractions
    my $sf = $self->flen;
    my $af = $addend->flen;
    if($sf<$af) {    # self fract is shorter
        $$self .= '0' x ($af-$sf);
        $sf = $self->flen;
    } elsif ($af<$sf) {
        $$addend .= '0' x ($sf-$af);
        $af = $addend->flen;
    }

    #print join "\n", 'add',
    #    'self  :'.$$self.':'.$sdot.':'.$sf,
    #    'addend:'.$$addend.':'.$adot.':'.$af;

    die "add: could not get the lengths the same:\n\t>>$$self<<\n\t>>$$addend<<\n"  if length($$self)!=length($$addend);
    my $l = length($$self);
    my $c = 0;
    for my $p (1 .. $l) {
        next if $sdot == ($l-$p);
        my $s = substr $$self,   -$p, 1;
        my $a = substr $$addend, -$p, 1;
        my $total = $s + $a + $c;
        substr($$self, -$p, 1) = $total % 10;
        $c = int($total/10);
    }
    substr $$self, 0, 0, $c if $c;

    #print join "\n", '', 'answer:'.$$self;
    return $self;
}

1;

package main;
# ULP: 2**-52
my $ulp = FixedPoint->new(1);
$ulp->halve for 1..52;
# make 1+2**-52
my $x = FixedPoint->new(1)->add($ulp);
print "\n";
printf "%s\n", '-'x(52+2);
printf "%-8s = %s\n%-8s = %s\n", 'ulp(1)', $$ulp, '1+ulp(1)', $$x;
printf "%s\n", '-'x(52+2);
# bring it down to the smallest normal double plus one ULP
#   = (1+2**-52) * 2**-1022 = 2**-1022 + 2**-1074
for(1..1022) {
    $x->halve;
    printf "%s\n", $$x          if 0;
}

# details on the FixedPoint:
print "\n", '-'x40, "\n";
printf "(1+ulp(1))*2**-1022 = \n%s\n", $$x;
print "\n";
printf "FixedPoint: %-16d = Digits before the fixed decimal point\n", $x->ilen();
printf "FixedPoint: %-16d = Total digits the fixed decimal point\n", $x->flen();
my $p = 0;
$$x =~ m/[1-9]/g ;
$p = pos($$x);
$p -= $x->ilen()+1;
printf "FixedPoint: %-16d = Location of first significant digit (nth digit after decimal point)\n", $p;
printf "FixedPoint: %-16d = Number of significant digits\n", $x->flen() - $p + 1;

use POSIX qw(floor ceil log10 log2);
use Data::IEEE754::Tools qw(:ulp :constants :floatingpoint);
my $tiny = nextUp( POS_NORM_SMALLEST );
print "\n";
printf "Double:    %30.16e\n", $tiny;
printf "HexFloat:  %30s\n", to_hex_floatingpoint($tiny);
printf "DecFloat:  %30s\n", to_dec_floatingpoint($tiny);
print "\n";
printf "Math:       %-16d = should start at ceil(-log10(x))-th digit after decimal point\n", my $d = ceil(-log10($tiny));
printf "Math:       %-16d = total number of digits after decimal should be the lowest power of two = ceil(-log2(ulp(x)))\n", my $e = ceil(-log2(ulp($tiny)));
printf "Math:       %-16d = Number of significant digits", $e - $d + 1;

##</code><code>##

__RESULTS__
------------------------------------------------------
ulp(1)   = 0.0000000000000002220446049250313080847263336181640625
1+ulp(1) = 1.0000000000000002220446049250313080847263336181640625
------------------------------------------------------

----------------------------------------
(1+ulp(1))*2**-1022 =
0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000022250738585072018771558785585789482407880088486837041956131300312119688603996006965297904292212628858639037013670281908017171296072711910355127227413175152199055740043138804567803233377539881639177387328959246074229270113078053813397081653361296447449529789521218979090783852583365901851789618799885150427514782636076021680436220311292700454832073964845713103912225963935608322440623896907276890186717054549275173986589324810401738228328251245795065655738191038008646911615828719989708647293221449796971546706720399791990809160347625980385995424739847678861180095072511543762389603716215171729816011544604359531284325406441938645324905389137795680915804792405099227413854274942620542640408839836919187418172987793340279242767544565229087538682506419718265533447265625

FixedPoint: 1                = Digits before the fixed decimal point
FixedPoint: 1074             = Total digits the fixed decimal point
FixedPoint: 308              = Location of first significant digit (nth digit after decimal point)
FixedPoint: 767              = Number of significant digits

Double:           2.2250738585072019e-308
HexFloat:        +0x1.0000000000001p-1022
DecFloat:     +0d1.0000000000000002p-1022

Math:       308              = should start at ceil(-log10(x))-th digit after decimal point
Math:       1074             = total number of digits after decimal should be the lowest power of two = ceil(-log2(ulp(x)))
Math:       767              = Number of significant digits