Exp series error term help (probably easy for people who know calc)

[Latin4567]

I'm extending j#'s built in BigInteger and BigDecimal classes to have proper operators and standard math functions like sqrt, log, exponentiation with decimal powers, trig functions, etc, in Vb.net... I already have logs and nth roots and nth powers done, and I implemented the super fast karatsuba multiplication algorithim, and a very very fast factorial algorithim....

I need a FAST Exp(x) (e^x) routine so I can define exponentiation with decimal/irrational powers.

I'm using the series deffinition of e^x

e^x = SUM TO INFINITY OF x^n/n!

It works, but I can't figure out the exact correlation between number of iterations and error... like I need to know, given an X value of relatively arbitrary magnitude (within reason), and a specified precision to a certain number of decimal points, how high of an N value I will need to go before the series converges on a correct sum to that degree of precision

please advise

[jemidiah]

I would be shocked if that series definition converges more quickly than some clever algorithm. Off the top of my head you could use exponentiation by squaring to quickly find e^n's, perhaps using the series for the fractional part of the exponent. For small x (< 1 in magnitude) the series would converge quite quickly.

Lemme think... for |x| < 1, take the ith partial sum Pi (sum from n=0 to i-1 of x^n/n!) which must itself be less than the sum from n=0 to i-1 of 1^n/n! = sum of 1/n!. The correction term is similarly bounded and gives...

e^x - Pi = sum from n=i to infinity of x^n/n! < sum from n=i to infinity of 1/n!

That last sum can be bounded nicely since

sum from n=0 to infinity of 1/n! = e^1 = e = sum from n=0 to i-1 of 1/n! + sum from n=i to infinity of 1/n!

which implies

sum from n=i to infinity of 1/n! = e - sum from n=0 to i-1 of 1/n!

Put together, we have this error term when |x| < 1:

e^x - Pi < e - sum from n=0 to i-1 of 1/n!


I know this doesn't give you the N you wanted, but it does the job--presuming you can compute e to whatever precision you desire. You could loop through computing the partial sum, and each iteration you could update the error term until it's close enough. There are more complicated improvements you could make for very small x, where the error term would be greatly overestimated. Yeah you could use a geometric series to bound it in those cases decently well. Combined the algorithm I'm imagining would do very well; log(x)*[order of squaring and multiplication algorithm] time for the integer part and in the long run well over a digit per iteration for the fractional part of x.

[Latin4567]

Thanks, that should work great -- I actually have a table lookup setup for e and pi with 2 million digits of e and pi built in (LOL). Basically I have an e(digits) function, and whenever a specific value of digits is called, it creates a cache of the BigDecimal variable containing e to that amount of digits instead of substringing the stored e to 2 million digits variable every time. It runs very fast... when I time it there is no detectable difference between using it and using a const



While I still might attempt getting this series to work, I noticed something on wikipedia that I think might point me in a faster direction... on the page for exponentiation, it briefly mentions that "Rational number exponents can be defined in terms of nth roots, and arbitrary nonzero exponents can then be defined by continuity," but it doesn't elaborate. My nth root function is very fast because it uses newton's method, so I would be very interested in defining exponentiation with rational powers. Any idea how this would be done?

thanks for your help, btw

[jemidiah]

That quote is referencing the fact that x^(1/n) = nth root of x. You could use this to find e^x by finding e^r for some rational number r=a/b very close to x by noting e^r = e^(a/b) = (e^(1/b))^a--that's what it means "by continuity". You can find r basically by binary search. Estimating the error would be somewhat difficult. I would imagine you would end up with a very large numerator and denominator for large, irrational x which may or may not be a problem.

[Latin4567]

I'm trying to define exponentiation. I only have the ability right now to do nth integer powers of integers or decimals, and exponents with decimal or integer bases and integer powers. I need to define exponentiation where the base and the power can both be decimal. The base and power shown are just two random decimals.

here it is in mathematical notation:


So is there a quick change I can make that will make this more computationally efficient (it's pathetic at the moment, unless the decimals are short in which case it takes like 10 ms), or have I reached a dead end? I suspect, if I have reached a dead end, that I will have to look into implementing Egyptian fractions (http://en.wikipedia.org/wiki/Egyptian_fraction)

Here is my current function:

Shared Function RationalPower(ByVal base As BigDecimal, ByVal power As BigDecimal) As BigDecimal
        Dim k As BigInteger = 10 ^ New BigInteger((DecimalMag(power) - 1))
        Dim kx As BigInteger = k * power
        'Console.WriteLine(k.ToString)
        'Console.WriteLine(kx.ToString)
 
        Dim h As BigInteger = Gcd(kx, k)
        ' Console.WriteLine(h.ToString)
 
        Dim m As BigInteger = kx / h
        'Console.WriteLine("m = " & m.ToString)
        Dim n As BigInteger = k / h
        ' Console.WriteLine(n.ToString)
        ' Console.WriteLine("nth root(" & base.ToString & " ^ " & m.ToString & ", " & n.ToString & ", " & BigDecimal.GlobalScale.ToString & ")")
 
        Return func_internal.nth_root_itr(base ^ m, n, BigDecimal.GlobalScale)
    End Function

Here is the nth_root function, which is very fast and isn't the problem (the problem is the moronic parameters I'm sending it):

