experchange > fortran

robin.vowels (03-04-19, 11:11 PM)
On Monday, March 4, 2019 at 10:33:02 PM UTC+11, john.coll...@gmail.com wrote:
> In the first case I would expect left to right evaluation,


The standard requires it; Dick confirmed it a few posts back.

> so:
> X*Y * I/J
> would evaluate as
> (X*Y*I) / J
> and the division would be real.


That's what Dick said.

[..]
> x,y: 0.1 0.2
> x*y*i/j: 0.400E-01
> i/j*x*y: 0.400E-01


Not a good choice ov values for I [4] and j [2]
if you are trying to illustrate something, because I/J gives
an exact result of 2.

> i/(j*x*y): 0.100E+03
> i,j: 2 4
> x,y: 0.1 0.2
> x*y*i/j: 0.100E-01
> i/j*x*y: 0.000E+00
> i/(j*x*y): 0.250E+02
> i/j*x*y is evaluated as (i/j)*x*y, apparently with left to right ordering of operations at the same level of precedance.


As per the standard.
robin.vowels (03-04-19, 11:16 PM)
On Tuesday, March 5, 2019 at 5:41:45 AM UTC+11, john.coll...@gmail.com wrote:
> No problem - that's what I thought you meant.
> The issue here is that the results of the expressions
> r1 = x*y*i/j
> r2 = i/j*x*y
> or any multiplication - division chain where there is at least one real number and there are integers both sides of any division, are undefined by the standard.


Would you care to show us an example?
Dick Hendrickson (03-05-19, 04:22 AM)
On 3/4/19 12:41 PM, john.collins.simcon wrote:
> No problem - that's what I thought you meant.
> The issue here is that the results of the expressions
> r1 = x*y*i/j
> r2 = i/j*x*y
> or any multiplication - division chain where there is at least one real number and there are integers both sides of any division, are undefined by the standard.


NO, NO, NO, 1000 NOs. ;) The rules in the start if chapter 7 clearly
(IMO) describe expression interpretation as a left-to-right process and
subsequent text describes type conversions.

The section 7 syntax rule interpretation of X*Y*I/J is
( (X*Y) * real(I)) / real(J)
the processor is allowed to evaluate this in any order (presumably using
an efficient one). An optimizing processor might choose to do

T1 = real(I)/real(J)
T2 = X*Y
result = T2*T1

but it can't compute I/J using integer divides (unless it has a
whiz-bang value propagator algorithm for I and J and knows for sure that
I/J and real(I)/real(J) are equal).

As a practical matter, integer division is a special case and the
processor can't do integer division without taking a peek to the left to
see what the context before "I" is.

X*I/J and I/J*X are fundamentally different expressions; X*Y/Z and Y/Z*X
are merely mathematically different (in the low order bits, usually).

A compiler may always choose to evaluate (i/j) first, causing an integer
division, or not to evaluate (i/j) avoiding an integer division.

Again NO! For r1 it cannot use integer division, for r2, it must use
integer division.

> I will try the test program on other compilers, but even if a compiler behaves as though it uses left-to-right ordering on a simple test program it doesn't mean that it will always do so. I will also write a test for this situation and test our code library to see how often it occurs (But this will take some time).


My newsreader (or my knowledge :( ) doesn't allow me to easily go
upthread; but my recollection is that you used something like 4/2 for
I/J. Try 4/6 or 9/4 where there is a significant different difference
between real and integer divides. I think all the compilers will get it
right.

This stuff isn't easy.

Dick Hendrickson
Arjen Markus (03-05-19, 07:39 AM)
On Thursday, February 28, 2019 at 8:37:40 AM UTC+1, ga...@u.washington.edu wrote:

> My favorites are to calculate e or pi. (Multiplied by the
> appropriate power of 10 to make an integer.)
> For both, you only need to divide a multiple precision value
> by a single precision integer, which is much easier than multiple
> precision divisor.


Computing e was simple, as it has a nice, easy-to-remember fast-converging series. My program to do so is:

use libtommath

integer :: decimals, len_approx
type(mp_int) :: googol, numerator, denominator, e_approx
character(len=:), allocatable :: approx_string

decimals = 100
googol = 10
googol = googol ** decimals

numerator = 1 ! 0!
denominator = 1
i = 0

do while ( denominator < googol )
i = i + 1
denominator = i * denominator
numerator = i * numerator + 1
enddo

e_approx = (googol * numerator) / denominator
approx_string = to_string(e_approx)
len_approx = len(approx_string)

write(*,*) 'e = ', approx_string(1:len_approx-decimals), '.', approx_string(len_approx-decimals+1:)

And the answer is:

e = 2.718281828459045235360287471352662497757247093699 95957496696762772407663035354759457138217852516642 74

Now I have to find a suitable series for pi :).

