avs devs: plz integrate nicefps() into changefps()/assumefps()

[raymod2]

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]

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]

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]

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]

@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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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]

@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]

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]

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]

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]

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]

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]

@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]

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]

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]

@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]

@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]

@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]

@IanB: Can you give me an example of input that causes reduce_frac2() to fail and what is the correct output?

[raymod2]

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]

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]

heya, with the recent avisynth update, what has changed regarding the initial issue?

[foxyshadis]

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]

yeah i saw this too, but what does it mean in practice

[Wilbert]

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]

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]

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]

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]

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]

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]

You mean asumefps(23976,1000)? Yours will give 0.023 fps... Assumefps(24000,1001) is the true framerate you need, anyway.

[bond]

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]

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]

24000/1001 = 23.976023976023976023976023976024.... Which is why it's "more accurate". :p

[pyrates]

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.