Shared Function nth_root_itr(ByVal x As BigDecimal, ByVal n As BigInteger, ByVal prec As Integer) As BigDecimal
        If n = 1 Then
            Return x
        End If
        'calculates x^(1/n)
        'f(x) = x^n - x
        'f'(x) = n*x^(n-1)-1
        'x1 = x0 - f(x0)/f'(x0)
        '(f^n-x)/(n*f^(n-1)-1)
        Dim oldscale As Integer = BigDecimal.GlobalScale
        Dim oldroundingmode As RoundingMode = BigDecimal.GlobalRoundingMode
        Dim oldscalingmode As ScalingMode = BigDecimal.GlobalScalingMode
        BigDecimal.SetGlobalScale(prec, ScalingMode.SCALE_ADAPTIVE)
        BigDecimal.GlobalRoundingMode = RoundingMode.ROUND_HALF_EVEN
        Dim f As BigDecimal
        f = 10 ^ New BigInteger(CType(x, BigInteger).ToString.Length / n - 1)
        Dim squared As BigDecimal
        squared = 0
        Dim oldf As BigDecimal = 0
        Dim oldvalues As New List(Of String)
        Do
            oldf = f
            f = f - (f ^ n - x) / (n * f ^ (n - 1) - 1)
            If oldvalues.Contains(f.ToString) Then
                Exit Do
            End If
            oldvalues.Add(f.ToString)
            ' Console.WriteLine(f.ToString)
 
            If oldf = f Then
                Exit Do
            End If
        Loop
        BigDecimal.SetGlobalScale(oldscale, oldscalingmode)
        BigDecimal.GlobalRoundingMode = oldroundingmode
        Return f
    End Function

Here is the Gcd function, which also shouldn't be posing a problem:

Shared Function Gcd(ByVal a As BigInteger, ByVal b As BigInteger) As BigInteger
        If b > a Then
            Dim tmp As BigInteger
            tmp = a
            a = b
            b = tmp
        End If
        If b = 0 Then
            Return a
        Else
            If a = 0 Then
                Return b
            Else
                'Console.WriteLine("Return Gcd(" & a.ToString & ", " & a.ToString & " Mod " & b.ToString & ")")
                'Console.ReadLine()
                Return Gcd(b, a Mod b)
            End If
        End If
    End Function

here is the DecimalMag function, which, for all intensive purposes, returns the number of numbers after the decimal point:

Shared Function DecimalMag(ByVal x As BigDecimal) As Integer
        Dim n As Integer = 0
        Do
            x = x * 10
            n = n + 1
            If x = Round(x) Then
                Exit Do
            End If
        Loop
        Return n
    End Function

The BigInteger and BigDecimal classes are custom , and act as wrappers for the J# BigInteger and BigDecimal classes. I wrote / am writing custom operators and functions so they integrate nicely into vb. I replaced the slow integer multiplication routine with a routine that let's the super fast karasuba multiplication algorithm kick in only when the numbers are of a high enough magnitude for there to be a performance gain.

I'll be willing to post the whole library if everything works out when it's done. My goal is to be able to treat BigInteger and BigDecimal exactly the way you would treat normal number types in .NET, and I'm well on the way. Exponentiation with rational powers is the last step -- it will allow me to do the Exp function (I could use the taylor series but that is very slow). I already have a very fast routine for the sin, cos, and tan functions. It's also supposed to be able to do exp, cosh, sinh and tanh, but it fails on higher magnitudes. The only reason sin, cos and tan work great is because I am shifting everything to the first period before I send anything to the function (which didn't write, and thus don't really understand at all).

Alternatively to fixing my rational power routine, if anyone has any tips that will make this work for exp with larger values, here is the function that I didn't write that I am wrapping:

Shared Function transcendental(ByVal p As TransType, ByVal n As Long, ByVal x As BigDecimal) As BigDecimal
        Dim k As Long
        Dim r As BigDecimal
        Dim s As BigDecimal
        If p <= 3 Then
            r = -x * x 'trig
        Else
            r = x * x 'hyperbolic
        End If
        s = 4 * n + 2
        k = n
        Do
            k = k - 1
            If k = 0 Then
                Exit Do
            Else
                s = 4 * k - 2 + r / s
            End If
        Loop
        Dim pmod As Integer
        pmod = p Mod 4
 
        Select Case pmod
            Case 0
                Return (s + x) / (s - x) 'exp
            Case 1
                Return 2 * x * s / (s * s - r) 'sin, sinh
            Case 2
                Return (s * s + r) / (s * s - r) 'cos, cosh
            Case 3
                Return 2 * x * s / (s * s + r) 'tan, tanh
            Case Else
                Return (s + x) / (s - x) 'shound't get called, return exp by default
        End Select
    End Function

