#!/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) {
÷_out($n, 2);
if ($n == 1) {
print "\n";
next;
}
}
# Handle 3 as a special case
if (($n % 3) == 0) {
÷_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) {
÷_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) {
÷_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
}
}