PSGnuUtils

0.1.0

Utils/factor.ps1

Register-GnuUtil @{
    Name                   = 'factor'
    Category               = 'Shell helpers'
    Synopsis               = 'Factor numbers.'
    Usage                  = @('factor [NUMBER]...', 'factor OPTION')
    Description            = @(
        'Print the prime factors of each specified integer NUMBER. If none',
        'are specified on the command line, read them from standard input.'
    )
    NegativeNumberOperands = $true
    Options                = @(
        @{ Short = 'h'; Long = 'exponents'; Help = 'print repeated factors in form p^e unless e is 1' }
    )
    Notes                  = @(
        'Numbers of any size are accepted (System.Numerics.BigInteger, Miller-Rabin, Pollard rho).',
        'PowerShell note: much slower than GNU; numbers whose two smallest prime factors both have',
        'more than about 12 digits can take very long.'
    )
}



function script:Test-GnuProbablePrime {
    <#
    .SYNOPSIS
    Miller-Rabin test; deterministic for numbers below 3.3e24 with these bases.
    #>
    param([Parameter(Mandatory)][System.Numerics.BigInteger]$Value)

    if ($Value -lt 2) { return $false }
    foreach ($small in @(2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41)) {
        if ($Value -eq $small) { return $true }
        if ($Value % $small -eq 0) { return $false }
    }
    $d = $Value - 1
    $r = 0
    while ($d % 2 -eq 0) { $d = $d / 2; $r++ }
    foreach ($base in @(2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41)) {
        $x = [System.Numerics.BigInteger]::ModPow($base, $d, $Value)
        if ($x -eq 1 -or $x -eq $Value - 1) { continue }
        $composite = $true
        for ($i = 1; $i -lt $r; $i++) {
            $x = [System.Numerics.BigInteger]::ModPow($x, 2, $Value)
            if ($x -eq $Value - 1) { $composite = $false; break }
        }
        if ($composite) { return $false }
    }
    $true
}

function script:Find-GnuFactorPollard {
    <#
    .SYNOPSIS
    Finds a non-trivial factor of a composite number with Pollard's rho (Floyd cycle detection).
    #>
    param([Parameter(Mandatory)][System.Numerics.BigInteger]$Value)

    for ($c = 1; ; $c++) {
        $x = [System.Numerics.BigInteger]2
        $y = [System.Numerics.BigInteger]2
        $d = [System.Numerics.BigInteger]::One
        while ($d -eq 1) {
            $x = ($x * $x + $c) % $Value
            $y = ($y * $y + $c) % $Value
            $y = ($y * $y + $c) % $Value
            $d = [System.Numerics.BigInteger]::GreatestCommonDivisor([System.Numerics.BigInteger]::Abs($x - $y), $Value)
        }
        if ($d -ne $Value) { return $d }
    }
}

function script:Get-GnuPrimeFactor {
    <#
    .SYNOPSIS
    The prime factors of a number in ascending order (with repetition).
    #>
    param([Parameter(Mandatory)][System.Numerics.BigInteger]$Value)

    $factors = New-Object 'System.Collections.Generic.List[System.Numerics.BigInteger]'
    $n = $Value
    foreach ($p in @(2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47)) {
        while ($n -gt 1 -and $n % $p -eq 0) { $factors.Add($p); $n = $n / $p }
    }
    $pending = New-Object 'System.Collections.Generic.Stack[System.Numerics.BigInteger]'
    if ($n -gt 1) { $pending.Push($n) }
    while ($pending.Count -gt 0) {
        $current = $pending.Pop()
        if (Test-GnuProbablePrime -Value $current) { $factors.Add($current); continue }
        $divisor = Find-GnuFactorPollard -Value $current
        $pending.Push($divisor)
        $pending.Push($current / $divisor)
    }
    $factors.Sort()
    foreach ($factor in $factors) { $factor }
}