#!/usr/bin/perl

use 5.038.2;

my $version = 0.03;

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 para salir";
    say "Pulse sólo 'm' para obtener la mediana";
    say "version $version\n";
    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 $ans;
my $x;
my $r = 0;
while(1) {
    
    print "Gimme an x : " unless defined $ans;
    chomp (my $ans = <STDIN>);

    if ( (0+$ans ne $ans) && $ans ne 'm' ) {

	last;  # salida (menos la m o retorno)
	
    } elsif ($ans eq 'm') {      
	# M de mediana
	my $c = 0;
	$x = 0;
	while ($c++ <= $precision) {
	    $x = $min + ($max - $min) * $c/$precision; 
	    $r= probabs();
	    $ans = $x;
	    last if ($r)
	}                     # REALMENTE CON LA MITAD BASTARIA + -

	$ans = undef;
	say "-+- La mediana es = $x\n";
	   
    }elsif (0+$ans eq $ans) {	
	$x = $ans;
	probabs();
	$ans = undef;
    
    }else{
	
	say "Input error. Try again!";
	$ans = undef;
    }
}

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";
    if ($i >= 0.5-$eps && $i <= 0.5+$eps){
	return 1;
    }else{
	return 0;
    }
}

sub probabs {
    ($B1, $B2, $B3) = (0, $A1, $A1+$A2);         # a modo de reinicio, porque hay cálculo nuevo 
    if ($x < $min){
	$r = 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     
	$r = 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 
	$r = report ( $B2 , 1- $B2 );
    }elsif ($x <= $max){
	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
	$r = report ($B3 , 1- $B3);                   # siempre el mismo truco de complementariedad
    }else{
	$r = report (1,0);  	
    }
    return $r;
}

__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;

v.0.02 -> tratar de buscar la mediana