I see this as being in the form of z = a / b * c / d; considering the brackets.
If I am correct then I should be able to reform this into z = (a * c) / (b * d), whereupon I should be able to cobine the brackets and reduce the variable occurences.
I agree with all of your suggestions. However, I wonder what exactly 1d1**bkap does in that expression? That term looks odd enough that I might refer back to the original derivation to see if it is correct.
BTW, I would probably never write the expression as a / b * c / d. I would always group terms together as in the examples shown so that the order of operations is clear to a human reader. Also, unless there is some issue with overflow or underflow, I would tend to use an expression that had a single division in it, such as (a*c)/(b*d). Then I might look for cancellations of factors in the numerator and denominator, or other simplifications. For example, (a*c) has a factor that can be computed as temp**4, and there are two constants that could be multiplied together.
Are all garynewport’s names of variables or constants double precision = kind(1d0) scalars? Does the reordering cause overflows or underflows? Why use rho**5d-1 instead of sqrt(rho), which would run faster?
Truthfully, most compilers would recognize that and substitute a sqrt call. For the things such as temp**3d0, the compiler can PROBABLY convert this to a series of multiplies - it would have a better chance if the exponent was integer.
The reason I mentioned this originally is that I often see constants written as 1.234*10.0**5 by new programmers rather than the simpler (to me) 1.234e5. In this case, there are constants in the a and c terms that could be combined with that exponent, provided bkap is another constant with an integer value and not an actual floating point variable. It might well be that a compiler will convert this at compile time, but it just looked odd to me.
Thank you for the replies; I am processing them at the moment.
The equation has already been tested whilst broken down (first a/b * c/d) then mixed. It was originally in the form a0 * a1, with a0 equating to a/b and a1 being c/d. I have recombined them to reduce some of the single-use variables and identify what each function is doing.
I thought I had read somewhere that the power was to be floating; hence why I have gone through adding the d0 (x**2d0). Since this is not the case, I am far happier going for the value only. I assume this is true even if the power is a decimal instead of a float.
bkap is a calculated value; it relates to the opacity of a section of a star.
I removed sqrt and replaced with 5d-1 (soon to be 0.5) to aid in my reading of the equations; it helps me identify what each one is doing. I have also done it to reduce the powers elsewhere(where the original code had something like sqrt(a * b**4 * c**2 * d**8). If sqrt is faster I will revert back eventually.
Only a few constants have been left in (PI, SIGMA, G, C, for example) and, again, these are in to aid readability. These will remain to help others understand the equations being used and, if necessary, improve the code further.
Thought clarity might help (and it seemed rude not to reply).
What do you mean by “decimal instead of a float”?
Personally, I would favour writing a power of 10 like 10^b as 10**b and not 1d1**b.
A quick inspection with compiler explorer shows that sqrt(a) * b**2 * c * d**4 gives you the smallest number of instructions. From the latter and the fact that the numbers are smaller I would also expect this form to incur smaller round-off errors.
In a C++ code base I’ve worked with in the past they used the following recursive template code for positive integer powers:
/** Compile time integer power, returns `base` raised to power `exponent`. @ingroup utils */
template <unsigned int exponent>
inline double ipow(double base) {
return ipow<exponent - 1>(base) * base;
}
/** Compile time integer power (base case 0) @ingroup utils */
template <>
inline double ipow<0>(double) { return 1; }
/**
* Compute non-negative integer power `b^e`. This function is usually faster
* than `std::pow` for `e < 50` and matches the `std::pow` performance at
* approximately `e = 100`. For compile time constant `e` this function
* use @ref ipow.
* @ingroup utils
*/
template <typename T>
inline T ipow(T b, int e) {
T r = 1.0;
while (e > 0) {
r *= b;
--e;
}
return r;
}
I’m glad that Fortran compilers take care of such things automatically. Besides avoiding the mess of different forms like ipow<4>(d) or ipow(d,4) or std::pow(d,4) along with their performance/accuracy tradeoffs most developers including myself aren’t familiar with, the Fortran ** is type generic and doesn’t require any explicit import.
(Here’s the same expression using the C++ ipow templated function given above… not pretty.)
Typing error. I meant “decimal instead of integer”; so x**0.3333 rather than x**3.333d-1
Again, I assumed that I needed to use the floating point throughout. This arose in part from using the warning all notice in the compiler at some point and it indicating I needed to use floats more often. Therefore, whereever a double precision was involved, I made sure all values were presented as such.
You live and learn!
In terms of the sqrt side, it is in my notes to revert back from x**0.5 once I have recognised all the parts of the code (or the relevant ones, anyway).
Just as a side note, bkap is simply kappa with a b stuck in front because the original programmer allowed implicits - so k was reserved for integers. I will be redefining them all as kappa shortly.
I don’t do a lot of math in my type of programs; mostly business logic and reporting. I’ve been looking at Fortran only because of co-arrays which seem to be a cute trick for allocating out threads for multiple tasks.
But one thing I have always discovered is that all languages fail to perform math correctly. Even the much touted math type programs make mistakes. It might take extra steps to limit ALL math in your program to a=bcd but multiple times it’s been the only thing that gets the right answers.
It’s the only way to ensure that the math gets done when the compiler hates you. I know that some people write code that works with equations that rap three lines of code, but it never works past 4 items for me. Write the full equation in comments and the break each piece into 1=234.
