六面体の頂点をたどる（10）

前回の分解を、もうちょっとちゃんと考えて、改良してみる。

a > 0として、

$a+b\sqrt{c} = (a_1+b_1\sqrt{c})(a_2+b_2\sqrt{c})$

と分解できるとすると、

$a = a_1a_2 + cb_1b_2$
$b = a_1b_2 + a_2b_1$

ここで、

$a_3 \equiv a_1a_2$
$b_3 \equiv b_1b_2$

と置くと、a₃,b₃は正だから、両方とも1〜a-1。
b₃をその範囲で決めて、b₁とb₂も決まったとすると、

$ab_1 = a_1a_2b_1 + cb_1b_3$
$a_1b = a_1^2b_2 + a_1a_2b_1$
$b_2a_1^2 - ba_1 + b_1a_3$

a₁の2次方程式と見て、判別式は、

$D = b^2 - 4a_3b_3 \ge 0$

a = a₃ + cb₃より、

$4cb_3^2 - 4ab_3 + b^2 \ge 0$

判別式より、a² ≧ cb² ならこれは常に成り立つ。
そうでなければ、

$b \le \frac{a - \sqrt{a^2-cb^2}}, \frac{a + \sqrt{a^2-cb^2}} \le b$



sub alg_resolve {

    my @result;

    my ($a, $b, $c) = @{$_[0]};

    if($a < 0) {

        $a = -$a;

        $b = -$b;

    }

    

    if($b * $b < 4 * ($a - $c)) {

        return \@result;

    }

    else {

        my $det = $a * $a - $c * $b * $b;

        if($det >= 0) {

            my $s = sqrt($det);

            my $alpha = ($a - $s) / (2 * $c);

            my $beta  = ($a + $s + 2 * $c - 1) / (2 * $c);

            for(my $b3 = 1; $b3 <= $alpha; $b3++) {

                alg_resolve_core($a, $b, $b3, $c, \@result);

            }

            for(my $b3 = $beta; $b3 <= ($a - 1) / $c; $b3++) {

                alg_resolve_core($a, $b, $b3, $c, \@result);

            }

        }

        else {

            for(my $b3 = 1; $b3 <= ($a - 1) / $c; $b3++) {

                alg_resolve_core($a, $b, $b3, $c, \@result);

            }

        }

    }

    

    return \@result;

}
sub alg_resolve_core {

    my ($a, $b, $b3, $c, $result) = @_;

    my $a3 = $a - $b3 * $c;

    my $det = $b * $b - 4 * $b3 * $a3;

    next if $det < 0;

    my $r = sqrt($det);

    next if $r * $r != $det;

    

    for(my $b1 = 1; $b1 * $b1 <= $b3; $b1++) {

        if($b3 % $b1 == 0) {

            my $b2 = $b3 / $b1;

            my ($a1, $a2);

            if(($b - $r) % (2 * $b2) == 0) {

                $a1 = ($b - $r) / (2 * $b2);

                if($a3 % $a1 == 0) {

                    $a2 = $a3 / $a1;

                    alg_push($result, [ $a1, $b1, $c ]);

                    alg_push($result, [ $a2, $b2, $c ]);

                }

            }

            

            if(($b + $r) % (2 * $b2) == 0) {

                $a1 = ($b + $r) / (2 * $b2);

                if($a3 % $a1 == 0) {

                    $a2 = $a3 / $a1;

                    alg_push($result, [ $a1, $b1, $c ]);

                    alg_push($result, [ $a2, $b2, $c ]);

                }

            }

        }

    }

}
sub alg_push {

    my ($ary, $a) = @_;

    for my $e(@{$ary}) {

        if($e->[0] == $a->[0] && $e->[1] == $a->[1]) {

            return;

        }

    }

    push @{$ary}, $a;

}