-
Notifications
You must be signed in to change notification settings - Fork 0
Possible Additions for 1.1.0
The following are possible additions for the next update. Not all entries are guaranteed, and the list is not exhaustive.
For my own sake, I'm going to put (completed) after each change I've created. It doesn't mean the item has already been added to the newest release, or even that it's 100% bug tested, but it does mean that that step has been completed, and might give you an indication of my progress with this update.
double sq(double) (completed)
So you can more easily square doubles.
double cub(double) (completed)
So you can more easily cube doubles.
long nCr(int, int) (almost completed)
So you can perform the combination function without worrying about overflow or divide by 0.
Complex[] fsincos() (completed)
Computes the sine and cosine simultaneously. This will be added because both cos and sin require computing sin(re), cos(re), sinh(im), and cosh(im). If you compute the cosine and sine of a complex number, you compute each of those values twice. So why not create a function that computes those values once, then uses them to compute the complex cosine and sine simultaneously?
double ulpMin() (completed)
double ulpMax() (completed)
Computes the ulp (unit in the last place) for a given complex number. For real numbers, the ulp is essentially the smallest positive number you can add to the inputted number and still get a different number. If you add anything smaller, the difference won't be recognized. Ulps are always powers of 2. For complex numbers, there are 2 ways to compute the ulp: either take the ulp of the maximum of the real and imaginary components, or take the ulp of the minimum of said components. If a number is at least as big as the former, adding it to the input is guaranteed to yield a different number. If a number is smaller than the latter, adding it to the input is guaranteed to yield the same number. Both the real and imaginary parts usually have different levels of precision, thus they usually have different ulps, thus there must be two different ulp functions. If you're wondering why I don't just have ulpReal and ulpImag, the answer is because that could be easily compressed in one line (Math.ulp(z.re), Math.ulp(z.im)), whereas the min and max ulp are much more complicated to compress (Math.ulp(Math.max(Math.abs(z.re), Math.abs(z.im))), Math.ulp(Math.min(Math.abs(z.re), Math.abs(z.im)))).
double scalb(int scaleFactor) (completed)
Returns the product with the instance and 2^scaleFactor.
The main additions to Cpx2 will be for functions related to the error function. To make a long story short: They're helpful functions for certain calculus problems. I was originally going to add some stuff about how stupid the definitions are of the error functions and all the illogical steps made in their conception, but I'm just going to save that for the article on the functions themselves, when this update is released. To avoid confusion, I'm going to lie and pretend that the erf function is the integral of e^-x², instead of 2/√(π)*e^-x².
Computes what's known as the error function. This function was originally conceptualized to compute the area under the bell curve. The bell curve (μ=0, σ=1) can be plotted out by y = e^(-x²/2)/√(2π). In order to find the area under the curve, we simply take the integral of that function. The error function can be used for this purpose, as it computes the integral of e^-x² dx, evaluated from 0 to z. If we scale the input by 1/√(2), and scale the output by 1/√(π), we get the area under the bell curve. It's called the error function because calculating that area calculates the chance that there isn't an error; the chance of being within that range.
And guess what? Integrals still work for complex numbers, so this function works for complex numbers as well! To see why this can be useful, read the entries below.
Complex erfi(Complex z) (completed)
Computes what's known as the imaginary error function. This function is to the error function what sinh is to sin, or what tanh is to tan. It's essentially just erf(z*i)/i. Just as the erf function computes the integral of e^-x², the erfi function computes the integral of e^+x². How is this useful?
...
I don't really know. It doesn't really have any specific applications, like with the error function, it's more or less just an integral that happens to come up every once in a while, and can't be computed using elementary functions. So it's good to have around, at least.
Complex fresnelC(Complex z) (completed)
Complex fresnelS(Complex z) (completed)
Computes what are known as the Fresnel integrals. Technically, the notations for these functions are C(z) and S(z), but that notation is a bit ambiguous, so we're putting frensel before it. Fresnel C computes the integral of cos(x²), while Fresnel S computes the integral of sin(x²). Due to Euler's identity, these are equivalent to the integrals of (e^(x²i)+e^(-x²i))/2 and (e^(x²i)-e^(-x²i))/(2i). As such, these can be computed using the error function, with inputs proportional to √(i) and √(-i). These integrals are useful, as plotting the parametric x(t)=C(t), y(t)=S(t) plots out what's known as the Euler spiral. Which, apart from just being pretty to look at, is also the theoretical path of motion for a car that's moving forward while turning at a constant angular acceleration (not to be confused with a circle, the path of a car turning at a constant angular velocity).
Complex erfc(Complex z) (completed)
What's known as the complementary error function. While the error function yields the area between 0 and z under the bell curve, the complementary error function yields the area between z and ∞. This exists primarily because it allows you to compute the error function more precisely for values far from 0.
Complex invErf(Complex z)
Complex invErfi(Complex z)
Complex invErfc(Complex z)
The inverse of the error function (as well as the imaginary and complementary counterparts). The inverse error function tends to be handy in probability theory, as it tells you how many standard deviations away from the mean you need to travel before you achieve a certain level of certainty. The inverse complementary error function is also useful, as it allows you to more precisely describe the probability of being outside a certain range. The inverse imaginary error function isn't useful in probability theory, but is there anyway because it can be easily calculated using the inverse error function. I can't make inverse Fresnel integrals, though, since those are linear combinations of error functions, and would need their own dedicated functions. I might end up having to make multiple versions of the inverse error function, since it's a multivalued function over the complex plane, so you would need to evaluate it along certain branches.
Other error-related functions
The erfcx is the erfc times e^x². It exists to avoid arithmetic underflow for sufficiently large inputs.
There are two Dawson functions. One that's the error function times e^x², the other that's the imaginary error function times e^-x².
The cumulative distribution function, denoted Φ(x), computes the area under the normal (μ=0, σ=1) bell curve from -∞ to x, without any stretching, scaling, or translating. (completed)
The Faddeeva function is equal to erfcx(-ix). The output is complex valued for most real inputs. I don't know what this function is used for, but it's here anyway.
The error functions have already been written, and appear to be very accurate for all inputs I throw at it. The inverse error functions haven't been written yet. I made a prototype of the inverse error function before, but it hasn't been RIGOROUSLY tested yet. Also, the current implementation of the error function may change. Because the current implementation I have takes 45 iterations to converge, and I'm looking into seeing if there are any more efficient algorithms.
If you feed exactly one input into the logSum function, the result will always be 0. It's supposed to be the logarithm of that one input, but instead it's 0. This is a minor oversight, and is pretty easy to fix. (completed)
If you plug in the Gudermannian function for values close to an odd multiple of πi/2 OTHER THAN ±πi/2, the answer will be wrong. This is due to an oversight, where I told the function to approximate the cosh as z±πi/2. Instead, it should do that to z2, the computed value for the adjusted value of z MOD 2πi restricted to the range (-πi,πi]. Just like the previously mentioned bug, this is pretty easy to fix. (completed)
The div function will need to be adjusted. As it is right now, when you divide a small number by a subnormal number, you get ∞. Sometimes this is right, sometimes it's wrong. What it's doing is computing the reciprocal of the subnormal number (Infinity) and multiplying by that small number. I didn't see this as a problem to get worked up over, at first. However, I currently want the Complex class to uphold a specific relationship with doubles. Namely, any equation with doubles as the inputs should yield the exact same result if real-valued Complex numbers are the inputs, except when the former yields undefined (i.e. sqrt(-1) vs. sqrt(-1+0i)). In order to preserve this property, I must fix this bug. (completed)
Possibly in an upcoming patch, I might have to change Bernoulli to either bernoulli or BERNOULLI. Depending on if I plan on making it a constant array. I mean, making an array final doesn't prevent people from editing it, it just prevents them from changing its size. For that reason, I might even have to change it to be a private data member and add a getter function to access them.
The equals method throws a NullPointerException when comparing to a null Object. This is because the method asks for the getClass and asks if it's Complex.class. This is very easy to fix, as it just requires instead asking if (obj instanceof Complex).
Cube roots are not as accurate as they should be. When you take the cube of 2+i, you get 2+11i. However, if you take the cube root of 2+11i, you get 2+1.0000000000000002i, or 2+(1+2^-52)i. IEEE 754 specification requires that the result of any non-transcendental function be off from the true result by 1/2 an ulp or less. This is because such methods, if they are just slightly off, can always be pushed towards the true result by iterative methods, while transcendental functions, due to the table-maker's dilemma, can never be guaranteed to be rounded to the nearest ulp, and there is no way of knowing how many extra digits it could take to make sure the result is correctly rounded. This property holds just as much true for real numbers as it does for complex numbers. The cube root is non-transcendental. If you have an okay approximation of cbrt(x) = c, you can get an even better approximation by replacing c with (2c+x/c²)/3. For instance, if we have x = 2+11i, and c = 2+1.0000000000000002i, we find that 2c = 4+2.0000000000000004i, x/c² = 2+0.9999999999999996i, 2c+x/c² = 6+3i, and that divided by 3 = 2+i. Of course, in practice, we'd actually multiply by 0.3333333... at the end, instead of dividing by 3.
Division also seems to somewhat break the IEEE 754 standard on occasion, often being an ulp off from the true answer. Unfortunately, nothing I do seems to fix this issue. Even using the Newton-Raphson method, replacing r with r+r(1-xr), doesn't seem to help that much. I fear I might have no choice but to simply lower my standards for precision on this library.
I might add "import Math" to the Mafs class. Or I might not. On one hand, it makes it way easier to call math functions of doubles. On the other hand, it could potentially override other Processing functions. I'll have to look into it.
I plan on making the code for the polygamma reflector function better organized. Simply because I wrote that code a long time ago, and can't really understand what's going on in it. The only reason I didn't organize it before submitting this library is because it requires computing the n-th derivative of the cotangent, which is quite a hefty bit of work. And to reorganize it would require looking through the notebooks of math work I have, find which one(s) have that particular calculation in it, and try to reverse engineer how the algorithm I'd see got me to that code result. Then make sure there aren't any exceptions...
The factorial(int) function might start throwing an ArithmeticException when the input is negative after the next update. Or it might not, since returning the maximum long integer guarantees that division by it will always result in zero, unless the numerator is either the maximum long integer, or it's a floating point.
The abs, inv, and log functions will be changed so that, in the case of overflow/underflow, instead of dividing by the lazy absolute value, performing the function, then readjusting based on that change, we instead divide/multiply by 2 to the power of 1022, then adjust for that change. It's more efficient this way, since the scalb function with ±1022 as a scale factor is more efficient than division, and it's more accurate since multiplying by an exact power of 2 is an operation guaranteed not to change the precision of the end result (except when doing so gives you ∞ or a subnormal number, which cannot happen since we're being taken closer to 1). (completed)
I'm trying to make the loggamma, digamma, and polygamma functions more accurate...somehow. If you try adding the Mascheroni constant to digamma(1), you don't get 0, even though you really, really, really should. This kinda diminishes most of the purpose of the digamma function if adding GAMMA to it doesn't return the harmonic sum up to that point, it just returns something kinda close to it. (VERY DIFFICULT, fixing this is almost certainly going to be the majority of the work I put into this update)
Part of the solution to the aforementioned problem might to be to add more Bernoulli numbers. That might make things easier for those of you using this to develop extreme math algorithms. All three of you. Or, the solution might actually be to loop through fewer Bernoulli numbers. Technically, the Stirling approximation diverges eventually, since the Bernoulli numbers B_n go up at a rate of O(n!/(2π)^n), it just takes them a while to actually start increasing. Namely, it takes up to B_14 before they start increasing.
Speaking of Bernoulli numbers, I might also add the Bernoulli polynomials to the library. You might end up using them, and also they will come in handy when computing the polylogarithms. Which should probably come around in update 1.3.0? After the error functions in 1.1, and the Riemann zeta functions in 1.2. I think. They're not anything big, but they do allow you to compute the K function and Barnes G function, which were both things I wanted to add with the gamma functions, but I couldn't, because I needed to add polylogarithms first! Polylogarithms are a requisite to calculate the Barnes G function.
I might also (possibly?) have to make the gamma function more accurate. It uses the Lanczos approximation, which is pretty good. But it appears to be way off for certain values. For instance, gamma(new Complex(9.9375)) is off by 13 ulps. If I were to fix this, I'd have to learn how to generate Lanczos coefficients (you caught me, I stole the coefficients currently in use from Wikipedia) to a suitable accuracy. For some reason, every time I try doing it, the coefficients are way off. Then again, that was a while ago, before I knew how to use BigDecimal numbers, so I might as well give it another shot.
As of right now, my plans are:
1.1.0: Add the error function and related functions
1.2.0: Add the Riemann zeta function and related functions
1.3.0: Add the polylogarithms, the Clausen functions, the K-Function, and the Barnes G-Function
1.4.0: Add exponential integrals
1.5.0: Add elliptic integrals (complete and incomplete, the incomplete Π probably not included)
1.6.0: Add Bessel functions? (Including Airy functions)
1.7.0: Add the incomplete gamma function? And the incomplete beta function?
1.8.0: Add the Hurwitz zeta function?
1.9.0: Add the Lambert W function?
1.10.0: Add inverse elliptic functions? And the elliptic Π function? I have no idea when that function will be ready. Possibly earlier than this, possibly later.
1.11.0: The hypergeometric functions? I guess? If I do make this, it'll only be 2F1, since none of the other hypergeometric functions are well documented.
1.12.0: The scorer's functions? I guess? I don't know.
Honestly, I have no idea when any of these will be ready past 1.5. To be completely honest, I have no idea how long I plan on releasing updates for this library, at some point I might just give up and decide I have enough functions to be whole. Or I might decide using what I have here to develop a Complex Matrix library might be much more worth my time. I already have the Matrix Library started and understand the algorithms. That said, I should at least be able to make it up to elliptic integrals, update 1.5.0.
You can find more info about these functions and others at the official Wikipedia list of functions with their own name.
Before someone asks, no, I will not make tetration. Trust me, tetration is extremely difficult. Extremely, extremely difficult. If I ever end up doing tetration, it'll probably be for only a few bases, unless someone pays me a bunch of money to do all the math no one else seems to care to do to make it work for all bases. Which, who can blame them? Tetration has pretty much no application whatsoever. And pentation is even more useless, and I'd only do that if someone paid me a LOT of money. Because while tetration is really only hard if the power is a non-integer, pentation is hard if the base OR power is a non-integer. When the base is a non-integer, you have to use all the hard math it took just to perform tetration on non-integers, and when the power is a non-integer, you have to invent even MORE math to make the derivatives continuous. Sextation is even harder, since you have to invent even more math piled onto the pentation math, and so on and so forth. And most bases have an overflow for nearly all powers. And, on top of all this, complex number inputs!!!