Add Pow, Ln, and Exp functions #12

sam-vdp · 2018-02-27T10:24:24Z

Add Exp and Ln functions and Pow based on the combination of both. Ln uses Newton's method to approximate the solution, this code was ported from libfixmath, which your license already covers. The Exp function uses a power series

asik

Thank you so much for this contribution! I have made some comments here.

asik · 2018-02-27T21:30:21Z

tests/Fix64Tests.cs

+            {
+                var b = Fix64.FromRaw(m_testCases[j]);
+
+                if (b <= Fix64.Zero)


We should test what we expect to happen in this case rather than just skip it.

asik · 2018-02-27T21:31:36Z

tests/Fix64Tests.cs

+                if (b <= Fix64.Zero)
+                    continue;
+
+                if (b > Fix64.MaxValue / (Fix64)100)


I commented it out and the test passes, what is this for?

Probably a leftover from a previous algorithm I tried that didn't handle large arguments well. Removed.

asik · 2018-02-27T21:39:07Z

tests/Fix64Tests.cs

+                    continue;
+
+                // Reduced precision requirements for very small values
+                double maxDelta = b < (Fix64)0.000001m ? 0.1 : (double)(Fix64.Precision * 10);


Looks like you can get this down to 0.0000003m. The apparent sudden leap in inaccuracy probably has more to do with the lack of test cases than the algorithm itself. However 0.1 is quite inaccurate for a type with almost 10 decimals of precision... Not a deal-breaker though.

The worst delta actually seen from the test data is 0.011, I've adjusted the test to reflect that. My guess would be that the errors of the argument scaling (repeated multiplication with E^4) add up and taint the result?

asik · 2018-02-27T21:50:19Z

tests/Fix64Tests.cs

+            {
+                var b = Fix64.FromRaw(m_testCases[i]);
+
+                if (b < Fix64.Zero)


Same as above, we should test the behavior even if it's an exception.

asik · 2018-02-27T21:51:34Z

tests/Fix64Tests.cs

+
+                if (b < Fix64.Zero)
+                    continue;
+                if (b > Fix64.MaxValue / (Fix64)100)


I don't like that the test simply avoids boundary conditions, here and in the inner loop. What are the issues? If the behavior is acceptable then let's document and test it.

For very large and very small bases the precision becomes bad or terrible, mostly because Ln isn't very precise anymore.
Example:
Pow(0.0000000002328306436538696289, -0.5) = expected 65536 but got 65901.3926025061
Or somewhat pathological cases like this where it really falls apart:
Pow(2147483647.9999999997671693563, 0.9999999997671693563461303711) = expected 2147483637.25622 but got 1426713511.50211

Finally, maybe the error should be measured relatively, instead of absolutely?
Pow(2147483647.9999999997671693563, 0.7911919844336807727813720703) = expected 24173977.2437025 but got 24173942.9738885

I don't see a way to fully fix these cases, but also describing the exact limitations is somewhat difficult. I have refactored the test somewhat and am looking forward to your comments.

So, the test skips half the range [0..MaxValue], and in the range ]1.. MaxValue/2] the delta is artificially reduced to a fraction of itself. I'm not sure I understand the rationale behind this? In any case I'm not comfortable with an operation that's effectively untested over a large range of inputs and boundary conditions (even though I'm guilty of that myself e.g. Tan - I do document that it's not tested well though).

If Ln is the culprit then let's find a more accurate Ln. You got me reading on the topic :) I think an accurate base-2 logarithm for a 64-bit integer would be a good start. https://graphics.stanford.edu/~seander/bithacks.html#IntegerLog has ideas, https://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers as well.

Once we have that, Fix64 is just a long divided by 2^32, and log(x/y) is log(x) - log(y), so log2 of a Fix64 is log2(RawValue) - 32 (right?)

Once we have an accurate log2, (ln x) is just (log2 x)/(log2 e).

asik · 2018-02-27T22:44:19Z

src/Fix64.cs

