ishka/manticore
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

More continued fraction numeric and more tests.

Browse Source
This commit is contained in:
Stephan Barth 2023-12-28 13:01:16 +01:00
parent 0f1c43023e
commit 7049ec6a11
3 changed files with 65 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
Manticore/Num/ContinuedFraction.pm
View File

@ -6,6 +6,7 @@ use strict;
use warnings;
use Math::BigRat;
use Data::Dumper;
### Methods and constructor
# Approximation object:
@ -19,9 +20,9 @@ sub new {
my ($pkg, $x) = @_;
die "ContinuedFraction can only be constructed of positive numbers!" if $x<0;
my $prop = int $x;
if(0 == $prop) {
return undef; # TODO
}
#if(0 == $prop) {
# return undef; # TODO
#}
my $rem = $x - $prop;
return bless {
a => $prop,
@ -85,9 +86,13 @@ sub reconstruct {
# round to bigRat with ContinuedFraction of $steps parts
sub roundSteps {
my ($steps, $x) = @_;
my $cf = Manticore::Num::ContinuedFraction->new($x);
#print "[[ $steps -- $x ]]\n";
return Math::BigRat->new(0) if 0==$x;
my $sgn = $x < 0 ? -1 : 1;
my $cf = Manticore::Num::ContinuedFraction->new($sgn*$x);
#print Data::Dumper::Dumper($cf);
my $l = $cf->firstk($steps);
return reconstruct($l);
return reconstruct($l)*$sgn;
}
1;

35
Manticore/Num/Trig.pm
View File

@ -3,6 +3,9 @@ package Manticore::Num::Trig;
use strict;
use warnings;
use Math::BigRat lib => 'GMP';
# Sin, cos, first terms
# 1/0! - x^2/2! + x^4/4! - x^6/6!
sub proxCos {
my ($steps, $x) = @_;
@ -31,4 +34,36 @@ sub proxSin {
return $sum
}
# sin,cos, first terms and round by continued fractions
# 1/0! - x^2/2! + x^4/4! - x^6/6!
sub proxCosCF {
my ($steps, $x) = @_;
my $round = sub { Manticore::Num::ContinuedFraction::roundSteps($steps, $_[0]) };
my $mxsq = $round->(-$x*$x);
my $sum = 0;
my $fa = 1;
for(1..$steps) {
my $step2 = $_ * 2;
$sum = $round->($sum+$fa);
$fa = $round->($fa*$mxsq/($step2*($step2-1)))
}
return $sum
}
# x/1! - x^3/3! + x^5/5! - x^7/7!
sub proxSinCF {
my ($steps, $x) = @_;
my $round = sub { Manticore::Num::ContinuedFraction::roundSteps($steps, $_[0]) };
my $mxsq = $round->(-$x*$x);
my $sum = 0;
my $fa = $x;
for(1..$steps) {
my $step2 = $_ * 2;
$sum = $round->($sum+$fa);
$fa = $round->($fa*$mxsq/($step2*($step2+1)))
}
return $sum
}
1;

24
tests/testCf.pl
View File

@ -16,10 +16,10 @@ use Manticore::Num::ContinuedFraction;
my $pi = Math::BigRat->new(3);
my $poly = 30;
my $poly = 10;
for(1..5) {
$pi = $pi + Manticore::Num::Trig::proxSin($poly,$pi);
$pi = $pi + Manticore::Num::Trig::proxSinCF($poly,$pi);
print "===\ns: $pi\n";
$pi = Manticore::Num::ContinuedFraction::roundSteps($poly,$pi);
print "r: $pi\n";
@ -27,7 +27,7 @@ print "r: $pi\n";
my $pi1 = Manticore::Num::ContinuedFraction->new($pi);
$pi = $pi + Manticore::Num::Trig::proxSin($poly+2,$pi);
$pi = $pi + Manticore::Num::Trig::proxSinCF($poly+2,$pi);
print "===\ns: $pi\n";
$pi = Manticore::Num::ContinuedFraction::roundSteps($poly,$pi);
print "r: $pi\n";
@ -38,7 +38,23 @@ my $pi2 = Manticore::Num::ContinuedFraction->new($pi);
#
#print Data::Dumper::Dumper($pi);
my $c = $pi1->commonk($pi2, 25);
my $c = $pi1->commonk($pi2, 2500);
print "[[ @$c ]]\n";
#my $max = 1;
#for(1..200000) {
# my $r = Manticore::Num::ContinuedFraction->new(rand 1);
# my $l = $r->firstk(50);
# for(@$l) { $max = $_ if $_>$max }
# print "[[ @$l ]]\n";
#}
#print "M: $max\n";
while(<>) {
my $cf = Manticore::Num::ContinuedFraction::roundSteps(10,$_);
my $f = Math::BigRat->new(0+$_);
my $r = Math::BigRat->new($_);
print "$cf $f $r\n";
}