update readme with sqrt explication
This commit is contained in:
163
README.md
163
README.md
@@ -63,19 +63,176 @@ this project uses submodules (maybe recursively), so either :
|
|||||||
-> Δ == 0 -> 1 solution : x = ( -b / 2a )
|
-> Δ == 0 -> 1 solution : x = ( -b / 2a )
|
||||||
-> Δ < 0 -> 2 solutions : x = ( -b / 2a ) +- i( √|Δ| / 2a )
|
-> Δ < 0 -> 2 solutions : x = ( -b / 2a ) +- i( √|Δ| / 2a )
|
||||||
-> solution :
|
-> solution :
|
||||||
- a; // double
|
|
||||||
- b; // double
|
|
||||||
- c; // double
|
|
||||||
- delta_sign; // + or -
|
- delta_sign; // + or -
|
||||||
- delta_absolute; // |Δ|
|
- delta_absolute; // |Δ|
|
||||||
|
- delta_sqrt; // √|Δ|
|
||||||
|
|
||||||
- first_term_gcd; // gcd(b, 2a)
|
- first_term_gcd; // gcd(b, 2a)
|
||||||
- first_term_numerator; // -b / gcd
|
- first_term_numerator; // -b / gcd
|
||||||
- first_term_denominator; // 2a / gcd
|
- first_term_denominator; // 2a / gcd
|
||||||
- first_term; // double (-b / 2a)
|
- first_term; // double (-b / 2a)
|
||||||
|
|
||||||
- second_term_gcd; // gcd(√|Δ|, 2a)
|
- second_term_gcd; // gcd(√|Δ|, 2a)
|
||||||
- second_term_numerator; // √|Δ| / gcd
|
- second_term_numerator; // √|Δ| / gcd
|
||||||
- second_term_denominator; // 2a / gcd
|
- second_term_denominator; // 2a / gcd
|
||||||
- second_term; // double (√|Δ| / 2a)
|
- second_term; // double (√|Δ| / 2a)
|
||||||
|
|
||||||
- double solution1; // first_term + second_term
|
- double solution1; // first_term + second_term
|
||||||
- double solution2; // first_term - second_term
|
- double solution2; // first_term - second_term
|
||||||
6. print solution
|
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`
|
||||||
|
|||||||
@@ -132,14 +132,18 @@ typedef struct
|
|||||||
{
|
{
|
||||||
e_delta_sign delta_sign; // DELTA_PLUS or DELTA_MINUS or DELTA_ZERO
|
e_delta_sign delta_sign; // DELTA_PLUS or DELTA_MINUS or DELTA_ZERO
|
||||||
double delta_absolute; // |Δ| == |b² - 4ac|
|
double delta_absolute; // |Δ| == |b² - 4ac|
|
||||||
|
double delta_sqrt; // √|Δ|
|
||||||
|
//
|
||||||
int first_term_gcd; // gcd(b, 2a)
|
int first_term_gcd; // gcd(b, 2a)
|
||||||
int first_term_numerator; // -b / gcd
|
int first_term_numerator; // -b / gcd
|
||||||
int first_term_denominator; // 2a / gcd
|
int first_term_denominator; // 2a / gcd
|
||||||
double first_term; // double (-b / 2a)
|
double first_term; // double (-b / 2a)
|
||||||
|
//
|
||||||
int second_term_gcd; // gcd(√|Δ|, 2a)
|
int second_term_gcd; // gcd(√|Δ|, 2a)
|
||||||
int second_term_numerator; // √|Δ| / gcd
|
int second_term_numerator; // √|Δ| / gcd
|
||||||
int second_term_denominator; // 2a / gcd
|
int second_term_denominator; // 2a / gcd
|
||||||
double second_term; // double (√|Δ| / 2a)
|
double second_term; // double (√|Δ| / 2a)
|
||||||
|
//
|
||||||
double solution1; // first_term + second_term
|
double solution1; // first_term + second_term
|
||||||
double solution2; // first_term - second_term (not if DELTA_ZERO)
|
double solution2; // first_term - second_term (not if DELTA_ZERO)
|
||||||
} s_solution_degree_2;
|
} s_solution_degree_2;
|
||||||
|
|||||||
2
libft
2
libft
Submodule libft updated: 8edf428fed...ae9787ba6d
16
src/solver.c
16
src/solver.c
@@ -71,8 +71,6 @@ static int find_gcd(double numerator, double denominator)
|
|||||||
/*
|
/*
|
||||||
typedef struct
|
typedef struct
|
||||||
{
|
{
|
||||||
double a; // a in "ax + b"
|
|
||||||
double b; // b in "ax + b"
|
|
||||||
int solution_gcd; // gcd(b, a)
|
int solution_gcd; // gcd(b, a)
|
||||||
int solution_numerator; // -b / gcd
|
int solution_numerator; // -b / gcd
|
||||||
int solution_denominator; // a / gcd
|
int solution_denominator; // a / gcd
|
||||||
@@ -82,19 +80,20 @@ typedef struct
|
|||||||
/*
|
/*
|
||||||
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
|
e_delta_sign delta_sign; // DELTA_PLUS or DELTA_MINUS or DELTA_ZERO
|
||||||
double delta_absolute; // |Δ| == |b² - 4ac|
|
double delta_absolute; // |Δ| == |b² - 4ac|
|
||||||
|
double delta_sqrt; // √|Δ|
|
||||||
|
|
||||||
int first_term_gcd; // gcd(b, 2a)
|
int first_term_gcd; // gcd(b, 2a)
|
||||||
int first_term_numerator; // -b / gcd
|
int first_term_numerator; // -b / gcd
|
||||||
int first_term_denominator; // 2a / gcd
|
int first_term_denominator; // 2a / gcd
|
||||||
double first_term; // double (-b / 2a)
|
double first_term; // double (-b / 2a)
|
||||||
|
|
||||||
int second_term_gcd; // gcd(√|Δ|, 2a)
|
int second_term_gcd; // gcd(√|Δ|, 2a)
|
||||||
int second_term_numerator; // √|Δ| / gcd
|
int second_term_numerator; // √|Δ| / gcd
|
||||||
int second_term_denominator; // 2a / gcd
|
int second_term_denominator; // 2a / gcd
|
||||||
double second_term; // double (√|Δ| / 2a)
|
double second_term; // double (√|Δ| / 2a)
|
||||||
|
|
||||||
double solution1; // first_term + second_term
|
double solution1; // first_term + second_term
|
||||||
double solution2; // first_term - second_term (not if DELTA_ZERO)
|
double solution2; // first_term - second_term (not if DELTA_ZERO)
|
||||||
} s_solution_degree_2;
|
} 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;
|
delta = b * b - 4 * a * c;
|
||||||
solution->delta_sign = ft_fsign(delta);
|
solution->delta_sign = ft_fsign(delta);
|
||||||
solution->delta_absolute = ft_fabs(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_gcd = find_gcd(b, 2 * a);
|
||||||
solution->first_term_numerator = -b / solution->first_term_gcd;
|
solution->first_term_numerator = -b / solution->first_term_gcd;
|
||||||
solution->first_term_denominator = 2 * a / solution->first_term_gcd;
|
solution->first_term_denominator = 2 * a / solution->first_term_gcd;
|
||||||
solution->first_term = -b / (2 * a);
|
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_numerator = ;
|
||||||
// solution->second_term_denominator = ;
|
// solution->second_term_denominator = ;
|
||||||
// solution->second_term = ;
|
// solution->second_term = ;
|
||||||
|
|||||||
@@ -80,11 +80,12 @@ static void print_context_solution()
|
|||||||
dprintf(STDERR_FILENO, " delta == 0 ( -b/2a )\n");
|
dprintf(STDERR_FILENO, " delta == 0 ( -b/2a )\n");
|
||||||
dprintf(STDERR_FILENO, " delta < 0 ( -b/2a +- i√|Δ|/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_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_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_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_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 : %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_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_numerator : %25g\n", solution_2.second_term_numerator);
|
||||||
// dprintf(STDERR_FILENO, "second_term_denominator: %25g\n", solution_2.second_term_denominator);
|
// dprintf(STDERR_FILENO, "second_term_denominator: %25g\n", solution_2.second_term_denominator);
|
||||||
|
|||||||
Reference in New Issue
Block a user