#!/usr/bin/perl

use 5.038.2;

my $precision = 100_000;   # para tirar para adelante, quizá pocos dígitos 

my $altura = $precision;

# distribución estadística contínua
# con zona central uniforme (plana)
# y residuos triangulares (triángulos rectángulos)
# (c) Jesús Lozano Mosterín, 2026

my ($min, $max);     # 0j0: debe poder haber negativos (testzone)
my ($mulow, $muhigh);     # definen la zona uniforme
# estos cuatro son los parámetros que marcan todo

# no importa mucho si me acuerdo de imprimir sólo 6 decimales 

unless (@ARGV && scalar(@ARGV) == 4){
    say "Usage: $0 <p1> <p2> <p3> <p4>";
    say "";      
    say "Parámetros extremos: p1 p4";    
    say "           normales: p2 p3";
    say "";
    say "Sólo (y nada menos) que 4 números\n";
    say "Se devuelven probabilidades para un x";
    say "interactivamente. p1 <= x <= p4";
    say "";
    say "Pulse retorno o no número para salir";
    say "";
    exit -1;
}
    
my ($p1, $p2, $p3, $p4) = sort { $a <=> $b } @ARGV;

for my $i ($p1, $p2, $p3, $p4){
    if (0+$i ne $i){
	die "Error: $i not numeric";
    }
}

$min = $p1;
$mulow = $p2;
$muhigh = $p3;
$max = $p4; 

###########################################################################


die "Error: parámetros incorrectos" unless 
  ($min < $mulow &&
      $mulow < $muhigh && 
      $muhigh < $max);

my $eps = 2;

$eps -= $eps * 0.1 while ( 1+$eps != 1 );
say "\nEl epsilon considerado es $eps\n";

my ($h1, $h2, $h3);
$h1 = $h2 = $h3 = $altura;
my $a1 = area1($min, $mulow);
my $a2 = area2($mulow, $muhigh);
my $a3 = area3($muhigh, $max);
my $total = $a1 + $a2 + $a3;                   # este $total es el total original, para $precision dada

my ($A1, $A2, $A3) = ( $a1 / $total, $a2 / $total, $a3 / $total );

die "Error: round error" unless ($A1 + $A2 + $A3 <= 1 + $eps && 
$A1 + $A2 + $A3 >= 1 - $eps);

my $H = $A2 / ($muhigh - $mulow);       # altura máxima efectiva -> aquí ya $H es ajustado al 100% estadístico, 1

my ($B1, $B2, $B3) = (0, $A1, $A1+$A2);               # las $Bn son acumuladas, no por zonas
# no sigue siendo $total, si usa $H y usa $A's que son base 1

my $x;
while(1) {
    
    ($B1, $B2, $B3) = (0, $A1, $A1+$A2);   # a modo de reinicio, porque hay cálculo nuevo 
    
    print "Gimme an x : ";
    chomp ($x = <STDIN>);
    goto end if ($x =~ /[a-zA-Z]/ || length($x) < 1);

    if ($x < $min){
	report (0, 1);
    }elsif ($x < $mulow){
	
	$h1 = ($x-$min) * $H / ($mulow-$min);    # ppio. de proporcionalidad de triángulos sin()/cos()
	$B1 += area1( $min, $x );                # induce a confusión usar B1 en vez de b1, porque ahora es sobre 1, por $h1     
	report ($B1, 1- $B1);                    # aquí es relevante la altura $h1. Se parte de 0 (menor que $min, Probab. = 0, POR CONVENIENCIA)
	
    }elsif ($x < $muhigh){

	$h2 = $H;
	$B2 += ($x - $mulow) * $h2;              # acumulado, $B2 parte con $B1 completo 
	report ( $B2 , 1- $B2 );
	
    }elsif ($x <= $max){
	
#	$h3 = $H;
	my $b3 = $H * ($max - $muhigh) / 2;
	$h3 = $H * ($max - $x) / ($max - $muhigh);    # la proporción preciosa de triángulos, sin trigonometría 
	$B3 += $b3 - $h3 * ($max-$x) /2;   # tramo de $B3, entero, que 'entra' menos tramo que 'no entra', de la izda.
        # idem. $B3 parte con $B1 + $B2 asegurado
	report ($B3 , 1- $B3);                  # siempre el mismo truco de complementariedad
	
    }else{
	report (1,0);  	
    }
}

end:
  say "\nQuit...\n";

exit 3;


sub area2 {
    # rectángulo
    my ($x0, $x1) = @_ ;
    return ($x1 - $x0) * $h2;
}

sub area1 {
    # $h1 <= $altura
    my ($x0, $x1) = @_;
    return ($x1 - $x0) * $h1 / 2;     # triángulo fácil rectángulo
}

sub area3 {
    # $h3 <= $altura, pero hay 2 fases: fase 1 -> cáculo de áreas y fase 2 -> cálculo de probabilidades
    my ($x0, $x1) = @_;
    return ($x1 - $x0) * $h3 / 2;
}

sub report {
    my ($i, $j) = @_;
    say "\n";
    say "P[ z < $x ] = $i";
    say "";
    say "P[ $x < z ] = $j";
    say "\n";
}



__END__

Para quitar empleo a las salas de control, mayormente.

Es broma. ;-) Se trata de devolver la probabilidad 
  
[0 , 1]
  
para un $x cualquiera. Hay varias:
  
p[ z < x ] = p[ min < x ];

p [ min < x < max ];

p [ z > x ] = p[ max > x ];

La pista habitual es que 
  
p[ min <= z <= max ] = 1;

p[ x < min ] = p[ x > max ] = 0;