+        {
+            if (x.RawValue == 0) return One;
+            if (x == One) return E;
+            if (x >= LnMax) return MaxValue;


What are the checks against LnMax and LnMin for? The tests still pass without. Optimization?

Yes. I figured if we know precisely where the argument range of the function ends, might as well test for it, but I guess if the early bailout test rarely succeeds it's a waste of cycles. I'm fine with removing it, your call.

It's a good optimization I think, just wanted to make sure.

asik · 2018-02-28T19:26:39Z

src/Fix64.cs

+            Fix64 result = x + One;
+            Fix64 term = x;
+
+            for (int i = 2; i < 40; i++)


Update: see my comment above on Pow(). We need a different implementation of Ln altogether I think.

Original comment:
It looks like this always converges (at least within the LnMin-LnMax range). I've tested with random inputs for a while. So you don't need to bound it to 40 and can replace with a while loop:

var i = 2; while (term.RawValue != 0) { term = x * term / (Fix64)i; result += term; ++i; }

This gains some accuracy with our test cases, down to Precision * 5; the test can be adjusted in consequence.
On that note, the actual accuracy of this algorithm (with my fix) seems to be 12 decimal digits (from experiments); ideally that's what we would test for rather than an absolute delta. Proper accuracy testing is really lacking in this library, I've should have had a function for doing this a long time ago. Doesn't need to be done in the context of this PR though.

WRT to your comment below: Yes, a better Ln implementation would be great! However, in your comment below, I think it's half of the solution. The code in the links calculates only the integer part of the log2, but we'll still need the fractional part, right?
Wikipedia has an algorithm for that, but I would assume comes with a similar speed and precision issues:
https://en.wikipedia.org/wiki/Binary_logarithm#Iterative_approximation

This one looks promising, I'll have a look at that:
https://github.com/dmoulding/log2fix

asik · 2018-02-28T20:01:54Z

tests/Fix64Tests.cs

+
+                if (b < Fix64.Zero)
+                    continue;
+                if (b > Fix64.MaxValue / (Fix64)100)


So, the test skips half the range [0..MaxValue], and in the range ]1.. MaxValue/2] the delta is artificially reduced to a fraction of itself. I'm not sure I understand the rationale behind this? In any case I'm not comfortable with an operation that's effectively untested over a large range of inputs and boundary conditions (even though I'm guilty of that myself e.g. Tan - I do document that it's not tested well though).

If Ln is the culprit then let's find a more accurate Ln. You got me reading on the topic :) I think an accurate base-2 logarithm for a 64-bit integer would be a good start. https://graphics.stanford.edu/~seander/bithacks.html#IntegerLog has ideas, https://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers as well.

Once we have that, Fix64 is just a long divided by 2^32, and log(x/y) is log(x) - log(y), so log2 of a Fix64 is log2(RawValue) - 32 (right?)

Once we have an accurate log2, (ln x) is just (log2 x)/(log2 e).

sam-vdp · 2018-03-01T11:28:20Z

Thank you very much for pushing back on the crappy previous solution and for your very good idea of using log2. I've added Log2 and Ln based on it, which don't have the previous precision issues.

Furthermore I've removed Exp and added Pow2, which splits the exponentiation into an integer and fractional part and handles large values much better, while maintaining precision for smaller values.

Please note that the Pow test still has some heuristics regarding the expected precision. Firstly, for large results the absolute errors become larger, which needs to be compensated. Secondly, some test cases are slightly pathological. E.g. Pow(1-Fix64.FromRaw(1), Fix64.MaxValue) is evaluated by doing Pow2(Fix64.MaxValue * Fix64.Log2(1-Fix64.FromRaw(1)), which becomes Pow2(Fix64.MaxValue * Fix64.Precision) which of course is not very precise.

asik · 2018-03-01T20:49:58Z

src/Fix64.cs

-        public static readonly Fix64 EPow4 = new Fix64(EPOW4);
-        public static readonly Fix64 LnMax = new Fix64(LNMAX);
-        public static readonly Fix64 LnMin = new Fix64(LNMIN);
+        public static readonly Fix64 Log2Max = new Fix64(LOG2MAX);


We shouldn't add public constants like this unless they serve a clear purpose for users of the library. On that note, PiTimes2, PiInv and PiOver2Inv are bad examples (I really should remove them). "Two" and "Three" don't seem to add any value either. Just declare what you need locally in the function you're using and the JIT will take care of optimizing constant variables away.

asik · 2018-03-01T20:56:29Z

src/Fix64.cs

+        /// <summary>
+        /// Returns the fractional part of the specified number.
+        /// </summary>
+        public static Fix64 FractionalPart(Fix64 value)


I'm not sure about adding that as a public API. It would also require tests. Let's keep this PR focused, 3 new operations is already a lot :)

asik · 2018-03-01T21:01:13Z

src/Fix64.cs

@@ -368,12 +375,18 @@ public partial struct Fix64 : IEquatable<Fix64>, IComparable<Fix64> {
            return x.m_rawValue <= y.m_rawValue;
        }

-        public static Fix64 Exp(Fix64 x)
+        public static Fix64 Pow2(Fix64 x)


I wouldn't make this a public API either, it doesn't seem particularly useful. Pow(2, x) does the trick.

The thing with Pow2 and Log2 is, that I would still like to have dedicated tests for them, not just their combined results. With the tests in the separate assembly, I don't see how that would be possible without making them public?

Good point. Make it internal and add a InternalsVisibleTo attribute to the assembly in an AssemblyInfo.cs file.

That's a neat trick. Done.

asik · 2018-03-01T21:03:04Z

src/Fix64.cs

