diff --git a/README.md b/README.md index c00f452..11abbcb 100644 --- a/README.md +++ b/README.md @@ -63,19 +63,176 @@ this project uses submodules (maybe recursively), so either : -> Δ == 0 -> 1 solution : x = ( -b / 2a ) -> Δ < 0 -> 2 solutions : x = ( -b / 2a ) +- i( √|Δ| / 2a ) -> solution : - - a; // double - - b; // double - - c; // double - delta_sign; // + or - - delta_absolute; // |Δ| + - delta_sqrt; // √|Δ| + - first_term_gcd; // gcd(b, 2a) - first_term_numerator; // -b / gcd - first_term_denominator; // 2a / gcd - first_term; // double (-b / 2a) + - second_term_gcd; // gcd(√|Δ|, 2a) - second_term_numerator; // √|Δ| / gcd - second_term_denominator; // 2a / gcd - second_term; // double (√|Δ| / 2a) + - double solution1; // first_term + second_term - double solution2; // first_term - second_term 6. print solution + +--- + +# sqrt implementation + +finding the square root of x + +## dichotomy method + +``` + solution + ↓ +|------------------------------ 1. define range bound_1 : 0 +|-----------------------------| 2. choose range bound_2 : x +|--------------+--------------| 3. take mid : (bound_1 + bound_2) / 2 +|--------------|-------ø------| 4. choose range bound : bound_1 or bound_2 + +---------------|--------------- 1. range bound_1 = mid +|--------------|--------------- 2. range bound_2 = chosen previous bound +|------+-------|--------------- 3. mid +|--ø---|-------|--------------- 4. choose range + +-------|----------------------- 1. range bound_1 = mid +-------|-------|--------------- 2. range bound_2 = chosen previous bound +-------|---+---|--------------- 3. mid +-------|---|-ø-|--------------- 4. choose range + +-----------|------------------- 1. range bound_1 = mid +-------|---|------------------- 2. range bound_2 = chosen previous bound +-------|-+-|------------------- 3. mid +-------|-|ø|------------------- 4. choose range + +---------|--------------------- 1. range bound_1 = mid +-------|-|--------------------- 2. range bound_2 = chosen previous bound +-------|+|--------------------- 3. mid +-------ø|---------------------- 4. choose range + +--------|---------------------- --> solution +``` + + +## Newton–Raphson method + +it's like a self-correcting binary search, we get rid of the step "choose range" : + +``` + solution + ↓ +------------------------------| 1. define range bound_1 : v +--|---------------------------| 2. choose range bound_2 : x/v +--|-------------+-------------| 3. take mid : (v + x/v) / 2 + +----------------|-------------- 1. range bound_1 : v = mid +------|---------|-------------- 2. range bound_2 : x/v +------|----+----|-------------- 3. mid + +-----------|------------------- 1. range bound_1 : v = mid +-|---------|------------------- 2. range bound_2 : x/v +-|----+----|------------------- 3. mid + +------|------------------------ 1. range bound_1 : v = mid +------|-------|---------------- 2. range bound_2 : x/v +------|---+---|---------------- 3. mid + +----------|-------------------- 1. range bound_1 : v = mid +------|---|-------------------- 2. range bound_2 : x/v +------|-+-|-------------------- 3. mid + +--------|---------------------- --> solution +``` + +### mathematical proof that each range is automatically in the right range : + - if the value was higher than the answer, then new value is below old value, and vice versa + - how ? : + - define `x`, solution `s = √x`, and value `v = (old_value + x / old_value) / 2` + - supposition : if `v < s` , then `new_v > v`, else `new_v < v` : + - demonstration : + 1. if `v < s` : + v < s + <=> v < √x + <=> v² < x (that's actually how we know that v < s) + <=> v²/v < x/v + <=> v < x/v + -> and is s < x/v ? : + v < s + <=> v < √x + <=> v² < x (as previously) + <=> v² < x²/x + <=> v² * x < x² + <=> (v * √x)² < x² + <=> v * √x < x + <=> v * √x < v * x/v + <=> √x < x/v + -> so indeed : if v < √x, then v < √x < x/v == v < s < x/v + -> conclusion, the new range < v , x/v > contains the solution + 2. the same demonstration works for `v > s` + +### here is a more intuitive demonstration, with x = 20 : + +1. **show that if `v² > x` (== `v > s`) then `v > s > x/v`, and if `v² < x` (== `v < s`) then `v < s < x/v` :** + 1.1. **for value too high `v > s` :** + 1.1.1 **why `v > x/v` :** + - let's take initial value v = 5 : + - is 5² the solution ? 5² == 25 -> so no, 5 is not the sqrt, it's too high + ``` + v v² + 0 (5) 10 15 20 25 + v : |----|----|----|----|----| + x/v: |---|---|---|---|---| <----- squiz it, so the previous 5 portions fit the x = 20 size + 0 (4) 8 12 16 20 + x/v x + ``` + - the value of the new portion is 4, and we can visually see that it's lower than the previous portion 5 + - so : `v > x/v` + 1.1.2 **why `s > x/v` :** + - let's take the value v = 5 : + - we already showed that it's too high, now we will see that x/v == 20/5 is too low : + ``` + v + 0 (5) 10 15 20 25 + v : |----|----|----|----|----| + x/v: |---|---|---|---|---| <----- squizz + 0 *1 *2 *3 *4 *5 -> number of portions + 01234 -> portion size + (x/v)²: |---|---|---|---| + 0 4 8 12 16 + ``` + - the portion size is smaller than the number of portions, so it's too small to be the sqrt, indeed we visually see that this portion size `x/v` is a root of a smaller number : `(x/v)² == 16` + - so : `s > x/v` + 1.1.3. **conclusion :** + - v > s + - and v > x/v (<- this proof is not essential) + - and s > x/v (<- we actually only need this proof) + - so `v > s > x/v` + + 1.2. **for value too high `v < s` :** + - this is the same demonstration but in other direction, let's just summarize it : + - let's take initial value v = 4 : + ``` + v v² + 0 (4) 8 12 16 + (1) v : |---|---|---|---| + (2) x/v: |----|----|----|----| -----> stretch + (3) 0 (5) 10 15 20 + x/v x + (4) 0 *1 *2 *3 *4 -> number of portions + 012345 -> portion size + (5) (x/v)²: |----|----|----|----|----| + 0 5 10 15 20 25 + ``` + - (1) : 4 is not the sqrt of x == 20, it's too smalle : (4² == 16) < 20 + - (2) : stretch it, so the previous 4 portions fit the x = 20 size + - (3) : the new portion x/v == 5, is more than v == 4, so `v < x/v` + - (4) : portion size is bigger than number of portions, so it's too big to be the root + - (5) : indeed, we see that the portion² == (x/v)² is bigger than x, so √x < x/v == `s < x/v` + - so `v < s < x/v` diff --git a/headers/computorv1.h b/headers/computorv1.h index baebff8..2cc9e44 100644 --- a/headers/computorv1.h +++ b/headers/computorv1.h @@ -132,14 +132,18 @@ typedef struct { e_delta_sign delta_sign; // DELTA_PLUS or DELTA_MINUS or DELTA_ZERO double delta_absolute; // |Δ| == |b² - 4ac| + double delta_sqrt; // √|Δ| + // int first_term_gcd; // gcd(b, 2a) int first_term_numerator; // -b / gcd int first_term_denominator; // 2a / gcd double first_term; // double (-b / 2a) + // int second_term_gcd; // gcd(√|Δ|, 2a) int second_term_numerator; // √|Δ| / gcd int second_term_denominator; // 2a / gcd double second_term; // double (√|Δ| / 2a) + // double solution1; // first_term + second_term double solution2; // first_term - second_term (not if DELTA_ZERO) } s_solution_degree_2; diff --git a/libft b/libft index 8edf428..ae9787b 160000 --- a/libft +++ b/libft @@ -1 +1 @@ -Subproject commit 8edf428fed1e5e292b4539a0b3e1bfab1c3bd392 +Subproject commit ae9787ba6db99412ff4c6d0bde357e060f350096 diff --git a/src/solver.c b/src/solver.c index f155b88..8c2f1a2 100644 --- a/src/solver.c +++ b/src/solver.c @@ -71,8 +71,6 @@ static int find_gcd(double numerator, double denominator) /* typedef struct { - double a; // a in "ax + b" - double b; // b in "ax + b" int solution_gcd; // gcd(b, a) int solution_numerator; // -b / gcd int solution_denominator; // a / gcd @@ -82,19 +80,20 @@ typedef struct /* typedef struct { - double a; // a in "ax² + bx + c" - double b; // b in "ax² + bx + c" - double c; // c in "ax² + bx + c" e_delta_sign delta_sign; // DELTA_PLUS or DELTA_MINUS or DELTA_ZERO double delta_absolute; // |Δ| == |b² - 4ac| + double delta_sqrt; // √|Δ| + int first_term_gcd; // gcd(b, 2a) int first_term_numerator; // -b / gcd int first_term_denominator; // 2a / gcd double first_term; // double (-b / 2a) + int second_term_gcd; // gcd(√|Δ|, 2a) int second_term_numerator; // √|Δ| / gcd int second_term_denominator; // 2a / gcd double second_term; // double (√|Δ| / 2a) + double solution1; // first_term + second_term double solution2; // first_term - second_term (not if DELTA_ZERO) } s_solution_degree_2; @@ -107,11 +106,16 @@ static void solve_degree_2(s_solution_degree_2 *solution, double a, double b, do delta = b * b - 4 * a * c; solution->delta_sign = ft_fsign(delta); solution->delta_absolute = ft_fabs(delta); + solution->delta_sqrt = ft_sqrt(solution->delta_absolute); + + // first term solution->first_term_gcd = find_gcd(b, 2 * a); solution->first_term_numerator = -b / solution->first_term_gcd; solution->first_term_denominator = 2 * a / solution->first_term_gcd; solution->first_term = -b / (2 * a); - // solution->second_term_gcd = ; + + // second term + // solution->second_term_gcd = find_gcd(b, 2 * a); // solution->second_term_numerator = ; // solution->second_term_denominator = ; // solution->second_term = ; diff --git a/src/utils/errors.c b/src/utils/errors.c index 6db650d..1014b62 100644 --- a/src/utils/errors.c +++ b/src/utils/errors.c @@ -76,15 +76,16 @@ static void print_context_solution() if (solution_g_err->degree == 2) { solution_2 = solution_g_err->solution_degree_2; - dprintf(STDERR_FILENO, "degree 2 : delta > 0 ( -b/2a +- √Δ/2a )\n"); - dprintf(STDERR_FILENO, " delta == 0 ( -b/2a )\n"); + dprintf(STDERR_FILENO, "degree 2 : delta > 0 ( -b/2a +- √Δ/2a )\n"); + dprintf(STDERR_FILENO, " delta == 0 ( -b/2a )\n"); dprintf(STDERR_FILENO, " delta < 0 ( -b/2a +- i√|Δ|/2a )\n"); dprintf(STDERR_FILENO, "delta_sign : %25s\n", delta_sign_to_str(solution_2.delta_sign)); - dprintf(STDERR_FILENO, "delta_absolute : %25g (|b² - 4ac|)\n", solution_2.delta_absolute); + dprintf(STDERR_FILENO, "delta_absolute : %25g ( |Δ| == |b² - 4ac| )\n", solution_2.delta_absolute); + dprintf(STDERR_FILENO, "delta_sqrt : %25g (√|Δ| )\n", solution_2.delta_sqrt); dprintf(STDERR_FILENO, "first_term_gcd : %25i\n", solution_2.first_term_gcd); - dprintf(STDERR_FILENO, "first_term_numerator : %25i (-b / gcd)\n", solution_2.first_term_numerator); - dprintf(STDERR_FILENO, "first_term_denominator : %25i (2a / gcd)\n", solution_2.first_term_denominator); - dprintf(STDERR_FILENO, "first_term : %25g (-b / 2a)\n", solution_2.first_term); + dprintf(STDERR_FILENO, "first_term_numerator : %25i (-b / gcd )\n", solution_2.first_term_numerator); + dprintf(STDERR_FILENO, "first_term_denominator : %25i (2a / gcd )\n", solution_2.first_term_denominator); + dprintf(STDERR_FILENO, "first_term : %25g (-b / 2a )\n", solution_2.first_term); // dprintf(STDERR_FILENO, "second_term_gcd : %25g\n", solution_2.second_term_gcd); // dprintf(STDERR_FILENO, "second_term_numerator : %25g\n", solution_2.second_term_numerator); // dprintf(STDERR_FILENO, "second_term_denominator: %25g\n", solution_2.second_term_denominator);