basic routines for ContFrac

contfrac.cc

@@ -24,12 +24,124 @@
 
 namespace LA {
 
+template <typename T>
+ContFrac<T>::ContFrac(double x, const int n, const T thres) : NRVec<T>(n+1)
+{
+for(int i=0; i<=n; ++i)
+	{
+	NRVec<T>::v[i]=floor(x);
+	x -= NRVec<T>::v[i];
+	double y= 1./x;
+	if(x==0. || (thres && fabs(y)>thres)) {resize(i,true); return;}
+	x=y;
+	}
+}
+
+
+//we have to recursively first determine length and then allocate and fill the values during recursion unwinding
+template <typename T>
+static void cf_helper(ContFrac<T> *me, T p, T q, int level)
+{
+T div=p/q;
+{
+T rem=p%q;
+if(rem) cf_helper(me,q,rem,level+1);
+else me->resize(level);
+}
+(*me)[level]=div;
+}
+
+template <typename T>
+ContFrac<T>::ContFrac(const T p, const T q) : NRVec<T>()
+{
+cf_helper<T>(this,p,q,0);
+}
+
+
+template <typename T>
+ContFrac<T> ContFrac<T>::reciprocal() const
+{
+int n=this->length();
+if((*this)[0] == 0)
+	{
+	ContFrac<T> r(n-1);
+	for(int i=1; i<=n; ++i) r[i-1] = (*this)[i];
+	return r;
+	}
+else
+	{
+	ContFrac<T> r(n+1);
+	r[0]=0;
+	for(int i=0; i<=n; ++i) r[i+1] = (*this)[i];
+	return r;
+	}
+}
+
+
+template <typename T>
+void ContFrac<T>::convergent(T *p, T*q, const int trunc) const
+{
+int top=this->length();
+if(trunc != -1) top=trunc;
+NRVec<T> hh(top+3),kk(top+3);
+T *h= &hh[2]; 
+T *k= &kk[2];
+//start for recurrent relations
+h[-2]=k[-1]=0;
+h[-1]=k[-2]=1;
+for(int i=0; i<=top; ++i)
+	{
+	if(i>0 && (*this)[i]==0) //terminate by 0 which means infinity if not canonically shortened
+		{
+		*p=h[i-1];
+		*q=k[i-1];
+		return;
+		}
+	h[i] = (*this)[i]*h[i-1] + h[i-2];
+	k[i] = (*this)[i]*k[i-1] + k[i-2];
+	}
+*p=h[top];
+*q=k[top];
+}
+
+template <typename T>
+double ContFrac<T>::value(const int trunc) const
+{
+T p,q;
+convergent(&p,&q,trunc);
+double x=p;
+x/=q;
+return x;
+}
+
+
+template <typename T>
+void ContFrac<T>::canonicalize()
+{
+int n=this->length();
+if(n==0) return;
+if((*this)[n]==1) {(*this)[n]=0; ++(*this)[n-1];} //avoid deepest 1/1
+for(int i=1; i<=n; ++i) //truncate if possible
+	{
+	if((*this)[i]==0) //convention for infinity
+		{
+		resize(i-1,true);
+		return;
+		}
+	}
+}
+
 
 /***************************************************************************//**
  * forced instantization in the corresponding object file
  ******************************************************************************/
 template class ContFrac<int>;
 template class ContFrac<unsigned int>;
+template class ContFrac<long>;
+template class ContFrac<unsigned long>;
+template class ContFrac<long long>;
+template class ContFrac<unsigned long long>;
+
 
 #define INSTANTIZE(T) \
 

contfrac.h

@@ -25,17 +25,45 @@
 
 namespace LA {
 
+//simple finite continued fraction class
+//NOTE: 0 on any position >0 means actually infinity; simplify()  shortens the vector
+//presently implements just conversion to/from rationals and floats
+//maybe implement arithmetic by Gosper's method cf. https://perl.plover.com/classes/cftalk/TALK
+
+template <typename T>
+class Rational {
+public:
+	T num;
+	T den;
+
+	Rational(const T p, const T q) : num(p),den(q) {};
+};
 
 template <typename T>
 class ContFrac : public NRVec<T> {
+private:
+	int size() const; //prevent confusion with vector size
 public:
 	ContFrac():  NRVec<T>() {};
 	template<int SIZE>  ContFrac(const T (&a)[SIZE]) : NRVec<T>(a) {};
 	ContFrac(const NRVec<T> &v) : NRVec<T>(v) {}; //allow implicit conversion from NRVec
 	ContFrac(const int n) : NRVec<T>(n+1) {};
+	ContFrac(double x, const int n, const T thres=0); //might yield a non-canonical form
+	ContFrac(const T p, const T q); //should yield a canonical form
+	ContFrac(const Rational<T> r) : ContFrac(r.num,r.den) {};
 
+	void canonicalize();
+	void convergent(T *p, T*q, const int trunc= -1) const;
+	Rational<T> rational(const int trunc= -1) const {T p,q; convergent(&p,&q,trunc); return Rational<T>(p,q);};
+	double value(const int trunc= -1) const;
+	ContFrac reciprocal() const;
 	int length() const {return NRVec<T>::size()-1;};
-	void resize(const int n, const bool preserve=true) {NRVec<T>::resize(n+1,preserve);}
+	void resize(const int n, const bool preserve=true) 
+		{
+		int nold=length();
+		NRVec<T>::resize(n+1,preserve);
+		if(preserve) for(int i=nold+1; i<=n;++i) (*this)[i]=0;
+		} 
 	
 };
 

polynomial.h

@@ -30,6 +30,8 @@ namespace LA {
 
 template <typename T>
 class Polynomial : public NRVec<T> {
+private:
+	int size() const; //prevent confusion with vector size
 public:
 	Polynomial():  NRVec<T>() {};
 	template<int SIZE>  Polynomial(const T (&a)[SIZE]) : NRVec<T>(a) {};

t.cc

@@ -2172,7 +2172,7 @@ p[0]=p[1]=1;
 CycleIndex<int> c(Sn);
 PERM_RANK_TYPE d;
 Polynomial<int> pp=c.substitute(p,&d);
-for(int i=0; i<=p.size(); ++i) if(d!=pp[i]) laerror("error in cycle index");
+for(int i=0; i<p.degree(); ++i) if(d!=pp[i]) laerror("error in cycle index");
 }
 
 if(0)
@@ -2414,7 +2414,7 @@ double det = fit.solve();
 cout <<"det= "<<det<<" fit "<<fit<<endl;
 }
 
-if(1)
+if(0)
 {
 NRMat<double> m;
 cin >>m;
@@ -2431,4 +2431,37 @@ NRMat<double> mm=m.permuted_rows(p);
 cout <<mm;
 }
 
+if(1)
+{
+double x;
+cin >>x;
+ContFrac<long long> xx(x,20,1000000);
+cout<<xx;
+double y=xx.value();
+cout <<y<<endl;
+cout << "CF roundoff error = "<<x-y<<endl;
+}
+
+if(0)
+{
+int p,q;
+cin>>p>>q;
+ContFrac<int> xx(p,q);
+double z=p; z/=q;
+ContFrac<int> zz(z,20,100000);
+ContFrac<int> yy=xx.reciprocal();
+cout<<xx;
+cout<<yy;
+cout<<zz;
+int pp,qq,rr,ss;
+xx.convergent(&pp,&qq);
+yy.convergent(&rr,&ss);
+cout << pp<<" "<<qq<<endl;
+cout << rr<<" "<<ss<<endl;
+if(p!=pp ||q!=qq||qq!=rr||pp!=ss) cout<<"ContFrac error\n";
+double zzz=zz.value();
+cout <<z<<" "<<zzz<<endl;
+}
+
+
 }