Exponentiation of integers, was: Fwd: Re: G77BUG9 [fnbroude at phys1104a dot weizmann dot ac dot il: g77-0.5.17.tar.gz]

[Toon Moene <toon at moene dot org>]

Perhaps it helps if I (re-)sent the following from the dim ... and 
distant past:

-------- Original Message --------
Subject: Re: G77BUG9  [fnbroude@phys1104a.weizmann.ac.il: g77-0.5.17.tar.gz]
Date: Mon, 22 Jan 96 08:24:09 +0100
From: Toon Moene <toon@moene.indiv.nluug.nl>
To: fnbroude@phys1104a.weizmann.ac.il
CC: g77-bug@gnu.ai.mit.edu

Hi Cyril (and the readers of g77-bug),

I wrote [about the treatment of i**j with j<0]:

>  What's clearly missing in the (exponent .LT. 0) case is
>  the loop.  Furthermore, on most machines it is probably
>  more efficient to first compute the result with -exponent
>  (using multiplies) followed by a final division, than to
>  have a loop with divisions.  Or the other way around -
>  first inverting the base and using -exponent, and then
>  doing the same as for the positive-exponent-case.  Perhaps
>  this was intended here, but didn't work out as expected
>  ...

This all sounds very reasonable, but it is plain wrong.  In Fortran
integers are _signed_ quantities in arithmetic context, so we have
the following 5 x 5 matrix of possibilities for i ** j:

i ** j		j < -1	j = -1	j = 0	j = 1	j > 1

i < -1		0	0	1	i	i**j

i = -1		i**(-j)	-1	1	-1	i**j	

i = 0		Error	Error	1 (!)	0	0

i = 1		1	1	1	1	1

i > 1		0	0	1	i	i**j

based on the truncation rules I think are valid for Fortran.

I hope I didn't make a mistake.  In that case, row 1 and 5 can be
folded, and row 2 can be dumbed down to -1 for j odd and +1 for j
even.  The remaining i**j case can be implemented with the loop as
generated by g77 now (see my previous mail).

Hope this helps,

--
Toon Moene (toon@moene.indiv.nluug.nl)
Saturnushof 14, 3738 XG  Maartensdijk, The Netherlands
Phone: +31 346 214290; Fax: +31 346 214286
URL: http://www.knmi.nl/hirlam