-            } while ((count++ < 10) && (delta.RawValue != 0));
+            for (int i = 0; i < FRACTIONAL_PLACES; i++)
+            {
+                z = z * z;


FastMul would be a big performance win here.

asik · 2018-03-01T21:05:02Z

src/Fix64.cs

            {
-                term = x * term / (Fix64)i;
+                term = x * term * Ln2 / (Fix64)i;


I haven't tested it but it looks like term should remain smaller than MaxValue, so FastMul could be a big perf improvement here.

asik · 2018-03-01T22:51:48Z

tests/Fix64Tests.cs

        {
-            double maxDelta = (double)(Fix64.Precision * 10);
+            var test = (double)Fix64.Pow2((Fix64)3);


asik · 2018-03-01T23:08:53Z

tests/Fix64Tests.cs

@@ -428,12 +445,12 @@ public void Pow()
                    {
                        var expected = e == Fix64.Zero ? 1 : b == Fix64.Zero ? 0 : Math.Min(Math.Pow((double)b, (double)e), (double)Fix64.MaxValue);


This is hiding at least one bug: it ignores the case Pow(0, negative power) which should give MaxValue (the closest we have to infinity).

Mathematically, 0^x (x<0) is undefined, so I would argue 0, MaxValue or Infinity are pretty much equally wrong. Windows Calculator for example returns 0. I've added an exception for this case, after all, the division operator also throws instead of return MaxValue.

asik · 2018-03-01T23:09:17Z

tests/Fix64Tests.cs

-                            continue;
+                        // Absolute precision deteriorates with large result values, take this into account
+                        // Similarly, large exponents reduce precision, even if result is small.
+                        double maxDelta = Math.Abs((double)e) > 100000000 ? 0.5 : expected > 100000000 ? 10 : expected > 1000 ? 0.5 : 0.00001;


Actually if the "real" Pow(x, y) > MaxValue then you expect exactly MaxValue with 0 delta.

That would be nice, but actually the Pow approximation is still subject to the limitations that increase the maxDelta to 0.5 before. I.e. even if the expected value is MaxValue, the actual result may be slightly below, but inside the 0.5 delta.

asik · 2018-03-01T23:11:36Z

src/Fix64.cs

            if (neg) result = One / result;

            return result;
        }

-        public static Fix64 Ln(Fix64 x)
+        public static Fix64 Log2(Fix64 x)


I wouldn't expose Log2 as public API. A more useful API (which can be a separate PR) would be Log(a, newBase) like System.Math offers.

asik · 2018-03-01T23:22:17Z

src/Fix64.cs


+        public static Fix64 Ln(Fix64 x)


All public APIs need a ///summary. I try to stay close to System.Math in descriptions. Also discussion of accuracy if it's not accurate down to Fix64.Precision (see for ex. Sin in the latest master - I just updated it yesterday).

…e appropriate, add documentation to public APIs, add DivideByZeroException to Pow(0, -x)

asik · 2018-03-11T18:43:42Z

I've made a few stylistic changes locally and merged. Thank you for this contribution!

Add Pow, Ln, and Exp functions

59972b0

sam-vdp mentioned this pull request Feb 27, 2018

Could you add Pow()? #10

Closed

asik suggested changes Feb 27, 2018

View reviewed changes

Fix range check on Ln(0), refactor tests

08a6e30

asik suggested changes Feb 28, 2018

View reviewed changes

Bartl added 2 commits March 1, 2018 12:17

Add Log2, Pow2 algorithms, base Pow on them

835dfdb

Include log2fix in LICENSE.txt

59b2473

asik suggested changes Mar 1, 2018

View reviewed changes

Bartl added 2 commits March 2, 2018 12:54

Removed public constants, remove FractionalPart API, use FastMul wher…

1863296

…e appropriate, add documentation to public APIs, add DivideByZeroException to Pow(0, -x)

Make Pow2 and Log2 internal

c222903

asik closed this Mar 11, 2018

		@@ -428,12 +445,12 @@ public void Pow()
		{
		var expected = e == Fix64.Zero ? 1 : b == Fix64.Zero ? 0 : Math.Min(Math.Pow((double)b, (double)e), (double)Fix64.MaxValue);

Add Pow, Ln, and Exp functions #12

Add Pow, Ln, and Exp functions #12

Conversation

sam-vdp commented Feb 27, 2018

asik left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

asik Feb 28, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

asik Feb 28, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

asik Feb 28, 2018 • edited Loading

Choose a reason for hiding this comment

sam-vdp commented Mar 1, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

asik Mar 1, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

sam-vdp Mar 2, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

asik commented Mar 11, 2018

asik Feb 28, 2018 •

edited

Loading

asik Feb 28, 2018 •

edited

Loading

asik Feb 28, 2018 •

edited

Loading

asik Mar 1, 2018 •

edited

Loading

sam-vdp Mar 2, 2018 •

edited

Loading