Problem with AMOD intrinsic function using SGI f77 compilers
Author 
Message 
Liam Gumle #1 / 25

Problem with AMOD intrinsic function using SGI f77 compilers
Quote:
>I am having a problem with the AMOD intrinsic function using f77 >compilers on IRIX 5.3 or IRIX 6.0.1. >If I run this on HPUX, I get: > amod(5.0,0.2) = .0 >If I run this on an IBM RS/6000 running xlf, I get: > amod(5.0,0.2) = 0.1999999285 >Now, when I run this on SGI, I get: > amod(5.0,0.2) = 7.4505806E08
The solution is to code your algorithm in a portable manner. As you show, MOD gives different answers on different platforms for your application. It would be better to recode the algorithm so that it gives consistent answers on all platforms. For example: implicit none double precision a,b,rmd a=5.0d0 b=0.2d0 write(*,*) rmd(a,b),mod(a,b) end double precision function rmd(a,b) implicit none double precision a,b rmd=b*(a/bdble(int(a/b))) end rmd gives a consistent result of 0.0d0 on sgi, hp, and ibm, whereas mod gives a different result on each. Try not to fix the problem on one platform, but rather make your code robust enough to give reliable answers on any platform. Cheers, Liam.

Sat, 07 Feb 1998 03:00:00 GMT 


Prof. Loren Meissne #2 / 25

Problem with AMOD intrinsic function using SGI f77 compilers
Quote:
> >I am having a problem with the AMOD intrinsic function using f77 > >compilers on IRIX 5.3 or IRIX 6.0.1. > rmd gives a consistent result of 0.0d0 on sgi, hp, and ibm, whereas mod gives > a different result on each. Try not to fix the problem on one platform, but > rather make your code robust enough to give reliable answers on any platform.
I think that going to double precision merely postpones the problem, it doesn't solve it. If Gumley's DP version happens to work for the example in question, it will probably fail for some other values. Unless you can recast your example in terms of integers, you will always have to worry about inexact results from real calculations. PS. Since Gumley is apparently writing in F90 [note "implicit none"] he could get higher precision from the generic INTRINSIC Mod by setting the "kind" type parameter of the arguments. If you don't want to use the intrinsic, you could make the DP Mod function an internal procedure. Loren Meissner

Sun, 08 Feb 1998 03:00:00 GMT 


Prof. Loren Meissne #3 / 25

Problem with AMOD intrinsic function using SGI f77 compilers
Quote:
> It is easy enough for me to correct this problem, since I only have to > check if the result is negative and add 0.2 if it is the case, but why > could not amod do this by itself ? Am I too picky ?
I guess I vote for "you are too picky". The MOD function on reals is naturally errorprone. Library function writers SHOULD try to get it as close as possible whenever they can, but users have to expect some error in almost any calculation on real numbers. MOD on integers is a different story. The intrinsic function should get the right answer in all cases for integers (except MOD(I, 0) which is undefined, and perhaps some extremely large or small values of I) Loren Meissner

Sun, 08 Feb 1998 03:00:00 GMT 


Michel OLAGN #4 / 25

Problem with AMOD intrinsic function using SGI f77 compilers
Quote: >I am having a problem with the AMOD intrinsic function using f77 >compilers on IRIX 5.3 or IRIX 6.0.1. Here is the code I am using to >show the problem: > print *, 'amod(5.0,0.2) = ', amod(5.0,0.2) > end >If I run this on HPUX, I get: > amod(5.0,0.2) = .0 >which is correct since this is the analytical answer. >If I run this on an IBM RS/6000 running xlf, I get: > amod(5.0,0.2) = 0.1999999285 >which is also correct if one considers that the binary representation of >5.0 and 0.2 is truncated and that if 5.0 is truncated to 4.999999... one >could get an answer like the above. >Now, when I run this on SGI, I get: > amod(5.0,0.2) = 7.4505806E08 >which I think is wrong, since the answer should be between 0.0 >(inclusive) and 0.2 (exclusive). >It is easy enough for me to correct this problem, since I only have to >check if the result is negative and add 0.2 if it is the case, but why >could not amod do this by itself ? Am I too picky ?
I guess that you are. But you are not the only one, and I think that this problem is one of the many reasons why MODULO was introduced in F90, with the constraint that the result should stand between 0.0 (inclusive) and 0.2 (exclusive). MOD(5.0, 0.2) is just a shortcut for 5.0  INT(5.0/0.2)*0.2, with no constraint on the result. Michel 
 IFREMER: Institut Francais de Recherches pour l'Exploitation de la Mer

