Design input for model state

@damianmoz : FWIW, Chapel converts exponentiation by small compile-time known integers into multiplication in order to take away some of the performance sting that you may have felt in other languages. As an example, here's the generated code for taking a variable x (which was an integer in my source code) to the 1..12 powers:

  writeln(x);
  writeln(((int64_t)((x * x))));
  writeln(((int64_t)((((int64_t)((x * x))) * x))));
  call_tmp38 = (x * x);
  writeln(((int64_t)((call_tmp38 * call_tmp38))));
  call_tmp39 = (x * x);
  writeln(((int64_t)((((int64_t)((call_tmp39 * call_tmp39))) * x))));
  call_tmp40 = (x * x);
  writeln(((int64_t)((((int64_t)((call_tmp40 * call_tmp40))) * call_tmp40))));
  writeln(pow(x, INT64(7)));
  call_tmp41 = (x * x);
  call_tmp42 = (call_tmp41 * call_tmp41);
  writeln(((int64_t)((call_tmp42 * call_tmp42))));
  writeln(pow(x, INT64(9)));
  writeln(pow(x, INT64(10)));
  writeln(pow(x, INT64(11)));
  writeln(pow(x, INT64(12)));

I'm curious about your comment about accuracy concerns, though.

In practice in my code, I think most of my uses of exponentiation are either 2**n or n**2.

-Brad

Ah, yeah, that gives much more insight. Think the pow function might be a good way to go.

As a programmer, exponentiation is ideal for clearly expressing various engineering relationships such as

stress = stress0 + K * strain ** n; // Holloman's law
// or
stress = K * ( strain + strain0 ) ** n; // Swift-Krupkowsky law 

where all those identifiers are real(w).

But as a numerical analyst, I find that some people get lazy with an exponentiation operator and I tend to discourage its use. Formula like polynomials should be written using Horner's rule or Estrin's scheme. especially if one has access to the fused multiply-add instruction.

Also, something like

t = x ** 2 - y ** 2

is more accurately written as

t = (x - y) * (x + y)

and even

t = x ** 2 + y **  2

will often need to be written more accurately unless you want it to bite you with accuracy loss. And ...

t = x ** 2 + y ** 2 + z ** 2

Well, that is a whole new can of worms.

And I have seen exponentiations with a negative exponent in the numerator which can hide from view some latent inaccuracy. Such problems can often be cured if the exponent were positive and the exponentiation pushed into the denominator.

Chapel solves some of the problems with exponentiation by expanding small exponents into simple multiplications but the manual is unclear of the details so it makes it hard to do an error analysis that is certain to be correct.

Also, some library implementations of pow() that sit behind exponentiation are pretty dodgey so its use should be minimized and made really obvious.

My 2c.

And if you think that my worrying about inaccuracies are overly paranoid, have I got a problem for you!

1 Like

I don't think they're overly paranoid at all, thanks very much for sharing your experience here @damianmoz. I also agree with your point about Chapel needing to document which cases it transforms and how is important for those who want to be aware of the performance or accuracy implications of using **.

Given a file in the --savec directory, what compiler options are needed to compile that code? That should allow me to stick a -S on the command line to see the assembler output nice and easily. Or is it not that simple?

@damianmoz - say you do chpl myprogram.chpl --savec=tmp and are using the C backend. Then you can do make -f tmp/Makefile to compile the generated C code. You could edit the code / Makefile if that helps you.

You can also supply --ccflags -save-temps but this kind of thing tends to put the output all together in one file. With the LLVM backend you can ask for an LLVM IR dump of a specific function, though.