File Coverage

blib/lib/Math/Polynomial/Cyclotomic.pm

Criterion	Covered	Total	%
statement	218	218	100.0
branch	82	82	100.0
condition	48	53	90.5
subroutine	29	29	100.0
pod	8	16	50.0
total	385	398	96.7

# Copyright (c) 2019-2021 by Martin Becker.  This package is free software,
# licensed under The Artistic License 2.0 (GPL compatible).

package Math::Polynomial::Cyclotomic;

use 5.006;
use strict;
use warnings;
use Carp qw(croak);
use Math::BigInt try => 'GMP';
use Math::Polynomial 1.019;
use Math::Prime::Util qw(
    divisors factor factor_exp euler_phi gcd is_square_free is_power
    moebius kronecker vecprod
);
require Exporter;

our @ISA         = qw(Exporter);
our @EXPORT_OK   = qw(
    cyclo_poly cyclo_factors cyclo_plusfactors cyclo_poly_iterate
    cyclo_lucas_cd cyclo_schinzel_cd cyclo_int_factors cyclo_int_plusfactors
);
our %EXPORT_TAGS = (all => \@EXPORT_OK);
our $VERSION     = '0.004';

# some polynomial with coefficient type Math::BigInt
my $default = Math::Polynomial->new(Math::BigInt->new('1'));
$default->string_config({fold_sign => 1, times => '*'});

# ----- private subroutines -----

sub _cyclo_poly {
    my ($poly, $table, $n, $divisors) = @_;
    return $table->{$n} if exists $table->{$n};
    my $m = vecprod(map { $_->[0] } factor_exp($n));    # $m = radix($n)
    my $o = $m & 1;
    $m >>= 1 if !$o;
    my $p;
    if (exists $table->{$m}) {
        $p = $table->{$m};
    }
    else {
        my $one = $poly->coeff_one;
        my @d = $divisors? (grep { not $m % $_ } @{$divisors}): divisors($m);
        $p = $poly->monomial(pop @d)->sub_const($one);
        if (@d) {
            my $b = $d[-1];
            $p /= $poly->monomial($b)->sub_const($one);
            for (my $i = 1; $i < $#d; ++$i) {
                my $r = $d[$i];
                if ($b % $r) {
                    $p /= $table->{$r} || _cyclo_poly($poly, $table, $r, \@d);
                }
            }
        }
        $table->{$m} = $p;
    }
    if (!$o) {
        $p = -$p if $m == 1;
        $m <<= 1;
        $p = $p->mirror;
        $table->{$m} = $p;
    }
    if ($m != $n) {
        $p = $p->inflate($n/$m);
        $table->{$n} = $p;
    }
    return $p;
}

# for a positive integer n, return r,s so that r is square-free and n=rs^2
sub _square_free_square {
    my ($n) = @_;
    my ($r, $s) = (1, 1);
    foreach my $be (factor_exp($n)) {
        my ($b, $e) = @{$be};
        $r *= $b if $e & 1;
        $e >>= 1;
        while ($e) {
            $s *= $b if $e & 1;
            $e >>= 1 and $b *= $b;
        }
    }
    return ($r, $s);
}