Sun, 08 Feb 1998 03:00:00 GMT 


Michel OLAGN #5 / 25

Problem with AMOD intrinsic function using SGI f77 compilers
Quote:
>>I am having a problem with the AMOD intrinsic function using f77 >>compilers on IRIX 5.3 or IRIX 6.0.1. >>If I run this on HPUX, I get: >> amod(5.0,0.2) = .0 >>If I run this on an IBM RS/6000 running xlf, I get: >> amod(5.0,0.2) = 0.1999999285 >>Now, when I run this on SGI, I get: >> amod(5.0,0.2) = 7.4505806E08 >The solution is to code your algorithm in a portable manner. As you show, >MOD gives different answers on different platforms for your application. >It would be better to recode the algorithm so that it gives consistent >answers on all platforms. For example: > implicit none > double precision a,b,rmd > a=5.0d0 > b=0.2d0 > write(*,*) rmd(a,b),mod(a,b) > end > double precision function rmd(a,b) > implicit none > double precision a,b > rmd=b*(a/bdble(int(a/b))) > end >rmd gives a consistent result of 0.0d0 on sgi, hp, and ibm, whereas mod gives >a different result on each. Try not to fix the problem on one platform, but >rather make your code robust enough to give reliable answers on any platform.
The fact that the result is consistent on 3 machines for a given pair of numbers is no evidence by itself that it is exact. Did you try it with a=1007999999999.999d0 and b=0.7d0 ? (With my old Sun f77, mod = 0.69908738634622, rmd = 1006496761447.10 with a warning about some arithmetic exceptions: Inexact; Invalid Operand; That is not what I call robust and portable. And I did not try on a Pentium ;)). Michel 
 IFREMER: Institut Francais de Recherches pour l'Exploitation de la Mer

Sun, 08 Feb 1998 03:00:00 GMT 


Liam Gumle #6 / 25

Problem with AMOD intrinsic function using SGI f77 compilers
Quote: >I guess I vote for "you are too picky". The MOD function on reals is >naturally errorprone. Library function writers SHOULD try to get it >as close as possible whenever they can, but users have to expect some >error in almost any calculation on real numbers.
Absolutely. My feeble demonstration with rmd was illdirected. MOD on reals is errorprone. Stick to integers when using MOD where possible.

Sun, 08 Feb 1998 03:00:00 GMT 


Prof. Loren Meissne #7 / 25

Problem with AMOD intrinsic function using SGI f77 compilers
Quote:
> ... MODULO was introduced in F90, > with the constraint that the result should stand between 0.0 (inclusive) > and 0.2 (exclusive).
I think that the only real difference between Mod [in F77 and F90] and Modulo [in F90] is that the latter handles negative numbers differently. As a clue to getting better results from Mod OR Modulo with real arguments, you could look at your numbers and see whether they can be expressed as ratios of integers; then multiply by the common denominator and change to integer type. For example, if you know that X and Y are SUPPOSED TO BE "exact" multiples of 0.1, you can write Mod (NInt (10.0 * X), NInt (10.0 * Y)) / 10.0  here NInt is the F90 rounding intrinsic; for F77 you have to add 0.5 (with the correct sign) and then apply Int. Since the args to Mod have integer type, the generic Mod function will compute with integers; the result after division by 10.0 will be coerced back to real. This is "guaranteed" to work IF your assumption about the denominators is correct. (You could even make a test comparing NInt (10.0 * X) to X within some tolerance.) Note that Mod is still useful. For example, in decomposing an integer you may want to combine integer division with Mod. Since integer division in fortran truncates the same way as Mod, it would be a disaster to replace Mod by Modulo here. [Modulo truncates like Floor, which is a F90 intrinsic; but Fortran does not use Floor automatically when it truncates integer division.]  Loren Meissner

Sun, 08 Feb 1998 03:00:00 GMT 


Joe Landma #8 / 25

