## Calculate the power of a floating point number

High Level Assembler(HLASM) for MVS & VM & VSE

### Re: Calculate the power of a floating point number

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

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
cheers
enrico
When I tell somebody to RTFM or STFW I usually have the page open in another tab/window of my browser,
so that I am sure that the information requested can be reached with a very small effort
enrico-sorichetti
Global moderator

Posts: 2995
Joined: Fri Apr 18, 2008 11:25 pm
Has thanked: 0 time
Been thanked: 164 times

### Re: Calculate the power of a floating point number

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.
steve-myers
Global moderator

Posts: 2105
Joined: Thu Jun 03, 2010 6:21 pm
Has thanked: 4 times
Been thanked: 243 times

Previous