Page 2 of 2

Re: Calculate the power of a floating point number

PostPosted: Sat Nov 24, 2018 12:09 am
by enrico-sorichetti
I got essentially the same result

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


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

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


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

Re: Calculate the power of a floating point number

PostPosted: Sat Nov 24, 2018 7:20 am
by steve-myers
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 -
         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.