max/ffmpeg
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

swscale/rational64: add 64-bit rational type

Browse Source 
This is needed by the ops code, to represent intermediate values for 32-bit
formats, which can exceed the value range of int32_t (especially for
intermediate products).

I copied the math almost 1:1 from rational.c, but adapted to use the 128-bit
integer wrappers defined by int128.h. I added a generous amount of tests in
any case.

Sponsored-by: Sovereign Tech Fund
Signed-off-by: Niklas Haas <git@haasn.dev>
This commit is contained in:
Niklas Haas
2026-06-10 14:58:39 +02:00
parent 6a853ce327
commit b568f20f32
4 changed files with 417 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
libavutil/tests/rational.c
+123
View File
@@ -20,12 +20,14 @@
*/
#include "libavutil/rational.c"
#include "libswscale/rational64.c"
#include "libavutil/integer.h"
#include "libavutil/intfloat.h"
int main(void)
{
AVRational a,b,r;
AVRational64 a64,b64,r64;
int i,j,k;
static const int64_t numlist[] = {
INT64_MIN, INT64_MIN+1, INT64_MAX, INT32_MIN, INT32_MAX, 1,0,-1,
@@ -53,6 +55,69 @@ int main(void)
}
}
for (a64.num = -2; a64.num <= 2; a64.num++) {
for (a64.den = -2; a64.den <= 2; a64.den++) {
for (b64.num = -2; b64.num <= 2; b64.num++) {
for (b64.den = -2; b64.den <= 2; b64.den++) {
const double adbl = av_q2d_64(a64);
const double bdbl = av_q2d_64(b64);
const int c = av_cmp_q64(a64,b64);
const int d = adbl == bdbl ? 0 :
adbl > bdbl ? 1 :
adbl < bdbl ? -1 : INT_MIN;
if (c != d)
av_log(NULL, AV_LOG_ERROR, "%lld/%lld %lld/%lld, %d != %d\n",
(long long) a64.num, (long long) a64.den,
(long long) b64.num, (long long) b64.den, c,d);
// Check arithmetic result
if (a64.den && b64.den) {
double rdbl;
r64 = av_add_q64(a64, b64);
rdbl = av_q2d_64(r64);
if (rdbl != adbl + bdbl) {
av_log(NULL, AV_LOG_ERROR, "%f + %f = %f != %f\n",
adbl, bdbl, rdbl, adbl + bdbl);
}
r64 = av_mul_q64(a64, b64);
rdbl = av_q2d_64(r64);
if (rdbl != adbl * bdbl) {
av_log(NULL, AV_LOG_ERROR, "%f * %f = %f != %f\n",
adbl, bdbl, rdbl, adbl * bdbl);
}
}
// Check addition round-trip
r64 = av_sub_q64(av_add_q64(a64, b64), b64);
if (b64.den && (r64.num*a64.den != a64.num*r64.den ||
!r64.num != !a64.num ||
!r64.den != !a64.den))
{
av_log(NULL, AV_LOG_ERROR, "%lld/%lld != %lld/%lld\n",
(long long) a64.num, (long long) a64.den,
(long long) r64.num, (long long) r64.den);
}
if (b64.num) {
// Check multiplication round-trip
r64 = av_div_q64(av_mul_q64(a64, b64), b64);
if (b64.den && (r64.num*a64.den != a64.num*r64.den ||
!r64.num != !a64.num ||
!r64.den != !a64.den))
{
av_log(NULL, AV_LOG_ERROR, "%lld/%lld != %lld/%lld\n",
(long long) a64.num, (long long) a64.den,
(long long) r64.num, (long long) r64.den);
}
}
}
}
}
}
/* Check overflow behavior and edge cases */
static const AVRational unit_mul_q[][3] = {
{{INT_MAX, 2}, { 2, 1}, { INT_MAX, 1}},
@@ -106,6 +171,64 @@ int main(void)
}
}
static const AVRational64 unit_mul_q64[][3] = {
{{INT64_MAX, 2}, { 2, 1}, { INT64_MAX, 1}},
{{INT64_MAX, 2}, {-2, 1}, {-INT64_MAX, 1}},
{{INT64_MAX, 2}, { 0, 1}, {0, 1}},
{{INT64_MIN, 2}, { 2, 1}, {-INT64_MAX, 1}}, /* not INT64_MIN */
{{INT64_MIN, 2}, {-2, 1}, { INT64_MAX, 1}},
{{INT64_MIN, 2}, { 0, 1}, {0, 1}},
{{INT64_MAX >> 8, 1}, {INT64_MAX >> 8, 1}, {INT64_MAX, 1}},
{{1, INT64_MAX >> 8}, {1, INT64_MAX >> 8}, {0, 1}},
{{1, 1}, {0, 0}, {0, 0}},
{{0, 1}, {0, 0}, {0, 0}},
};
for (i = 0; i < FF_ARRAY_ELEMS(unit_mul_q64); i++) {
for (int c = 0; c < 2; c++) { /* test commutativity */
AVRational64 a = unit_mul_q64[i][c ? 1 : 0];
AVRational64 b = unit_mul_q64[i][c ? 0 : 1];
AVRational64 c = unit_mul_q64[i][2];
AVRational64 r = av_mul_q64(a, b);
if (r.num != c.num || r.den != c.den) {
av_log(NULL, AV_LOG_ERROR, "%lld/%lld * %lld/%lld = %lld/%lld, expected %lld/%lld\n",
(long long) a.num, (long long) a.den,
(long long) b.num, (long long) b.den,
(long long) r.num, (long long) r.den,
(long long) c.num, (long long) c.den);
}
}
}
static const AVRational64 unit_add_q64[][3] = {
{{INT64_MAX, 1}, { 2, 2}, { INT64_MAX, 1}},
{{INT64_MAX, 1}, {-2, 2}, { INT64_MAX - 1, 1}},
{{INT64_MAX, 1}, { 0, 2}, { INT64_MAX, 1}},
{{INT64_MIN, 1}, { 2, 2}, {-INT64_MAX, 1}},
{{INT64_MIN, 1}, {-2, 2}, {-INT64_MAX, 1}},
{{INT64_MIN, 1}, { 0, 2}, {-INT64_MAX, 1}},
{{INT64_MAX - 10, 1}, {20, 1}, { INT64_MAX, 1}},
{{2, INT64_MAX}, {2, INT64_MAX}, {4, INT64_MAX}},
{{1, 1}, {0, 0}, {0, 0}},
{{0, 1}, {0, 0}, {0, 0}},
};
for (i = 0; i < FF_ARRAY_ELEMS(unit_add_q64); i++) {
for (int c = 0; c < 2; c++) { /* test commutativity */
AVRational64 a = unit_add_q64[i][c ? 1 : 0];
AVRational64 b = unit_add_q64[i][c ? 0 : 1];
AVRational64 c = unit_add_q64[i][2];
AVRational64 r = av_add_q64(a, b);
if (r.num != c.num || r.den != c.den) {
av_log(NULL, AV_LOG_ERROR, "%lld/%lld + %lld/%lld = %lld/%lld, expected %lld/%lld\n",
(long long) a.num, (long long) a.den,
(long long) b.num, (long long) b.den,
(long long) r.num, (long long) r.den,
(long long) c.num, (long long) c.den);
}
}
}
for (i = 0; i < FF_ARRAY_ELEMS(numlist); i++) {
int64_t a = numlist[i];
libswscale/Makefile
+1
View File
@@ -34,6 +34,7 @@ OBJS-$(CONFIG_UNSTABLE) += \
ops_dispatch.o \
ops_memcpy.o \
ops_optimizer.o \
rational64.o \
uops.o \
uops_backend.o \
libswscale/rational64.c
+146
View File
@@ -0,0 +1,146 @@
/*
* 64-bit rational numbers
* Copyright (c) 2025 Niklas Haas
* Copyright (c) 2003 Michael Niedermayer <michaelni@gmx.at>
*
* This file is part of FFmpeg.
*
* FFmpeg is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* FFmpeg is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with FFmpeg; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
/**
* @file
* 64-bit rational numbers
* @author Niklas Haas
*/
#include <limits.h>
#include "libavutil/int128.h"
#include "rational64.h"
static av_int128 gcd128(av_int128 a, av_int128 b)
{
while (av_test128(b)) {
av_int128 tmp = b;
b = av_mod128(a, b);
a = tmp;
}
return a;
}
static AVRational64 reduce64(av_int128 num, av_int128 den)
{
const av_int128 zero = av_to128i(0);
const av_int128 max = av_to128i(INT64_MAX);
const int num_sign = av_cmp128(num, zero) < 0;
const int den_sign = av_cmp128(den, zero) < 0;
if (num_sign)
num = av_sub128(zero, num);
if (den_sign)
den = av_sub128(zero, den);
const av_int128 gcd = gcd128(num, den);
if (av_test128(gcd)) {
num = av_div128(num, gcd);
den = av_div128(den, gcd);
}
av_uint128 a0n = av_to128u(0), a0d = av_to128u(1);
av_uint128 a1n = av_to128u(1), a1d = av_to128u(0);
if (av_cmp128(num, max) <= 0 && av_cmp128(den, max) <= 0) {
a1n = num;
a1d = den;
goto done;
}
while (av_test128(den)) {
av_int128 x = av_div128(num, den);
av_int128 next_den = av_sub128(num, av_mul128(den, x));
av_uint128 a2n = av_add128(av_mul128(x, a1n), a0n);
av_uint128 a2d = av_add128(av_mul128(x, a1d), a0d);
if (av_cmp128(a2n, max) > 0 || av_cmp128(a2d, max) > 0) {
if (av_test128(a1n))
x = av_div128(av_sub128(max, a0n), a1n);
if (av_test128(a1d)) {
av_uint128 tmp = av_div128(av_sub128(max, a0d), a1d);
x = av_min128(x, tmp);
}
av_uint128 x1d = av_mul128(x, a1d);
av_uint128 a = av_mul128(den, av_add128(av_add128(x1d, x1d), a0d));
av_uint128 b = av_mul128(num, a1d);
if (av_cmp128(a, b) > 0) {
a1n = av_add128(av_mul128(x, a1n), a0n);
a1d = av_add128(x1d, a0d);
}
break;
}
a0n = a1n;
a0d = a1d;
a1n = a2n;
a1d = a2d;
num = den;
den = next_den;
}
done:;
AVRational64 res = { av_from128i(a1n), av_from128i(a1d) };
if (num_sign ^ den_sign)
res.num = -res.num;
return res;
}
int av_cmp_q64(AVRational64 a, AVRational64 b)
{
const av_int128 p = av_mul128(av_to128i(a.num), av_to128i(b.den));
const av_int128 q = av_mul128(av_to128i(b.num), av_to128i(a.den));
const int test = av_cmp128(p, q);
if (test)
return (a.den < 0) ^ (b.den < 0) ? -test : test;
else if (b.den && a.den)
return 0;
else if (a.num && b.num)
return (a.num >> 63) - (b.num >> 63);
else
return INT_MIN;
}
AVRational64 av_mul_q64(AVRational64 b, AVRational64 c)
{
return reduce64(av_mul128(av_to128i(b.num), av_to128i(c.num)),
av_mul128(av_to128i(b.den), av_to128i(c.den)));
}
AVRational64 av_div_q64(AVRational64 b, AVRational64 c)
{
return av_mul_q64(b, av_inv_q64(c));
}
AVRational64 av_add_q64(AVRational64 b, AVRational64 c) {
return reduce64(av_add128(av_mul128(av_to128i(b.num), av_to128i(c.den)),
av_mul128(av_to128i(c.num), av_to128i(b.den))),
av_mul128(av_to128i(b.den), av_to128i(c.den)));
}
AVRational64 av_sub_q64(AVRational64 b, AVRational64 c)
{
return reduce64(av_sub128(av_mul128(av_to128i(b.num), av_to128i(c.den)),
av_mul128(av_to128i(c.num), av_to128i(b.den))),
av_mul128(av_to128i(b.den), av_to128i(c.den)));
}
libswscale/rational64.h
+147
View File
@@ -0,0 +1,147 @@
/*
* 64-bit rational numbers
* Copyright (c) 2025 Niklas Haas
* Copyright (c) 2003 Michael Niedermayer <michaelni@gmx.at>
*
* This file is part of FFmpeg.
*
* FFmpeg is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* FFmpeg is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with FFmpeg; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
/**
* @file
* @ingroup lavu_math_rational
* 64-bit extension of AVRational.
* @author Niklas Haas
*/
#ifndef SWSCALE_RATIONAL64_H
#define SWSCALE_RATIONAL64_H
#include <stdint.h>
#include <limits.h>
#include "libavutil/attributes.h"
/**
* @defgroup AVRational64
* 64-bit extension of AVRational
*
* Offers a 64-bit extended version of AVRational. This is less efficient (and
* may revolve around emulated 128-bit multiplications internally), but allows
* to represent a much larger range of rational numbers without overflow.
*
* @{
*/
/**
* 64-bit Rational number (pair of numerator and denominator).
*/
typedef struct AVRational64 {
int64_t num; ///< Numerator
int64_t den; ///< Denominator
} AVRational64;
/**
* Create an AVRational64.
*
* Useful for compilers that do not support compound literals.
*
* @note The return value is not reduced.
*/
static inline AVRational64 av_make_q64(int64_t num, int64_t den)
{
AVRational64 r = { num, den };
return r;
}
/**
* Compare two 64-bit rationals.
*
* @param a First rational
* @param b Second rational
*
* @return One of the following values:
* - 0 if `a == b`
* - 1 if `a > b`
* - -1 if `a < b`
* - `INT_MIN` if one of the values is of the form `0 / 0`
*/
int av_cmp_q64(AVRational64 a, AVRational64 b);
/**
* Convert an AVRational64 to a `double`.
* @param a AVRational64 to convert
* @return `a` in floating-point form
* @see av_d2q()
*/
static inline double av_q2d_64(AVRational64 a){
return a.num / (double) a.den;
}
/**
* Multiply two 64-bit rationals.
* @param b First multiplicant
* @param c Second multiplicant
* @return b*c
*/
AVRational64 av_mul_q64(AVRational64 b, AVRational64 c) av_const;
/**
* Divide one 64-bit rational by another.
* @param b Dividend
* @param c Divisor
* @return b/c
*/
AVRational64 av_div_q64(AVRational64 b, AVRational64 c) av_const;
/**
* Add two 64-bit rationals.
* @param b First addend
* @param c Second addend
* @return b+c
*/
AVRational64 av_add_q64(AVRational64 b, AVRational64 c) av_const;
/**
* Subtract one 64-bit rational from another.
* @param b Minuend
* @param c Subtrahend
* @return b-c
*/
AVRational64 av_sub_q64(AVRational64 b, AVRational64 c) av_const;
/**
* Invert a 64-bit rational.
* @param q value
* @return 1 / q
*/
static av_always_inline AVRational64 av_inv_q64(AVRational64 q)
{
AVRational64 r = { q.den, q.num };
return r;
}
/**
* Return the best rational so that a and b are multiple of it.
* If the resulting denominator is larger than max_den, return def.
*/
AVRational64 av_gcd_q64(AVRational64 a, AVRational64 b, int max_den, AVRational64 def);
/**
* @}
*/
#endif /* SWSCALE_RATIONAL64_H */