

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.400E01 > i/j*x*y: 0.400E01 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.100E01 > 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. 


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? 


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 lefttoright 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 whizbang 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 lefttoright 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 


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, easytoremember fastconverging 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_approxdecimals), '.', approx_string(len_approxdecimals+1:) And the answer is: e = 2.718281828459045235360287471352662497757247093699 95957496696762772407663035354759457138217852516642 74 Now I have to find a suitable series for pi :). Regards, Arjen 


On Monday, March 4, 2019 at 9:39:13 PM UTC8, 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_approxdecimals), '.', approx_string(len_approxdecimals+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. 


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)+mk*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 


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 


On Friday, 15 February 2019 00:30:00 UTC, spectrum wrote:
[..] 


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 finegrain 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 


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. 


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. 


On Monday, February 18, 2019 at 12:31:28 AM UTC+9, FortranFan wrote:
> On Sunday, February 17, 2019 at 12:33:42 AM UTC5, spectrum wrote: > @spectrum, > Do you inform your colleague, who recently started using Fortran, to a) work only with DEFINED kinds of floatingpoint types and b) to ensure the literal constants of floatingpoint 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 floatingpoint exponential form). Also, my second code snippet above is just for my interest (to see the result of gfortran8), 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.) 


# 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 gfortran8.2 (with no option) on my Mac gives 159383552 6.01999981E+23 and "gfortran8 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 standardconforming (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?) 


Just for comparison, Go language seems to be using yet another approach fornumerical
constants. According to the following pages, both integer and floatingpoint 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, 


On Tuesday, March 5, 2019 at 4:16:53 AM UTC8, 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. 