Double kind coefficients as double kind rationals

Help
1

Hello, Fortran Community! As always, thank you very much for the amazing job you’re doing!

Suppose I have a very long polynomial expression with double kind coefficients (for example, Legendre Polynomials on the unit sphere). I can write the coefficients as something like 21.381638320872348_dk, where dk = kind(1.0d0), or I could write them as 24185702762325.0_dk / 524288.0_dk. Moreover, I can put the common factor multiplying in front as in 1.0_dk / 2.0_dk * (3.0_dk * x + 1.0_dk), or I could write it as 3.0_dk / 2.0_dk * x + 1.0_dk / 2.0_dk, only with much longer numbers. Just for the sake of giving an example, here is a piece of my code:

a(1) = 325.0_dk / 4194304.0_dk * (52003.0_dk * r ** 24 &
       - 16848972.0_dk * r ** 22 * x(3) ** 2 &
       + 895803678.0_dk * r ** 20 * x(3) ** 4 &
       - 18513276012.0_dk * r ** 18 * x(3) ** 6 &
       + 196372963413.0_dk * r ** 16 * x(3) ** 8 &
       - 1221876216792.0_dk * r ** 14 * x(3) ** 10 &
       + 4794938487108.0_dk * r ** 12 * x(3) ** 12 &
       - 12329841823992.0_dk * r ** 10 * x(3) ** 14 &
       + 21063479782653.0_dk * r ** 8 * x(3) ** 16 &
       - 23679206030172.0_dk * r ** 6 * x(3) ** 18 &
       + 16824699021438.0_dk * r ** 4 * x(3) ** 20 &
       - 6846414320412.0_dk * r ** 2 * x(3) ** 22 &
       + 1215486600363.0_dk * x(3) ** 24) * x(1)

I am wondering if there is a “best-practice” way to do this, specially from the numerical point of view. If there is any other advice you could give me on this type of code, it is much appreciated.

Thank you very much in advance!

3

Maybe,
for coefficients (fixed), you could use an array with the attribute parameter
example:
real(dk), parameter :: coeff(*) = [0.1_dk, 0.2_dk, 0.3_dk]
or
real(dk), parameter :: coeff(3) = (/0.1_dk, 0.2_dk, 0.3_dk/)

a(1) = coeff(1) + coeff(2) * x + coeff(3) * x**2

you can also consult the post: Fortran Best Practice Minibook - #2 by Beliavsky

4

Yes, definitely use Horner’s rule. I think the only other question is whether to write out the polynomials explicitly or to evaluate them in a loop. This is a question of what the compiler can recognize and optimize. I suspect that writing out small-order polynomials (with Horner’s rule) is more efficient than the loop, but the loop is more efficient for large-order polynomials. The switchover is probably around order 3 or 4, but actual testing with your compiler and with your production optimization options is probably best.