mirror of
https://github.com/jankae/LibreVNA.git
synced 2025-12-06 07:12:10 +01:00
90 lines
1.9 KiB
C++
90 lines
1.9 KiB
C++
#include "algorithm.hpp"
|
|
|
|
#include "stm.hpp"
|
|
#include <cmath>
|
|
#include <algorithm>
|
|
|
|
Algorithm::RationalApproximation Algorithm::BestRationalApproximation(float ratio, uint32_t max_denom) {
|
|
RationalApproximation result;
|
|
uint32_t a = 0, b = 1, c = 1, d = 1;
|
|
while (b + d <= max_denom) {
|
|
auto mediant = (float) (a + c) / (b + d);
|
|
if (ratio == mediant) {
|
|
if (b + d <= max_denom) {
|
|
result.num = a + c;
|
|
result.denom = b + d;
|
|
return result;
|
|
} else if (d > b) {
|
|
result.num = c;
|
|
result.denom = d;
|
|
return result;
|
|
} else {
|
|
result.num = a;
|
|
result.denom = b;
|
|
return result;
|
|
}
|
|
} else if (ratio > mediant) {
|
|
a = a + c;
|
|
b = b + d;
|
|
} else {
|
|
c = a + c;
|
|
d = b + d;
|
|
}
|
|
}
|
|
// check which of the two is the better solution
|
|
float dev_ab = (float) a / b - ratio;
|
|
float dev_cd = (float) c / d - ratio;
|
|
if(fabs(dev_cd) < fabs(dev_ab)) {
|
|
result.num = c;
|
|
result.denom = d;
|
|
} else {
|
|
result.num = a;
|
|
result.denom = b;
|
|
}
|
|
return result;
|
|
}
|
|
|
|
uint32_t gcd(uint32_t u, uint32_t v) {
|
|
if(u==0) {
|
|
return v;
|
|
} else if(v==0) {
|
|
return u;
|
|
}
|
|
|
|
uint8_t i = __builtin_ctz(u);
|
|
u >>= i;
|
|
uint8_t j = __builtin_ctz(v);
|
|
v >>= j;
|
|
|
|
uint8_t k = i < j ? i : j;
|
|
while(true) {
|
|
if(u > v) {
|
|
std::swap(u, v);
|
|
}
|
|
v -= u;
|
|
if(v==0) {
|
|
return u << k;
|
|
}
|
|
v >>= __builtin_ctz(v);
|
|
}
|
|
}
|
|
|
|
Algorithm::RationalApproximation Algorithm::BestRationalApproximation(
|
|
RationalApproximation ratio, uint32_t max_denom) {
|
|
if(ratio.denom <= max_denom) {
|
|
// nothing to do, we can just return the ratio as it is
|
|
return ratio;
|
|
}
|
|
// Try to simplify the ratio.
|
|
// Find greatest common divider
|
|
uint32_t div = gcd(ratio.num, ratio.denom);
|
|
ratio.num /= div;
|
|
ratio.denom /= div;
|
|
if(ratio.denom <= max_denom) {
|
|
// small enough now, can use as is
|
|
return ratio;
|
|
}
|
|
// no good, we need to approximate
|
|
return Algorithm::BestRationalApproximation((float) ratio.num / ratio.denom, max_denom);
|
|
}
|