diff options
author | Michael Niedermayer <michaelni@gmx.at> | 2012-02-22 01:30:35 +0400 |
---|---|---|
committer | Michael Niedermayer <michaelni@gmx.at> | 2012-02-27 03:00:55 +0400 |
commit | 59affed23c3a30b85c006c45f9e5d3aa8ce50741 (patch) | |
tree | 67175bea0de833a2d1055b0041438e9d72e3024e | |
parent | 8c5ee45d13b7a0d14ff75b7db72d5400b7750571 (diff) |
eval: add root() to solve f(x)=0
Signed-off-by: Michael Niedermayer <michaelni@gmx.at>
-rw-r--r-- | doc/eval.texi | 4 | ||||
-rw-r--r-- | libavutil/eval.c | 48 | ||||
-rw-r--r-- | tests/ref/fate/eval | 6 |
3 files changed, 57 insertions, 1 deletions
diff --git a/doc/eval.texi b/doc/eval.texi index c4a8ac564f..c4b669fbcb 100644 --- a/doc/eval.texi +++ b/doc/eval.texi @@ -114,6 +114,10 @@ then 0 is assumed. note, when you have the derivatives at y instead of 0 taylor(expr, x-y) can be used When the series does not converge the results are undefined. + +@item root(expr, max) +Finds x where f(x)=0 in the interval 0..max. +f() must be continuous or the result is undefined. @end table The following constants are available: diff --git a/libavutil/eval.c b/libavutil/eval.c index 5147f7f3f2..7c3e182509 100644 --- a/libavutil/eval.c +++ b/libavutil/eval.c @@ -26,6 +26,7 @@ * see http://joe.hotchkiss.com/programming/eval/eval.html */ +#include <float.h> #include "avutil.h" #include "eval.h" #include "log.h" @@ -134,7 +135,7 @@ struct AVExpr { e_squish, e_gauss, e_ld, e_isnan, e_mod, e_max, e_min, e_eq, e_gt, e_gte, e_pow, e_mul, e_div, e_add, - e_last, e_st, e_while, e_taylor, e_floor, e_ceil, e_trunc, + e_last, e_st, e_while, e_taylor, e_root, e_floor, e_ceil, e_trunc, e_sqrt, e_not, e_random, e_hypot, e_gcd, e_if, e_ifnot, } type; @@ -199,6 +200,48 @@ static double eval_expr(Parser *p, AVExpr *e) p->var[id] = var0; return d; } + case e_root: { + int i; + double low = -1, high = -1, v, low_v = -DBL_MAX, high_v = DBL_MAX; + double var0 = p->var[0]; + double x_max = eval_expr(p, e->param[1]); + for(i=-1; i<1024; i++) { + if(i<255) { + p->var[0] = av_reverse[i&255]*x_max/255; + } else { + p->var[0] = x_max*pow(0.9, i-255); + if (i&1) p->var[0] *= -1; + if (i&2) p->var[0] += low; + else p->var[0] += high; + } + v = eval_expr(p, e->param[0]); + if (v<=0 && v>low_v) { + low = p->var[0]; + low_v = v; + } + if (v>=0 && v<high_v) { + high = p->var[0]; + high_v = v; + } + if (low>=0 && high>=0){ + while (1) { + p->var[0] = (low+high)*0.5; + if (low == p->var[0] || high == p->var[0]) + break; + v = eval_expr(p, e->param[0]); + if (v<=0) low = p->var[0]; + if (v>=0) high= p->var[0]; + if (isnan(v)) { + low = high = v; + break; + } + } + break; + } + } + p->var[0] = var0; + return -low_v<high_v ? low : high; + } default: { double d = eval_expr(p, e->param[0]); double d2 = eval_expr(p, e->param[1]); @@ -342,6 +385,7 @@ static int parse_primary(AVExpr **e, Parser *p) else if (strmatch(next, "st" )) d->type = e_st; else if (strmatch(next, "while" )) d->type = e_while; else if (strmatch(next, "taylor")) d->type = e_taylor; + else if (strmatch(next, "root" )) d->type = e_root; else if (strmatch(next, "floor" )) d->type = e_floor; else if (strmatch(next, "ceil" )) d->type = e_ceil; else if (strmatch(next, "trunc" )) d->type = e_trunc; @@ -727,6 +771,8 @@ int main(int argc, char **argv) "ifnot(1, NaN) + if(0, 1)", "taylor(1, 1)", "taylor(eq(mod(ld(1),4),1)-eq(mod(ld(1),4),3), PI/2, 1)", + "root(sin(ld(0))-1, 2)", + "root(sin(ld(0))+6+sin(ld(0)/12)-log(ld(0)), 100)", NULL }; diff --git a/tests/ref/fate/eval b/tests/ref/fate/eval index ed0d8bc124..67b107a930 100644 --- a/tests/ref/fate/eval +++ b/tests/ref/fate/eval @@ -175,5 +175,11 @@ Evaluating 'taylor(1, 1)' Evaluating 'taylor(eq(mod(ld(1),4),1)-eq(mod(ld(1),4),3), PI/2, 1)' 'taylor(eq(mod(ld(1),4),1)-eq(mod(ld(1),4),3), PI/2, 1)' -> 1.000000 +Evaluating 'root(sin(ld(0))-1, 2)' +'root(sin(ld(0))-1, 2)' -> 1.570796 + +Evaluating 'root(sin(ld(0))+6+sin(ld(0)/12)-log(ld(0)), 100)' +'root(sin(ld(0))+6+sin(ld(0)/12)-log(ld(0)), 100)' -> 60.965601 + 12.700000 == 12.7 0.931323 == 0.931322575 |