#!/usr/bin/perl -w
# Copyright 2010, 2011, 2012 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.006;
use strict;
use warnings;
use POSIX 'floor', 'fmod';
use Math::Trig 'pi', 'atan';
use Math::BigFloat try => 'GMP';
use Math::Libm 'hypot';
use Math::PlanePath::TheodorusSpiral;
use Smart::Comments;
{
# Euler summation
# k<n n n
# sum f(k) = integ f(x) dx - (f(n)-f(1))/2 + integ B1(frac(x))*f'(x) dx
# k=1 1 1
#
# B1(x) = x-1/2
#
# k<n n
# sum f(k) = integ f(x) dx
# k=1 1
#
# - (f(n)-f(1))/2
#
# + B2/2! (f'(n) - f'(1))
# ...
# (-1)^m*Bm (m-1) (m-1)
# + --------- * (f (n) - f (1))
# m!
# ...
# B0=1
# B1=-1/2
my $B2 = 1/6;
my $B3 = 0;
my $B4 = -1/30;
# f(x) = arctan 1/sqrt(x)
# f'(x) = 1/(1+x^2)
# f'2(x) = (-1 * (2 * x)) / (((x ^ 2) + 1) ^ 2)
# = -2x / (x^2 + 1)^2
foreach my $x (1 .. 40) {
my $sum = fsum($x);
my $ifx = ifx($x) - ifx(1) - (fx($x)-fx(1))/2;
my $t1 = $B2/2 * (dfx($x) - dfx(1));
my $if2 = $ifx + $t1;
printf "%.4f %.4f[%.4f] %.4f[%.4f]\n",
$sum, $ifx, $ifx-$sum, $if2, $if2-$sum;
}
sub ifx {
my ($x) = @_;
return sqrt($x) + $x*fx($x) - atan2($x,1);
}
sub fx {
my ($x) = @_;
return atan2(1, sqrt($x));
}
sub dfx {
my ($x) = @_;
return 1/($x*$x+1);
}
sub ddfx {
my ($x) = @_;
return -2*$x/(($x*$x+1)**2);
}
sub fsum {
my ($x) = @_;
my $ret = 0;
foreach my $i (1 .. $x-1) {
$ret += fx($x);
}
return $ret;
}
exit 0;
}
{
{
package Math::Symbolic::Custom::MySimplification;
use base 'Math::Symbolic::Custom::Simplification';
# use Math::Symbolic::Custom::Pattern;
# my $formula = Math::Symbolic->parse_from_string("TREE_a * (TREE_b / TREE_c)");
# my $pattern = Math::Symbolic::Custom::Pattern->new($formula);
use Math::Symbolic::Custom::Transformation;
my $trafo = Math::Symbolic::Custom::Transformation::Group->new
(',',
'TREE_a * (TREE_b / TREE_c)' => '(TREE_a * TREE_b) / TREE_c',
'TREE_a * (TREE_b + TREE_c)' => 'TREE_a * TREE_b + TREE_a * TREE_c',
'(TREE_b + TREE_c) * TREE_a' => 'TREE_b * TREE_a + TREE_c * TREE_a',
# '(TREE_a / TREE_b) / TREE_c' => 'TREE_a / (TREE_b * TREE_c)',
'(TREE_a / TREE_b) / (TREE_c / TREE_d)'
=> '(TREE_a * TREE_d) / (TREE_b * TREE_c)',
'1 - TREE_a / TREE_b' => '(TREE_b - TREE_a) / TREE_b',
'TREE_a / TREE_b + TREE_c' => '(TREE_a + TREE_b * TREE_c) / TREE_b',
'(TREE_a / TREE_b) * TREE_c' => '(TREE_a * TREE_c) / TREE_b',
'TREE_a - (TREE_b + TREE_c)' => 'TREE_a - TREE_b - TREE_c',
'(TREE_a - TREE_b) - TREE_c' => 'TREE_a - TREE_b - TREE_c',
);
sub simplify {
my $tree = shift;
### simplify(): "$tree"
### traf: ($trafo->apply_recursive($tree)//'').''
return $trafo->apply_recursive($tree) || $tree;
# if (my $m = $pattern->match($tree)) {
# $m = $m->{'trees'};
# ### trees: $m
# ### return: ($m->{'a'} * $m->{'b'}) / $m->{'c'}
# return ($m->{'a'} * $m->{'b'}) / $m->{'c'};
# } else {
# ### no match
# return $tree;
# }
}
__PACKAGE__->register();
}
require Math::Symbolic;
require Math::Symbolic::Derivative;
{
my $t = Math::Symbolic->parse_from_string('1/(x^2+1)');
$t = Math::Symbolic::Derivative::total_derivative($t, 'x');
$t = $t->simplify;
print "$t\n";
exit 0;
}
{
my $a = Math::Symbolic->parse_from_string(
'(x+y)/(1-x*y)'
);
my $z = Math::Symbolic->parse_from_string(
'z'
);
my $t = ($a + $z) / (1 - $a*$z);
$t = $t->simplify;
print $t;
exit 0;
}
}
{
my $path = Math::PlanePath::TheodorusSpiral->new;
my $prev_x = 0;
my $prev_y = 0;
#for (my $n = 10; $n < 100000000; $n = int($n * 1.2)) {
foreach my $n (2000, 2010, 2020, 2010, 2000, 2010, 2000, 2010) {
my ($x,$y) = $path->n_to_xy($n);
my $rsq = $x*$x+$y*$y;
my $dx = $x - $prev_x;
my $dy = $y - $prev_y;
my $dxy_dist = hypot($dx,$dy);
printf "%d %.2f,%.2f %.2f %.4f\n", $n, $x,$y, $rsq, $dxy_dist;
($prev_x, $prev_y) = ($x,$y);
}
exit 0;
}
sub integral {
my ($x) = @_;
print "log ", log(1+$x*$x), " at x=$x\n";
return $x * atan($x) - 0.5 * log (1 + $x*$x);
}
print "integral 0 = ", integral(0), "\n";
print "integral 1 = ", integral(1)/(2*pi()), "\n";
print "atan 1 = ", atan(1)/(2*pi()), "\n";
sub est {
my ($n) = @_;
my $k = $n-1;
if ($k == 0) { return 0; }
my $K = 2.1577829966;
my $root = sqrt($k);
my $a = 2*pi()*pi();
my $radians;
$radians = integral(1/$root); # - integral(0);
# $radians = ($k+1)*atan(1/$root) + $root - 1/($root*$k);
return $radians / (2*pi());
# $radians = 2*$root;
# return $radians / (2*pi());
#
# $radians = $root - atan($root) + $k*atan(1/$root);
# return $radians / (2*pi());
#
# return $k / $a; # revolutions
# return $k / pi();
#
# return 2*$root / $a;
# $radians = 2*sqrt($k+1) + $K + 1/(6*sqrt($k+1)); # plus O(n^(-3/2))
# return 0.5 * $a * ($k * sqrt(1+$k*$k) + log($k + sqrt(1+$k*$k))) / $k;
# return $root + ($k+1)*atan(1/$root);
}
print "est 1 = ", est(1), "\n";
print "est 2 = ", est(2), "\n";
{
require Math::Polynomial;
open OUT, '>', '/tmp/theodorus.data' or die;
my @n;
my @theta;
my $total = 0;
foreach my $n (2 .. 120) {
my $inc = Math::Trig::atan(1/sqrt($n-1)) / (2*pi()); # revs
$total += $inc;
my $est = est($n);
my $diff = $total - $est;
# $diff = 1/$diff;
if ($n > 50) {
push @n, $n-51;
push @theta, $diff;
print OUT "$n $diff\n";
}
print "$n $inc $total $est $diff\n";
}
print "\n";
Math::BigFloat->accuracy(500);
my $p = Math::Polynomial->new; # (Math::BigFloat->new(0));
$p = $p->interpolate(\@n, \@theta);
foreach my $i (0 .. $p->degree) {
print "$i ",$p->coeff($i),"\n";
}
# $p->string_config({ fold_sign => 1,
# variable => 'n' });
# print "theta = $p\n";
close OUT or die;
system "xterm -e 'gnuplot <devel/theodorus.gnuplot; read'";
exit 0;
}
{
my $next = 1;
my $total = 0;
my $n = 1;
my $prev_n = 0;
my $prev_diff = 0;
my $total_diff_diff;
my $count_diff_diff;
for (;;) {
my $inc = Math::Trig::atan2(1,sqrt($n++)) / (2*pi());
$total += $inc;
if ($total >= $next) {
$next++;
my $diff = $n - $prev_n;
my $diff_diff = $diff - $prev_diff;
$total_diff_diff += $diff_diff;
$count_diff_diff++;
print "$n +$diff +$diff_diff $total\n";
if ($next >= 1000) {
last;
}
$prev_n = $n;
$prev_diff = $diff;
}
}
my $avg = $total_diff_diff / $count_diff_diff;
print "average $avg\n";
print "\n";
exit 0;
}
{
my $c2 = 2.15778;
my $t1 = 1.8600250;
my $t2 = 0.43916457;
my $z32 = 2.6123753486;
my $tn1 = 2*$t1 - 2*$t2 - $z32;
my $n = 1;
my $x = 1;
my $y = 0;
while ($n < 10000) {
my $r = sqrt($n); # before increment
($x, $y) = ($x - $y/$r, $y + $x/$r);
$n++;
$r = sqrt($n); # after increment
my $theta = atan2($y,$x);
if ($theta < 0) { $theta += 2*pi(); }
my $root;
$root = 2*sqrt($n) - $c2;
# $root += .01/$r;
# $root = -atan(sqrt($n)) + $n*atan(1/sqrt($n)) + sqrt($n);
# $root = atan(1/sqrt($n)) - pi()/2 + $n*atan(1/sqrt($n)) + sqrt($n);
$root = 2*sqrt($n)
+ 1/sqrt($n)
- $c2
# - 1/($n*sqrt($n))/3
# + 1/($n*$n*sqrt($n))/5
# - 1/($n*$n*sqrt($n))/7
# + 1/($n*$n*$n*sqrt($n))/9
;
# $root = -pi()/4 + Arctan($r);
# foreach my $k (2 .. 1000000) {
# $root += atan(1/sqrt($k)) - atan(1/sqrt($k + $r*$r - 1));
# # $root += atan( ($r*$r - 1) / ( ($k + $r*$r)*sqrt($k) + ($k+1)*sqrt($k+$r*$r-1)));
# }
# $root = -pi()/2 + Arctan($r) + $t1 *$r*$r/2 + ($tn1 - $t1)*$r**2/8;
$root = fmod ($root, 2*pi());
my $d = $root - $theta;
$d = fmod ($d + pi(), 2*pi()) - pi();
# printf "%10.6f %10.6f %23.20f\n", $theta, $root, $d;
printf "%23.20f\n", $d;
}
exit 0;
}
{
my $t1 = 0;
foreach my $k (1 .. 100) {
$t1 += 1 / (sqrt($k) * ($k+1));
printf "%10.6f\n", $t1;
}
exit 0;
}
sub Arctan {
my ($r) = @_;
return pi()/2 - atan(1/$r);
}
{
Math::BigFloat->accuracy(200);
my $bx = Math::BigFloat->new(1);
my $by = Math::BigFloat->new(0);
my $x = 1;
my $y = 0;
my $n = 1;
my @n = ($n);
my @x = ($x);
my @y = ($y);
my $count = 0;
my $prev_n = 0;
my $prev_d = 0;
my @dd;
while ($n++ < 10000000) {
my $r = hypot($x,$y);
my $py = $y;
($x, $y) = ($x - $y/$r, $y + $x/$r);
if ($py < 0 && $y >= 0) {
my $d = $n-$prev_n;
my $dd = $d-$prev_d;
push @dd, $dd;
printf "%5d +%4d +%3d %7.3f %10.6f %10.6f\n",
$n,
$d,
$dd,
# (sqrt($n)-1.07)/pi(),
sqrt($n),
$x, $y;
$prev_n = $n;
$prev_d = $d;
if (++$count >= 10) {
push @n, $n;
push @x, $x;
push @y, $y;
$count = 0;
}
}
}
print "average dd ", List::Util::sum(@dd)/scalar(@dd),"\n";
# require Data::Dumper;
# print Data::Dumper->new([\@n],['n'])->Indent(1)->Dump;
# print Data::Dumper->new([\@x],['x'])->Indent(1)->Dump;
# print Data::Dumper->new([\@y],['y'])->Indent(1)->Dump;
# require Math::Polynomial;
# my $p = Math::Polynomial->new(0);
# $p = $p->interpolate([ 1 .. @nc ], \@nc);
# $p->string_config({ fold_sign => 1,
# variable => 'd' });
# print "N = $p\n";
exit 0;
}
{
Math::BigFloat->accuracy(200);
my $bx = Math::BigFloat->new(1);
my $by = Math::BigFloat->new(0);
my $x = 1;
my $y = 0;
my $n = 1;
while ($n++ < 10000) {
my $r = hypot($x,$y);
($x, $y) = ($x - $y/$r, $y + $x/$r);
my $br = sqrt($bx*$bx + $by*$by);
($bx, $by) = ($bx - $by/$br, $by + $bx/$br);
}
my $ex = "$bx" + 0;
my $ey = "$by" + 0;
printf "%10.6f %10.6f %23.20f\n", $ex, $x, $ex - $x;
exit 0;
}