# generate polynomials C, D, satisfying the identity of Aurifeuille,
# Le Lasseur and Lucas: for n > 1, square-free,
#   n === 1 (mod 4): Phi_n(x)  = C^2 - n*x*D^2
#   else:            Phi_2n(x) = C^2 - n*x*D^2
sub _lucas_cd {
    my ($poly, $table, $n) = @_;
    my $c1 = 1 == ($n & 3);
    my $n1 = $c1? $n: $n + $n;
    if (exists $table->{"$n1 $n"}) {
        return (@{$table->{"$n1 $n"}})
    }
    my $c2 = !($n & 1);
    my $ee = euler_phi($n1);
    my $e  = $ee >> 1;
    my $one = $poly->coeff_one;
    my $nul = $poly->coeff_zero;

    my $E = $e | 1;
    my @q = ($one);
    my $cos = $one;
    my %MP = ();
    for (my $k = 2; $k <= $E; ++$k) {
        if ($k & 1) {
            push @q, $nul + kronecker($n, $k);
            next;
        }
        if ($c2 && $k & 2) {
            push @q, $nul;
            next;
        }
        $cos = -$cos if !$c1;           # $cos is cos((n-1)*k*pi/4)
        my $gk = gcd($n1, $k);
        my $mp = $MP{$gk} ||= moebius($n1 / $gk) * euler_phi($gk);
        push @q, $cos * $mp;
    }

    my @cs = ($one);
    my @ds = ($one);
    $cs[$e] = $ds[$e-1] = $one;
    for(my $k = 1, my $kk = $e - 1; $k <= $kk; ++$k, --$kk) {
        my ($c, $d) = ($nul, $nul);
        for (my $i = 0, my $j = $k-1; $j >= 0; $i += 2, --$j) {
            my $cc = $cs[$j];
            my $dd = $ds[$j];
            my ($q1, $q2) = @q[$i, $i+1];
            $c += $q1 * $n * $dd - $q2 * $cc;
            $d += $q[$i+2] * $cc - $q2 * $dd;
        }
        $c /= ($k + $k);
        $cs[$k] = $cs[$kk] = $c;
        $ds[$k] = $ds[$kk-1] = ($d + $c) / ($k + $k + 1) if $k < $kk;
    }

    my $cp = $poly->new(@cs);
    my $dp = $poly->new(@ds);
    $table->{"$n1 $n"} = [$cp, $dp];
    return ($cp, $dp);
}

# generate polynomials C, D, satisfying the identity of Beeger and Schinzel:
# for k > 1, square-free, k1 = k (if k === 1 (mod 4)), k1 = 2*k (else),
# n = odd multiple of k1: Phi_n(x) = C^2 - k*x*D^2
sub _schinzel_cd {
    my ($poly, $table, $n, $k) = @_;
    my $K = ($k & 3) == 1? $k: $k << 1;
    my $f = $n / $K;
    return () if not $f & 1;
    my $r = vecprod(map { $_->[0] } factor_exp($f));
    my $m = $f / $r * gcd($f, $k);
    if ($m > 1) {
        $n /= $m;
        $f /= $m;
    }
    my ($C, $D) = ();
    if (exists $table->{"$n $k"}) {
        ($C, $D) = @{$table->{"$n $k"}};
    }
    else {
        ($C, $D) = _lucas_cd($poly, $table, $k);
        my ($L, $M, $Q) = ();
        foreach my $p (factor($f)) {
            my $p2 = $p >> 1;
            my $kp2 = Math::BigInt->new($k)->bpow($p2);
            if (!defined $M) {
                $Q = $poly->monomial(2, $k);
                my $CC = $C->nest($Q);
                my $DD = $D->nest($Q)->shift_up(1)->mul_const($k);
                $L = $CC - $DD;
                $M = $CC + $DD;
            }
            if (kronecker($k, $p) < 0) {
                ($L, $M) = (
                    $L->nest($poly->monomial($p, $kp2)) / $M,
                    $M->nest($poly->monomial($p, $kp2)) / $L,
                );
            }
            else {
                ($L, $M) = (
                    $M->nest($poly->monomial($p, $kp2)) / $M,
                    $L->nest($poly->monomial($p, $kp2)) / $L,
                );
            }
        }
        if (defined $M) {
            $C = ($M + $L)->div_const(2)->unnest($Q);
            $D = (($M - $L) / $poly->monomial(1, $k << 1))->unnest($Q);
            $table->{"$n $k"} = [$C, $D];
        }
    }
    if ($m > 1) {
        $C = $C->inflate($m);
        $D = $D->inflate($m)->shift_up($m >> 1);
    }
    return ($C, $D);
}

# for r > 1, integer, square-free, s > 0, integer, x = r*s^2,
# r1 = r (if r === 1 (mod 4)), r1 = 2*r (else), n = odd multiple of r1,
# generate integer factors l < m such that l * m = Phi_n(x).
# if l = 1, return (m) else (l, m).
sub _schinzel_lm {
    my ($poly, $table, $n, $r, $s, $x) = @_;
    my ($c, $d) = _schinzel_cd($poly, $table, $n, $r);
    my $cx = $c->evaluate($x);
    my $dx = $d->evaluate($x) * $r * $s;
    my ($l, $m) = ($cx - $dx, $cx + $dx);
    return $m if $l == $poly->coeff_one;
    return ($l, $m);
}