Regards,

Arjen
gah4 (03-05-19, 09:58 AM)
On Monday, March 4, 2019 at 9:39:13 PM UTC-8, Arjen Markus wrote:

(snip, I wrote)

> > My favorites are to calculate e or pi. (Multiplied by the
> > appropriate power of 10 to make an integer.)


[..]
> approx_string = to_string(e_approx)
> len_approx = len(approx_string)
> write(*,*) 'e = ', approx_string(1:len_approx-decimals), '.', approx_string(len_approx-decimals+1:)


Looks more complicated than it should be.

Well, if you want to be able to test mp_int divisor, then
it is probably a good way.

But you should only need to divide mp_int by a regular int.

I did it a few years ago on an IBM 360/20 in decimal.
(The machine has 16 bit binary integers, but can do 15 digit BCD arithmetic.
Only add and subtract in binary, but add, subtract, multiply, and
divide for BCD.) And only need two mp_int values.

For pi, it is slightly more complicated, as you need add and subtract,
along with divide by single integer, and you need three mp_int values.
gah4 (03-05-19, 11:08 AM)
Here is 400 digits of e, without using any mp package.

integer x(101), y(101), base
base=10000
x=0
y=0
x(1)=1000
y(1)=1000
do i=1,210
! divide x by i
m=0
do j=1,ubound(x,1)
k=(x(j)+m)/i
m=(x(j)+m-k*i)*base
x(j)=k
enddo
! add x to y
m=0
do j=ubound(x,1),1,-1
y(j)=y(j)+x(j)+m
m=y(j)/base
y(j)=y(j)-m*base
enddo
enddo
write(*,'(9999I4.4)') y(1:ubound(y,1)-1)
end
Arjen Markus (03-05-19, 12:45 PM)
On Tuesday, March 5, 2019 at 10:08:58 AM UTC+1, ga...@u.washington.edu wrote:
[..]
> enddo
> write(*,'(9999I4.4)') y(1:ubound(y,1)-1)
> end


Intriguing :)

Regards,

Arjen
john.collins.simcon (03-05-19, 01:17 PM)
On Friday, 15 February 2019 00:30:00 UTC, spectrum wrote:
[..]
john.collins.simcon (03-05-19, 01:19 PM)
In the case of
r = x*y*i/j

Steve Kargl replied:
"That is one possible choice. A conforming processor can
evaluate the original expression as if it were written as
(x*y)*(i/j) or x*(y*i)/j or x * i * (y / j) or ..."

I am happy with Dick Hendrickson's reading of the standard, that there cannot be an integer divide in this case.

The preprocessor we wrote to convert all expressions to triples for a fine-grain parallel system would always convert

r = x*y*i/j

to

rtemp_1 = x*y
rtemp_2 = rtemp_1*i
r = rtemp_2/j

where the rtemp_* variables are real.

This has converted quote a lot of code and we have never seen a difference (except in precision) in the behaviour of the converted code.

Similarly, in
r = i/j*x*y

I assume that the divide must always be an integer divide?

Best wishes,

John
robin.vowels (03-05-19, 02:08 PM)
On Tuesday, March 5, 2019 at 10:19:21 PM UTC+11, john.coll...@gmail.com wrote:
[..]
> Similarly, in
> r = i/j*x*y
> I assume that the divide must always be an integer divide?


The standard is quite clear and unambiguous.
Integer division is required, as the interpretation for
equal priority operators is L to R.
robin.vowels (03-05-19, 02:16 PM)
On Tuesday, March 5, 2019 at 4:39:13 PM UTC+11, Arjen Markus wrote:
[..]
> And the answer is:
> e = 2.718281828459045235360287471352662497757247093699 95957496696762772407663035354759457138217852516642 74
> Now I have to find a suitable series for pi :).


There is an algorithm for computing PI to an unlimited
number of digits, without using multiple integer arithmetic.
See Rosetta code.

