#!/usr/bin/perl
# bigfactor - calculate prime factors
use strict;
use Math::BigInt;
use integer;
use vars qw{ $opt_s };
use Getopt::Std;

@ARGV && getopts('s')         or die "usage: $0 [-s] number ...";

foreach my $orig (@ARGV) {
    my ($quo, $rem, $n, $p);
    $n = $orig;
    $n =~ s/\D//g;
    $n =~ s/^0+//;
    unless (length($n)) {
        warn "Invalid positive integer $orig\n" ;
        next;
    }
    if (length ($n) == 1) {
        # Primality test fails on primes of length 1
        if ($n == 1) {
            print "1: Neither prime nor composite\n";
        }
        elsif ($n == 2) {
            print "2: prime\n";
            next;
        }
        elsif ($n == 3) {
            print "3: prime\n";
            next;
        }
        elsif ($n == 5) {
            print "5: prime\n";
            next;
        }
        elsif ($n == 7) {
            print "7: prime\n";
            next;
        }
    }
    my $ptest = $n;
    my $n = Math::BigInt->new($orig);
    if ($n < 2**4000000000) {
        $n =~ s/\+//; # Make it a normal number
    }
    print "$ptest:";
    # Handle 2 as a special case
    if (($n % 2) == 0) {
        &divide_out($n, 2);
        if ($n == 1) {
            print "\n";
            next;
        }
    }
    # Handle 3 as a special case
    if (($n % 3) == 0) {
        &divide_out($n, 3);
        if ($n == 1) {
            print "\n";
            next;
        }
    }
    # Now the general case
    # Start with 1, and move up 4, try, 2, try, 4, try, 2, try etc
    # keeping the quot and rem up to date
    # This skips things divisible by 2 and 3 without calculation.
    my $p = 1;
    while (1) { # Endless loop, will leave with last, enter at 1
        # Move to the next quotient and remainder (ie +4)
        $p += 4;
        my ($quo, $rem);
        if (ref($n) =~ /Math/) {
               no integer;
            ($quo, $rem) = $n->bdiv($p);
               last if $quo < $p;
           }
           else {
               last if  $n / $p < $p;
               $rem = $n % $p;
           }
        if (0 == $rem) {
            &divide_out($n, $p); # Got one, divide and print
        }
        # Move to the next quotient and remainder (ie +2)
        $p += 2;
        $rem = $n % $p;
        if (0 == $rem) {
            &divide_out($n, $p); # Got one, divide and print
        }
    }
    if ($n == 1) {
        print "\n";
    }
    else {
        $n =~ s/\+//;
        if ($n eq $ptest) {
            print " prime\n";
        }
     else {
            print " $n\n";
        }
    }
}

sub divide_out {
    # Efficiency does not matter here
    my $factor = pop;
    my @factors;
    while (($_[0]%$factor) == 0) {
        push @factors, $factor;
        $_[0] = $_[0]/$factor;
    }
    s/\+// foreach @factors;
    # Do the print here so that the user gets feedback...
    if ($opt_s and (1 < scalar @factors)) {
        print " $factors[0]**", (scalar @factors);
    }
    else {
        print map {" $_"}  @factors;
    }
    if ($_[0] < 2**32) {
        $_[0] =~ s/\+//; # Make it a normal number
    }
}