! test that single lnegth is accurate. Wrong answers arise if
! this routine is worked at double length accuracy
%begin
   ! Modified     10/January/91   19:00                                  SPOWER2

!Code for f_powrr and f_powri

!---History:
    !
    !SPOWER2    in POWRR, allow a negative real to be raised by a
          !     integral valued real (such as 123.0)

    !SPOWER1    includes spheader.inc rather than spheader#.i (the former
          !     additionally defines IMPROUNDS which was otherise declared
          !     in-line within this file)



%CONSTANTINTEGER INTPTSEL= 0  { if 0 use TRUNC, if 1 then ok INTPT}
                              { if 2 use INTPT with correction of neg vals}

%CONSTANTINTEGER IMPROUNDS= 1 {Rounding on conversion from dp to sp = 1 }
%ROUTINE PRHEX(%INTEGER VALUE, PLACES)
%CONSTSTRING(1)%ARRAY HEX(0:15)="0","1","2","3","4",
               "5","6","7","8","9","A","B","C","D","E","F"
%INTEGER I
      %CYCLE I=PLACES<<2-4, -4, 0
         PRINT STRING(HEX(VALUE>>I&15))
      %REPEAT
%END
%routine error(%integer val)
printstring("error")

%end


                !* Procedures Declared *!

!       %realfnspec POWRI %alias "f_powri"(%real arg , %integer narg)
       %realfnspec POWRR (%real argx,    %real argy)


!SINGLE PRECISION ERROR NUMBERS
%OWNINTEGER %C
 POWERNEG    = 23, { NEGATIVE S.P. VALUE RAISED TO A NON-INTEGER POWER}
 POWERZERO   = 22, { S.P. ZERO RAISED TO NON-POSITIVE POWER   }
 POWERBIG    = 24, { S.P. VALUES RAISED TO S.P. POWER TOO BIG}
 POWERSMALL  = 25  { S.P. VALUE TO S.P. POWER UNDERFLOWS. NOT USED}

%CONSTANTINTEGER  EXPEXCESS= 126, MAXEXP= 254, MINEXP= 1

%CONSTANTREAL     GREATEST = 3.4028234664 @ 38,
                  ONE      =  1.0,
                  SIXTEENTH=  0.0625,
                  ZERO     =  0.0


%CONSTANTINTEGER %C
IEXP= X'7F800000',
IMANT= X'007FFFFF',
EXPEXT23= X'3F000000', {126*2^23}
EXPADD1= X'00800000',  {ADD THIS TO MULTIPLY BY 2, SUBTRACT TO DIVIDE}
EXPADD3= X'02000000'   {ADD THIS TO MULTIPLY BY 16, SUB TO MULT BY 0.0625}




!******************************************************************************
%REALFN POWRR %ALIAS "f_powrr"(%REAL ARG1,ARG2)

!  ARRAY OF 2**((1-J)/16), J= 1,2,...,17 ROUNDED TO WORKING PRECISION

%CONSTANTREALARRAY A1(1:17)= %C
R'3F800000',
R'3F75257D',
R'3F6AC0C7',
R'3F60CCDF',
R'3F5744FD',
R'3F4E248C',
R'3F45672A',
R'3F3D08A4',
R'3F3504F3',
R'3F2D583F',
R'3F25FED7',
R'3F1EF532',
R'3F1837F0',
R'3F11C3D3',
R'3F0B95C2',
R'3F05AAC3',
R'3F000000'

!  ARRAY 2**((1-2*J)/16) - A1(2*J), J= 1,2,...,8. I.E. CORRECTION TO A1.

%CONSTANTREALARRAY A2(1:8)= %C
R'31A92430',
R'B19EAB58',
R'31A8FC20',
R'B2C14FE6',
R'B1ADEAF0',
R'32C12342',
R'32E75622',
R'32CF9890'

%CONSTANTREAL K= 0.44269504088896341 { LOG2 - 1 }

%CONSTANTREAL  %C
P1= 0.83357 541 @-1,
Q1= 0.69314 675 @ 0,
Q2= 0.24018 510 @ 0,
Q3= 0.54360 383 @-1


%CONSTANTINTEGER IW1MAX= 4030, IW1MIN= -4159

%INTEGER PP,P,N,M,M16,IW1,IY {the integer part of ARG2}
    %REAL D,Z,R,V,G,X,Y,U1,U2,W,W1,W2
%LONGREAL UU1,UU2,WW

X= ARG1; Y= ARG2
%IF X<=ZERO %START
   %IF X=ZERO %AND Y>ZERO %THENRESULT= ZERO
   %IF INTPTSEL# 0 %THENSTART                   ; !get the
      IY= INTPT(Y)                              ; !    integer
      %IF INTPTSEL=2 %AND IY< 0 %THEN IY= IY + 1; !       part of ARG2
   %FINISHELSE IY= TRUNC(Y)

   %IF Y# FLOAT(IY) %THEN ERROR (POWERNEG)      ; !error if ARG2 not an integal