and here is my wrapper for it. I have it cut off at 31 because beyond that it begins to fail. The function didn't come with any pointers about how high of an n value you need for different orders of magnitude, so what you see below is my attempt to get it to work (this works the best from what I've tried, but there is probably a much more specific correlation between n and the magnitude of the inputs):

Shared Function Exp(ByVal x As BigDecimal) As BigDecimal
        If x > 31 Then
            Return Round(e(BigDecimal.GlobalScale) ^ Round(x))
        Else
            Return func_internal.transcendental(func_internal.TransType.Exp, BigDecimal.GlobalScale * 2, x)
        End If
    End Function

anyway thanks to anyone who can help lead me to properly defining exponentiation with decimal powers!

FYI: I just graduated from high school, and I have only taken pre-cal. While I can figure a lot of stuff out myself (I got newton's method to work for nth powers using only the wikipedia page on newton's method and mathematica's built in derivative calculator), assume that I will require extra explanation to understand anything beyond pre-calculus


by the way, as an interesting curiosity, I believe I have developed an extremely efficient factorial algorithm. You can even do it on paper for small-ish values of n with relative ease... it's based on the principal of keeping the magnitude of all multiplications as small as possible for as long as possible. Most of the magnitude gets saved for the last 3-4 multiplications.

it's easier to demonstrate visually than to explain... here is the process for calculating 10!

start with a list of 1 through 10:
1 2 3 4 5 6 7 8 9 10
instead of multiplying through stupidly, combine the first term with the last, second with second to last, third with third to last, etc,
10 18 24 28 30
Since there are an odd number of items, add one to the beginning
1 10 18 24 28 30
now repeat the previous steps again and again until you are left with 1 number

so in the end it looks like this:
1 2 3 4 5 6 7 8 9 10
1 10 18 24 28 30
1 30 280 432
432 8400
3628800


Which is extremely efficient compared to the naive routine, which would look like this

1
2
6
24
120
720
5040
40320
362880
3628800

As an example benchmark, my routine calculates 1000!, which has 2568 digits, in 35 milliseconds, whereas the traditional method takes 20 or so minutes (even with karasuba multiplication), and Stirling's approximation takes 433 ms. Of course, it's possible my implementation of Stirling's approximation is horribly inefficient. Suprisingly, my routine is more efficient, even though it involves system.List(Of BigInteger)'s with thousands of items. I'll provide the function if/when I finish / release my library. I'll probably post it on codeplex or something (of course I'll have to eliminate / get permission for some of the code that I didn't write, but most of it I did write)

[jemidiah]

I would have to think through that transcendental routine for a long time to figure it out without any comments. At first glance it looks like some iterative approximation method. I don't think it would be worth the time, though; if it fails at higher magnitudes it's probably only meant to do a good approximation near zero.

How are you doing integer exponents? That is, how is the "base ^ m" call computed? Exponentiation by squaring has fantastic performance gains for gigantic numbers, perhaps you could try that. Since you're also taking f^n and similar in the nth root routine, a faster exponentiation algorithm might be the thing you need.

If that ends up not working you could try my original idea, which I'll flesh out slightly more:

Compute x^y for x and y decimals as follows. Let y = int(y) + frac(y). Compute x^int(y) through exponentiation by squaring. Compute x^frac(y) through the Taylor series, which converges very quickly for y near zero--the only somewhat slow case is for frac(y) ~= 1, though even that should converge quite quickly since the factorial would "take off" after a few iterations. Multiply your numbers to get x^int(y) * x^frac(y) = x^(int(y) + frac(y)) = x^y. Done.


Also, I liked that factorial computation. I'm surprised the performance gains are as large as you've gotten but I suppose that's the magic of orders of algorithm complexity.

[Latin4567]

ah, yes this would explain why my *ahem* implementation of normal exponentiation is slow.... I'm literally multiplying by the base n times.....

I'll take a look at that article!


edit: wow that's a major improvement. I feel like such an idiot now. This is also definitely the reason why my Stirling approximation was so slow

[Latin4567]

my brain hurts trying to think how I would express my factorial algorithm in mathematical notation

[Latin4567]

ok now everything is fast but my log routine (which was horribly inefficent and can't really be salvaged) is slowing everything down. I'm about to begin the search for a more efficient algorithm, but any suggestions on what would be the fastest? I'd imagine log2 would be fastest to compute, and once I have that I have all bases

[jemidiah]

Notation for your factorial algorithm would be a little cumbersome to write up. You'd need two indicies, one tracking which level you're at and one tracking which number in that level. But, the recursive definition from there isn't that bad.

Either log2 or ln would probably be best to implement, depending on the algorithm you find. There's of course a Taylor series for them, and the one for ln happens to work out nicely. As for algorithms, I don't actually know. I can imagine my own algorithm where I would repeatedly divide by 2, 2^(1/2), 2^(1/4), ... and get the logarithm that way but it wouldn't be efficient for incredibly small numbers. If you wanted 1000 binary bits of precision you'd need to evaluate 2^(1/2^1000). I'm sure you can find something better than that online, though.