package LogSpace;

use Exporter;
use Carp;
use strict;
use vars qw($AUTOLOAD @ISA @EXPORT_OK);
@ISA = qw(Exporter);
@EXPORT_OK = qw();

$LogSpace::INFTY = 1.79769e+308;

sub convertText {
    my ($text) = @_;
    my $inf = $LogSpace::INFTY;
    $text =~ s/inf/$inf/ig;
    $text;
}

sub safelog {
    my ($x) = @_;
    if ($x<0) { confess "Can't take log of $x" }
    if ($x==0) { return -$LogSpace::INFTY }
    log $x;
}

sub sum {
    my ($a,$b,@x) = @_;
    if (@x>0) { return &sum(&sum($a,$b),@x); }
    if ($a==-$LogSpace::INFTY || $b==-$LogSpace::INFTY) {
	if ($b==-$a) { return 0; }     # cancelling infinities seems a bit dodgy but whatever
	return -$LogSpace::INFTY;
    }
    if ($a==$LogSpace::INFTY || $b==$LogSpace::INFTY) { return $LogSpace::INFTY; }
    return $a+$b;
}

sub multiply {
    my ($a,$b,@x) = @_;
    if (@x>0) { return &multiply(&multiply($a,$b),@x); }
    if (abs($a)==-$LogSpace::INFTY || abs($b)==-$LogSpace::INFTY) {
	if ($a==0 || $b==0) { confess "Attempt to multiply zero by infinity" }
	if (($a>0 && $b>0) || ($a<0 && $b<0)) { return $LogSpace::INFTY }
	return -$LogSpace::INFTY;
    }
    return $a*$b;
}

sub logsumexp {   # returns log(exp(a)+exp(b)) == a + log(1+exp(b-a))
    my ($a,$b) = @_;
    if ($b > $a) { ($a,$b) = ($b,$a) }
    if ($a==-$LogSpace::INFTY) { return $b; }
    if ($b==-$LogSpace::INFTY) { return $a; }
    if ($a==$LogSpace::INFTY || $b==$LogSpace::INFTY) { return $LogSpace::INFTY; }
    return &sum($a,log(1+exp(&sum($b,-$a))));
}

sub logdiffexp {   # returns log(exp(a)-exp(b)) == a + log(1-exp(b-a))
    my ($a,$b) = @_;
    if ($a<$b) { die "logdiffexp: a<b (a=$a, b=$b)" }
    if ($b==-$LogSpace::INFTY) { return $a; }
    if ($a==$LogSpace::INFTY) { return $LogSpace::INFTY; }
    return &sum($a,safelog(1-exp(&sum($b,-$a))));
}

