// This is my implementation of H.Cohen's Algorithm 5.7.2 "Fundamental Unit using Continued Fractions"
//Given a Fundamental discriminate D > 0
//The result of Algorithm is the Fundamental Unit Q(sqrt D).
//The main aim is to show that for CFRAC "Although expansion is ultimately periodic
//, the expansion never computes on a whole period or even a half period.
//It simply computes " from Pg 478 "A course in computational algebraic number theory" by Henri Cohen
//So why I included 5.7.2. Well worth including, as Modulo decision practise makes perfect ! ie (//if d = D(mod 2))?
//I altered this somewhat, as the Algorithm as "I interpreted it" Didn't work for Modulo as can be seen ! Sorry
//My approach to Modulo is via 'Prime Ends' and Continued Fraction PQa. The idea of a conditional check Modulo, alien to me.
//If this was an alien, then it's 'visible form' could reduce to something much smaller. Not C++ Code, I would copy !
//
//This is better code than Algorithm 5.7.1 Check all solutions as implementation of Modulo differs.
#include
#include
using namespace std; //Some say using namespace is not prefered to manually entering std
static double cfrac5_7_2(double D1) {
double au1, au2, b1, v1, v2, p1, p2, dD1, q1, q2, A1, t1;
int result1 = 0;
bool correct1{}; // So is False
//step 1 [Initialise]
dD1 = sqrt(D1);
if (((dD1 - D1) / 2) == (floor((dD1 - D1) / 2))) { //if d = D(mod 2)
b1 = dD1;
}
else{
b1 = (dD1 - 1);
}
au1 = -b1;
au2 = 2;
v1 = 1;
v2 = 0;
p1 = b1;
q1 = 2;
//q2 = (D1 - (pow(p1, 2))) / q1;
//step 2 [Euclidean step]
do{
A1 = ((p1 + dD1) / q1);
t1 = p1;
p1 = ((A1 * q1) - p1);
if ((t1 == p1) && (v2 != 0)) {
au1 = ((pow(au2, 2)) + (pow(v2, 2)))/q1;
v1 = ((2 * au2 * v2) / q1);
result1 = ((au1 + (v1 * sqrt(D1))) / 2);
return result1;
correct1 == true;
break;
}
else {
t1 = (A1 * au2) + au1;
au1 = au2;
au2 = t1;
t1 = (A1 * v2) + v1;
v1 = v2;
v2 = t1;
}
if ((q1 == t1) && (v2 != 0)) {
au1 = (((au1 * au2) + (D1 * v1 * v2)) / q1);
v1 = (((au1 * v2) + (au2 * v1)) / q1);
result1 = ((au1 + (v1 * sqrt(D1))) / 2);
return result1;
correct1 == true;
break;
}
} while (correct1 == false);
}
int main() {
int D1 = 0;
cout << "Enter D " << endl;
cin >> D1;
cout << "Check Fundamental Unit of Q(sqrt D) is = " << cfrac5_7_2(D1) << endl;
}