help implementing frexp

[srvaldez]

I want to implement the frexp function, here's what I have

Code: Pascal  [Select][+][-]

1. program test_frexp;
2. uses math;
3.  
4. type dblptr=^double;
5.  
6. var x,y:double;
7.    m, e:int64;
8.    z:dblptr;
9.    n:int32;
10. begin
11.     x:=3.1415926535897932;
12.     writeln('x = ', x);
13.     m:=int64((@x)^);
14.     e:=int64(m shr 52); // shift-out the mantissa, leaving the exponent
15.     e:=e-$3FE; //exponent - offset
16.     writeln('exponent - offset = ',e);
17.     writeln('hex representation of x = ',hexstr(m,16));
18.     frexp(x,y,n);
19.     m:=int64((@y)^);
20.     writeln('frexp(x,y,n)');
21.     writeln('hex representation of y = ',hexstr(m,16));
22.     writeln('decimal representation of y = ',y);
23.     writeln('n = ',n);
24. end.
25.  

I think that m:=((m shl 11) shr 11)+$3FE0000000000000; should give the same value as that of frexp for y in this example, my problem is how to inject that value into y
<note> x is assumed to be positive

[MathMan]

Just to confirm I got it right - you want to extract the binary exponent from 'xd' and then transfer that to 'yd' leaving sign and mantissa of 'yd' unchanged?

[srvaldez]

hello MathMan, moments after posting this I found the answer so I deleted this thread, why it's still here?
what I wanted was to implement the frexp procedure and as I don't do tricky pointer stuff I was stumped for a while
for example, if x=Pi then frexp[(x, y, e) sets y to .78539816339744828 and the e to 2
here's the procedure

Code: Pascal  [Select][+][-]

1. procedure myFrexp(xd:double; var yd:double; var ne:int32);
2.  
3. type dblptr=^double;
4.  
5. var x:double;
6.    m, e,p:int64;
7.    z:dblptr;
8. begin
9.     x:=abs(xd);
10.     m:=int64((@x)^);
11.     e:=int64(m shr 52); // shift-out the mantissa, leaving the exponent
12.     ne:=e-$3FE; //exponent - offset
13.     p:=((m shl 11) shr 11)+$3FE0000000000000;
14.     z:=@p;
15.     yd:=z^;
16.     if xd<0 then yd:=-yd;
17. end;
18.  

what threw me off was the expression @x^ which the compiler didn't accept so by trial and error I wound up with a convoluted expression that worked but all that was needed was (@x)^

[MathMan]

Hello srvaldez,

Your function seems to be handling 0.0, de-normals, NaNs and Infs incorrect. My take (untested!) would be

Code: Pascal  [Select][+][-]

1. procedure MyFrExp(
2.   x: Double;
3.   var y: Double;
4.   var e: NativeInt
5. );
6.  
7. var
8.   xUInt: UInt64;
9.   ExpBin: NativeInt;
10.  
11. begin
12.   xUInt := pUInt64( @x )^;             // transfer Double to UInt64;
13.   ExpBin := ( xUint shr 52 ) and $7ff; // get biased exponent
14.  
15.   if( ( ExpBin<>0 ) and ( ExpBin<>$7ff ) ) then begin // sort out zero, de-normal, NaN and Inf cases
16.     xUInt := ( xUInt and UInt64( $800fffffffffffff ) ) or $3fe0000000000000;
17.     y := pDouble( @xUInt )^;
18.     e := ExpBin - $3fe;
19.   end else begin
20.     if( ( ExpBin=0 ) and ( ( xUInt and $fffffffffffff )=0 ) ) then begin
21.       e := 0;
22.       y := pDouble( @xUInt )^;   
23.     end else begin
24.       // none of the remaining values - de-normal, NaN or Inf - can begin
25.       // 'normalised' in a meaningfull way. Add whatever you like in error
26.       // handling here
27.     end;
28.   end;
29. end;

[srvaldez]

thanks a million MathMan 
that's just what I needed