*        Nth root - x**(1/n)       - 29/07/2018
NTHROOT  CSECT
         USING  NTHROOT,R13        base register
         B      72(R15)            skip savearea
         DC     17F'0'             savearea
         SAVE   (14,12)            save previous context
         ST     R13,4(R15)         link backward
         ST     R15,8(R13)         link forward
         LR     R13,R15            set addressability
         BAL    R14,ROOTN          call rootn(x,n)
         LE     F0,XN              xn=rootn(x,n)
         LA     R0,6               decimals=6
         BAL    R14,FORMATF        edit xn
         MVC    PG(13),0(R1)       output xn
         XPRNT  PG,L'PG            print buffer
         L      R13,4(0,R13)       restore previous savearea pointer
         RETURN (14,12),RC=0       restore registers from calling sav
ROOTN    MVC    ZN,=E'0'           zn=0  ----------------------------
         MVC    ZN,N               n
         MVI    ZN,X'46'           zn=unnormalize(n)
         LE     F0,ZN              zn
         AE     F0,=E'0'           normalized
         STE    F0,ZN              zn=normalize(n)
         LE     F6,=E'0'           xm=0
         LE     F0,X               x
         DE     F0,ZN              /zn
         STE    F0,XN              xn=x/zn
WHILEA   LE     F0,XN              xn
         SER    F0,F6              xn-xm
         LPER   F0,F0              abs((xn-xm)
         DE     F0,XN              /xn
         CE     F0,EPSILON         while abs((xn-xm)/xn)>epsilon
         BNH    EWHILEA            ~
         LE     F6,XN                xm=xn
         LE     F0,ZN                zn
         SE     F0,=E'1'             zn-1
         MER    F0,F6                f0=(zn-1)*xm
         L      R2,N                 n
         BCTR   R2,0                 n-1
         LE     F2,=E'1'             xm
POW      MER    F2,F6                *xm
         BCT    R2,POW               f2=xm**(n-1)
         LE     F4,X                 x
         DER    F4,F2                x/xm**(n-1)
         AER    F0,F4                (zn-1)*xm+x/xm**(n-1)
         DE     F0,ZN                /zn
         STE    F0,XN                xn=((zn-1)*xm+x/xm**(n-1))/zn
         B      WHILEA             endwhile
EWHILEA  LE     F0,XN              xn
         BR     R14                return ---------------------------
         COPY   FORMATF            format a float
X        DC     E'2'               x  <== input
N        DC     F'2'               n  <== input
EPSILON  DC     E'1E-6'            imprecision
XN       DS     E                  xn :: output
ZN       DS     E                  zn=float(n)
PG       DC     CL80' '            buffer
         REGEQU
         END    NTHROOT