Problem with AMOD intrinsic function using SGI f77 compilers
Quote:
>I am having a problem with the AMOD intrinsic function using f77 >compilers on IRIX 5.3 or IRIX 6.0.1. Here is the code I am using to >show the problem: > print *, 'amod(5.0,0.2) = ', amod(5.0,0.2) > end >If I run this on HPUX, I get: > amod(5.0,0.2) = .0 >which is correct since this is the analytical answer. >If I run this on an IBM RS/6000 running xlf, I get: > amod(5.0,0.2) = 0.1999999285 >which is also correct if one considers that the binary representation of >5.0 and 0.2 is truncated and that if 5.0 is truncated to 4.999999... one >could get an answer like the above. >Now, when I run this on SGI, I get: > amod(5.0,0.2) = 7.4505806E08 >which I think is wrong, since the answer should be between 0.0 >(inclusive) and 0.2 (exclusive).
Did I miss something here? THe IBM result differs from 0.2 in the 8th digit after the decimal (where one expects to see single precision type errors) and our result differs from 0.0 (which you accept as correct) in the 8th digit... Actually the 7.45E08 is very close to 0. And amod does the following. MOD(3F) Silicon Graphics MOD(3F) NAME mod, imod, jmod, amod, dmod  FORTRAN remaindering intrinsic functions SYNOPSIS r3 = amod(r1, r2) r3 = mod(r1, r2) DESCRIPTION mod returns the integer remainder of its first argument divided by its second argument. imod returns the integer*2 remainder of its two integer*2 arguments. jmod returns the integer*4 remainder of its two integer*4 arguments. amod and dmod return, respectively, the real and doubleprecision whole number remainder of the integer division of their two arguments. The generic version mod will return the data type of its arguments. The result of these intrinsics is undefined when the value of the second argument is zero. So my question is, if you will not accept our result, why do you accept the IBM result (which suffers from exactly the same problem) ? But I wish to point out something with regard to the IBM. Quote: >If I run this on an IBM RS/6000 running xlf, I get: > amod(5.0,0.2) = 0.1999999285
You should not get this result. You should get a small positive or negative number very close to 0 which is the correct (analytical) answer, which both the SGI and HP machine give you. I suspect that you may need to specify one of the q options to xlf to set the precision appropriately, or check to program for bugs. Quote: >It is easy enough for me to correct this problem, since I only have to >check if the result is negative and add 0.2 if it is the case, but why >could not amod do this by itself ? Am I too picky ?
It is not an SGI problem, though I would suggest checking the xlf compiler flags. If it really bothers you try using the r8 compiler flag on the SGI which will do the following. r8 Use REAL*8 and COMPLEX*16 as the defaults for real and complex variables respectively when they are not explicitly declared with a length. This option is usually needed when porting program from 64bit machines which may result in convergence problem (and very long execution time) if the floating point accuracy is inadequate. Quote:
>professionnel de recherche tel: (514)3695223 fax: (514)3693880 >CERCA (CEntre de Recherche en Calcul Applique) >5160, boul. Decarie, bureau 400, Montreal (Quebec), Canada, H3X 2H9
 Joe Landman Systems Engineer  Silicon Graphics, 24155 Drake Rd, Farmington MI 48335  {}{}{}{} voice: (810) 615 2169  \ O / fax: (810) 478 3181   
+ / \ _ _

Sun, 08 Feb 1998 03:00:00 GMT 


Liam Gumle #9 / 25

