#!/usr/bin/perl

# Se trata de una variación de la idea de 'dist_std_res2'
# La masa central, que es la que puede ser discutible como uniforme,
# sea simétrica y semicircular entre puntos regulares;
# Los residuos triangulares y asimétricos finitos creo que están
# bien en la formulación anterior

# Los inputs serían los mismos: $min, $regular_low, $regular_high, $max
# y, como siempre, para que esté chupado, densidad($regular_low) =
# densidad($regular_high), siendo la distribución diferente en ambos
# puntos, claro. 

use 5.038.2;
use Math::Trig;      # for functions asin() and acos() and atan(y/x)

my $DEBUG = 1;

my $version = 0.3141592;

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

my $altura = $precision;

my $pi = 4 * atan2(1,1);
# my $radiande90 = ($pi/2);

# distribución estadística contínua
# con zona central semicircular sobre rectángulo
# 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";
    }
}

my ($centro, $radio, $semicirculo, $semicirculoregularizado);

$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";
say "Los parámetros son $p1 $p2 $p3 $p4" if $DEBUG;

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
$semicirculoregularizado = $semicirculo / $total;

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

say "Las áreas correspondientes a las zonas son:";
say "Residuo inferior  = $A1";
say "Zona central      = $A2";
say "Residuo superior  = $A3";
if ($DEBUG){
    say "\nDe la zona central:";
    say "Zona uniforme     = ", $A2 - $semicirculoregularizado;
    say "Zona semicircular = $semicirculoregularizado";     
}
say "";

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

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

# jue 29 ene 2026 11:05:17 CET   Hasta aquí no veo nada

my ($ans, $x, $r, $c);

$c = 0;

while(1) {
    
    print "Gimme an x [return to exit]: " unless defined $ans;
    $ans = <STDIN>;
    chomp $ans;
#    trim($ans);
    
    if ( $ans !~ /[\d|\.|\-]/ && $ans !~ /m/i ) {

	last;  # salida (menos la m o retorno)
	
    } elsif ( $ans !~ /[\d|\.|\-]/ && $ans =~ /m/i ) {              # MEDIANA
	
	while ($r <= 0.5-$eps || $r <= 0.5+$eps) {                 # $r va a saltos
	    $x = $min + ($max - $min) * $c/$precision; 
	    $r = probabs();
	    $c++;
	    $ans = $x;
	}                     # REALMENTE CON LA MITAD de $c aprox. BASTARIA + -

	say "-+- La mediana es = $x\n";
	undef $ans;
	$c = 0;
	
    }elsif ($ans =~ /[\d|\.|\-]/ && $ans !~ /m/i ) {	                              # NUMÉRICO (sólo, y no formato general /\de\d/)
	$x = $ans;
	probabs();
	undef $ans;
	
    }else{ 
	# do nothing -> loop
    }
}

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

exit 3;

sub trim {
    # poor man trim
    my $s = shift;
    if ($s =~ s/^\s*(\$*)\s*$/$1/){
	return $s;
    }
}

sub area2 {
    # rectángulo + semicírculo
    my ($x0, $x1) = @_ ;
    my $rectangulo = ($x1 - $x0) * $h2;
    $centro = ($x1 - $x0) /2 ;
    $radio = $centro - $x0;
    $semicirculo = $pi * $radio ** 2;
    return $rectangulo + $semicirculo; 
}

sub semisphere {
    my $x = shift;
    my $lado = $muhigh - $mulow;
#    my $alturatotal = $H + $centro;
    my $area;
    my $hh;
    if ($x < $mulow){
	return 0;
    }
    if ($x >= $muhigh){
	return $semicirculoregularizado;
    }
    if ($x < $centro){
	# $mulow, $x, $centro    $radio = $centro - $mulow
	# ($centro-$x)**2 + $hh**2 = $radio**2 
	$hh = sqrt($radio**2 - ($centro-$x)**2);
	# atan(y/x) devuelve [-1, 1] en 180 grados [0,1] en 90
	$area = abs( atan( $hh/($centro-$x) ) * $semicirculoregularizado/2 );
    }else{
	$area = $semicirculoregularizado/2;
	# radio ** 2 = (x-centro)**2 + hh**2 
	if ($radio**2 < ($x - $centro)**2){
	    $area = $semicirculoregularizado;
	    return $area;
	}
	$hh = sqrt($radio**2 - ($x - $centro)**2);
	$area += atan( $hh/ ($x - $centro) ) * $semicirculoregularizado / 2;	
    }
    return $area;
}

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 "";
    say "P[ z < $x ] = $i\n";
    say "P[ $x < z ] = $j\n";
    return $i;
}

sub probabs {
    ($B1, $B2, $B3) = (0, 0, 0);         # 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 = $A1 + ($x - $mulow) * $h2 + semisphere($x);              # acumulado, $B2 parte con $B1 completo 
	$r = report ( $B2 , 1- $B2 );
    }elsif ($x <= $max){
	my $b3 = $H;
	$h3 = $H * ($max - $x) / ($max - $muhigh);    # la proporción preciosa de triángulos, sin trigonometría 
	$B3 = 1 - $A1 - $A2 - $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