// This is an 'interpretation' of H.Cohen's Algorithm 5.7.1 "Fundamental Unit using Continued Fractions"
//The result of Algorithm is the Fundamental Unit, independently of the initial reduced Form
//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.1. Well worth a look, as Modulo decision practise makes perfect !
//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 !
// Condition being 4a | (D - b^2) and a > 0;
//So 4 * 3 = 12 | 144 = 153 - 9.
//So a = 3, D = 153, b = 3
//Reduced Form == (a,b, ((b^2 - D) / 4a)
//Reduced Form == 12
#include
#include
using namespace std; //Some say using namespace is not prefered to manually entering std
static double cfrac(double an, double bx, double D1) {
double au1, au2, v1, v2, p1, q1, q2, d1, t1, A1;
int result1 = 0;
bool correct1{}; // So is False
//step 1 [Initialise]
au1 = -bx;
au2 = (2 * an);
v1 = 1;
v2 = 0;
p1 = bx;
q1 = (2 * an);
d1 = sqrt(D1);
//step 2 [Euclidean step]
do {
A1 = ((p1 + d1) / q1);
p1 = floor((A1 * q1) - p1); //altered Algorithm doesn't say to put in 'floor'
q1 = ((D1 - (pow(p1, 2))) / q1);
// Then
t1 = ((A1 * au2) + au1);
au1 = au2;
au2 = t1;
t1 = ((A1 * v2) + v1);
v1 = v2;
v2 = t1;
//Step 3 [End of Period ?]
int q2 = (2 * an); /*Mod variable */
int bx1 = (bx); /*Mod variable */
// cout << " Mod " << (bx1 % q2) << endl; The correct "if" ((q1 == q2) && (p1 == (bx % (2 * an)) ie altered Algorithm
if ((q1 == q2) && (p1 == abs(((pow(bx, 2)) - D1 )/(4 * an)))){ //altered Algorithm
cout << "Reduced Form == " << p1 << endl;
au1 = (au2 / an); //exact division
v1 = (v2 / an); //exact division
result1 = ((au1 + (v1 * sqrt(D1))) / 2);
return result1;
correct1 == true;
break;
}
} while (correct1 == false);
}
int main() {
int an = 0;
int bx = 0;
int D1 = 0;
cout << "Enter a " << endl;
cin >> an;
cout << "Enter b " << endl;
cin >> bx;
cout << "Enter D " << endl;
cin >> D1;
cout << "Check if Fundamental Unit ? is = "<< cfrac(an, bx, D1) << endl;
}