ContFrac arithmetics

This commit is contained in:
Jiri Pittner 2022-02-21 16:45:44 +01:00
parent b4aaa77da4
commit a032344c66
3 changed files with 229 additions and 67 deletions

View File

@ -135,68 +135,206 @@ for(int i=1; i<=n; ++i) //truncate if possible
} }
template <typename T>
Homographic<T> Homographic<T>::input(const T &z, const bool inf) const
{
Homographic<T> hnew;
if(inf) //effective infinity, end of input
{
hnew.v[0][0]=hnew.v[0][1]=v[0][1];
hnew.v[1][0]=hnew.v[1][1]=v[1][1];
}
else
{
hnew.v[0][0]=v[0][1];
hnew.v[1][0]=v[1][1];
hnew.v[0][1]=v[0][0]+v[0][1]* z;
hnew.v[1][1]=v[1][0]+v[1][1]* z;
}
return hnew;
}
template <typename T>
Homographic<T> Homographic<T>::output(const T &z) const
{
Homographic<T> hnew;
hnew.v[0][0]=v[1][0];
hnew.v[0][1]=v[1][1];
hnew.v[1][0]=v[0][0]-v[1][0]*z;
hnew.v[1][1]=v[0][1]-v[1][1]*z;
return hnew;
}
template <typename T>
bool Homographic<T>::outputready(T &z) const
{
bool inf=0;
T q0,q1;
if(v[1][0]==0) inf=1; else q0=v[0][0]/v[1][0];
if(v[1][1]==0) inf=1; else q1=v[0][1]/v[1][1];
if(!inf && q0==q1) {z=q0; return true;}
return false;
}
template <typename T>
bool Homographic<T>::terminate() const
{
return v[1][0]==0&&v[1][1]==0;
}
template <typename T> template <typename T>
ContFrac<T> Homographic<T>::value(const ContFrac<T>&x) const ContFrac<T> Homographic<T>::value(const ContFrac<T>&x) const
{ {
Homographic<T> h(*this); Homographic<T> h(*this);
std::list<T> l; std::list<T> l;
typename ContFrac<T>::iterator px=x.begin(); for(typename ContFrac<T>::iterator px=x.begin(); px!=x.beyondend(); ++px)
do //scan all input
{ {
//digest next input term //digest next input term
Homographic<T> hnew; h=h.input(*px,px==x.end()|| px!=x.begin()&& *px==0);
if(px==x.end()|| px!=x.begin()&& *px==0) //input is infinity
//output as much as possible
T out;
while(h.outputready(out))
{ {
hnew.x[0][0]=hnew.x[0][1]=h.x[0][1]; l.push_back(out);
hnew.x[1][0]=hnew.x[1][1]=h.x[1][1]; h=h.output(out);
}
else
{
hnew.x[0][0]=h.x[0][1];
hnew.x[1][0]=h.x[1][1];
hnew.x[0][1]=h.x[0][0]+h.x[0][1]* *px;
hnew.x[1][1]=h.x[1][0]+h.x[1][1]* *px;
} }
//std::cout<<"hnew\n"<< hnew.x[0][0]<<" "<< hnew.x[0][1]<<"\n"<< hnew.x[1][0]<<" "<< hnew.x[1][1]<<"\n"; //terminate if exhausted
if(h.terminate())
//check if we are ready to output
bool inf=0;
T q0,q1;
if(hnew.x[1][0]==0) inf=1; else q0=hnew.x[0][0]/hnew.x[1][0];
if(hnew.x[1][1]==0) inf=1; else q1=hnew.x[0][1]/hnew.x[1][1];
while(!inf && q0==q1) //ready to output
{ {
l.push_back(q0);
//std::cout <<"out "<<q0<<std::endl;
Homographic<T> hnew2;
hnew2.x[0][0]=hnew.x[1][0];
hnew2.x[0][1]=hnew.x[1][1];
hnew2.x[1][0]=hnew.x[0][0]-hnew.x[1][0]*q0;
hnew2.x[1][1]=hnew.x[0][1]-hnew.x[1][1]*q0;
//std::cout<<"hnew2\n"<< hnew2.x[0][0]<<" "<< hnew2.x[0][1]<<"\n"<< hnew2.x[1][0]<<" "<< hnew2.x[1][1]<<"\n";
hnew=hnew2;
inf=0;
if(hnew.x[1][0]==0) inf=1; else q0=hnew.x[0][0]/hnew.x[1][0];
if(hnew.x[1][1]==0) inf=1; else q1=hnew.x[0][1]/hnew.x[1][1];
}
//termination
if(hnew.x[1][0]==0&&hnew.x[1][1]==0) //terminate
{
//std::cout<<"terminate at "<<i<<"\n";
if(px!=x.end()) laerror("unexpected termination in Homographic::value"); if(px!=x.end()) laerror("unexpected termination in Homographic::value");
break; break;
} }
h=hnew; }
++px;
} while(px!=x.beyondend());
//std::cout <<"size= "<<l.size()<<std::endl;
return ContFrac<T>(l); return ContFrac<T>(l);
} }
template <typename T>
BiHomographic<T> BiHomographic<T>::inputx(const T &x, const bool inf) const
{
BiHomographic<T> hnew;
for(int i=0; i<2; ++i)
{
hnew.v[i][0][0]= v[i][0][1];
hnew.v[i][0][1]= inf?v[i][0][1] : v[i][0][0]+v[i][0][1]*x;
hnew.v[i][1][0]= v[i][1][1];
hnew.v[i][1][1]= inf?v[i][1][1] : v[i][1][0]+v[i][1][1]*x;
}
return hnew;
}
template <typename T>
BiHomographic<T> BiHomographic<T>::inputy(const T &y, const bool inf) const
{
BiHomographic<T> hnew;
for(int i=0; i<2; ++i)
{
hnew.v[i][0][0]= v[i][1][0];
hnew.v[i][0][1]= v[i][1][1];
hnew.v[i][1][0]= inf?v[i][1][0] : v[i][0][0]+v[i][1][0]*y;
hnew.v[i][1][1]= inf?v[i][1][1] : v[i][0][1]+v[i][1][1]*y;
}
return hnew;
}
template <typename T>
BiHomographic<T> BiHomographic<T>::output(const T &z) const
{
BiHomographic<T> hnew;
for(int i=0; i<2; ++i) for(int j=0; j<2; ++j)
{
hnew.v[0][i][j]= v[1][i][j];
hnew.v[1][i][j]= v[0][i][j] - v[1][i][j]*z;
}
return hnew;
}
template <typename T>
int BiHomographic<T>::inputselect() const
{
if(v[1][0][0]==0)
{
if(v[1][0][1]==0) return 1;
else return 0;
}
if(v[1][0][1]==0) return 0;
if(v[1][1][0]==0) return 1;
if(MYABS(v[0][0][1]/v[1][0][1] - v[0][0][0]/v[1][0][0]) > MYABS(v[0][1][0]/v[1][1][0] - v[0][0][0]/v[1][0][0])) return 0;
return 1;
}
template <typename T>
bool BiHomographic<T>::outputready(T &z) const
{
T q[2][2];
for(int i=0; i<2; ++i) for(int j=0; j<2; ++j)
{
if(v[1][i][j]==0) return false;
else q[i][j]=v[0][i][j]/v[1][i][j];
if(q[i][j]!=q[0][0]) return false;
}
z=q[0][0];
return true;
}
template <typename T>
bool BiHomographic<T>::terminate() const
{
return v[1][0][0]==0&&v[1][0][1]==0&&v[1][1][0]==0&&v[1][1][1]==0;
}
template <typename T>
ContFrac<T> BiHomographic<T>::value(const ContFrac<T>&x, const ContFrac<T>&y) const
{
BiHomographic<T> h(*this);
std::list<T> l;
typename ContFrac<T>::iterator px=x.begin();
typename ContFrac<T>::iterator py=y.begin();
do
{
//select next input term
int which;
if(px==x.beyondend()) which=1;
else if(py==y.beyondend()) which=0;
else which = h.inputselect();
if(which) {h=h.inputy(*py,py==y.end()|| py!=y.begin()&& *py==0); ++py;}
else {h=h.inputx(*px,px==x.end()|| px!=x.begin()&& *px==0); ++px;}
//output as much as possible
T out;
while(h.outputready(out))
{
l.push_back(out);
h=h.output(out);
}
//terminate if exhausted
if(h.terminate())
{
if(px!=x.end()&&px!=x.beyondend() || py!=y.end()&&py!=y.beyondend()) laerror("unexpected termination in Homographic::value");
break;
}
}
while(px!=x.beyondend() || py!=y.beyondend());
return ContFrac<T>(l);
}
template <typename T> template <typename T>
void Rational<T>::simplify() void Rational<T>::simplify()
{ {
@ -206,7 +344,7 @@ if(den<0)
den= -den; den= -den;
} }
T g=gcd(num,den); T g=gcd(num,den);
if(g>1) if(MYABS(g)>1)
{ {
num/=g; num/=g;
den/=g; den/=g;
@ -218,7 +356,7 @@ Rational<T> & Rational<T>::operator*=(const T &rhs)
{ {
T r=rhs; T r=rhs;
T g=gcd(r,den); T g=gcd(r,den);
if(g>1) if(MYABS(g)>1)
{ {
r/=g; r/=g;
den/=g; den/=g;
@ -233,7 +371,7 @@ Rational<T> & Rational<T>::operator/=(const T &rhs)
{ {
T r=rhs; T r=rhs;
T g=gcd(r,num); T g=gcd(r,num);
if(g>1) if(MYABS(g)>1)
{ {
r/=g; r/=g;
num/=g; num/=g;
@ -270,19 +408,20 @@ Rational<T> & Rational<T>::operator*=(const Rational &rhs)
Rational r(rhs); Rational r(rhs);
T g; T g;
g=gcd(num,r.den); g=gcd(num,r.den);
if(g>1) if(MYABS(g)>1)
{ {
num/=g; num/=g;
r.den/=g; r.den/=g;
} }
g=gcd(den,r.num); g=gcd(den,r.num);
if(g>1) if(MYABS(g)>1)
{ {
den/=g; den/=g;
r.num/=g; r.num/=g;
} }
num*=r.num; num*=r.num;
den*=r.den; den*=r.den;
if(den<0) {den= -den; num= -num;}
return *this; return *this;
} }
@ -293,32 +432,20 @@ return *this;
* forced instantization in the corresponding object file * forced instantization in the corresponding object file
******************************************************************************/ ******************************************************************************/
template class Rational<int>; template class Rational<int>;
template class Rational<unsigned int>;
template class Rational<long>; template class Rational<long>;
template class Rational<unsigned long>;
template class Rational<long long>; template class Rational<long long>;
template class Rational<unsigned long long>;
template class ContFrac<int>; template class ContFrac<int>;
template class ContFrac<unsigned int>;
template class ContFrac<long>; template class ContFrac<long>;
template class ContFrac<unsigned long>;
template class ContFrac<long long>; template class ContFrac<long long>;
template class ContFrac<unsigned long long>;
template class Homographic<int>; template class Homographic<int>;
template class Homographic<unsigned int>;
template class Homographic<long>; template class Homographic<long>;
template class Homographic<unsigned long>;
template class Homographic<long long>; template class Homographic<long long>;
template class Homographic<unsigned long long>;
template class BiHomographic<int>; template class BiHomographic<int>;
template class BiHomographic<unsigned int>;
template class BiHomographic<long>; template class BiHomographic<long>;
template class BiHomographic<unsigned long>;
template class BiHomographic<long long>; template class BiHomographic<long long>;
template class BiHomographic<unsigned long long>;

View File

@ -34,6 +34,8 @@ namespace LA {
template <typename T> template <typename T>
class ContFrac; class ContFrac;
//@@@operator== > >= etc.
template <typename T> template <typename T>
class Rational { class Rational {
public: public:
@ -93,6 +95,8 @@ class BiHomographic;
//@@@operator== > >= etc.
template <typename T> template <typename T>
class ContFrac : public NRVec<T> { class ContFrac : public NRVec<T> {
private: private:
@ -133,7 +137,10 @@ public:
ContFrac operator/(const T &rhs) const {Homographic<T> h({{0,1},{rhs,0}}); return h.value(*this);}; ContFrac operator/(const T &rhs) const {Homographic<T> h({{0,1},{rhs,0}}); return h.value(*this);};
//arithmetics with two ContFrac operands //arithmetics with two ContFrac operands
//@@@ ContFrac operator+(const ContFrac &rhs) const {BiHomographic<T> h({{{0,1},{1,0}},{{1,0},{0,0}}}); return h.value(*this,rhs);};
ContFrac operator-(const ContFrac &rhs) const {BiHomographic<T> h({{{0,1},{-1,0}},{{1,0},{0,0}}}); return h.value(*this,rhs);};
ContFrac operator*(const ContFrac &rhs) const {BiHomographic<T> h({{{0,0},{0,1}},{{1,0},{0,0}}}); return h.value(*this,rhs);};
ContFrac operator/(const ContFrac &rhs) const {BiHomographic<T> h({{{0,1},{0,0}},{{0,0},{1,0}}}); return h.value(*this,rhs);};
//iterator //iterator
class iterator { class iterator {
@ -161,11 +168,16 @@ public:
template <typename T> template <typename T>
class Homographic { class Homographic {
public: public:
T x[2][2]; //{{a,b},{c,d}} for (a+b.z)/(c+d.z) T v[2][2]; //{{a,b},{c,d}} for (a+b.x)/(c+d.x) i.e. [denominator][power_x]
Homographic(){}; Homographic(){};
explicit Homographic(const T (&a)[2][2]) {memcpy(x,a,2*2*sizeof(T));}; explicit Homographic(const T (&a)[2][2]) {memcpy(v,a,2*2*sizeof(T));};
ContFrac<T> value(const ContFrac<T>&x) const; ContFrac<T> value(const ContFrac<T>&z) const;
Homographic input(const T &x, const bool inf) const;
Homographic output(const T &x) const;
bool outputready(T &x) const;
bool terminate() const;
}; };
@ -173,11 +185,19 @@ T x[2][2]; //{{a,b},{c,d}} for (a+b.z)/(c+d.z)
template <typename T> template <typename T>
class BiHomographic { class BiHomographic {
public: public:
T x[2][2][2]; T v[2][2][2]; //{{{a,b},{c,d}},{{e,f},{g,h}}} i.e.[denominator][power_y][power_x]
BiHomographic(){}; BiHomographic(){};
explicit BiHomographic(const T (&a)[2][2][2]) {memcpy(x,a,2*2*2*sizeof(T));}; explicit BiHomographic(const T (&a)[2][2][2]) {memcpy(v,a,2*2*2*sizeof(T));};
ContFrac<T> value(const ContFrac<T>&x, const ContFrac<T>&y) const; ContFrac<T> value(const ContFrac<T>&x, const ContFrac<T>&y) const;
BiHomographic inputx(const T &x, const bool inf) const;
BiHomographic inputy(const T &y, const bool inf) const;
BiHomographic output(const T &z) const;
int inputselect() const;
bool outputready(T &x) const;
bool terminate() const;
}; };

17
t.cc
View File

@ -2431,7 +2431,7 @@ NRMat<double> mm=m.permuted_rows(p);
cout <<mm; cout <<mm;
} }
if(1) if(0)
{ {
Rational<int> p,q; Rational<int> p,q;
cin>>p>>q; cin>>p>>q;
@ -2483,5 +2483,20 @@ ContFrac<int> z= x*Rational<int>({2,3});
cout<<Rational<int>(z)<<endl; cout<<Rational<int>(z)<<endl;
} }
if(1)
{
ContFrac<int> x(11,101);
ContFrac<int> v(3,7);
ContFrac<int> y= x+v;
cout<<Rational<int>(y)<<endl;
ContFrac<int> z= x*v;
cout<<Rational<int>(z)<<endl;
BiHomographic<int> h({{{12,4},{3,1}},{{0,1},{-1,0}}});
ContFrac<int> zz=h.value(x,v);
cout<<Rational<int>(zz)<<endl;
cout<<(Rational<int>(x)+3)*(Rational<int>(v)+4)/(Rational<int>(x)-Rational<int>(v))<<endl;
}
} }