#!/usr/bin/perl -w
# Copyright 2010, 2011 Kevin Ryde
# This file is part of Math-PlanePath.
#
# Math-PlanePath is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the Free
# Software Foundation; either version 3, or (at your option) any later
# version.
#
# Math-PlanePath is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
# for more details.
#
# You should have received a copy of the GNU General Public License along
# with Math-PlanePath. If not, see <http://www.gnu.org/licenses/>.
use 5.010;
use strict;
use warnings;
use POSIX 'fmod';
use Math::BigRat;
use Math::Prime::XS;
#use Smart::Comments;
use constant PHI => (1 + sqrt(5)) / 2;
# (3n-1)*n/2 pentagonal
# (3n+1)*n/2 second pentagonal
# http://www.research.att.com/~njas/sequences/A005449
# sum of n consecutive numbers >= n (n+1)+(n+2)+...+(n+n)
# triangular+square (n+1)*n/2 + n*n
# (3n+1)*n/2-2 = offset (3n+7)*n/2
# http://www.research.att.com/~njas/sequences/A140090
# sum n+1 to n+n-3 or some such
# (3n+1)*n/2
# (3n+1)*n/2 - 1
# (3n+1)*n/2 - 2
sub three {
my ($i) = @_;
return (3*$i+1)*$i/2 - 2;
}
sub is_perfect_square {
my ($n) = @_;
$n = sqrt($n);
return ($n == int($n));
}
{
my $prev_k = 0;
foreach my $k (0 .. 1000) {
my $sq = 24*$k+1;
if (is_perfect_square($sq)) {
printf "%4d %+4d %4d %4d\n", $k, $k-$prev_k, $k%24, $sq;
$prev_k = $k;
}
}
exit 0;
}
{
# i==0mod4 or 1mod4 always even
#
foreach my $k (0 .. 100) {
my $i = 4*$k + 2;
my $n = three($i);
my $factors = factorize($n);
printf "%4d %4d %s\n", $i,$n,$factors;
# unless ($factors =~ /\Q*/) {
# die;
# }
}
exit 0;
}
{
local $, = ',';
print map {three($_)} 0..20;
exit 0;
}
{
my $a = Math::BigRat->new('3/2');
my $b = Math::BigRat->new('1/2');
my $c = Math::BigRat->new('-2');
my $x = -$b;
my $sq = ($b*$b-4*$a*$c);
my $y = $sq;
$y->bsqrt;
print "$x $sq $y\n";
my $r1 = ($x + $y)/(2*$a);
my $r2 = ($x - $y)/(2*$a);
print "$r1 $r2\n";
exit 0;
}
{
foreach my $i (5 .. 500) {
my $n = three($i);
if (Math::Prime::XS::is_prime($n)) {
say "$i $n";
last;
}
}
exit 0;
}
sub factorize {
my ($n) = @_;
my @factors;
foreach my $f (2 .. int(sqrt($n)+1)) {
while (($n % $f) == 0) {
push @factors, $f;
### $n
$n /= $f;
}
}
if ($n != 1) {
push @factors, $n;
}
return join ('*',@factors);
}
exit 0;