!   %RESULT= POWRI (X,IY)
%FINISH

M= ((INTEGER(ADDR(X))&IEXP)>>23) - EXPEXCESS
G= X
INTEGER(ADDR(G))= EXPEXT23 ! (IMANT & INTEGER(ADDR(G)))
!prhex(integer(addr(g)),8); newline
P= 1
%IF G<=A1(9) %THEN P= 9
%IF G<=A1(P+4) %THEN P= P + 4
%IF G<=A1(P+2) %THEN P= P + 2
PP= P + 1
Z= (G - A1(PP)); z=z - A2(PP>>1)
!prhex(integer(addr(z)),8); newline
D= (G + A1(PP))
!prhex(integer(addr(d)),8); newline
INTEGER(ADDR(D))= INTEGER(ADDR(D)) - EXPADD1
!prhex(integer(addr(d)),8); newline
Z= Z/D
V= Z*Z
R= P1*V*Z
!printfl(r,8); newline
R= R + K*R
!printfl(v,8); printfl(p1,8); printfl(z,8); newline
!printfl(v,8); printfl(p1,8); printfl(k,8); newline
!prhex(integer(addr(r)),8); space; prhex(integer(addr(Z)),8); newline
U2= (R + Z*K) + Z
%IF M>0 %THEN M16=M<<4 %ELSE M16= -(-M)<<4
U1= FLOAT(M16 - P)
!prhex(integer(addr(u1)),8); space; prhex(integer(addr(u2)),8); newline
%IF M16#P %THEN %C
      INTEGER(ADDR(U1))= INTEGER(ADDR(U1)) - EXPADD3

UU1= U1
UU2= U2
WW= Y*(UU1 + UU2)
W= WW
%IF IMPROUNDS=0 %START
  WW= WW + (WW - W)
  W= WW
%FINISH
!W1= W + TRICK3
!prhex(integer(addr(w)),8);space;  prhex(integer(addr(ww)),8);
!prhex(integer(addr(ww)+4),8);newline
%IF W#ZERO %START
  INTEGER(ADDR(W))= INTEGER(ADDR(W)) + EXPADD3
  W1= FLOAT(INT(W))
  W= W1
  %IF W1#ZERO %THEN INTEGER(ADDR(W1))= INTEGER(ADDR(W1)) - EXPADD3
  W2= WW - W1

! correction of problems with INTPT

  %IF INTPTSEL=0 %THEN IW1= TRUNC(W) %ELSE IW1= INTPT(W)
  %IF INTPTSEL=2 %START
    %IF IW1<0 %THEN IW1= IW1 + 1
  %FINISH
%FINISHELSESTART
  W1= ZERO
  W2= ZERO
  IW1= 0
%FINISH

!ALLOW OVERFLOW TO CAUSE TERMINATION HERE.
!write(iw1,5); write(integer(addr(w)),8); write(integer(addr(w1)),8); 

!%IF IW1>IW1MAX %START
!PRINTSTRING(" IW1 IS >4030.") %AND NEWLINES(2)
!%IF IW1<IW1MIN %START
!PRINTSTRING(" IW1 IS <-4159.") %AND NEWLINES(2)

%IF W2>ZERO %START
  W2= W2 - SIXTEENTH
  IW1= IW1 + 1
%FINISH

!I= 1
!%IF IW1<0 %THEN I= 0
!M= IW1/16 + I
!P= M*16 - IW1 + 1

%IF IW1<0 %START
  P= ((-IW1)>>4)
  M= -P
  P= -(P<<4) - IW1 + 1
%FINISHELSESTART
  M= IW1>>4 + 1
  P= M<<4 - IW1 + 1
%FINISH

Z= (((Q3*W2 + Q2)*W2 + Q1)*W2)*A1(P)+A1(P)
!Z= A1(P) + A1(P)*Z
%IF Z=ZERO %THENSTART
   %RESULT= ZERO
%FINISHELSESTART
   N= (INTEGER(ADDR(Z))&IEXP)>>23 + M
   %IF N>MAXEXP %OR N<MINEXP %START
      %IF N<MINEXP %START
         {PRINTSTRING(" EXPONENTIATION UNDERFLOW."); NEWLINES(2)}
         {ERROR ( POWERSMALL )}
         %RESULT= ZERO
      %FINISH
      ERROR(  POWERBIG )
      %RESULT= GREATEST
   %FINISH
   INTEGER(ADDR(Z))= N<<23 ! (IMANT & INTEGER(ADDR(Z)))
   %RESULT= Z
%FINISH
%END {OF POWRR}

!******************************************************************************
%real x,y,z
x=2; y=4
z=powrr(x,y)
%if mod(z-16)<0.000001 %then printstring("Test ok") %else printfl(16-z,9)
%ENDOFPROGRAM