Problem with AMOD intrinsic function using SGI f77 compilers
Quote: >The fact that the result is consistent on 3 machines for a given pair of numbers >is no evidence by itself that it is exact. >Did you try it with a=1007999999999.999d0 and b=0.7d0 ? >(With my old Sun f77, mod = 0.69908738634622, rmd = 1006496761447.10 with a warning >about some arithmetic exceptions
Duly noted. Try replacing rmd=b*(a/bdble(int(a/b))) with rmd=b*mod(a/b,1.0d0) and see what you get. Of course division by zero is not checked. Cheers, Liam.

Sun, 08 Feb 1998 03:00:00 GMT 


Ron Sverdlove x25 #10 / 25

Problem with AMOD intrinsic function using SGI f77 compilers
Quote: > ... >denominator and change to integer type. For example, if you know >that X and Y are SUPPOSED TO BE "exact" multiples of 0.1, you can >write Mod (NInt (10.0 * X), NInt (10.0 * Y)) / 10.0  here NInt is >the F90 rounding intrinsic; for F77 you have to add 0.5 (with the >correct sign) and then apply Int. Since the args to Mod have integer > ... > Loren Meissner
The NINT intrinsic function also existed in Fortran 77. It is defined in the same way in both standards: For A > 0, NINT(A) = INT(A+0.5); for A <= 0, NINT(A) = INT(A0.5).  Ronald Sverdlove Computational Science Research
Tel. 6097342517 CN 5300 FAX 6097342662 Princeton, NJ 085435300

Mon, 09 Feb 1998 03:00:00 GMT 


Ron Shepa #11 / 25

Problem with AMOD intrinsic function using SGI f77 compilers
Quote:
>I am having a problem with the AMOD intrinsic function using f77 >compilers on IRIX 5.3 or IRIX 6.0.1. Here is the code I am using to >show the problem: > print *, 'amod(5.0,0.2) = ', amod(5.0,0.2) > end >If I run this on HPUX, I get: > amod(5.0,0.2) = .0 >which is correct since this is the analytical answer.
Just because it is the analytical answer, doesn't mean that it is correct. For P.ne.0, mod(A,P) is defined to be mod(A,P) = A  int(A / P) * P So, for floating point numbers there are two issues: (1) how accurate are A and P (compared to what you think they are), and (2) how accurate are the division, multiplication and substraction to arrive at the final result. The number 5 can be represented exactly in floating point (this isn't required by the standard, but it is true). The number 2/10 cannot be represented exactly, so that is the first part of the issue  the arguments you are using in your mod() are not what you think they are. If "A/P" is equal to or slightly larger than 25, int(a/p) will reduce it to exactly 25. This is then converted to a floating point 25.0 (which is exact because integers are represented in fp exactly) and then multiplied by P. This result can be either slightly smaller or slightly larger than 5.0, so the final result R will be a small floating point number close to zero. I'm looking now on page 223 of the ISO/IEC F90 document, and there are no further restrictions (e.g. nonnegative) that are imposed on the result. Perhaps elsewhere such restrictions are imposed? But it appears that the result in this case will be a small number close to zero of arbitrary sign. But if "A/P" is slightly smaller than 25, int(a/p) will be converted to 24. This will be converted to a fp 24.0, then multiplied by P, to end up being a number close to 4.8. The final subtraction will be a number close to 0.2. Again, there are no further restrictions on the final result, it could be less than 0.2 or greater than 0.2. You would think that the result should be less than "P", but apparently even this is not required for a floating point mod(). The MODULO() function with integer arguments does require the result R to satisfy 0<=R<P for positive P, but I don't see any such restrictions for floating point arguments. (BTW, some arguably lesser languages do not even define the results of some of these operations for integer arguments. :) The lesson, I think, is not to expect more of a floating point MOD() (or MODULO() for that matter) than is guaranteed from its definition. The fortran standards have always stayed clear of specifying floating point accuracy requirements, and this is the heart of the issue. Many programmers avoid using floating point MOD() for this reason. Quote: >If I run this on an IBM RS/6000 running xlf, I get: > amod(5.0,0.2) = 0.1999999285 >which is also correct if one considers that the binary representation of >5.0 and 0.2 is truncated and that if 5.0 is truncated to 4.999999... one >could get an answer like the above.
In practice, the 5.0 is exact and it is the 0.2 that is approximated. However, even if the numbers are represented exactly, e.g. mod(5.0, 0.5), there is still no guarantee of getting the "analytical" result because of the three floating point operations. As far as the language standard is concerned, errors of essentially arbitrary size can be introduced in each of these three operations. Quote: >Now, when I run this on SGI, I get: > amod(5.0,0.2) = 7.4505806E08 >which I think is wrong, since the answer should be between 0.0 >(inclusive) and 0.2 (exclusive).
I'm not so sure. It still looks like it could be right to me. Compute the intermediate results yourself on the SGI and see. Quote: >It is easy enough for me to correct this problem, since I only have to >check if the result is negative and add 0.2 if it is the case, but why >could not amod do this by itself ? Am I too picky ?
You're not too picky, you just want a different function. :) The function you want does not satisfy the above definition. $.02 Ron Shepard

