View Full Version : avs devs: plz integrate nicefps() into changefps()/assumefps()
bond
27th December 2005, 11:50
changefps/assumefps() allow to input float fps and people make extensively use of this
now the problem is that it doesnt seem to output for eg 23.976 the standard 24000/1001 numerator,denominator but some other high values
this now leads to that for example the x264 encoder is forced to use 64bit precision times, which again causes players like mplayer or quicktime not handling these 64bit, but expecting 32bit, to bork
now we had lots of people reporting problems in the mpeg-4 avc forum reporting these problems and they could be tracked down to the useage of change/assumefps() with float fps input or with input .avi files that have been created from such scripts
now there are three solutions for fixing this:
1) make the whole world use change/assumefps(24000/1001) instead of changefps(23.976) or similar for other framerates
2) make the whole world use the nicefps() filter, as discussed here (http://forum.doom9.org/showthread.php?p=717181#post717181), which scales the framerate in a "nice" way
3) change the behaviour of change/assumefps() to not output these big values anymore by maybe integrating the nicefps()'s function into change/assumefps() or making change/assumefps process the string itself instead of using float math
now i would be really happy if 3) could be done (as 1) and 2) is surely impossible), cause it would remove like "50%" of the x264 error reports ;)
stickboy
27th December 2005, 12:50
If only AviSynth used doubles instead of floats, then it more easily could rationalize 23.976 to 2997/125.
Overloading AssumeFPS to take a string might be okay (please, don't add any more *FPS functions), but then you still have the problem of educating people to use it.
Bidoche
27th December 2005, 13:49
If only AviSynth used doubles instead of floats, then it more easily could rationalize 23.976 to 2997/125.I don't see how that would help...
Users generally mean 24000/1001 when they put 23.976 and not 2997/125.
Besides FPS is rational, not float or even double, every api with some sense use it like that.
VFW certainly does, avs should do the same (3.0 does)
but then you still have the problem of educating people to use itEasy : add overloaded *FPS functions and change all the old ones to throw a message reporting them deprecated and asking to use the new ones.
tritical
27th December 2005, 15:43
If only AviSynth used doubles instead of floats, then it more easily could rationalize 23.976 to 2997/125. I think the main thing is that because it uses floats one has to realize that the precision of the input value is limited to roughly 7 decimal digits. However, that is still enough for values like 23.976 or 29.970. For example, if the current code in assumefps, which is:double n = args[1].AsFloat();
int d = 1;
while (n < 16777216 && d < 16777216) { n*=2; d*=2; }
return new AssumeFPS(args[0].AsClip(), int(n+0.5), d, args[2].AsBool(false), env);were changed to: double n = args[1].AsFloat();
int d = 1;
while (n < 1000000.0 && d < 1000000) { n*=10.0; d*=10; }
int x = int(n+0.5), y = d, t;
while (y) { t = x%y; x = y; y = t; }
return new AssumeFPS(args[0].AsClip(), int(n+0.5)/x, d/x, args[2].AsBool(false), env);then Assumefps(23.976) would result in 2997/125, Assumefps(29.970) would result in 2997/100, Assumefps(119.88) would result in 2997/25, etc... As opposed to what you currently get which is 12570329/524288, 15712911/524288, and 982057/8192. The problem with 30000/1001 or 24000/1001 is that those are repeating decimals and cannot be accurately expressed using a limited precision decimal number. Thus, I don't see any way avisynth could deduce 24000/1001 from a decimal input. To me it seems that any rational fps that requires more than 7 decimal digits of precision (or 15 if it was changed to double) would need to be set using assumefps(numerator,denominator).
sh0dan
27th December 2005, 16:27
Can anyone predict problems with the patch above, otherwise it'll go into 2.6.
stickboy
27th December 2005, 20:08
I don't see how that would help...
Users generally mean 24000/1001 when they put 23.976 and not 2997/125.Well, a conventional algorithm of repeatedly adding the numerators and denominators of low/high bounds rationalizes 23.976 fine for doubles but not for floats.
And rationalizing to 2997/125 is still an improvement over 12570329/524288 that AviSynth currently generates, yes? Would 2997/125 not address bond's "high-values lead to unnecessary 64-bit computation" issue?Easy : add overloaded *FPS functions and change all the old ones to throw a message reporting them deprecated and asking to use the new ones.Sure, if you want to break backwards compatibility.
foxyshadis
27th December 2005, 21:36
The drift between 2997/125 and 24000/1001 is only one-hundredth of a second after three hours. If you never told anyone they'd probably never find out.
Fizick
28th December 2005, 02:07
IMO, assumefps(numerator, denominator) is most correct way.
So, if somebody use it explicitly, Avisynth must not change these values.
Mug Funky
28th December 2005, 02:54
another possibility is to simply pick the /1001 framerate that fits the input the best.
AFAIK the only rational framerates are those of (colour) NTSC. when someone types 29.97 it is plainly obvious that they mean 30000/1001, so really this number can be assumed without any drawbacks (as stickboy says, who's going to notice 1/100th of a second after 3 hours? considering that 1-frame desync is quite common on DVDs anyway, particularly if they've been standards converted, and nobody seems to notice that).
presets would be interesting too - something like "assumefps("film")" to give 24000/1001, or "assumefps("NTSC") gives 30000/1001, etc. this would certainly be easier to type (and understand...) than 30000/1001 etc.
Bidoche
28th December 2005, 09:44
Easy : add overloaded *FPS functions and change all the old ones to throw a message reporting them deprecated and asking to use the new ones.
Sure, if you want to break backwards compatibility.I was answering on how educating users.
That method certainly works
presets would be interesting too - something like "assumefps("film")" to give 24000/1001, or "assumefps("NTSC") gives 30000/1001, etc. this would certainly be easier to type (and understand...) than 30000/1001 etc.Looks interesting
d'Oursse
28th December 2005, 17:15
afaik, 24000/1001 is not the fps of the "films". it's the fps used to get the NTSC framerate from a "photographic material produced for the cinema". the framerate of the films is 24 fps (http://www.doom9.org/ivtc-tut.htm). Maybe the word "film" is not correct.
tritical
28th December 2005, 17:38
another possibility is to simply pick the /1001 framerate that fits the input the best.
AFAIK the only rational framerates are those of (colour) NTSC. when someone types 29.97 it is plainly obvious that they mean 30000/1001, so really this number can be assumed without any drawbacks (as stickboy says, who's going to notice 1/100th of a second after 3 hours? considering that 1-frame desync is quite common on DVDs anyway, particularly if they've been standards converted, and nobody seems to notice that).
While this would work for getting the 24000/1001 and 30000/1001 ratios, I personally don't think this is a good idea. When someone uses assumefps(23.976) they should, IMHO, get what they asked for. Adding work arounds to have it adjust to another value is just going to add to the confusion when someone really does want 23.976 and not 24000/1001 for whatever reason.
presets would be interesting too - something like "assumefps("film")" to give 24000/1001, or "assumefps("NTSC") gives 30000/1001, etc. this would certainly be easier to type (and understand...) than 30000/1001 etc. Sounds like a good idea.
bond
31st December 2005, 12:28
so will my proposal be followed? :)
tritical
31st December 2005, 20:51
If anyone is interested in testing, I have merged the code I posted before and Mug Funky's preset idea (though I used the names "ntsc_film", "ntsc_video", "ntsc_double" for 24000/1001, 30000/1001, 60000/1001 respectively) into the changefps/convertfps/assumefps filters in the custom version of avisynth I keep on my website.
IanB
2nd January 2006, 14:34
Avery's latest blog, Wikipedia and AVI fractions (http://www.virtualdub.org/blog/pivot/entry.php?id=81#body) refers to this article on continued fractions (http://en.wikipedia.org/wiki/Continued_fraction). This might give the narrow minded something to broaden their horizons.
As for guessing what the true FPS the user really intended from a single precision float approximations of FPS, I think I would have to class that with vodoo and other paranomal superstition. Doing so is very fraught. The current code is able to exactly return to the same single precision floating point value originally set, none of the above suggested code mods provide this.
If we change the definition of 29.97 and 23.976 FPS then there will be some lazy users who may be transiently happy but in the long run smart alec hacks, like the proposed, always return to bite you big time. And what do we do with the well intensioned 29.97002997 or 23.976023976 values (which don't fit is a float anyway).
For a case in point consider :- A current user builds a script and set the fps to 23.976 (not 24000/1001), they add an audio track edited in an external editor that has explicit knowledge of single precision 23.976 fps framerate. It works as expected. 9 months later we have a new AVS version, with hacked code, the user is not happy his previous working script inexplicably is now getting audio desync....
The idea for easy to spec string equivalences has merit. Of course it can be implemented immediatly as a simple .avsi script, no need to polute the core.
Bottom line is if you want 30000/1001 fps or 24000/1001 fps then that is what you should spec. So maybe we need to provide text presets to make this easy to understand for the innumerate masses.
P.S. For the uninformed this codewhile (n < 16777216 && d < 16777216) { n*=2; d*=2; }is crudely extracting the mantisa and exponent from the float and storinging them loselessly (well for values less than 2^31) as a rational fraction.
stickboy
2nd January 2006, 23:39
BTW, would it be useful to modify Info() so that instead of:
Frames per second: 23.9760
it showed:
Frames per second: 23.9760 (24000/1001)
?
sh0dan
3rd January 2006, 00:13
@stickboy: Info fractions added to CVS.
Wilbert
3rd January 2006, 00:42
Some discussion on #x264 with respect to this problem
<funmain> bond: high in your opinion is above 23.976 right?
<|bond|> (funmain) bond: high in your opinion is above 23.976 right? <- no
<funmain> |bond|: define 'other high values'
<|bond|> funmain: well it was more meant as a joke that you post there, but it would definitely dont hurt if you would post there
<funmain> 23.976 is based on really really old standards (they should stick to it), i will do bond
<Haali> hehe you seriously propose changing avisynth because of broken qt player?
<|bond|> funmain: 23.976 can be represented as 24000/1001, but thats not 100% correct. avs uses some very high values which make use of 64bit. quicktime doesnt handle that
<|bond|> (Haali) hehe you seriously propose changing avisynth because of broken qt player? <- and mplayer and ffmpeg
<pengvado> we propose to change avisynth because avisynth is _broken_
<pengvado> just because some players happen to play the resulting file doesn't make it correct
<|bond|> (@pengvado) we propose to change avisynth because avisynth is _broken_ <- post your opinion here, cause they dont seem to change it: http://forum.doom9.org/showthread.php?t=104681
<Haali> avisynth is definitely not broken
<|bond|> at least they could add the feature of nicefps()
<|bond|> it wouldnt hurt and help newbies
<Haali> it either passes through what it got from source file or what you input in the script
<Haali> if you do mean 24000/1001 then say so in the script
<|bond|> newbies have no clue about that
<pengvado> except that it doesn't
<pengvado> assumefps(23.976) doesn't set fps=23976/1000
<Haali> because you'd use assumefps(23976,1000) for that, no?
<Haali> anyway, this could be nicely improved by treating the argument as a string, not a double
<pengvado> it seems perfectly reasonable to me to expect assumefps(23.976) and assumefps(23976,1000) to do the same thing
<pengvado> hell, even a double would work. just single precision isn't enough.
<Haali> i agree it can be handled better
<Haali> but in this instance it's the players that can't read valid files, not broken avisynth
<Haali> and pretty obscure players
<Haali> and at least ffmpeg shouldn't be too hard to fix
<funmain> Haali: have u been reading the topic on doom9?
<Haali> yes
<|bond|> some tried to fix ffmpeg to handle 64bit but its still not working correctly and he doesnt answer anymore
<|bond|> noone on the mplayer maillist responded
<Haali> send them some sample they'd definitely want to watch
<|bond|> mplayer ftp isnt working
<Haali> pengvado: btw avisynth uses doubles for args already
<pengvado> I didn't look at the code. but the fps you get from assumefps(23.976) is _exactly_ the float representation on 23.976, not the double representation.
<Haali> yeah, because assumefps only extracts 24 bits from it
<Haali> should have done all 32 :P
<pengvado> ok, so they just need to do the double->rational using nicefps's algo instead of truncating
<Haali> i'd think avoiding a round trip through double can be a better idea
<funmain> pengvado: you dont know.. ?
<pengvado> I know how to fix it, but the patch has already been posted. we just need some avisynth dev to agree to apply it, and I'm not one.
<algern0n> pengvado: time for you to become an avisynth dev then :D
<Koepi> bond is right. you should make an announcement
tritical
3rd January 2006, 05:09
The current code is able to exactly return to the same single precision floating point value originally set, none of the above suggested code mods provide this. This is definitely not true for all cases... why do you think avisynth currently returns a num/den of 982057/8192 for assumefps(119.88) when the actual binary single precision fp approximation of 119.88 is equal to 15712911/131072? The two reasons it doesn't work for all cases are:
1.) d is limited to < 2^24, which means various values in the range 0 < x < 0.5 (granted it's not a common fps value) will not be translated into a num/den combo that when converted back to float will equal the original. Also, in the current code n only needs to be limited to < 2^23 for the loop condition since single precision floating point only has 23 bits for the mantissa.
2.) It will not work for other various values (119.88 is one example), due to how avisynth's script parser converts to float. Take a look at the floating point part of the Tokenizer::GetNumber() routine, it should be doing all the intermediate stuff in double (otherwise too much error accumulates), but it currently does it in float. For example, if you do assumefps(119.88) avisynth spits out 982057/8192 for a num/den. The 982057/8192 comes from assumefps() calling vi.setfps(num,den) which finds the gcd and and divides it out, in this case it is 32. That means original numbers that get passed are 31425824/262144 while the actual binary fp representation of 119.88 is equal to 31425822/262144. The error comes from the script parser conversion, which, as written, doesn't convert 119.88 to 0x42efc28f (hex value of the 4 bytes at the float memory location), but 0x42efc290. If that routine is changed to use doubles for intermediate steps then the correct value (0x42efc28f) is produced.
Anyways, whether or not the code is changed doesn't really matter to me (that's why I keep a custom version :)).
bond
6th January 2006, 15:22
another problem with parsers on such 64bit requiring framerates:
http://forum.doom9.org/showthread.php?t=105154
edit: is the behaviour the same than with assume/changefps when setting the fps via directshowsource?
can we plz change the handling of 29.97, 23.976 aso in avisynth :>
IanB
9th January 2006, 13:55
Okay, I reworked fps.cpp in CVS along the lines from Triticals mods. I added a massive stack of presets to seed the idea mill, I don't intended to leave all those options in, I just wanted them written down so we can make a thought out final selection.
I also reworked the number parser, it now uses doubles for the internals.
One approach I experimented with (sorry not in CVS) is bit stuffing from the float to directly pick out the mantissa and exponent and build the rational values directly from there, no surprise here it gives same result as the original *=2 code.
Bottom line is floats just do not have enough precision to accurately represent the exact rational values required. In cases where exactness is required then you really should be specifying the rational values, not a half baked single precision floating point approximation.
Also for the 0<x<0.5 case if it really matters the following can be done to maximise the precission.double d = 1.0/double(args[1].AsFloat());
int n = 1;
while (n < 16777216 && d < 16777216) { n*=2; d*=2; }
return new AssumeFPS(args[0].AsClip(), int(n+0.5), d, ...
bond
9th January 2006, 14:14
IanB, will this now change the behaviour avoiding triggering the problems i described in the first post if someone uses changefps(23.976) or assumefps(29.97) aso...?
Richard Berg
9th January 2006, 22:06
Aside from code churn, would it break anything if the Avisynth script type 'float' were re-implemented to be a double under the hood?
raymod2
10th January 2006, 09:04
I'd like to weigh in on this issue. I don't agree with Tritical's modification. It ignores the real problem which is the (lossy) string to float conversion. It is easy enough to convert directly from the string to the rational representation (losslessly). I have included an algorithm called string_to_frac() which shows how this can be done.
I have also included an algorithm called reduce_frac() which performs the same function as NiceFPS() but much more efficiently. NiceFPS() uses a brute force method that uses floating point math to evaluate every denominator in the interval [1,limit]. reduce_frac() uses continued fractions and best rational approximations to perform the same task using far fewer calculations without using any floating point math.
I am NOT suggesting that reduce_frac() be incorporated into assumefps() as default behavior. It might be useful, however, to add it to the core so that people who want to use it (and accept the resultant loss in accuracy) don't need to download an external plugin.
#include <stdio.h>
typedef unsigned long UINT32;
//
// This function converts a decimal string (ie. "23.976") into a fraction
// (ie. 2997 / 125). By avoiding a round trip through floating point
// representation, this process is lossless.
//
static char string_to_frac(char *str, UINT32 *num, UINT32 *denom)
{
char ch, digits = 0, dp_found = 0;
UINT32 n = 0, d = 1, x, y, t;
while ( (ch = *str++) )
{
if ( ch >= '0' && ch <= '9' )
{
if (++digits > 9) { return(1); } // too many digits
n *= 10; n += (ch - '0');
if (dp_found) { d *= 10; }
}
else if ( ch == '.' )
{
if (dp_found) { return(1); } // more than one decimal point
dp_found = 1;
}
else
{
return(1); // invalid character
}
}
x = (n > d) ? (n) : (d); // set x to the larger of n,d
y = (n > d) ? (d) : (n); // set y to the smaller of n,d
while (y) { t = x%y; x = y; y = t; } // compute the GCD (result in x)
*num = n/x; *denom = d/x;
return(0);
}
//
// This function uses continued fractions to find the best rational
// approximation that satisfies (denominator <= limit). The algorithm
// is from Wikipedia, Continued Fractions.
//
static void reduce_frac(UINT32 *num, UINT32 *denom, UINT32 limit)
{
UINT32 n0 = 0, n1 = 1, n2, nx = *num; // numerators
UINT32 d0 = 1, d1 = 0, d2, dx = *denom; // denominators
UINT32 a2, ax, amin; // integer parts of quotients
UINT32 f1, f2; // fractional parts of quotients
UINT32 i = 0; // number of loop iterations
while (1) // calculate convergents
{
a2 = nx / dx;
f2 = nx % dx;
n2 = n0 + n1 * a2;
d2 = d0 + d1 * a2;
if (f2 == 0) { break; }
if ((i++) && (d2 >= limit)) { break; }
n0 = n1; n1 = n2; d0 = d1; d1 = d2; f1 = f2;
nx = dx; dx = f2;
}
if (d2 <= limit)
{
*num = n2; *denom = d2; // use last convergent
}
else // (d2 > limit)
{
ax = (limit - d0) / d1; // set d2 = limit and solve for a2
if ((a2 % 2 == 0) && (d0 * f1 > f2 * d1))
{
amin = a2 / 2; // passed 1/2 a_k admissibility test
}
else
{
amin = a2 / 2 + 1;
}
if (ax < amin)
{
*num = n1; *denom = d1; // use previous convergent
}
else
{
// calculate best semiconvergent
*num = n0 + n1 * ax;
*denom = d0 + d1 * ax;
}
}
}
int main(int argc, char *argv[])
{
UINT32 n, d;
string_to_frac("29.97003", &n, &d);
reduce_frac(&n, &d, 5000);
printf("%u / %u\n", n, d);
}
Bidoche
10th January 2006, 09:32
Isn't it already double, though it calls it float ?
stickboy
10th January 2006, 10:28
Nope. AVSValue internally stores a float value. The tokenizer also retrieves floating-point values as floats.
IanB
10th January 2006, 13:22
will this now change the behaviour avoiding triggering the problems i described in the first post if someone uses changefps(23.976) or assumefps(29.97)Well correctly parsing the decimal number string into the single precision floating point value should help some. I am still investigating the effect of the *=10 style code versus the original *=2 code.
From my bit stuffing experiments, I better appreciate this statementd is limited to < 2^24, which means various values in the range 0 < x < 0.5Well to get 0.015625 < x < 0.5 to behave better...while (n < 16777216 && d < 0x3FFFFFFF) { n*=2; d*=2; }
would it break anything if the Avisynth script type 'float' were re-implemented to be a double under the hoodProbably every existing plugin that used AVSValue.AsFloat() due to the included (baked in :devil: ) code from avisynth.h
@raymod2, Yep thats the stuff Avery Lee is so excited about ;)
raymod2
11th January 2006, 20:34
@AviSynth devs: Have you considered incorporating string_to_frac() into changefps()/assumefps()?
1) It is fast and simple (and uses no floating point math).
2) It allows the user to specify up to 9 digits (previous methods allow 7 or 8).
3) It satisfies IanB's request that the user gets exactly what he specifies.
4) It satisfies bond's original request (common fps choices such as 29.97 and 23.976 result in small dwScale values).
All the other proposed solutions suffer the same flaw: The conversion from string to float adds a small error to the specified value. Any attempt to remove that error will be unreliable since the error is unknown.
IanB
11th January 2006, 23:53
@raymod2,
I need to play some more with your routines, but initial thoughts.
To clarify, the string_to_frac() routine effectively always returns NUM as an integer of the string with the "." removed and DENOM as 10**P where P is the number of decimal places. NUM and DENOM then have the GCD stripped out. Thus "29.97"->[2997, 100] and "23.976"->[23976, 1000]->[2997, 125], GCD=8. In fact here the GCD is always going to be (2**A + 5**B) as 2 and 5 are the only prime factors of 10 (and also 10**P). Triticals mods effectively force this in the existing code.
Now the reduce_frac() routine is the interesting bit. Your example "29.97003"->[2997003, 100000]. Then you want to reduce the denominator of this to <=5000. And I seem to get get [60000, 2002], (I must have translated the algo incorrectly, I am away from my C compiler and have to use VB inside Excel, YUCK! More playing at the weekend required)
So okay asuming the code actually produced [30000, 1001] how do we input/choose the "5000" value. This is becoming a little tangled. String presets I can forcibly swallow but I think I might choke if I put a special number string parser into the xxxFPS() routines, it's just to plain hacky.
More thoughts latter,
IanB
raymod2
12th January 2006, 00:21
And I seem to get get [60000, 2002], (I must have translated the algo incorrectly, I am away from my C compiler and have to use VB inside Excel, YUCK! More playing at the weekend required)
I compiled the example program using MinGW and it outputs 30000 / 1001 as expected.
So okay asuming the code actually produced [30000, 1001] how do we input/choose the "5000" value. This is becoming a little tangled.
If reduce_frac() were included I would vote for making its use optional (maybe by providing an optional argument to assumefps() that chooses the limit). The user should be the one to make this decision because it introduces error into the originally specified FPS.
String presets I can forcibly swallow but I think I might choke if I put a special number string parser into the xxxFPS() routines, it's just to plain hacky.
A string->float parser already exists. What is hacky about adding a string->ratio parser? It is the right tool for the job and it is very simple.
IanB
12th January 2006, 01:11
What is hacky about adding a string->ratio parser? We don't have an AVSValue subtype for ratio, so xxxFPS() exclusively would grow a String parameter that is actually a number. I am not sure it is necessary now I have repaired the string to float routine. --total rubbish deleted--
I am still playing with reduce_frac and other repeated fraction algorithms, be patient.
stickboy
12th January 2006, 07:13
What is hacky about adding a string->ratio parser? It is the right tool for the job and it is very simple.If you allow AssumeFPS("29.97"), then you either can't have the preset system proposed earlier or the function needs to decide whether the string argument represents a number or a preset. And sure, it's doable, but I agree it's hacky.
I am not sure it is necessary now I have repaired the string to float routine.So what does AviSynth generate for 23.976 and 29.97 now? If they're reasonable, then I think this is pretty moot.
raymod2
12th January 2006, 08:22
@stickboy: I know a hack when I see one and I despise them. A hack is making unnecessary conversions that introduce error. A hack is trying to fix that error with algorithms that work for some values but not others. A hack is making guesses about what your user might have meant. A hack is saying "the error is small and the user probably won't notice if you don't tell him".
IanB
12th January 2006, 12:32
@stickboy,
There must be something in the water, I have obviously been halocinating. The error only occasionally effected the very bottom bit of the mantissa, so the net result is no (significant) change.
mg262
12th January 2006, 13:29
Ian,
Whatever you end up with, would you overload the relevant SetFPS method in AVISynth.h with SaidFPS(float) so filters have access to the mechanism?
(I guess you've played with replacing { n*=2; d*=2; } with { n*=10; d*=10; } before reduction? That's off the top of my head; hope it makes sense...
Edit: Oops... thiss was already in tritical`s code fragment above...)
bond
27th January 2006, 12:26
so did my feature request die slowly or will actually some change be made to the way avisynth currently acts fixing the initial problem?
IanB
28th January 2006, 13:19
@bond,
The Named Presets have made it into the 2.5.7 CVS, they just need to be selected, organized and reduced. Triticals mod 10 idea has been parked in the code as a compile option so it won't get forgotten. I have finally tested raymod2's reduce_frac code on a C compiler and I like it, however I am not sure that chosing a denominator limit is the most useful mode of operation, as doing so is really just tantemount to just straight specing the numerator and denominator. Off the top of my head I would expect an epsilon accuracy limit would be more useful, and I also would like to see what Avery does in VDub. There will be some form of continued fraction functionality, I am just still researching the best presentation.
bond
28th January 2006, 13:54
ok this means the initial problem will not be solved, as getting people to use the presets is like getting them to set the numerator/denominator. also it would mean to get gui developers to change things
this all will not happen. :(
what about the proposals to solve this from haali and pengvado, wilbert posted?
IanB
28th January 2006, 16:31
@Bond,
As to just forcing AssumeFPS(29.97), AVSValue=0x41efc28, to return 2997/100 or 30000/1001 and AssumeFPS(23.976), AVSValue=0x41bfced9, to return 2997/125 or 24000/1001. I have a real problem with doing this, as it would make discontinuities in the number space and materially effect scripts that calculate FPS values near these 2 magic values.
Unfortunatly Triticals mod 10 suggestion and NiceFPS limits the overall precision to a shade less than 20 bits, even adding an extra iteration with an appropriate factor of 2, 4 or 5 still only brings precision back to approx. 22.25 bits.
I am still experimenting with continued fractions as a way to keep both full 24 bit precision and convienient denominator values. The current conversion only uses 24 distinct denominators, i.e 2**n, 0<=n<24 out of a possible 2**24 so there is a lots of hope.
In the meantime to get around the GUI generator restrictions put an .AVSI in your pluggin directory that overloads the AssumeFPS(float) signature.Function AssumeFPS(Clip clip, Float fps) {
# Triticals algorithm (approx)
Return AssumeFPS(clip, Round(fps*1000000.0), 1000000)
}orFunction AssumeFPS(Clip clip, Float fps) {
# Force explicit values
num=(fps == 29.97)?30000:(fps == 23.976)?24000:Round(fps*16777216)
den=(fps == 29.97)?1001:(fps == 23.976)?1001:16777216
Return AssumeFPS(clip, num, den)
}The GCD code in the AssumeFPS(clip, num, den) signature will normalize the final numerator and denominator values.
raymod2
28th January 2006, 21:20
As to just forcing AssumeFPS(29.97) to return 2997/100... I have a real problem with doing this
I guess I must be missing something here. Since 29.97 is EXACTLY 2997/100 why would you have a problem with that? When a user specifies 29.97 he should get 29.97, not 29.97 plus the small error introduced by converting it to the 32-bit IEEE floating point representation of 29.97. It is exactly that small error that is causing the large denominators.
IanB
29th January 2006, 06:00
I guess I must be missing something here. Since 29.97 is EXACTLY 2997/100...Ding! I can never know the text in the script was exactly 29.97. Inside the code I only ever get to see an AVSValue containing a IEEE single precision float, the value of which is the result the parser turning the text "29.97" into that float. In fact a range of values all result in the same bit pattern for the float representation, as shown in the table belowValue Low High Double Hex
29.97 [ 29.96999836 - 29.97000026 ] = 29.969999313354492, 0x41efc28f
30/1.001 [ 29.97002888 - 29.97003078 ] = 29.970029830932617, 0x41efc29f
23.976 [ 23.97599888 - 23.97600078 ] = 23.975999832153320, 0x41bfced9
24/1.001 [ 23.97602368 - 23.97602558 ] = 23.976024627685547, 0x41bfcee6And I guess this is why you proposed the string input with string_to_frac. However this and the string presets do not address Bond's core problem which he has finaly made me understand. He is using a GUI, it emits a script with the text AssumeFPS(29.97) and he has no hope of getting the Author(s) to mend their ways and emit a more correct text of AssumeFPS(30000, 1001) or AssumeFPS(2997, 100).
Your repeated fraction implementation, reduce_frac, or something similar is probably the ultimate solution, I just need some time to research and experiment with how best to tune it. (remember whatever we do Bond wants exactly AssumeFPS(29.97) to parse nicely.
Manao
29th January 2006, 08:28
and he has no hope of getting the Author(s) to mend their ways and emit a more correct text of AssumeFPS(30000, 1001) or AssumeFPS(2997, 100).It's imho a shame. It's up to them to properly emit such functions, especially since _they_ know what the user want.
stickboy
29th January 2006, 09:59
However this and the string presets do not address Bond's core problem which he has finaly made me understand. He is using a GUI, it emits a script with the text AssumeFPS(29.97) and he has no hope of getting the Author(s) to mend their ways and emit a more correct text of AssumeFPS(30000, 1001) or AssumeFPS(2997, 100).This seems totally backwards. Why is asking the GUI authors to adjust their output more hopeless than getting AviSynth developers to adjust the parser?
IanB
29th January 2006, 11:00
@Manao,
I agree, it is a shame and I am not sure they do actually know. Most people are not mathemagicians and simple don't known the issues about dealing with binary floating point representations and the incompatibilities between base 2 and base 10 numbering systems on the right hand side of the decimal/binary point.
@stickboy,
You are probably right, it is a bit backwards, but fixing it in AVS provides the greatest good. (It's not really the parser, it's the crappy way we currently generate rational FPS from a Float.)
Manao
29th January 2006, 11:04
There's no "uncrappy way" to generate fps from a float. That's the whole issue. A float is unprecise, while a fps must be completely precise.
Imho, the best way is to remove the float parameter from those filters ( but i guess a lot of people wouldn't like that ).
IanB
29th January 2006, 12:05
:D :D :D
Damn vBulletin lower casing my smilie
IanB
30th January 2006, 13:58
Okay, here is a test program I have been exercising for extracting an exact rational pair from a float. The reduce_float routine (hacked from raymod2's code) uses repeated fractions to quickly extract a reasonable rational pair, however it is not always the pair with the smallest denominator. The other routine brute_float is very slow but does find the rational pair with the smallest denominator. The program uses both to cross check.
From inspection, it seems that the smallest denominator for an exact pair is between one half the denominator and the denominator from the reduce_float routine. Any budding mathemagicians got a proof or justification for this, and/or can see a way to get the repeated fraction approach to produce the smallest denominator directly or iteratively faster than brute force.
Personally I am very impressed with the performance of the algorithm and I think it will address Bond's concernes, it produces the expected rational pair for most of the popular fps values, there is just not quite enough precision in a float to correctly reverse 120000/1001 = 119.8801198801, it comes back as 40999/342 = 119.8801169591 a 0.00000243665% error but on the plus side 119.88 does yield 2997/25
For implementation in Avisynth I am undecided whether to just use the repeated fraction algorithm or to supplement it with a range limited brute force to get the absolute smallest denominator. The brute force aproach is very slow try it on a range around 1.0
I await your thoughts, suggestions and comments.
#include <stdlib.h>
#include <stdio.h>
#define UINT32 unsigned long
static void ratio_float(UINT32 *nx, UINT32 *dx, float value)
{
union { float f; UINT32 i; } conv;
if (value < 0) value = -value;
conv.f = value;
*nx = (conv.i & 0x7FFFFF) + 0x800000;
*dx = 0x800000;
int exp = (conv.i >> 23) - 127;
if (exp < 0)
*dx <<= -exp;
else
*dx >>= exp;
}
static bool brute_float(UINT32 lo, UINT32 hi,
UINT32 *num, UINT32 *denom, float value)
{
UINT32 nt, dt;
for (dt=lo; dt<hi; dt++) // Brute force BEST ratio!
{
nt = (UINT32)((double)value * dt); // no rounding
if ((float)((double)nt/(double)dt) == value) // floor
{
*num = nt; *denom = dt;
return true;
}
if ((float)((double)(nt+1)/(double)dt) == value) // ceil
{
*num = nt+1; *denom = dt;
return true;
}
}
*num = nt; *denom = dt;
return false;
}
static void reduce_float(UINT32 *num, UINT32 *denom, float value)
{
UINT32 n0 = 0, n1 = 1, n2, nx; // numerators
UINT32 d0 = 1, d1 = 0, d2, dx; // denominators
UINT32 a2; // integer parts of quotients
UINT32 f2; // fractional parts of quotients
ratio_float(&nx, &dx, value);
for (;;) // calculate convergents
{
a2 = nx / dx;
f2 = nx % dx;
n2 = n0 + n1 * a2;
d2 = d0 + d1 * a2;
if ((f2 == 0) || ((float)((double)n2/(double)d2) == value))
{
*num = n2; *denom = d2; // use convergent
return;
}
n0 = n1; n1 = n2;
d0 = d1; d1 = d2;
nx = dx; dx = f2;
}
}
int main(int argc, char *argv[])
{
int count = 2;
UINT32 n0, d0, n1, d1;
union { float f; UINT32 i; } u;
if (argc > 2) count = atoi(argv[2]);
u.f = 29.97;
if (argc > 1) u.f = atof(argv[1]);
u.i = u.i-count;
for (int i=0;i<=count*2;i++){
reduce_float(&n0, &d0, u.f);
if (brute_float(2, d0/2, &n1, &d1, u.f))
printf("%8u /%7u %.9f (%8u /%7u) 1st Half!\n",
n1, d1, (float)((double)n1/(double)d1), n0, d0);
else if (brute_float(d0/2, d0, &n1, &d1, u.f))
printf("%8u /%7u %.9f (%8u /%7u)\n",
n1, d1, (float)((double)n1/(double)d1), n0, d0);
else
printf("%8u /%7u %.9f\n",
n0, d0, (float)((double)n0/(double)d0));
u.i+=1;
}
}
squid_80
30th January 2006, 14:48
Is speed really an issue? Wouldn't the calculation only be done once for a clip?
mg262
31st January 2006, 00:59
Ian,
I think you can get the best of both worlds... if you expand the continued fraction expansions for the 2 endpoints of the interval you are interested in (I think you call these lo and hi) until they disagree, I believe that
1) The continued fraction expansion for the optimal result agrees with the above expansions (as far as they agree),
and if we truncate the expansion at this point so that one remaining fraction a/b needs to be found
2) The interval within which a/b is allowed to lie is much larger than the original interval (probably on the order of 1/2)
3) The denominator of the final product is a linear function of a and b (with coefficients that can be calculated)
[2 and 3 together imply that finding a/b by brute force is relatively fast.]
Please take this with a pinch of salt as I don't presently have a chance to check this on paper... I will try and do that soon.
IanB
31st January 2006, 08:36
@squid_80,
Well trying to find a rational pair for the float (1.0-epsilon), and then brute forcing for a better pair on my 3GHz PIV takes just on 1 second, on my 400Mhz PII it take over 20 seconds. True, this is one of the worst values, and I can tune and control this, but it pays to be aware.
@mg262,
lo and hi apply to the brute force routine and give it a range to search. In a final version I might use the results of the continued fraction to contrive to scan a much smaller range.
@all,
The routine I am using is not the same as Raymod2's, in his you specify a limit for the denominator and it returns the most accurate rational pair within that constraint. In mine I iteratate for a rational pair that exactly converts back to the original float.
raymod2
31st January 2006, 09:20
IanB, you are not generating EXACT ratios here. The only way to do that is to take the output of ratio_float() and optionally perform GCD on it.
What it appears you are trying to do is come up with a best rational APPROXIMATION that satisfies the condition (error < epsilon) but your testing methodology is flawed because you haven't defined epsilon. Converting from float to rational to float and testing for equivalence is a meaningless test because epsilon changes depending on the magnitude of the input!
I have modified my original function reduce_frac() to do what you want. I have replaced the third parameter 'limit' with 'epsilon' and made the appropriate changes:
typedef unsigned long UINT32;
typedef long long INT64;
//
// This function uses continued fractions to find the best rational
// approximation that satisfies (error <= epsilon). The algorithm
// is from Wikipedia, Continued Fractions.
//
static void reduce_frac2(UINT32 *num, UINT32 *denom, float epsilon)
{
UINT32 n0 = 0, n1 = 1, n2, nx = *num; // numerators
UINT32 d0 = 1, d1 = 0, d2, dx = *denom; // denominators
UINT32 a2, ax, amin; // integer parts of quotients
UINT32 f1, f2; // fractional parts of quotients
float error;
if ((dx == 0) || (epsilon < 0)) { return; } // ERROR: invalid input
while (1) // calculate convergents
{
a2 = nx / dx;
f2 = nx % dx;
n2 = n0 + n1 * a2;
d2 = d0 + d1 * a2;
if (f2 == 0) { break; } // no more convergents (n2 / d2 == input)
error = (INT64)(*num) * d2 - (INT64)(*denom) * n2;
error /= (INT64)(*denom) * d2;
error = (error < 0) ? (-error) : (error);
if (error <= epsilon) { break; }
n0 = n1; n1 = n2; d0 = d1; d1 = d2; f1 = f2;
nx = dx; dx = f2;
}
if (d2 != 1) // we have been through the loop at least twice
{
if ((a2 % 2 == 0) && (d0 * f1 > f2 * d1))
{
amin = a2 / 2; // passed 1/2 a_k admissibility test
}
else
{
amin = a2 / 2 + 1;
}
for (ax = amin; ax <= a2; ax++) // calculate semiconvergents
{
n2 = n0 + n1 * ax;
d2 = d0 + d1 * ax;
error = (INT64)(*num) * d2 - (INT64)(*denom) * n2;
error /= (INT64)(*denom) * d2;
error = (error < 0) ? (-error) : (error);
if (error <= epsilon) { break; }
}
}
*num = n2; *denom = d2;
}
P.S. I have an academic nitpick about ratio_float(). The variable 'exp' could fall anywhere in the interval [128,-127] but this function will return erroneous results if 'exp' falls outside the interval [8,-23].
IanB
31st January 2006, 14:16
IanB, you are not generating EXACT ratios here. ...No indeed I am not, as I said, and I chose those words very carefully...a rational pair that exactly converts back to the original float.To be absolutely explicit :- I want the one rational pair with the smallest denominator value that when converted back is bit identical to the original float value.What it appears you are trying to do is come up with a best rational APPROXIMATION that satisfies the condition (error < epsilon) but your testing methodology is flawed because you haven't defined epsilon.By best, you would mean the one with the smallest denominator. And by implication in the way the code works I am using [magnitude of the least significant bit of the floating point number] as my definition for epsilon. Well actually +/-0.5 of the LSB.Converting from float to rational to float and testing for equivalence is a meaningless test because epsilon changes depending on the magnitude of the input!And because it changes, is why I short circuited the whole issue and tested directly for equivalence of the floats, so it is hardly meaninless. I am just being sneaky.I have modified my original function reduce_frac() to do what you want. I have replaced the third parameter 'limit' with 'epsilon' and made the appropriate changes:Seems to work well. Defining epsilon is tricky, I did it by bit bashing the target float, i.e. zero the mantisa and subtract 24 from the exponent.
I further adapted you code by chucking the error and epsilon calcs and substituted directly comparing the target float value. It is simpler and works the same.
Only problem is it is still linearly searching the possible semi-convergents for the first one that converts back exactly. In a round about way this is what I was already doing, I just didn't have the justification for starting at d2/2, which you have now provided.
I have a thought to use the ax=(limit-d0)/d1 from your original by setting limit=d2-1 iterating on d2 until the float no longer matches. If it works I'll post the results.
P.S. I have an academic nitpick about ratio_float(). The variable 'exp' could fall anywhere in the interval [128,-127] but this function will return erroneous results if 'exp' falls outside the interval [8,-23].Lighten up, it's a cheap and nasty test harness. The outrange values correspond to impossible FPS values i.e outside the range 1/(2^32-1) to (2^32-1)/1 so they don't matter here.
raymod2
31st January 2006, 18:16
I am using [magnitude of the least significant bit of the floating point number] as my definition for epsilon. Well actually +/-0.5 of the LSB.And because it changes, is why I short circuited the whole issue and tested directly for equivalence of the floats, so it is hardly meaninless.
I guess I don't understand why you want a moving epsilon. Do you think that any precision which is lost by calculating (float)((double)n/d) is not worth keeping? You should never perform that calculation when actually using the ratio. For example, if you want to find the frame number at x seconds you would calculate x*((double)n/d) not x*(float)((double)n/d).
Only problem is it is still linearly searching the possible semi-convergents for the first one that converts back exactly. In a round about way this is what I was already doing, I just didn't have the justification for starting at d2/2, which you have now provided.
Well you must search through the convergents linearly because each convergent can only be calculated using the previous two. I suppose you could use a binary search algorithm to search through the semiconvergents but why bother? The list of semiconvergents is generally quite small. The largest I saw while playing with it was 16.
Lighten up, it's a cheap and nasty test harness. The outrange values correspond to impossible FPS values i.e outside the range 1/(2^32-1) to (2^32-1)/1 so they don't matter here.
I said it was a nitpick. I just thought a comment should be added to inform the user that the function only accepted a limited range of floats.
raymod2
1st February 2006, 05:20
Ian, here is an example of a linear search through the semiconvergents with reduce_frac2(2997003, 100000, .00000001). Do you really think it needs to be optimized?
Evaluating semiconvergents over the range (17 <= ax <= 32)
until error is less than or equal to epsilon (1.00e-008)
ax = 17, ratio = (529001/17651), error = 2.66e-008
ax = 18, ratio = (559001/18652), error = 2.36e-008
ax = 19, ratio = (589001/19653), error = 2.09e-008
ax = 20, ratio = (619001/20654), error = 1.84e-008
ax = 21, ratio = (649001/21655), error = 1.62e-008
ax = 22, ratio = (679001/22656), error = 1.41e-008
ax = 23, ratio = (709001/23657), error = 1.23e-008
ax = 24, ratio = (739001/24658), error = 1.05e-008
ax = 25, ratio = (769001/25659), error = 8.96e-009
Found best semiconvergent!
IanB
1st February 2006, 10:58
@Raymod2,
We seem to be continually talking at cross purposes.
1. The aim in this thread (as I think I understand it) is to return a "Nice" rational pair when some gui user's generated script uses the AVS construct AssumeFPS(float).
Churning thru all the chaff and to form a defineable goal, I have assigned the meaning to "Nice" of The one rational pair with the smallest denominator value that when converted back is bit identical to the original float value. I chose "smallest denominator" value because a problem with large denominators seem to be what prompted this thread. And "is bit identical" so that the resulting FPS is as close as possible to the float value specified in the script.
[b]2. The term "epsilon" I am using comes from include/float.h
#define FLT_EPSILON 1.192092896e-07F /* smallest such that 1.0+FLT_EPSILON != 1.0 */
I know use of this term here is not strictly correct, but I hope context provides some additional meaning.
3. The problem with linear searching is the algorithm can be extremely slow. Yes values around 29.97 resolve extremely quickly. Try solving values of 1.0+FLT_EPSILON*N for moderate N.
4. I am not disputing that the algorithm as you have implemented it with either, a denominator limit or an error limit, is very useful. I will probably expose this functionality in version 2.6, but it does not quite meet what is required in this thread, so I keep pushing the envelope.
raymod2
2nd February 2006, 07:50
FLT_EPSILON defines the granularity (distance between adjacent values) in a 32-bit IEEE float when the exponent is 0. To find the granularity for an arbitrary number you must multiply FLT_EPSILON by the magnitude of the exponent part. For example, the granularity at 29.97 is FLT_EPSILON * 16 and the granularity at 119.88 is FLT_EPSILON * 64. Your definition of "nice" is tantamount to choosing epsilon to be FLT_EPSILON * 2^exp / 2. In other words, you are specifying the maximum error to be half the float granularity in the vicinity of the value you are converting. Fair enough.
I see now that the list of semiconvergents can be quite large (such as when there are a lot of zeroes to the right of the decimal point). I have rewritten reduce_frac2() to avoid the search altogether and immediately calculate the best semiconvergent.
It also occurred to me that ratio_float() could be improved to accept a larger range of inputs so I wrote a replacement called float_to_frac().
The code below implements my improvements and will perform float->ratio conversions that are "nice" as you have defined it.
#include <stdio.h>
#include <math.h>
#define FLT_EPSILON 1.192092896e-07F /* smallest such that 1.0+FLT_EPSILON != 1.0 */
typedef unsigned long UINT32;
typedef long long INT64;
//
// This function converts a 32-bit IEEE float into a fraction. This
// process is lossless but it may fail for very large or very small
// floats. It also discards the sign bit. Since the denominator will
// always be a power of 2 and the numerator will always be odd (except
// when the denominator is 1) the fraction will already be reduced
// to its lowest terms.
//
static char float_to_frac(float input, UINT32 *num, UINT32 *denom)
{
union { float f; UINT32 i; } value;
UINT32 mantissa;
int exponent;
value.f = input;
mantissa = (value.i & 0x7FFFFF) + 0x800000; // add implicit bit on the left
exponent = ((value.i & 0x7F800000) >> 23) - 127 - 23; // remove decimal pt
// minimize the mantissa by removing trailing 0's
while (!(mantissa & 1)) { mantissa >>= 1; exponent += 1; }
// if exponent is too large try removing leading 0's of mantissa
while( (exponent > 0) && !(mantissa & 0x80000000) )
{
mantissa <<= 1; exponent -= 1;
}
if (exponent < -31) { return(1); } // number too small
if (exponent > 0) { return(1); } // number too large
*num = mantissa;
*denom = 1 << (-exponent);
return(0);
}
//
// This function uses continued fractions to find the best rational
// approximation that satisfies (error <= epsilon). The algorithm
// is from Wikipedia, Continued Fractions.
//
static void reduce_frac2(UINT32 *num, UINT32 *denom, float epsilon)
{
UINT32 n0 = 0, n1 = 1, n2, nx = *num; // numerators
UINT32 d0 = 1, d1 = 0, d2, dx = *denom; // denominators
UINT32 a2, ax, amin; // integer parts of quotients
UINT32 f1, f2; // fractional parts of quotients
float error, x;
char sign;
if ((dx == 0) || (epsilon < 0)) { return; } // ERROR: invalid input
while (1) // calculate convergents
{
a2 = nx / dx;
f2 = nx % dx;
n2 = n0 + n1 * a2;
d2 = d0 + d1 * a2;
if (f2 == 0) { break; } // no more convergents (n2 / d2 == input)
error = n2 * (INT64)(*denom) - (INT64)(*num) * d2;
error /= d2 * (INT64)(*denom);
if (fabs(error) <= epsilon) { break; }
n0 = n1; n1 = n2; d0 = d1; d1 = d2; f1 = f2;
nx = dx; dx = f2;
}
if (d2 == 1)
{
*num = n2; *denom = d2;
}
else // we have been through the loop at least twice
{
if ((a2 % 2 == 0) && (d0 * f1 > f2 * d1))
{
amin = a2 / 2; // passed 1/2 a_k admissibility test
}
else
{
amin = a2 / 2 + 1;
}
// find the sign of the error (actual + error = n2/d2) and then
// set (n2/d2) = (*num/*denom + sign * epsilon) and solve for a2
sign = (error < 0) ? (-1) : (1);
x = ((INT64)(*denom) * n0 - (*num + (INT64)(*denom) * sign * epsilon) * d0);
x /= ((*num + (INT64)(*denom) * sign * epsilon) * d1 - (INT64)(*denom) * n1);
ax = ceil(x); // round up to nearest integer
if (ax < amin) { ax = amin; }
// calculate best semiconvergent
*num = n0 + n1 * ax;
*denom = d0 + d1 * ax;
}
}
int main(int argc, char *argv[])
{
UINT32 n, d;
union { float f; UINT32 i; } value;
int exp;
float epsilon;
value.f = 29.97003;
exp = ((value.i & 0x7F800000) >> 23) - 127;
epsilon = FLT_EPSILON * pow(2, exp) / 2;
float_to_frac(value.f, &n, &d);
reduce_frac2(&n, &d, epsilon);
printf("(%d/%d) epsilon = %e\n", n, d, epsilon);
}
stickboy
2nd February 2006, 11:00
You are probably right, it is a bit backwards, but fixing it in AVS provides the greatest good. (It's not really the parser, it's the crappy way we currently generate rational FPS from a Float.)Well, the parser currently generates 32-bit floats instead of a 64-bit doubles when a script (although maybe not with some of your recent modifications, I'm not sure), so yeah, I do think some of the fault lies in the parser. If you store 29.97 into a double instead of a float, even some naive rationalizing methods produce better results.
Manao
2nd February 2006, 15:25
Am I still the only one to think that since the user means 24000 / 1001 when he writes 23.976, writing 23.976 shouldn't been allowed ? It's convenient, it's "user friendly", but lets face it, it's not wanted, since it lures the user into using a fps he believes correct, while it isn't.
I was only half joking when i proposed to remove the float argument from the xxxxxxfps() filters.
Also, it would allow devs to spend more time on things that matters more, such as yuy2 -> yv12 chroma conversion.
bond
2nd February 2006, 19:31
well if the user specifies 23.976 avisynth should return 23.976, as its not up to avisynth to guess what the user wants
if the user wants ntsc (24000/1001) or whatever, he should use the presets or set 24000/1001 right away
mg262
2nd February 2006, 19:44
Q. When does the ~one part in a million difference between 23.976 and 24000/1001 actually cause a problem? :confused: Over a three-hour movie it would come to about a 10 ms difference, which is well below the threshold of perception. (It's about the same synchronisation error as that caused by the time taken for the sound from your set to reach you...)
well if the user specifies 23.976 avisynth should return 23.976, as its not up to avisynth to guess what the user wantsIf that's all that's wanted then it seems to me that this issue can be resolved without Ian having to mess around with continued fractions, etc, as in tritical's code fragment above...
Richard Berg
2nd February 2006, 20:17
well if the user specifies 23.976 avisynth should return 23.976
What do you mean exactly? We have to turn it into a fraction at some point, if nothing else when we generate the fake-AVI header. Are you saying we should turn 23.976 into 23976/1000?
Well, the parser currently generates 32-bit floats instead of a 64-bit doubles when a script (although maybe not with some of your recent modifications, I'm not sure), so yeah, I do think some of the fault lies in the parser. If you store 29.97 into a double instead of a float, even some naive rationalizing methods produce better results.
IanB correctly remarked on the previous page that we can't change the script semantics without modifying the AVSValue class. The script parser will have turned the numeral into a C++ float before the AssumeFPS filter gets access to it, so I don't think we can make a special case here, either :(
raymod2
3rd February 2006, 02:43
Here is a summary and analysis of the methods proposed to date:
epsilon = FLT_EPSILON * 2^exp / 2
string_to_frac():
-----------------
- Losslessly converts string to ratio
- Ratio is reduced to its lowest terms
Original AviSynth method:
-------------------------
- Lossy conversion of string to float (error <= epsilon)
- Losslessly converts float to ratio
- Limited range of input which could be improved
- Does not reduce ratio to its lowest terms
float_to_frac():
----------------
- Lossy conversion of string to float (error <= epsilon)
- Losslessly converts float to ratio
- Handles full range of input that can be expressed as a ratio
- Ratio is reduced to its lowest terms
float_to_frac() + reduce_frac2():
---------------------------------
- Lossy conversion of string to float (error <= epsilon)
- Losslessly converts float to ratio
- Lossy reduction of ratio using best rational approximation
(error limit can be chosen - IanB proposes using epsilon)
Tritical's modification:
------------------------
- Lossy conversion of string to float (error <= epsilon)
- Lossy conversion of float to ratio
(error <= 1/((10^(6-int(log10(input))))*2))
- Ratio is reduced to its lowest terms
The following is an example of the different methods using a randomly chosen input:
input = 120.89032, epsilon = 3.815e-006
-----------------------------------------
string_to_frac() ( 1511129/12500 )
AviSynth original ( 31690672/262144 ), error = -1.758e-007
float_to_frac() ( 1980667/16384 ), error = -1.758e-007
IanB's nice fraction ( 18738/155 ), error = +2.581e-006
Tritical's mod ( 1208903/10000 ), error = -2.000e-005
Obviously the magnitude of the error can and will change in a random manner depending on the input. But the first method is always lossless, the next two are bounded by epsilon, the fourth is bounded by 2*epsilon, and Tritical's mod is bounded by epsilon + x (where x in this example is 5.000e-005 which is over 13 times larger than epsilon).
mg262
3rd February 2006, 02:55
The error in tritical's version (IMO the right thing to do in this case) is governed by this line:
while (n < 1000000.0 && d < 1000000) { n*=10.0; d*=10; }
Edit: the remainder of this post would no longer make sense since the above post has been edited. Note that the '13 times larger than epsilon' is misleading and can be substantially reduced by increasing 1000000.0 in the above code.
raymod2
3rd February 2006, 09:16
tritical's version (IMO the right thing to do in this case)
@mg262: How can you read my analysis and then come to that conclusion? Tritical's modification performed the worst of all. It's error bound is very large and hard to define. Here is a step-by-step breakdown of his algorithm operating on the input 120.89032:
operation machine representation value error
-----------------------------------------------------------------------------------
1.208903200e+002
string->float 42f1c7d8 1.208903198e+002 +1.758e-007
float->(double/int) (405e38fb00000000/00000001) 1.208903198e+002 +0.000e+000
(double/int) *= (10/10) (4092e39ce0000000/0000000a) 1.208903198e+002 +0.000e+000
(double/int) *= (10/10) (40c79c8418000000/00000064) 1.208903198e+002 +0.000e+000
(double/int) *= (10/10) (40fd83a51e000000/000003e8) 1.208903198e+002 +0.000e+000
(double/int) *= (10/10) (4132724732c00000/00002710) 1.208903198e+002 +0.000e+000
(double/int) += 0.5 (41327247b2c00000/00002710) 1.208903698e+002 -5.000e-005
(double/int)->(int/int) (00127247/00002710) 1.208903000e+002 +6.982e-005
-----------------------------------------------------------------------------------
-2.000e-005
As you can see, using a multiplier of 10 changes the mantissa and the exponent. The mantissa inevitably grows and when it comes time to convert it to an integer you have to chop off part of the mantissa. That is where the error comes from (the step which adds 0.5 has no effect on the end result since it goes into the part of the mantissa that gets truncated).
Contrast this to the original AviSynth algorithm which uses a multiplier of 2 and doesn't touch the mantissa (just increments the exponent) so the double->int conversion is lossless.
Doom9
3rd February 2006, 10:00
Q. When does the ~one part in a million difference between 23.976 and 24000/1001 actually cause a problem?I think I should weigh in my POV seeing that I work on a GUI. The change in DGIndex to no longer report 23976 / 29970 (in current betas only.. or we'd have gotten a lot more reports all of a sudden) actually broke my application, and it is yet unclear whether there's going to be any influence on other tools. If you have a FPS dropdown with your standard values and you try to match the FPS reported from a script to it, if it's no longer 23976/29970, you're either screwed and go to the default (and then have people asking you why their "23.976" script comes out as 25 fps or whatever you set as the default), or you have to implement a fuzzy detection. It also remains to be seen what effect this will have on additional software, like muxers that need an FPS specified.. if you take the value from the source and you get 23.97869023743whatnot, they might just run into trouble themselves. Since so far, not a single person has reported any problem with the current way things are done (23976 / 1000), and that thing has been going on even longer than my site has been around (which would be 6 years next month), I am very sceptical towards making any change now that risks to break entire toolschains that have been stable for years, just because two players cannot handle 64 bit timestamps in MP4. It will be much less effort in the global scheme of things to fix these two players and leave things in a working state elsewhere.
Being very active in the AVC forum myself, I also have to violently disagree with cause it would remove like "50%" of the x264 error reportsIf quicktime or mplayer can't handle a file, that's a problem of those players, not of x264 and the percentage is wildly exagerated.
Richard Berg
3rd February 2006, 10:25
If you have a FPS dropdown with your standard values and you try to match the FPS reported from a script to it, if it's no longer 23976/29970, you're either screwed and go to the default (and then have people asking you why their "23.976" script comes out as 25 fps or whatever you set as the default), or you have to implement a fuzzy detection.
Any GUI that chooses a dropdown value from a script simply must use fuzzy detection. There are far far more scripts out there than what DGIndex produces; IMO we have to be able to handle any valid script the user wants to input.
mg262
3rd February 2006, 13:23
raymod2,
Look at your numerical error again: -2.000e-005 ... that's not caused by imprecision in floating point arithmetic. That's clearly caused by something decimal. In particular:
while (n < 1000000.0 && d < 1000000) { n*=10.0; d*=10; }
120.89032/1 -> 1208903.2/10000 -> 1208903/10000
The error has happened in the second step, where the 0.2 is truncated. Obviously the error can't be removed entirely, but it can be substantially reduced (and the magnitude of the error can be bounded above by a small multiple of epsilon). To test this quickly, try adding an extra 0 (or two) to the 1000000.0.
There is of course an issue with multiplying by too large a power of 10, which is presumably what the (n < 1000000.0 && d < 1000000) check is for... but I don't think that n < 1000000.0 is the appropriate bound. n needs to be bounded to prevent it overflowing an integer (unless we use __int64), and d needs to be bounded to prevent d*epsilon from becoming greater than one.
____________________
Doom9, Richard Berg,
To clarify my stance, it's that AVISynth should return 23976/1000, not mess with continued fractions,* and that I can see no case where the 'exact' frame rate needs to be represented with more precision than by 23976/1000.
*Are we all arguing on the same side? I've seen several people, including Bond, supporting 23976/1000 but no one asking for continued fractions.
I take the point about fuzzy detection; cf. this:
http://forum.doom9.org/showthread.php?t=79888
... that was actually a serious headache at the time.
raymod2
4th February 2006, 04:20
@mg262,
Here is a quantitative error analysis of Tritical's mod:
Analysis of error bound in Tritical's mod for input = 120.89032
---------------------------------------------------------------
exp = int(log2(input)) = 6
EXP = int(log10(input)) = 2
LIM = int(log10(numerator limit)) = 6
1
(string->float) error bound = ------------------ = 3.815e-006
2^(23 - exp) * 2
1
(float->ratio) error bound = -------------------- = 5.000e-005
10^(LIM - EXP) * 2
The first error is the same error I have called epsilon in previous analyses. The second (additional) error results from Tritical's mod because he uses a lossy method to extract the mantissa from the float. Note that changing the numerator limit (as you suggested) will decrease the magnitude of this error but it will never eliminate it completely. Hence, Tritical's mod will always perform worse than the original AviSynth method.
mg262
4th February 2006, 10:50
1. tritical clearly meant that code to convey an idea, not to be fine-tuned down to the last dot. Since you are editing your summary, why not actually try the change?
2. Obviously the original AVISynth method has a smaller error than anything else (excluding red herrings about strings). But no one is aiming for the method that has the smallest numerical error. Read Ian's comments again.
3. If one method had a much higher error than all the others, it might be an issue... but this simply isn't the case.
Doom9
4th February 2006, 11:44
To clarify my stance, it's that AVISynth should return 23976/1000, not mess with continued fractions,* and that I can see no case where the 'exact' frame rate needs to be represented with more precision than by 23976/1000.My point exactly.. your average 23.976 or 29.97 fps AVI reports 23976/1000 and 29970/1000 respectively (dwRate/dwScale in the AVIFile API).. If opening an Avisynth script that contains assumefps(23.976) were suddenly to return a 24000/1001 dwRate/dwScale fraction, there's no knowing what chain reaction that would have an existing tools relying on the knowledge of ntsc film = 23.976 and ntsc = 29.97. Most tools and input methods only use a float value, just like avisynth does in these commands.
raymod2
5th February 2006, 03:01
1. tritical clearly meant that code to convey an idea, not to be fine-tuned down to the last dot. Since you are editing your summary, why not actually try the change?
What is there to try? There is more to numerical analysis than plugging numbers. By providing a mathematical formula for the error limit I have shown that it can be minimized but never eliminated completely. Also, don't you realize that increasing the numerator limit also increases the denominator? That defeats the purpose behind Tritical's mod because now you will be getting more error AND a larger denominator than the original AviSynth method.
In defense of Tritical's mod it does work perfectly when the user specifies 7 or less significant decimal digits. By forcing the denominator to be a power of 10, the algorithm "guesses" what the user specified and the error it adds will cancel out the error from the string->float conversion. My only other criticism is that it unnecessarily limits the denominator (denominator limit should be 10^9 instead of 10^6) and I just noticed that GCD is done in SetFPS() so it is unnecessary to do it in AssumeFPS().
IanB
5th February 2006, 08:57
As said above, Tritical's mod assumes the input was originally decimal and uses the bottom bits of the mantissa as guard bits. This is trading precision for a specific functionality. The Continued Fractions based method returns results as good as or better than this. It also fits the bill for not affecting the output for existing scripts, in fact the smaller denominator values should improve matters.
@Raymod2,
I am currently testing your latest reduce_frac2, seems to get the "nicest" pair every time. :D
@[Those who still think typing 23.976 means 23976/1000],
1. There is a conversion to single precision floating point format. It cannot be avoided.
2. Most decimal numbers do NOT have an exact representation in floating point format (even using doubles).
3. Although the Continued Fractions algorithm does actually return the expected trivial result for 23.976, 29.97, 59.94 & 119.88, it should not be expected in the general case. For example 119.88012 does not return 120000/1001.
4. If the AVI framerate ratio really matters in your application, then explicitly specify the numerator and denominator. The construct to do so is Avisynth exactly for this reason.
mg262
5th February 2006, 10:47
raymod2,
Also, don't you realize that increasing the numerator limit also increases the denominator?Obviously. But the upper bound on the denominator (which is IMO what we care about) is unchanged.
What I was originally objecting to was the inference from a particular numerical example that was in your post before you edited it. Since you have edited that post substantially, I'm not going to argue with you.
IanB,
Edit: scratch that... I misunderstood what you meant by "expected trivial result"!
IanB
5th February 2006, 11:31
@mg262,
You cannot increase the 1000000 in Tritical's code. Doing so would exceed the maximum available precission and no longer achieve the desired decimal implied conversion. (No guard bits!)
@All,
Expected trivial results are 23976/1000=2997/125, 2997/100, 2997/50 & 2997/25
mg262
5th February 2006, 16:19
Ian,
i.e. there is an issue that e.g. 120.0 might turn into 12,000,001/100,000? Yes, you're completely right... I had my head in maths-mode not programming-mode and I was thinking of epsilon as a constant. :o
raymod2
5th February 2006, 23:52
Here is a plot of the error bounds using the different string->ratio methods. Keep in mind that the actual error for Tritical's mod is 0 when the user specifies 7 or less significant decimal digits.
IanB
7th February 2006, 12:22
@Raymod2,
Small problem with your graph, the way you join the points together is invalid. With all the methods the error does dip to zero as the returned ratio and the original number just happen to coincide. It's a bit like a broken clock being exactly right twice a day :D
raymod2
7th February 2006, 20:19
@Ian, it's a graph of the error bound not the actual error. It would be impractical to plot the actual error unless you chose a very small interval for x.
IanB
14th February 2006, 12:40
@All,
I have been thinking about the K*30000/1001 and K*24000/1001 rational pairs failing to be such for K>2 due to insufficient precision. In the latest test harness (below) I have prototype code to check if the float is within range of a desired dropframe multiple and print a message. It also tests for 30000/(K*1001) and 24000/(K*1001) sub-multiples.
Thorts? :)
@Raymod2,
I have finished testing your latest reduce_frac2, it only gets the "nicest" pair most of the time. It seems to calculate the AX value for the nicest semiconvergent a little to high, mostly by 1 or 2, occasionally by 3. I have tried to tune out the error, but was only mildly successfull. (See hackings below)
Also your epsilon value calculation seems borderline, occasionally the returned rational pair are off by 1 epsilon. Setting epsilon 1 step lower seems to fix the problem, i.e ((1.-1.>>24)>>24) instead of (1.>>24)
My latest test harness code:-#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#define FLT_EPSILON 1.192092896e-07F /* smallest such that 1.0+FLT_EPSILON != 1.0 */
typedef unsigned long UINT32;
typedef long long INT64;
//
// This function converts a 32-bit IEEE float into a fraction. This
// process is lossless but it may fail for very large or very small
// floats. It also discards the sign bit. Since the denominator will
// always be a power of 2 and the numerator will always be odd (except
// when the denominator is 1) the fraction will already be reduced
// to its lowest terms. output range (2^32-2^(32-24))/1 -- (1/2^31)
// i.e. 4294967040/1 -- 1/2147483648(4.65661287307e-10)
//
static char float_to_frac(float input, UINT32 *num, UINT32 *denom)
{
union { float f; UINT32 i; } value;
UINT32 mantissa;
int exponent;
value.f = input;
mantissa = (value.i & 0x7FFFFF) + 0x800000; // add implicit bit on the left
exponent = ((value.i & 0x7F800000) >> 23) - 127 - 23; // remove decimal pt
// minimize the mantissa by removing trailing 0's
while (!(mantissa & 1)) { mantissa >>= 1; exponent += 1; }
// if exponent is too large try removing leading 0's of mantissa
while( (exponent > 0) && !(mantissa & 0x80000000) )
{
mantissa <<= 1; exponent -= 1;
}
if (exponent < -31) { return(1); } // number too small
if (exponent > 0) { return(1); } // number too large
*num = mantissa;
*denom = 1 << (-exponent);
return(0);
}
//
// This function uses continued fractions to find the best rational
// approximation that satisfies (denom <= limit). The algorithm
// is from Wikipedia, Continued Fractions.
//
static void reduce_frac(UINT32 *num, UINT32 *denom, UINT32 limit)
{
UINT32 n0 = 0, n1 = 1, n2, nx = *num; // numerators
UINT32 d0 = 1, d1 = 0, d2, dx = *denom; // denominators
UINT32 a2, ax, amin; // integer parts of quotients
UINT32 f1, f2; // fractional parts of quotients
UINT32 i = 0; // number of loop iterations
while (1) // calculate convergents
{
a2 = nx / dx;
f2 = nx % dx;
n2 = n0 + n1 * a2;
d2 = d0 + d1 * a2;
if (f2 == 0) { break; }
if ((i++) && (d2 >= limit)) { break; }
n0 = n1; n1 = n2; d0 = d1; d1 = d2; f1 = f2;
nx = dx; dx = f2;
}
if (d2 <= limit)
{
*num = n2; *denom = d2; // use last convergent
}
else // (d2 > limit)
{
ax = (limit - d0) / d1; // set d2 = limit and solve for a2
if ((a2 % 2 == 0) && (d0 * f1 > f2 * d1))
{
amin = a2 / 2; // passed 1/2 a_k admissibility test
}
else
{
amin = a2 / 2 + 1;
}
if (ax < amin)
{
*num = n1; *denom = d1; // use previous convergent
}
else
{
// calculate best semiconvergent
*num = n0 + n1 * ax;
*denom = d0 + d1 * ax;
}
}
}
//
// This function uses continued fractions to find the best rational
// approximation that satisfies (error <= epsilon). The algorithm
// is from Wikipedia, Continued Fractions.
//
static void reduce_frac2(UINT32 *num, UINT32 *denom, float epsilon)
{
UINT32 n0 = 0, n1 = 1, n2, nx = *num; // numerators
UINT32 d0 = 1, d1 = 0, d2, dx = *denom; // denominators
UINT32 a2, ax, amin; // integer parts of quotients
UINT32 f1, f2; // fractional parts of quotients
float error, x;
int sign;
if ((dx == 0) || (epsilon < 0)) { return; } // ERROR: invalid input
while (1) // calculate convergents
{
a2 = nx / dx;
f2 = nx % dx;
n2 = n0 + n1 * a2;
d2 = d0 + d1 * a2;
error = n2 * (INT64)(*denom) - (INT64)(*num) * d2;
error /= d2 * (INT64)(*denom);
if (f2 == 0) { break; } // no more convergents (n2 / d2 == input)
if (fabs(error) <= epsilon) { break; }
n0 = n1; n1 = n2; d0 = d1; d1 = d2; f1 = f2;
nx = dx; dx = f2;
}
if (d2 == 1)
{
*num = n2; *denom = d2;
}
else // we have been through the loop at least twice
{
if ((a2 % 2 == 0) && (d0 * f1 > f2 * d1))
{
amin = a2 / 2; // passed 1/2 a_k admissibility test
}
else
{
amin = a2 / 2 + 1;
}
// find the sign of the error (actual + error = n2/d2) and then
// set (n2/d2) = (*num/*denom + sign * epsilon) and solve for a2
INT64 nE = (error < 0) ? -(INT64)(*denom) : (*denom);
x = ((INT64)(*denom)*n0) - ((INT64)(*num)*d0+llroundf((nE*d0)*epsilon));
x /= ((INT64)(*num)*d1+llroundf((nE*d1)*epsilon)) - ((INT64)(*denom)*n1);
ax = (UINT32)lroundf(x); // round up to nearest integer
// sign = (error < 0) ? (-1) : (1);
// x = ((INT64)(*denom) * n0 - (*num + (INT64)(*denom) * sign * epsilon) * d0);
// x /= ((*num + (INT64)(*denom) * sign * epsilon) * d1 - (INT64)(*denom) * n1);
// ax = (UINT32)ceilf(x); // round up to nearest integer
if (ax < amin) { ax = amin; }
// calculate best semiconvergent
*num = n0 + n1 * ax;
*denom = d0 + d1 * ax;
}
}
////////////////////////////////////////////////////////////
// Extract ratio from float using bit bashing
// No error checking!
static void ratio_float(float value, UINT32 *nx, UINT32 *dx)
{
union { float f; UINT32 i; } conv;
if (value < 0) value = -value;
conv.f = value;
*nx = (conv.i & 0x7FFFFF) + 0x800000;
*dx = 0x800000;
int exp = (conv.i >> 23) - 127;
if (exp < 0)
*dx <<= -exp;
else
*dx >>= exp;
}
// Find by brute force the smallest rational pair between lo and
// hi such that the result truncated to a float is equal to value.
static bool brute_float(float value, UINT32 lo, UINT32 hi, UINT32 *num, UINT32 *denom)
{
UINT32 nt, dt;
for (dt=lo; dt<hi; dt++) // Brute force BEST ratio!
{
nt = (UINT32)((double)value * dt); // no rounding
if ((float)((double)nt/(double)dt) == value) // floor
{
*num = nt; *denom = dt;
return true;
}
if ((float)((double)(nt+1)/(double)dt) == value) // ceil
{
*num = nt+1; *denom = dt;
return true;
}
}
*num = nt; *denom = dt;
return false;
}
// Find by continued fractions the rational pair such
// that the result truncated to a float is equal to
// value. Then binary search the semiconvergents for
// the smallest such rational pair.
static void reduce_float(float value, UINT32 *num, UINT32 *denom)
{
UINT32 n0 = 0, n1 = 1, n2, nx; // numerators
UINT32 d0 = 1, d1 = 0, d2, dx; // denominators
UINT32 a2, ax, amin; // integer parts of quotients
UINT32 f1, f2; // fractional parts of quotients
ratio_float(value, &nx, &dx);
while (1) // calculate convergents
{
a2 = nx / dx;
f2 = nx % dx;
n2 = n0 + n1 * a2;
d2 = d0 + d1 * a2;
if (f2 == 0) { break; } // no more convergents (n2 / d2 == input)
if ((float)((double)n2/(double)d2) == value) break;
n0 = n1; n1 = n2; d0 = d1; d1 = d2; f1 = f2;
nx = dx; dx = f2;
}
*num = n2; *denom = d2;
if (d2 != 1) { // we have been through the loop at least twice
if ((a2 % 2 == 0) && (d0 * f1 > f2 * d1))
amin = a2 / 2; // passed 1/2 a_k admissibility test
else
amin = a2 / 2 + 1;
UINT32 ahi = (d2-1 - d0) / d1; // set d2 = limit and solve for a2
while (amin <= ahi) { // search for nicest semiconvergent - binary chop
ax = (ahi + amin)/2;
n2 = n0 + n1 * ax;
d2 = d0 + d1 * ax;
if ((float)((double)n2/(double)d2) == value) {
*num = n2; *denom = d2;
ahi = ax-1;
}
else
amin = ax+1;
}
}
}
int main(int argc, char *argv[])
{
int count;
UINT32 n0, d0, n1, d1, n2, d2, n3, d3;
union { float f; UINT32 i; } u, eps;
if (argc == 4) {
count = atoi(argv[3]);
u.f = atof(argv[1]) / atof(argv[2]);
}
else {
count = (argc > 2) ? atoi(argv[2]) : 8;
u.f = (argc > 1) ? atof(argv[1]) : 29.970014572;
}
u.i = u.i-count;
for (int i=0;i<=count*2;i++){
reduce_float(u.f, &n0, &d0);
int i = (int)(u.f*1001+0.5);
float x = (i/24000*24000)/1001.0;
if (x == u.f) printf("%8u / 1001 ---- 24\n",(n0*1001+d0/2)/d0, i);
x = (i/30000*30000)/1001.0;
if (x == u.f) printf("%8u / 1001 ---- 30\n",(n0*1001+d0/2)/d0, i);
if (u.f < 15) {
i = (int)(30000/u.f+0.5);
x = 30000.0/(i/1001*1001);
if (x == u.f) printf(" 30000 /%7u ---- 1/30\n", i);
i = (int)(24000/u.f+0.5);
x = 24000.0/(i/1001*1001);
if (x == u.f) printf(" 24000 /%7u ---- 1/24\n", i);
}
/* if (brute_float(u.f, d0/2, d0, &n1, &d1))
printf("%8u /%7u %.9f (%8u /%7u) %.9f\n",
n0, d0, (float)((double)n0/(double)d0), n1, d1, u.f);
else
*/
printf("%8u /%7u %.9f\n",
n0, d0, (float)((double)n0/(double)d0));
// eps.i = (u.i & 0x7F800000) - (24 << 23); // epsilon=(1>>24)*2^K
eps.i = (u.i & 0x7F800000) - (25 << 23) + 0x4fffff; // epsilon=((1-1>>24)>>24)*2^K
float_to_frac(u.f, &n2, &d2);
reduce_frac2(&n2, &d2, eps.f);
if (brute_float(u.f, d2/2, d2, &n3, &d3))
printf("%8u /%7u %.9f (%8u /%7u) %.9f\n\n",
n2, d2, (float)((double)n2/(double)d2), n3, d3, u.f);
else
printf("%8u /%7u %.9f\n\n",
n2, d2, (float)((double)n2/(double)d2));
u.i+=1; // u.f += epsilon!
}
}
raymod2
14th February 2006, 18:15
@IanB: Can you give me an example of input that causes reduce_frac2() to fail and what is the correct output?
raymod2
16th February 2006, 06:57
All right Ian, I used a random number generator to generate (num,denom) pairs and fed them to reduce_frac2() to see if I could find any failures. I discovered that there were some cases where "sign" was assigned incorrectly (ie. when the current convergent is also the last convergent). I also discovered that there were some cases where float math was not precise enough when calculating ax so I switched to doubles. I have attached the corrected version. I tested this version with 5,000,000 random inputs and it always produced the correct output.
//
// This function uses continued fractions to find the best rational
// approximation that satisfies (error <= epsilon). The algorithm
// is from Wikipedia, Continued Fractions.
//
static void reduce_frac2(UINT32 *num, UINT32 *denom, float epsilon)
{
UINT32 n0 = 0, n1 = 1, n2, nx = *num; // numerators
UINT32 d0 = 1, d1 = 0, d2, dx = *denom; // denominators
UINT32 a2, ax, amin; // integer parts of quotients
UINT32 f1, f2; // fractional parts of quotients
float error;
double x;
char sign;
if ((dx == 0) || (epsilon < 0)) { return; } // ERROR: invalid input
while (1) // calculate convergents
{
a2 = nx / dx;
f2 = nx % dx;
n2 = n0 + n1 * a2;
d2 = d0 + d1 * a2;
if (f2 == 0) { break; } // no more convergents (n2 / d2 == input)
error = n2 * (INT64)(*denom) - (INT64)(*num) * d2;
error /= d2 * (INT64)(*denom);
if (fabs(error) <= epsilon) { break; }
n0 = n1; n1 = n2; d0 = d1; d1 = d2; f1 = f2;
nx = dx; dx = f2;
}
if (d2 == 1)
{
*num = n2; *denom = d2;
}
else // we have been through the loop at least twice
{
if ((a2 % 2 == 0) && (d0 * f1 > f2 * d1))
{
amin = a2 / 2; // passed 1/2 a_k admissibility test
}
else
{
amin = a2 / 2 + 1;
}
// We need the sign of the error for the semiconvergents. This will be
// opposite to the sign of the error for the previous convergent (n1/d1).
sign = (n1 * (INT64)(*denom) > (INT64)(*num) * d1) ? (-1) : (1);
// Set (n2/d2) = (*num/*denom + sign * epsilon) and solve for a2:
//
// (*denom)(n0) - (*num + (*denom)(sign)(epsilon))(d0)
// ax = ---------------------------------------------------
// (*num + (*denom)(sign)(epsilon))(d1) - (*denom)(n1)
x = ((double)(*denom) * n0 - (*num + (double)(*denom) * sign * epsilon) * d0);
x /= ((*num + (double)(*denom) * sign * epsilon) * d1 - (double)(*denom) * n1);
ax = ceil(x); // round up to nearest integer
if (ax < amin) { ax = amin; }
// calculate best semiconvergent
*num = n0 + n1 * ax;
*denom = d0 + d1 * ax;
}
}
IanB
16th February 2006, 08:26
Sorry, below applies to your previous post. (last one of the prev page, damn I can really hate vbulletin sometimes :mad: )
And yes this version seems to be spot on ;)
@Raymod2,
24.999996185 = 5242899 / 209716
(8738149/349526) 24.999998093 epsilon=9.536743e-07
24.999998093 = 8738149 / 349526
(26214399/1048576) 25.000000000 epsilon=9.536743e-07
25.000001907 = 8738151 / 349526
(26214401/1048576) 25.000000000 epsilon=9.536743e-07
25.000003815 = 5242901 / 209716
(8738151/349526) 25.000001907 epsilon=9.536743e-07
Reducing epsilon to 7.748603e-07, i.e. the next lower float seem to fix this problem, but then it does not alway return the nicest rational pair, i.e. ax comes out a little high.
Can you please explain how you intend you ax=... code to work. Please do NOT reread your code first, just explain how you intended it to work. I suspect there is something subtle missing from the code as written.
And yes error was not current.
bond
27th March 2006, 10:59
heya, with the recent avisynth update, what has changed regarding the initial issue?
foxyshadis
27th March 2006, 11:31
Tritical added presets and Raymond's code:
* xxxFPS("preset") string preset FPS values. (Tritical)
* xxxFPS(float) now uses continued fraction to generate a minimal rational pair. (Raymod2)
bond
27th March 2006, 11:38
yeah i saw this too, but what does it mean in practice
Wilbert
27th March 2006, 11:45
yeah i saw this too, but what does it mean in practice
It means that floats like 29.97 and 23.976 are approximated accurately by ratios of small integers, like i explained in the off-line documentation (the exact meaning of this is not explained there). I wonder why i put much effort in updating the documentation :)
bond
27th March 2006, 13:13
It means that floats like 29.97 and 23.976 are approximated accurately by ratios of small integers, like i explained in the off-line documentation (the exact meaning of this is not explained there). I wonder why i put much effort in updating the documentation :)thx, i checked the online documentation but didnt find a change mentioned
i checked this new feature and it seems to work indeed fine, fixing the problems i started the thread about
thx a lot :)
IanB
27th March 2006, 15:51
Okay the nitty gritty ;)
For the set of floats matching the condition 30000*N/1001 or 24000*N/1001, the ratio is set to {30000*N, 1001} or {24000*N, 1001}
For the set of floats matching the condition 30000/(1001*N) or 24000/(1001*N), the ratio is set to {30000, 1001*N} or {24000, 1001*N} Normalized by the GCD.
For all other float values the rational pair {Num, Den} with the smallest possible value Den such that
(float)((double)Num/(double)Dem) == value
is returned.
All the popular values mentioned in this thread, their multiple and their factors return a nice value. :D
pyrates
8th September 2006, 07:37
Okay the nitty gritty ;)
For the set of floats matching the condition 30000*N/1001 or 24000*N/1001, the ratio is set to {30000*N, 1001} or {24000*N, 1001}
For the set of floats matching the condition 30000/(1001*N) or 24000/(1001*N), the ratio is set to {30000, 1001*N} or {24000, 1001*N} Normalized by the GCD.
For all other float values the rational pair {Num, Den} with the smallest possible value Den such that
(float)((double)Num/(double)Dem) == value
is returned.
All the popular values mentioned in this thread, their multiple and their factors return a nice value. :D
Has this been integrated yet or would avisynth 2.5.6a not have it yet?
foxyshadis
8th September 2006, 08:02
2.5.6 will never have it, but 2.5.7 does have an implementation of it. Ian's enough of a perfectionist that he might still rewrite it all though. ;)
pyrates
8th September 2006, 21:11
2.5.6 will never have it, but 2.5.7 does have an implementation of it. Ian's enough of a perfectionist that me might still rewrite it all though. ;)
I thought so. Just added assumefps(23.976,1000) at the end of my avisynth script when doing my third x264 encoding of this movie. Maybe this time it will actually seek. The previous 2 times it didn't.
foxyshadis
8th September 2006, 21:52
You mean asumefps(23976,1000)? Yours will give 0.023 fps... Assumefps(24000,1001) is the true framerate you need, anyway.
bond
9th September 2006, 14:29
Has this been integrated yet or would avisynth 2.5.6a not have it yet?2.5.7 has it and it works great :)
pyrates
12th September 2006, 21:31
You mean asumefps(23976,1000)? Yours will give 0.023 fps... Assumefps(24000,1001) is the true framerate you need, anyway.
Ah yes you're right, just misspelled it in that post, not in the avisynth script. I did put assumefps(23976,1000). As for assumefps(24000,1001), doesn't that then amount to 23.98 instead of the native fps which is 23.976?
foxyshadis
12th September 2006, 22:51
24000/1001 = 23.976023976023976023976023976024.... Which is why it's "more accurate". :p
pyrates
15th September 2006, 08:49
24000/1001 = 23.976023976023976023976023976024.... Which is why it's "more accurate". :p
Oh that's right. Guess I'll use assumefps(24000,1001) from now on. Thanks.
vBulletin® v3.8.5, Copyright ©2000-2012, Jelsoft Enterprises Ltd.