mathematics.c 5.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211
  1. /*
  2. * Copyright (c) 2005-2012 Michael Niedermayer <michaelni@gmx.at>
  3. *
  4. * This file is part of FFmpeg.
  5. *
  6. * FFmpeg is free software; you can redistribute it and/or
  7. * modify it under the terms of the GNU Lesser General Public
  8. * License as published by the Free Software Foundation; either
  9. * version 2.1 of the License, or (at your option) any later version.
  10. *
  11. * FFmpeg is distributed in the hope that it will be useful,
  12. * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  14. * Lesser General Public License for more details.
  15. *
  16. * You should have received a copy of the GNU Lesser General Public
  17. * License along with FFmpeg; if not, write to the Free Software
  18. * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  19. */
  20. /**
  21. * @file
  22. * miscellaneous math routines and tables
  23. */
  24. #include <stdint.h>
  25. #include <limits.h>
  26. #include "mathematics.h"
  27. #include "libavutil/intmath.h"
  28. #include "libavutil/common.h"
  29. #include "avassert.h"
  30. #include "version.h"
  31. /* Stein's binary GCD algorithm:
  32. * https://en.wikipedia.org/wiki/Binary_GCD_algorithm */
  33. int64_t av_gcd(int64_t a, int64_t b) {
  34. int za, zb, k;
  35. int64_t u, v;
  36. if (a == 0)
  37. return b;
  38. if (b == 0)
  39. return a;
  40. za = ff_ctzll(a);
  41. zb = ff_ctzll(b);
  42. k = FFMIN(za, zb);
  43. u = llabs(a >> za);
  44. v = llabs(b >> zb);
  45. while (u != v) {
  46. if (u > v)
  47. FFSWAP(int64_t, v, u);
  48. v -= u;
  49. v >>= ff_ctzll(v);
  50. }
  51. return (uint64_t)u << k;
  52. }
  53. int64_t av_rescale_rnd(int64_t a, int64_t b, int64_t c, enum AVRounding rnd)
  54. {
  55. int64_t r = 0;
  56. av_assert2(c > 0);
  57. av_assert2(b >=0);
  58. av_assert2((unsigned)(rnd&~AV_ROUND_PASS_MINMAX)<=5 && (rnd&~AV_ROUND_PASS_MINMAX)!=4);
  59. if (c <= 0 || b < 0 || !((unsigned)(rnd&~AV_ROUND_PASS_MINMAX)<=5 && (rnd&~AV_ROUND_PASS_MINMAX)!=4))
  60. return INT64_MIN;
  61. if (rnd & AV_ROUND_PASS_MINMAX) {
  62. if (a == INT64_MIN || a == INT64_MAX)
  63. return a;
  64. rnd -= AV_ROUND_PASS_MINMAX;
  65. }
  66. if (a < 0)
  67. return -(uint64_t)av_rescale_rnd(-FFMAX(a, -INT64_MAX), b, c, rnd ^ ((rnd >> 1) & 1));
  68. if (rnd == AV_ROUND_NEAR_INF)
  69. r = c / 2;
  70. else if (rnd & 1)
  71. r = c - 1;
  72. if (b <= INT_MAX && c <= INT_MAX) {
  73. if (a <= INT_MAX)
  74. return (a * b + r) / c;
  75. else {
  76. int64_t ad = a / c;
  77. int64_t a2 = (a % c * b + r) / c;
  78. if (ad >= INT32_MAX && b && ad > (INT64_MAX - a2) / b)
  79. return INT64_MIN;
  80. return ad * b + a2;
  81. }
  82. } else {
  83. #if 1
  84. uint64_t a0 = a & 0xFFFFFFFF;
  85. uint64_t a1 = a >> 32;
  86. uint64_t b0 = b & 0xFFFFFFFF;
  87. uint64_t b1 = b >> 32;
  88. uint64_t t1 = a0 * b1 + a1 * b0;
  89. uint64_t t1a = t1 << 32;
  90. int i;
  91. a0 = a0 * b0 + t1a;
  92. a1 = a1 * b1 + (t1 >> 32) + (a0 < t1a);
  93. a0 += r;
  94. a1 += a0 < r;
  95. for (i = 63; i >= 0; i--) {
  96. a1 += a1 + ((a0 >> i) & 1);
  97. t1 += t1;
  98. if (c <= a1) {
  99. a1 -= c;
  100. t1++;
  101. }
  102. }
  103. if (t1 > INT64_MAX)
  104. return INT64_MIN;
  105. return t1;
  106. #else
  107. /* reference code doing (a*b + r) / c, requires libavutil/integer.h */
  108. AVInteger ai;
  109. ai = av_mul_i(av_int2i(a), av_int2i(b));
  110. ai = av_add_i(ai, av_int2i(r));
  111. return av_i2int(av_div_i(ai, av_int2i(c)));
  112. #endif
  113. }
  114. }
  115. int64_t av_rescale(int64_t a, int64_t b, int64_t c)
  116. {
  117. return av_rescale_rnd(a, b, c, AV_ROUND_NEAR_INF);
  118. }
  119. int64_t av_rescale_q_rnd(int64_t a, AVRational bq, AVRational cq,
  120. enum AVRounding rnd)
  121. {
  122. int64_t b = bq.num * (int64_t)cq.den;
  123. int64_t c = cq.num * (int64_t)bq.den;
  124. return av_rescale_rnd(a, b, c, rnd);
  125. }
  126. int64_t av_rescale_q(int64_t a, AVRational bq, AVRational cq)
  127. {
  128. return av_rescale_q_rnd(a, bq, cq, AV_ROUND_NEAR_INF);
  129. }
  130. int av_compare_ts(int64_t ts_a, AVRational tb_a, int64_t ts_b, AVRational tb_b)
  131. {
  132. int64_t a = tb_a.num * (int64_t)tb_b.den;
  133. int64_t b = tb_b.num * (int64_t)tb_a.den;
  134. if ((FFABS(ts_a)|a|FFABS(ts_b)|b) <= INT_MAX)
  135. return (ts_a*a > ts_b*b) - (ts_a*a < ts_b*b);
  136. if (av_rescale_rnd(ts_a, a, b, AV_ROUND_DOWN) < ts_b)
  137. return -1;
  138. if (av_rescale_rnd(ts_b, b, a, AV_ROUND_DOWN) < ts_a)
  139. return 1;
  140. return 0;
  141. }
  142. int64_t av_compare_mod(uint64_t a, uint64_t b, uint64_t mod)
  143. {
  144. int64_t c = (a - b) & (mod - 1);
  145. if (c > (mod >> 1))
  146. c -= mod;
  147. return c;
  148. }
  149. int64_t av_rescale_delta(AVRational in_tb, int64_t in_ts, AVRational fs_tb, int duration, int64_t *last, AVRational out_tb){
  150. int64_t a, b, this;
  151. av_assert0(in_ts != AV_NOPTS_VALUE);
  152. av_assert0(duration >= 0);
  153. if (*last == AV_NOPTS_VALUE || !duration || in_tb.num*(int64_t)out_tb.den <= out_tb.num*(int64_t)in_tb.den) {
  154. simple_round:
  155. *last = av_rescale_q(in_ts, in_tb, fs_tb) + duration;
  156. return av_rescale_q(in_ts, in_tb, out_tb);
  157. }
  158. a = av_rescale_q_rnd(2*in_ts-1, in_tb, fs_tb, AV_ROUND_DOWN) >>1;
  159. b = (av_rescale_q_rnd(2*in_ts+1, in_tb, fs_tb, AV_ROUND_UP )+1)>>1;
  160. if (*last < 2*a - b || *last > 2*b - a)
  161. goto simple_round;
  162. this = av_clip64(*last, a, b);
  163. *last = this + duration;
  164. return av_rescale_q(this, fs_tb, out_tb);
  165. }
  166. int64_t av_add_stable(AVRational ts_tb, int64_t ts, AVRational inc_tb, int64_t inc)
  167. {
  168. int64_t m, d;
  169. if (inc != 1)
  170. inc_tb = av_mul_q(inc_tb, (AVRational) {inc, 1});
  171. m = inc_tb.num * (int64_t)ts_tb.den;
  172. d = inc_tb.den * (int64_t)ts_tb.num;
  173. if (m % d == 0)
  174. return ts + m / d;
  175. if (m < d)
  176. return ts;
  177. {
  178. int64_t old = av_rescale_q(ts, ts_tb, inc_tb);
  179. int64_t old_ts = av_rescale_q(old, inc_tb, ts_tb);
  180. return av_rescale_q(old + 1, inc_tb, ts_tb) + (ts - old_ts);
  181. }
  182. }