IBM Mainframe Forum

by **enrico-sorichetti** » Sat Nov 24, 2018 12:09 am

I got essentially the same result

I did not know that a new comparison operator had been introduced

a quick and dirty test using C ( float,double,long double)
the snippet

 ****** ***************************** Top of Data ******************************
 000001 #include <math.h>                                                       
 000002 #include <stdio.h>                                                      
 000003 #include <stdlib.h>                                                     
 000004 int main()                                                              
 000005 {                                                                       
 000006     float  y  = 57.2;                                                   
 000007     float  z  = 2.23;                                                   
 000008     float  r  ;                                                         
 000009                                                                         
 000010     double yd = 57.2;                                                   
 000011     double zd = 2.23;                                                   
 000012     double rd ;                                                         
 000013                                                                         
 000014     long double yl = 57.2;                                              
 000015     long double zl = 2.23;                                              
 000016     long double rl ;                                                    
 000017                                                                         
 000018     printf("Single     y = %.4f\n", y);                                 
 000019     printf("Single     z = %.4f\n", z);                                 
 000020     printf("Single     ^ = %.4f\n", powf( y,  z) );                     
 000021                                                                         
 000022     printf("Double     y = %.4lf\n", yd);                               
 000023     printf("Double     z = %.4lf\n", zd);                               
 000024     printf("Double     ^ = %.4lf\n", pow(yd, zd) );                     
 000025                                                                         
 000026     printf("Long Doubl y = %.4Lf\n", yl);                               
 000027     printf("Long Doubl z = %.4Lf\n", zl);                               
 000028     printf("Long Doubl ^ = %.4Lf\n", powl(yl, zl) );                    
 000029                                                                         
 000030     return (0);                                                         
 000031 }                                                                       
 ****** **************************** Bottom of Data ****************************

the results on zOS

Select all

Single     y = 57.2000     
Single     z = 2.2300      
Single     ^ = 8298.3594   
Double     y = 57.2000     
Double     z = 2.2300      
Double     ^ = 8298.3767   
Long Doubl y = 57.2000     
Long Doubl z = 2.2300      
Long Doubl ^ = 8298.3767  

the results on apple

Select all

Single     y = 57.2000
Single     z = 2.2300
Single     ^ = 8298.3779
Double     y = 57.2000
Double     z = 2.2300
Double     ^ = 8298.3767
Long Doubl y = 57.2000
Long Doubl z = 2.2300
Long Doubl ^ = 8298.3767
 

while for double and long the results match, quite a bit of difference for a plain float

by **steve-myers** » Sat Nov 24, 2018 7:20 am

One doesn't expect accurate results in single precision hexadecimal floating point, though what we got here was disappointing. Downright embarrassing seems to be a better choice of words.

I thought about trying to disassemble the library module /360 Fortran used for Y^Z. After preliminary research I determined the library module calls another library module - ALOG10 I think it is - and gave up that idea. I had no idea how deep that call chain would go.

X * (1 + Y) ^ Z becomes -

Select all

          LE    0,=E'1.'
         AE    0,Y
         STE   0,ARG1
         CALL  Y^Z,(ARG1,Y)
         ME    0,X
         STE   0,RESULT
         ...
ARG1     DS    E
X        DS    E
RESULT   DS    E

The Y^Z library module returns its result in floating point register 0. This is pretty typical of Fortran library functions of that era.

As you can see, excerpt for the Y^Z issue, this was pretty simple. I'm pretty certain the Fortran compiler would have generated something like this had I told it to compile X * ( 1 + Y) ** Z.

IBM Mainframe Forum

Calculate the power of a floating point number

Re: Calculate the power of a floating point number

Re: Calculate the power of a floating point number