Mon, 09 Feb 1998 03:00:00 GMT 


Robin Vowe #12 / 25

Problem with AMOD intrinsic function using SGI f77 compilers
>I am having a problem with the AMOD intrinsic function using f77 >compilers on IRIX 5.3 or IRIX 6.0.1. >If I run this on HPUX, I get: > amod(5.0,0.2) = .0 >If I run this on an IBM RS/6000 running xlf, I get: > amod(5.0,0.2) = 0.1999999285 >Now, when I run this on SGI, I get: > amod(5.0,0.2) = 7.4505806E08 > It is easy enough for me to correct this problem, since I only have to > check if the result is negative and add 0.2 if it is the case, but why > could not amod do this by itself ? Am I too picky ? 0.2 cannot be represented exactly, thus mod(5.0,0.2) won't yield 0 exactly. You could try multiplying everything by 10 (and rounding up & converting to integer in the process) to ensure that you get 50 and 2, respectively. Then your mod (50, 2) will give you an integer result of 0 .

Tue, 10 Feb 1998 03:00:00 GMT 


Gordon Slishm #13 / 25

Problem with AMOD intrinsic function using SGI f77 compilers
It may be of interest to note that mod(x,y) can be computed exactly whenever x and y are model numbers, y nonzero. To see why, consider two cases: (1) If x < y, then mod(x,y) = x exactly. (2) If x >= y, then depict x and y in fixed point: XXX...XXX ( aligned radix points not shown ) YYY...YYY, and repeatedly add/subtract the y to/from x exactly until the remaining x is less than y. What remains of x has no more bits than y and therefore has an exact model representation (possibly unnormalizable). Gordon 
phone: 88622439, 9149452439

Wed, 11 Feb 1998 03:00:00 GMT 


Prof. Loren Meissne #14 / 25

Problem with AMOD intrinsic function using SGI f77 compilers
Quote:
> It may be of interest to note that mod(x,y) can be computed exactly > whenever x and y are model numbers, y nonzero. .. > (2) If x >= y, then depict x and y in fixed point: > XXX...XXX ( aligned radix points not shown ) > YYY...YYY, > and repeatedly add/subtract the y to/from x exactly until the > remaining x is less than y.
Maybe this is theoretically interesting, but is it practical? Or am I missing something? For example, let X = 2.0 ** 20 Y = 1.0  2.0 ** ( 20) I think that both X and Y are model numbers. [I am using the Fortran "real" model as described on P. 51 of F90 Handbook, with b = 2, p = 24] It takes about 2**20 "exact" subtractions to get a remainder less than Y. How do you "exactly" subtract a very small model number from a very large one, as X  Y in this case? Seems the result of the exact subtraction would not be a model number, as it would require about 40 binary digits in the "mantissa". Loren Meissner

Wed, 11 Feb 1998 03:00:00 GMT 


Gordon Slishm #15 / 25

Problem with AMOD intrinsic function using SGI f77 compilers
: Maybe this is theoretically interesting, but is it practical? Or am : I missing something? For example, let : X = 2.0 ** 20 : Y = 1.0  2.0 ** ( 20) : I think that both X and Y are model numbers. [I am using the Fortran : "real" model as described on P. 51 of F90 Handbook, with b = 2, p = 24] : It takes about 2**20 "exact" subtractions to get a remainder less than Y. : How do you "exactly" subtract a very small model number from a very : large one, as X  Y in this case? Seems the result of the exact : subtraction would not be a model number, as it would require about 40 : binary digits in the "mantissa". : Loren Meissner  In practice, one does not compute mod(x,y) by repeatedly subtracting y. Your apropos example, mod( 2.0**20, 1.02.0**(20) ) equals 0.953674316406250000D06 exactly  with 0.0 roundoff error  as my append predicted since both arguments are model numbers. Regards, Gordon
phone: 88622439, 9149452439

Fri, 13 Feb 1998 03:00:00 GMT 


Page 1 of 2

[ 25 post ] 

Go to page:
[1]
[2] 