There's a similar one for e, also without using MIA.
spectrum (03-05-19, 04:16 PM)
On Monday, February 18, 2019 at 12:31:28 AM UTC+9, FortranFan wrote:
> On Sunday, February 17, 2019 at 12:33:42 AM UTC-5, spectrum wrote:
> @spectrum,
> Do you inform your colleague, who recently started using Fortran, to a) work only with DEFINED kinds of floating-point types and b) to ensure the literal constants of floating-point types are coded consistently with such kinds?
> Meaning, given the "fixed numerical value of the Avogadro constant" per IUPAC, about the only approach worth considering in scientific (e.g., chemistry) code written in Fortran is now:
> integer, parameter :: XX = selected_real_kind( .. )
> .
> real(kind=XX), parameter :: NA = 6.02214076E23_xx
> .
> That is, aren't all the other variations you are bringing up in this thread needlessly confusing, if not entirely irrelevant, for a "newcomer"?


Hi @FortranFan,
I'm sorry for a late reply, and no worry about 10^23 (for the newcomer). I'm not recommending its use but opposite (i.e., recommend the use of floating-point exponential
form). Also, my second code snippet above is just for my interest (to see
the result of gfortran-8), so no worry also here.

Because some languages (like Python3) give a correct result even with
10^23 etc, I guess it might cause some misunderstanding for newcomers
about other languages (like Fortran or C++ with machine integer arithmetic).

# Indeed, I had come across this page some time ago

and the OP mentioned something about Fortran, Python, Matlab, so
I worried that the same misunderstanding might occur for new people...
(FYI, I had read some more replies in the above page, but the story seems to have
deviated much to different topics, and only a few initial replies seem to be related here.)
spectrum (03-05-19, 10:59 PM)
# Just to be clear, my sentence above ("...but the story seems to have
deviated much...") refers to the linked page, not this CLF thread.
(The linked page seems to be talking about rational numbers or
and some speed things later.)

--------
Going back to Fortran,

On Saturday, March 2, 2019 at 2:45:00 AM UTC+9, steve kargl wrote:
[..]
> 2 | x = 6.02 * 10**23
> | 1
> Error: Result of exponentiation at (1) exceeds range of INTEGER(4)


For comparison, I've just tried the following code (essentially the same as above),

program main
implicit none
integer i
real x
i = 10**23
x = 6.02 * 10**23
print *, i, x
end

gfortran-8.2 (with no option) on my Mac gives

-159383552 6.01999981E+23

and "gfortran-8 -Wall" gives

x = 6.02 * 10**23
1
Warning: Change of value in conversion from 'INTEGER(4)' to 'REAL(4)' at (1)
[-Wconversion]

Here, although my current understanding (from the above thread) is that
the use of 10**23 itself is not standard-conforming (correct?) and so the compiler can do
anything, I feel that the response by "gfcx" may be more helpful for users who
are potentially not aware of the limitation of default integer
(because of the error stop and a more explanatory message).
I expected that -Wall may give such a warning, but it seems not... (maybe need to use
other warning options?)
spectrum (03-06-19, 01:16 AM)
Just for comparison, Go language seems to be using yet another approach fornumerical
constants. According to the following pages,





both integer and floating-point literals are evaluated *exactly* (using
an arbitrary precision library) and then assigned to variables (of finite precision)
afterward. Overflow occurs only when actual assignment to a variable is performed.
So, a code like this

const myhuge = 1.0e10000
var x = 10.0
fmt.Println( x + myhuge / 1.0e9999 )

gives 20 (which can be tested with the following link (in the TryItOnline site)).



# A related Q/A here:


But I'm not sure whether this approach is good or not because variables areafter all
evaluated with finite precision. Python3 seems to be somewhere in between (or outside??)
by using arbitrary precision for integers (but not for floating points)... So,
gah4 (03-06-19, 02:12 AM)
On Tuesday, March 5, 2019 at 4:16:53 AM UTC-8, robin....@gmail.com wrote:

(snip)

> There is an algorithm for computing PI to an unlimited
> number of digits, without using multiple integer arithmetic.
> See Rosetta code.


For those not following such things, there is an algorithm for
computing individual hexadecimal digits of pi, without computing
all the digits in between. (None for decimal, as far as I know.)

This is used to help validate the results of binary calculations.

Then you only need to validate the binary to decimal conversion.

Similar Threads