r/C_Programming • u/onecable5781 • 1d ago
Converting fractions to integers
I have a double* which has the following entries:
0.3333333333333334
0.6666666666666668
0.1249999999999999
1
Here, the last entry, 1, can be considered the right hand side of an inequality:
0.3333333333333334 x + 0.6666666666666668 y + 0.1249999999999999 z >= 1
These numbers come from a numerical linear algebra library over which I don't have any control. What is the easiest way to "convert" this to the following equivalent inequality (subject to a user provided tolerance of what counts as an epsilon so that epsilon within an integer is to be counted as an integer)?
8 x + 16 y + 3 z >= 24
Is there a package that does such conversion, if reasonably possible? I consider it unreasonably possible by multiplying everything in the original equation by a large enough power of 10. But I do not want that.
4
u/MyTinyHappyPlace 1d ago
It’s tough that you don’t get the actual fractions but their nearest double values. Otherwise the GNU MP (https://gmplib.org/manual/Rational-Number-Functions) library could carry you a bit on your way
2
u/Educational-Paper-75 1d ago
By brute force multiplication of your 1/3, 2/3 and 1/8 by successive integers (2, 3, ...) until all are within epsilon of their floor values.
1
u/Forever_DM5 1d ago
Multiply by 10^n where n is the number of significant digits you want(your tolerance).
Apply Euclid’s algorithm to find the LCD then multiply the original equation by that.
1
u/dstroy0 23h ago
/**
* Rounds the 128-bit significand to a 64-bit integer, a half going up.
*
* [in] args The significand and its exponent [BORROWS].
* The rounded integer, 0 when the value rounds below one, or all ones when it needs
* more than 64 bits.
* k is the negated exponent, so it says how far right the significand must move to become an integer.
* The three branches split on word_shift, the part past 64: none of it, under 64 more, and
* 64 or more, in that order.
* Only the halfway bit decides. Each branch reads a shift and a mask and nothing else, which is
* less for the compiler to carry than gathering the bits below the half and then testing them.
* u/note Reads neither args->rest nor args->above. A half is taken up whatever either holds.
*/
EMBED_INLINE embed_u64 muto_to_u64(const MutoCtx *args)
{
if ((args->hi | args->lo) == 0u)
{
return 0u;
}
const embed_iword total_shift = -(args->fe2);
if (total_shift > 128)
{
return 0u;
}
if (total_shift < 64)
{
// Explicit cast pins the all-ones saturation value at the embed_u64 this returns
return ~(embed_u64)0;
}
// Explicit cast holds the shift width at the embed_word the shifts below take. The tests above
// established total_shift is between 64 and 128, so the subtraction cannot wrap
const embed_word word_shift = (embed_word)(total_shift - 64);
embed_u64 whole;
embed_u64 half;
if (word_shift == 0u)
{
whole = args->hi;
half = (args->lo >> 63) & 1u;
}
else if (word_shift < 64u)
{
whole = args->hi >> word_shift;
half = (args->hi >> (word_shift - 1u)) & 1u;
}
else
{
whole = 0u;
half = args->hi >> 63;
}
// A half rounds up. Nothing below the halfway bit can change that: a value past the half is going
// up already, and one under it has no half bit to find. So the bits below are neither gathered nor
// tested, and each arm is left at a shift and a mask the compiler can take an early out of.
// Measured against the tie-to-even form this replaced, at no decimals: 140 cycles to 81 on an
// ESP32-S3 and 159 to 107 on an ESP32-C6. At six decimals, 378 to 322 and 413 to 368.
if (half != 0u)
{
whole += 1u;
}
return whole;
}
/**
* Turns *args->mant times ten raised to args->ex into a rounded 64-bit integer.
*
* [in,out] args The mantissa, its binary exponent, the decimal exponent and the tie bias [BORROWS].
* The rounded integer, or 0 when *args->mant is zero.
* Always takes the 128-bit path, with no exact double shortcut, unlike muto_scale.
* Reads args->e2, which muto_scale leaves at zero.
* u/warning Does not bound args->ex against MMGR_POW5_MAX the way muto_scale does, so a larger one loses its high bits.
*/
EMBED_INLINE embed_u64 muto_scale_to_u64(MutoCtx *args)
{
if (*args->mant == 0u)
{
return 0u;
}
muto_seat(args);
muto_apply_pow10(args);
return muto_to_u64(args);
}
/**
* Loads *args->mant into the 128-bit significand and normalizes it.
*
* [in,out] args The mantissa, its binary exponent, and the bits already dropped [BORROWS].
* Puts the mantissa in args->hi with args->lo zero, so args->fe2 starts 64 below args->e2.
* u/note Carries args->dropped into args->rest, so a truncation the caller already made still reaches the rounding.
*/
EMBED_INLINE void muto_seat(MutoCtx *args)
{
args->hi = *args->mant;
args->lo = 0u;
args->fe2 = args->e2 - 64;
args->rest = args->dropped;
muto_norm(args);
}
EMBED_INLINE void muto_apply_pow10(MutoCtx *args)
{
// A positive exponent is applied as exact powers of ten rather than by walking the bits of the
// wide tables. Ten to the eighteenth is the widest that fits 64 bits, so a larger exponent goes
// on in chunks of it. Each chunk is a 128 by 64 multiply, where every set bit of the walk is a
// 128 by 128 one, and the walk needs one per bit rather than one per eighteen.
if (args->ex > 0)
{
// Below the table's reach one power covers the whole exponent, so it needs no loop
if (args->ex <= MMGR_MUTO_EXACT_U64_POW10)
{
args->right = mmgr_muto_pow10[args->ex];
muto_mul_pow10(args);
return;
}
embed_iword left = args->ex;
while (left > 0)
{
const embed_iword take = (left > MMGR_MUTO_EXACT_U64_POW10) ? MMGR_MUTO_EXACT_U64_POW10 : left;
args->right = mmgr_muto_pow10[take];
muto_mul_pow10(args);
left -= take;
}
return;
}
// Explicit cast holds the negated exponent at the embed_iword the walk counts down in
embed_iword magnitude = (args->ex < 0) ? (embed_iword)(-args->ex) : args->ex;
// Two bounds: the step count keeps an exponent past the tables from reading off the end, and the
// emptied magnitude ends the walk once no bits are left to apply
for (embed_iword step = 0; (step < MMGR_POW5_STEPS) && (magnitude != 0); ++step)
{
if ((magnitude & 1) != 0)
{
args->pow = (args->ex < 0) ? &mmgr_pow5_down[step] : &mmgr_pow5_up[step];
muto_mul_pow5(args);
}
magnitude >>= 1;
}
args->fe2 += args->ex;
}
/**
* Ten raised to 0 through 18, as exact 64-bit integers.
*
* muto_apply_pow10 applies these for a positive decimal exponent, one outright when the
* exponent is within the table and otherwise in chunks of the largest, rather than walking
* the pow5 tables a set bit at a time.
*/
static const embed_u64 mmgr_muto_pow10[MMGR_MUTO_EXACT_U64_POW10 + 1] = {1ull,
10ull,
100ull,
1000ull,
10000ull,
100000ull,
1000000ull,
10000000ull,
100000000ull,
1000000000ull,
10000000000ull,
100000000000ull,
1000000000000ull,
10000000000000ull,
100000000000000ull,
1000000000000000ull,
10000000000000000ull,
100000000000000000ull,
1000000000000000000ull};
/**
* Multiplies the 128-bit significand by one exact 64-bit power of ten, then renormalizes.
*
* [in,out] args The significand, its exponent, and the power to apply as args->right [BORROWS].
* Two multiplies where muto_mul_pow5 takes four, and one column sum where it takes three.
* The power is applied whole, both its five and its two, so nothing is added to args->fe2 for it.
* Sets args->rest from the discarded low 64 bits, as muto_mul_pow5 does from its columns.
* u/note Adds 64 to args->fe2, for the bits the 192-bit product was taken down by.
*/
EMBED_INLINE void muto_mul_pow10(MutoCtx *args)
{
const embed_u64 fhi = args->hi;
const embed_u64 flo = args->lo;
const embed_u64 power = args->right;
args->left = fhi;
args->right = power;
muto_mul(args);
const embed_u64 hh_h = args->phi;
const embed_u64 hh_l = args->plo;
args->left = flo;
args->right = power;
muto_mul(args);
const embed_u64 lh_h = args->phi;
const embed_u64 lh_l = args->plo;
const embed_u64 col = hh_l + lh_h;
const embed_u64 carry = (col < hh_l) ? 1u : 0u;
if (lh_l != 0u)
{
args->rest = 1;
}
args->hi = hh_h + carry;
args->lo = col;
// Explicit cast holds the summed exponent at the embed_iword fe2 carries
args->fe2 = (embed_iword)(args->fe2 + 64);
muto_norm(args);
}
Sorry for pasting a ton of stuff, but this is a fun problem, I solve it this way, please ignore my dispatch and entry macros. What you want is to follow muto scale to u64, that will give you the fast pow10 we know how large and small these numbers can be represented path, which is exact from smallest normal to largest normal. (-1022 to +1023) someone else mentioned using pow10 already, this is taken to the extreme using language semantics forcing compiler behavior for opt, you don't need to reproduce this exactly for it to work the way you want, but this is one of the fastest methods available to perform the mutation/transform.
1
u/dstroy0 23h ago
It's from here if you wanted to read the whole thing: https://github.com/dstroy0/MMgr/tree/main/src/transformo
it's meant as a replacement for libc/newlib <string> and mem functions plus some other transforms for embedded targets but is in general way faster than what is out there because of the primitives used, per target/function disassembly+optimization and the order of the conditionals, and it requires aligned memory to be constructed for it to assume dumb load/store by default (machine word width at a time) there's also a full test/ suite that uses outside sources as the oracles.
13
u/aioeu 1d ago edited 23h ago
I would convert each of the coefficients to a rational number by using continued fractions. You can choose what accuracy you would like by cutting off this process just before the denominator exceeds a certain threshold. A continued fraction will always yield "best" rational approximations to a real number. With these coefficients, you're going to hit 1/3, 2/3 and 1/8 pretty quickly, with the rational approximations following these all having very large denominators.
Once you've got rational approximations, it's only a small amount of extra work to find the LCM of the denominators so you can convert everything to integers.