# ----- Math::Polynomial extension -----

sub Math::Polynomial::cyclotomic {
    my ($this, $n, $table) = @_;
    return _cyclo_poly($this, $table || {}, $n);
}

sub Math::Polynomial::cyclo_factors {
    my ($this, $n, $table) = @_;
    my @d = divisors($n);
    $table ||= {};
    return map { _cyclo_poly($this, $table, $_, \@d) } @d;
}

sub Math::Polynomial::cyclo_int_factors {
    my ($this, $x, $n, $table) = @_;
    return ($this->coeff_zero) if !$n || $x == $this->coeff_one;
    return ($x - $this->coeff_one) if $n == 1 || !$x;
    $table ||= {};
    if (my $exp = is_power($x, 0, \my $root)) {
        $x = $root;
        $n *= $exp;
    }
    my ($r, $s) = _square_free_square($x);
    my $r1 = (1 == ($r & 3))? $r: $r+$r;
    my @d = divisors($n);
    my $i = $x == 2 || 0;
    return
        map {
            !($_ % $r1) && ($_ / $r1) & 1?
            _schinzel_lm($this, $table, $_, $r, $s, $x):
            (_cyclo_poly($this, $table, $_, \@d))->evaluate($x)
        } @d[$i .. $#d];
}

sub Math::Polynomial::cyclo_plusfactors {
    my ($this, $n, $table) = @_;
    my @d = divisors($n << 1);
    my $m = $n ^ ($n - 1);
    $table ||= {};
    return
        map { _cyclo_poly($this, $table, $_, \@d) }
        grep { !($_ & $m) } @d;
}

sub Math::Polynomial::cyclo_int_plusfactors {
    my ($this, $x, $n, $table) = @_;
    my $u = $this->coeff_one;
    return ($u + $u) if !$n;
    return ($x + $u) if $n == 1 || !$x || $x == $u;
    if (my $exp = is_power($x, 0, \my $root)) {
        $x = $root;
        $n *= $exp;
    }
    my ($r, $s) = _square_free_square($x);
    my $r1 = (1 == ($r & 3))? $r: $r+$r;
    my @d = divisors($n << 1);
    my $m = $n ^ ($n - 1);
    $table ||= {};
    return
        map {
            !($_ % $r1) && ($_ / $r1) & 1?
            _schinzel_lm($this, $table, $_, $r, $s, $x):
            (_cyclo_poly($this, $table, $_, \@d))->evaluate($x)
        } grep { !($_ & $m) } @d;
}

sub Math::Polynomial::cyclo_lucas_cd {
    my ($this, $n, $table) = @_;
    if ($n <= 1 || !is_square_free($n)) {
        croak "$n: not a square-free integer greater than one";
    }
    return _lucas_cd($this, $table || {}, $n);
}

sub Math::Polynomial::cyclo_schinzel_cd {
    my ($this, $n, $k, $table) = @_;
    if ($k <= 1 || !is_square_free($k)) {
        croak "$k: not a square-free integer greater than one";
    }
    my @cd = _schinzel_cd($this, $table || {}, $n, $k);
    if (!@cd) {
        my $k1 = ($k & 3) == 1? 'k': '2*k';
        croak "$n: n is not an odd multiple of $k1";
    }
    return @cd;
}

sub Math::Polynomial::cyclo_poly_iterate {
    my ($this, $n, $table) = @_;
    $n ||= 1;
    --$n;
    $table ||= {};
    return
        sub {
            ++$n;
            _cyclo_poly($this, $table, $n);
        };
}

# ----- public subroutines -----

sub cyclo_poly            { $default->cyclotomic(@_)            }
sub cyclo_factors         { $default->cyclo_factors(@_)         }
sub cyclo_plusfactors     { $default->cyclo_plusfactors(@_)     }
sub cyclo_poly_iterate    { $default->cyclo_poly_iterate(@_)    }
sub cyclo_lucas_cd        { $default->cyclo_lucas_cd(@_)        }
sub cyclo_schinzel_cd     { $default->cyclo_schinzel_cd(@_)     }
sub cyclo_int_factors     { $default->cyclo_int_factors(@_)     }
sub cyclo_int_plusfactors { $default->cyclo_int_plusfactors(@_) }

1;

__END__