From 7049ec6a110345c02518db1b8b2d9d53c0eb2ef1 Mon Sep 17 00:00:00 2001
From: Stephan Barth <stephan.barth@protonmail.com>
Date: Thu, 28 Dec 2023 13:01:16 +0100
Subject: [PATCH] More continued fraction numeric and more tests.

---
 Manticore/Num/ContinuedFraction.pm | 15 ++++++++-----
 Manticore/Num/Trig.pm              | 35 ++++++++++++++++++++++++++++++
 tests/testCf.pl                    | 24 ++++++++++++++++----
 3 files changed, 65 insertions(+), 9 deletions(-)

diff --git a/Manticore/Num/ContinuedFraction.pm b/Manticore/Num/ContinuedFraction.pm
index cd55c15..1d559e2 100644
--- a/Manticore/Num/ContinuedFraction.pm
+++ b/Manticore/Num/ContinuedFraction.pm
@@ -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;
diff --git a/Manticore/Num/Trig.pm b/Manticore/Num/Trig.pm
index 1135f00..4b0be7b 100644
--- a/Manticore/Num/Trig.pm
+++ b/Manticore/Num/Trig.pm
@@ -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;
diff --git a/tests/testCf.pl b/tests/testCf.pl
index 5a2b01f..88fcd98 100755
--- a/tests/testCf.pl
+++ b/tests/testCf.pl
@@ -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";
+}
+