Problem with AMOD intrinsic function using SGI f77 compilers 
Author Message
 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 HP-UX, 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.4505806E-08

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 re-code 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/b-dble(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  
 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  
 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 error-prone. 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  
 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 HP-UX, 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.4505806E-08

>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  
 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 HP-UX, 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.4505806E-08

>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 re-code 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/b-dble(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  
 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 error-prone. 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 ill-directed. MOD on reals
is error-prone. Stick to integers when using MOD where possible.


Sun, 08 Feb 1998 03:00:00 GMT  
 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  
 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 HP-UX, 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.4505806E-08

>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.45E-08 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
     double-precision 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
          64-bit 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)369-5223  fax: (514)369-3880
>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  
 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/b-dble(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  
 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(A-0.5).

--
Ronald Sverdlove               Computational Science Research

Tel. 609-734-2517              CN 5300
FAX  609-734-2662              Princeton, NJ 08543-5300



Mon, 09 Feb 1998 03:00:00 GMT  
 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 HP-UX, 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.4505806E-08

>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  
 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 HP-UX, 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.4505806E-08

        > 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  
 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 non-zero. 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:  8-862-2439,  914-945-2439



Wed, 11 Feb 1998 03:00:00 GMT  
 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 non-zero.
..
> (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  
 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.0-2.0**(-20) ) equals
0.953674316406250000D-06 exactly -- with 0.0 roundoff error -- as
my append predicted since both arguments are model numbers.
Regards, Gordon


phone:  8-862-2439,  914-945-2439



Fri, 13 Feb 1998 03:00:00 GMT  
 
 [ 25 post ]  Go to page: [1] [2]

 Relevant Pages 

1. Internal compiler error when both using and passing an intrinsic function

2. Intel MMX intrinsic function, compiler problem

3. HELP ON INTRINSIC FUNCTIONS CONVERSION F77-90 ?

4. F77 Intrinsic Functions

5. F77 intrinsic functions

6. MOD intrinsic with the sun f77 compiler

7. Possible bug in SGI Irix 6.x MIPSpro F77 compiler

8. SGI f77 compiler

9. SGI f77 compiler bug?

10. Q: SGI f77 v6.4 compiler switch needed

11. SGI f77 compiler options

12. bug in SGI f77 compiler linker ?

 

 
Powered by phpBB® Forum Software