/* REXX - polynom analysis (simplified Newton interpolation):  */

/* p(n) = pq * (q ** n) + pk * (n ** k) + ... p1 * n + p0 for  */
/* k < 10 = DMAX - 2 can be determined by the k-th difference  */
/* d(k,n) = d(k-1,n+1) - d(k-1,n) using d(0,n) = p(n) or by    */
/* the k+1st difference d(k+1,n) if q is neither 0 nor -1.     */

/* If q is 1, then the form c0' = pq + c0 and q' = 0 is found  */
/* in the k-th difference, therefore q = 1 can be ignored.     */

/* The obvious case cq = ck = ... = c0 = 0 is not handled, by  */
/* definition p(n) = 0 has NO degree, i.e. maximal k with non- */
/* zero ck.                                                    */

/* The k+1st difference determines split sequences in the form */
/* n=2i-1: p(n) = o(2i-1) = ok * (i ** k) + ... + o1 * i + o0, */
/* n=2j  : p(n) = e(2j)   = ek * (j ** k) + ... + e1 * j + e0, */
/* for constant d(k+1,n) - d(k+1,n-1) = d(k+1,n) - d(k+1,n+1). */

/* Of course k+3 equations for d(0,1), ..., d(k,1), d(k+1,1),  */
/* and d(k+1,2) cannot determine 2*k+2 variables ok, ek, ...., */
/* o1, e1, o0, e0, and so oj = ej is required for 1 < j <= k.  */

/* The alternating case d(k+1,n) + d(k+1,n+1) = 0 is shown in  */
/* closed form by adding pq * (-1 ** n).  This results in k+2  */
/* pq, pk, ...,, p0 instead of k+3 ok=ek, ..., o2=e2, o1, e1,  */
/* o0, e0 in an equivalent split form.  The k+3rd variable of  */
/* the closed form is obviously q = -1 in (q ** n).            */

   signal on novalue ;  numeric form scientific
   numeric digits 40 ;  numeric fuzz 8
   DMAX = 12         ;  STOP = '(+-*/|&%<=>)!'  /* added ! op. */
   TEST = 28         ;  INIT = 'STOP DMAX E P. D. F.'
   call INIT         ;  arg FORM ;  E = X( 1 )

   select
      when verify( FORM, '(N', 'M' ) > 0 then   /* got formula */
         call DIFF DMAX, WELL( FORM )
      when words( FORM ) > 1 then               /* or sequence */
         call DIFF words( FORM ), SEQU( FORM, 'N' )
      when words( FORM ) = 0 then do forever    /* interactive */
         call charout /**/, 'enter formula (e.g. 3 ** N % 2)'
         say ' or upto' DMAX 'numbers of unknown sequence:'
         pull FORM
         select
            when FORM = '' then exit 0
            when verify( FORM, '(N', 'M' ) > 0 then
               call DIFF DMAX, WELL( FORM )
            when words( FORM ) > 1 then
               call DIFF words( FORM ), SEQU( FORM, 'N' )
            otherwise
               say 'unknown argument:' FORM
         end
      end
      when FORM = '-' then do FORM = -1 to TEST /* test suite: */
         call DEMO FORM ;  say left( '', 79, '-' )
      end
      when FORM = '+' then do                   /* added info  */
         say 'added functions simplifying input of test formulas:'
         say 'x n   = x * (n)                after !0123456789)'
         say 'x (y) = x * (y)                implied product "*"'
         say 'E(x)  = x - O(x)               i.e. select Even n'
         say 'O(x)  = x * (n//2)             i.e. select Odd  n'
         say 'L(x)  = ln(x), ld(x)=L(x)/L(2) logarithm, x > 0'
         say 'X(y)  = exp(y), x = X(L(x))    e**y, e = X(1)'
         say 'G(x,y)= Greatest common divisor of x and y <> 0'
         say 'M(x,y)= (((x//y) + y ) // y)   x Modulo y'
         say 'R(x,y)= x ** (1/y)             Root, y <> 0 <= x'
         say 'D(x)  = x + sum of previous D  inverse Difference'
         say 'F(n-1)= predecessor, F(0)=0    sum: x + F(n-1)'
         say 'x!!, (x)1, L(0), x/0, R(-1,2)  some invalid terms'
      end
      when datatype( FORM, 'W' ) = 0 then do    /* usage info  */
         parse upper source . . FORM
         FORM = substr( FORM, 1+ lastpos( '\', FORM ))
         FORM = substr( FORM, 1+ lastpos( '/', FORM ))
         if pos( '.', FORM ) > 0
            then FORM = left( FORM, lastpos( '.', FORM ) -1 )
         say FORM 'sequence'
         say '     analyzes a sequence of 2 upto' DMAX 'numbers:'
         say '     at least k+3 numbers are required to find'
         say '     a formula of k-th degree for any sequence'
         say FORM 'formula'   ;  DEGR = DMAX || ', where'
         say '     analyzes a formula for N = 1 to' DEGR
         say '     x! = x * (x-1)! and 0! = 1 are supported'
         say FORM
         say '     asks for sequence or formula'
         say FORM '+'
         say '     additional informations'
         say FORM 'k'
         say '     analyzes test sequence k = -1,0,1,2,..'
         say '     use - to run all tests from -1 upto' TEST
      end
      otherwise call DEMO FORM   ;  pull
   end
   exit 0

DEMO: procedure expose (INIT)
   arg FORM ;  TEST = 80 % 8

   if arg() = 1 then select
      when -1 = FORM then do
         say '0 has NO degree by definition'
         FORM = '(0)'               ;  call DEMO 3, FORM
      end
      when 0 <= FORM & FORM < 23 then do
         FORM = FORM - 1   ;  DEGR = FORM // 6 + 1
         TYPE = FORM % 6   ;  FORM = '(n**' || DEGR || ')'
         if DEGR > 2
            then FORM = FORM '+4(n**' || DEGR % 2 + 1 || ')'
         select
            when TYPE = 0 then   /* use polynom of degree 1..6 */
               FORM = FORM '-3'
            when TYPE = 1 then   /* (x ** n) + polynom degrees */
               FORM = DEGR % 2 + 2 || '**n -(' || FORM || '-2)'
            when TYPE = 3 then   /* type 3 not yet implemented */
               FORM = '(n%2)**' || DEGR
            otherwise            /* split odd and even polynom */
               FORM = space( FORM )
               FORM = 'O(' || FORM '+20n +15) + E(' || FORM '-2)'
         end
         if TYPE < 3 | DEGR < 3
            then say 'degree' DEGR 'example (type' TYPE || ')'
            else say 'degree' DEGR '(type 3) not yet recognized'
         call DEMO DEGR + 3 + TYPE, FORM
      end
      when FORM = 23 then do     /* test output format & DMAX: */
         FORM = '(n**9)-8(n**8)'    ;  call DEMO DMAX, FORM
      end
      when FORM = 24 then do     /* test rounding (2**8)*3*3*5 */
         say '1st rounding test (2**8)*45'
         FORM = '16807 * n / 11520' ;  call DEMO 5, FORM
         say left( '', 79, '-' )
         say '2nd rounding test (-1**n)/4'
         FORM = '3 1 1 1 3 5 9 13'
         call DEMO words( FORM ) - 2, SEQU( FORM, 'N+1' )
         say left( '', 79, '-' )
         say '3rd rounding test, degree 6'
         FORM = '2 0 0 0 0 0 0 0 2 12 44 124'
         call DEMO words( FORM ) - 2, SEQU( FORM, 'N+1' )
      end
      when FORM = 25 then do     /* test WELL formed formulas: */
         say 'sigN( N ) test showing limits of extrapolation'
         FORM = 'sign(n)'           ;  call DEMO    3, FORM
         say left( '', 79, '-' )
         say 'N+3 over 3 test'
         FORM = '(n+3)!/(3!n!)'     ;  call DEMO TEST-1, FORM
         say left( '', 79, '-' )
         say 'N+2 over 3 test: bad F(0) caused by 0!'
         FORM = '(n+2)!/(3!(n-1)!)' ;  call DEMO TEST, FORM
      end
      when FORM = 26 then do     /* unrecognized diff. pattern */
         say 'miscellaneous unrecognized patterns'
         FORM = 'O(n) + E(n**2)'    ;  call DEMO TEST, FORM
         say left( '', 79, '-' )
         FORM = '(-1**n)*(n**2)'    ;  call DEMO TEST, FORM
         say left( '', 79, '-' )
         FORM = '(3**n)-(2**n)'     ;  call DEMO TEST, FORM
         say left( '', 79, '-' )
         FORM = '2*N -(N//4 == 2)'  ;  call DEMO TEST, FORM
      end
      when FORM = 27 then do     /* test inverse diff. formula */
         say 'constant inverse diff., ignore undefined F(0)'
         FORM = 'D( 1 )'            ;  call DEMO TEST, FORM
         say left( '', 79, '-' )
         FORM = 'D( 5 )'            ;  call DEMO TEST, FORM
         say left( '', 79, '-' )
         say 'variable inverse diff. not yet recognized'
         FORM = 'D( n )'            ;  call DEMO TEST, FORM
         say left( '', 79, '-' )
         FORM = 'D( n%2 )'          ;  call DEMO TEST, FORM
      end
      when FORM = 28 then do     /* special interest formulas: */
         say 'unrecognized: F(n) = 2**n + 2 * F(n-1) = D( 8*n )'
         FORM = 'D( 8*N )'          ;  call DEMO 8, FORM
         say left( '', 79, '-' )
         FORM = '2**N + 2*F(N-1)'   ;  call DEMO 9, FORM
         say left( '', 79, '-' )
         say 'unrecognized: fibonacci F(n) = F(n-1) + F(n-2)'
         FORM = 'format((( 1 + R( 5 )) / 2 ) ** N / R( 5 ),, 0 )'
         call DEMO TEST, FORM
         say left( '', 79, '-' )
         FORM = 'F(N-1)+F(N-2)+(N<2)'
         call DEMO TEST, FORM
      end
      otherwise say 'unknown argument:' FORM
   end

   if arg() = 1 then return
   arg TEST, FORM ;  TYPE = ''   ;  FORM = WELL( FORM )

   do N = 0 to 1 + TEST
      interpret 'K =' FORM ;  TYPE = TYPE K  ;  F.N = K
   end N
   parse var TYPE K TYPE   ;  N = word( TYPE, 1 + TEST )
   TYPE = subword( TYPE, 1, TEST )  ;  N = QUOT( N )
   K = 'F(' || 1 + TEST || ') =' N || ', F(0) =' QUOT( K )

   if abbrev( translate( FORM ), 'WORD' )
      then say 'test sequence:' TYPE
      else say 'test formula:'  FORM
   pull
   if abbrev( translate( FORM ), 'WORD' )
      then call DIFF TEST, SEQU( TYPE, 'N' )
      else call DIFF TEST, FORM
   say K 'expected'  ;  return

DIFF: procedure expose (INIT)
   arg LAST, FORM ;  LAST = min( LAST, DMAX )
   drop F. D.     ;  T = 7 ;  Z = ''

   do N = 1 to LAST              /* T = 7 digits for N = 1: */
      interpret 'F.N =' FORM        ;  Z = Z ||  F.N
      if 80 < N * 8 then iterate N  ;  A = FIXX( F.N )
      if length( A ) > T - ( A > 0 ) then do
         A = format( A, 2, T - 6,, 0 )
         if length( A ) > T then A = format( A, 2, T - 7,, 0 )
      end                        /* T = 8 digits for N > 1: */
      call charout /**/, right( A, T ) ;  T = 8
   end N

   T = 0 ;  A = 0 ;  D = 0

   if Z <> 0 then do D = 1 to LAST - 1
      K = sign( F.1 ) * sign( F.2 ) ;  D.D = F.1   ;  say
      A = ( K < 0  )             /* alternating: +, -, +, etc. */
      T = ( K <> 0 )             /* geometrical: Q, Q**2, etc. */
                                 /* ignore const 0, 0, 0, etc. */
      do N = 3 to LAST - D + 1 while T | A
         L = N - 1
         T = T & ( F.2 * F.L = F.1 * F.N )
         A = A & ( F.1 + F.2 = ( F.L + F.N ) * ( -1 ** N ))
         A = A & sign( abs( F.N ))
      end N

      T = T & ( 3 < N )          /* & \A: optional split -1**n */
      A = A & ( 3 < N ) & ( F.1 <> F.3 )
      /* F.1 <> F.3: don't match c d c d (found in next diff.) */

      if T | A then leave D      /* difference pattern matched */
      call charout /**/, left( '', 3 * D - 1 )

      do N = 1 to LAST - D
         K = N + 1   ;  F.N = F.K - F.N   ;  K = FIXX( F.N )
         if 80 <= N * 8 + 3 * D - 1 then iterate N /* too long */
         if length( K ) > 8 - ( K > 0 ) then do    /* 7 digits */
            K = format( K, 2, 2,, 0 )
            if length( K ) > 8 then K = format( K, 2, 1,, 0 )
         end
         call charout /**/, right( K, 8 )
      end N
   end D
   else say

   LAST = LAST + 1   ;  say

   if T | A then do N = D + 1 to DMAX
      D.N = 0
   end N
   else do
      /* NOT YET IMPLEMENTED match known D.N patterns */
      /* D.N = K**N => F.N = (K+1)**N already handled */
      /* D.N = C => F.N = 2**(N-1)*C  already handled */
   end

   select
      when T & F.1 = F.2 then          /* const diff. */
         if D < DMAX then do
            FORM = POLY( D, 'n', 0, 0 )
            say 'F(n) =' FORM
            N = 0    ;  interpret 'K =' FORM ;  K = QUOT( K )
            N = LAST ;  interpret 'N =' FORM ;  N = QUOT( N )
            say 'F(' || LAST || ') =' N || ', F(0) =' K
         end
         else say 'n **' D - 2 'formula (unexpected)'
      when T then do                /* geometrical    */
         F.1 = F.2 / F.1
         F.2 = ( D.D / F.1 ) / (( F.1 - 1 ) ** ( D - 1 ))
         FORM = '(' || QUOT( F.1 ) || '**n)'
         if F.2 <> 1 then FORM = space( FACT( F.2 ) || FORM )
         if D <= DMAX then do
            FORM = FORM POLY( D - 1, 'n', F.1 * F.2, F.1 - 1 )
            say 'F(n) =' FORM
            N = 0    ;  interpret 'K =' FORM ;  K = QUOT( K )
            N = LAST ;  interpret 'N =' FORM ;  N = QUOT( N )
            say 'F(' || LAST || ') =' N || ', F(0) =' K
         end
         else say 'F(n) =' FORM '+ formula of' D-2 || 'th degree'
      end
      when A then do                   /* alter diff. */
         do N = D + 1 to 3             /* needs D > 2 */
            F.1 = F.2 - F.1   ;  D.N = F.1
            F.2 = F.3 - F.2   ;  F.3 = F.4 - F.3
         end
         D = max( 3, D )               /* min. degree */
         do N = D + 1 to 2 * D - 2     /* force diff. */
            L = N - 1   ;  D.N = F.2 - D.L
            F.2 = -( D.L + 3 * F.2)    /* alternating */
         end N
         do N = 1 to 2 * D - 2         /* reset input */
            F.N = D( D.N )             /* add D+1 ... */
         end N                         /* ... 2*(D-1) */
         if D <= DMAX then do
            FORM = ODDS( 1, D, 'i') ;  say 'odd  F(2*i-1) =' FORM
            I = LAST % 2 + 1  ;  interpret 'I=' FORM
            FORM = ODDS( 0, D, 'j') ;  say 'even F(2*j  ) =' FORM
            J = 0             ;  interpret 'K=' FORM
            J = LAST % 2      ;  interpret 'J=' FORM
            call charout /**/, 'e.g. F(' right( LAST, 3 ) ') = '
            if LAST // 2
               then say QUOT( I ) || ', F(0) =' QUOT( K )
               else say QUOT( J ) || ', F(0) =' QUOT( K )
         end
         else say 'n **' D - 2 'formula (unexpected)'
      end            /* selection for altern. diff.   */
      when Z = 0 then do               /* const. zero */
         say 'F(n) = 0'
         say 'F(' || LAST || ') = 0, F(0) =' 0
      end
      otherwise say  /* ignore unknown diff. patterns */
   end

   return

ODDS: procedure expose (INIT)
   parse arg O, D, V

   do N = 1 to D - 1
      L = 2 * N - O  ;  D.N = F.L
   end N
   do N = 2 to D - 1
      do K = D - 1 to N by - 1
         L = K - 1   ;  D.K = D.K - D.L
      end K
   end N
   return POLY( D - 1, V, 0, 0 )

POLY: procedure expose (INIT)       /* Newton interpolation */
   parse arg N, V, C, Q             /* x-distance 1 and ... */
   FORM = ''
   if DMAX <= N then return abs()   /* assert known diff.   */

   do D = N to 1 by -1        /* D <= DMAX - 1  */
      do L = D + 1 to N       /* M <= DMAX - 2  */
         M = L - 1   ;  D.D = D.D - P.M.D * D.L
      end L                   /* P.3.4 = 3! = 6 */
      M = D - 1      ;  D.D = ( D.D - C * ( Q ** M )) / P.M.D
   end D                            /* ... Q**N differences */

   do N = N to 3 by -1              /* trunc. fuzz() digits */
      if FIXX( D.N ) = 0 then iterate N
      FORM = FORM FACT( D.N ) || '(' || V || '**' || N - 1 || ')'
   end N                            /* trunc. fuzz() digits */
   if N > 1 & FIXX( D.2 ) <> 0 then FORM = FORM FACT( D.2 ) || V
   if N > 0 & FIXX( D.1 )  > 0 then FORM = FORM '+' || QUOT( D.1 )
   if N > 0 & FIXX( D.1 ) <  0 then FORM = FORM        QUOT( D.1 )
   if FORM > '' | C <> 0
      then return space( FORM )     /* void if Q * C**N + 0 */
      else return 0                 /* (odd/even) polynom 0 */

INIT: procedure expose (INIT)       /* Newton interpolation */
   do L = 0 to DMAX - 2             /* difference constants */
      do D = 1 to L + 1
         F.D = D**L
      end D
      do D = 1 to L + 1       /* D <= DMAX - 1  */
         P.L.D = F.1          /* L <= DMAX - 2  */
         do N = 1 to L + 1 - D
            K = N + 1   ;  F.N = F.K - F.N
         end N
      end D
   end L
   return

FACT: procedure expose (INIT)
   select                           /* add non-zero factor: */
      when arg( 1 ) = +1 then return '+'
      when arg( 1 ) = -1 then return '-'
      when arg( 1 ) <  0 then return  QUOT( arg( 1 )) || '*'
      otherwise         return '+' || QUOT( arg( 1 )) || '*'
   end

QUOT: procedure expose (INIT)
   arg N ;  S = sign( N )  ;  N = abs( N )   ;  R = FIXX( N )
   if R = R % 1 then return S * R            ;  R = 1

   do I = 2 until R > 10 ** fuzz()
      R = R * I            /* 11!=39916800: 8 digits fuzz() */
      if R * N = trunc( FIXX( R * N )) then do
         N = FIXX( R * N ) ;  I = G( R, N )  ;  N = N % I
         return QUOT( S * N ) || '/' || ( R % I )
      end                  /* try to show rational numbers: */
   end I                   /* e.g. 0.0027777777778 => 1/360 */
   return arg( 1 )         /* e.g. 56.888888888889 => 512/9 */

FIXX: procedure            /* round to significant digits   */
   N = digits() - fuzz() - length( trunc( abs( arg( 1 ))))
   N = trunc( abs( arg( 1 )) + 0.5 / (10 ** N), max( 0, N ))
   return sign( arg( 1 )) * N / 1

WELL: procedure expose (INIT)
   arg FORM ;  F.. = ''

   do while FORM > ''                  /* ... N ==> ... (N) */
      TERM = FORM ;  FORM = ''
      if sign( pos( 'N', TERM )) then do
         parse var TERM TERM 'N' FORM
         N = pos( right( strip( TERM ), 1 ), '0123456789)!' )
         if N > 0 then do
            N = left( strip( FORM ), 1 )
            N = pos( N, STOP ) + ( N = '' )
         end
         if N = 0
            then TERM = TERM || 'N'    /* preserving SIGN() */
            else TERM = TERM || '(N)'  /* 5N () ==> 5(N)()  */
      end
      F.. = F.. || TERM
   end
   do while F.. > ''                   /* () () ==> () * () */
      TERM = F..  ;  F.. = ''          /*  5 () ==>  5 * () */
      if sign( pos( '(', TERM )) then do
         parse var TERM TERM '(' F..
         if pos( right( strip( TERM ), 1 ), '0123456789)!' ) > 0
            then TERM = TERM || '*('   /* ...() ==> ...* () */
            else TERM = TERM || '('    /* preserving SIGN() */
      end
      FORM = FORM || TERM
   end
   FORM = reverse( FORM )              /* right to left all */
   do while FORM > ''                  /* x!   ==> @(x)     */
      TERM = FORM ;  FORM = ''         /* (x)! ==> @(x)     */
      if sign( pos( '!', TERM )) then do
         parse var TERM TERM '!' FORM
         if abbrev( strip( FORM ), ')' ) then do
            I = 1 ;  N = pos( ')', FORM )
            do while I > 0          /* decrement ')' count  */
               I = I - 1   ;  J = N + 1
               N = pos( '(', FORM, J )
               K = pos( ')', substr( FORM, J, N - J ))
               do while K > 0
                  I = I + 1   ;  J = J + K
                  K = pos( ')', substr( FORM, J, N - J ))
               end
            end                     /* (...)! ==> @(...)    */
            if N = 0 then N = length( FORM )
            J = substr( FORM, N + 1 )
            FORM = left( FORM, N ) || '@'
         end
         else do                    /* ... x! ==> ... @(x)  */
            N = verify( FORM, STOP, 'M' ) - 1
            if N < 0 then N = length( FORM )
            J = substr( FORM, N + 1 )
            FORM = ')' || left( FORM, N ) || '(@'
         end
         if abbrev( J, '!' )
            then FORM = FORM '*' J  /* x!y! ==> @(x) * @(y) */
            else FORM = FORM J      /* but x!! won't work   */
      end
      F.. = reverse( TERM ) || F..
   end
   F.. = translate( F.., '!', '@' ) /* replace @(x) by !(x) */
   signal on syntax name ZERO       /* syntax test for N=1: */
   N = 1 ;  interpret 'N =' F..  ;  return F..     /* okay  */
ZERO:
   say 'syntax error in' F..        /* else replace dummy 0 */
   say errortext( rc )  ;  pull  ;  F.. = 0  ;  return F..

SEQU: procedure                     /* formula of sequence: */
   return 'word( "' || arg( 1 ) || '",' arg( 2 ) ')'

!: procedure expose (INIT)          /* !(x): WELL formed x! */
   F = 1                            /* assume 0 <= whole x  */
   do N = arg( 1 ) to 1 by -1 ;  F = F * N   ;  end N
   return F

E: procedure expose N               /* E(x) = 0 for  odd N: */
   return arg( 1 ) - O( arg( 1 ))   /* simplify split input */

O: procedure expose N               /* O(x) = 0 for even N: */
   return arg( 1 ) * (N // 2)       /* simplify split input */

M: procedure expose (INIT)          /* remainder <> modulo: */
   arg X, Y ;  return (( X // Y ) + Y ) // Y

D: procedure expose N (INIT)        /* input sequence given */
   if N < 1 then return 0           /* by N-th differences: */
   arg D.N                          /* also used to rebuild */
   do K = N to 2 by -1              /* extended split seq.  */
      L = K - 1   ;  D.L = D.L + D.K
   end K                            /* D(x) and F(x) assume */
   return D.1                       /* N = 1,2,3,... or N=0 */

F: procedure expose N (INIT)        /* primitive recursion: */
   arg K                            /* forced F(N)=0 if N<1 */
   if 0 < K & K < N  then return F.K   ;  else return 0

R: procedure expose (INIT)          /* approximate X**(1/Y) */
   arg X, Y ;  if arg( 2, 'o' ) then Y = 2
   numeric fuzz 0 ;  P = 0 ;  R = 1

   select                           /* force error if Y = 0 */
      when datatype( '1E' || (1/Y), 'n' ) then R = X ** ( 1 / Y )
      when X = 0 then R = 0
      when Y < 0 then R = 1 / R( X, -Y )
      when X < 0 then do            /* X < 0 okay for odd Y */
         do N = 1 to 2 + digits() % 2 while \ datatype( Y, 'w' )
            Y = Y * 10  ;  R = R * 10
         end N /* test Y*(10**N), other primes not reliable */
         if datatype( Y, 'w' ) then do
            N = G( R, Y )  ;  Y = Y / N   ;  R = R / N
         end                        /* R < 10** digits() -2 */
         else Y = R + 2             /* fake big square root */
         if Y // 2                  /* odd Y no square root */
            then if R // 2 then  R = -R( -X ** R, Y )
                           else  R = +R( -X ** R, Y )
            else R = X ** ( R / Y ) /* error if square root */
      end                           /* limit ** 999,999,999 */
      when datatype( '1E' || (Y-1), 'n' ) then do
         Z = Y - 1   ;  Q = 0       /* 500 iterations limit */
         do N = 1 to 500 until O = Q & P = R
            O = P ;  P = Q ;  Q = R /* stop if oscillating  */
            R = ( Z * R + X / ( R ** Z )) / Y
         end N
      end
      otherwise P = -1              /* else: e ** (ln(X)/Y) */
   end
   if P <> R & P <> 0 then R = X( L( X ) / Y )
   return R / 1

X: procedure expose (INIT)          /* e ** X, e = X( 1 )   */
   numeric fuzz 0 ;  arg X ;  S = 1 ;  P = 1
   select                           /* e**(N+x)=e**N * e**x */
      when X < 0 then S = 1 / X(-X) /* e**(-X) = 1 / (e**X) */
      when X > 2 then S = ( E ** ( X % 1 )) * X( X // 1 )
      otherwise do N = 1 until R = S
         R = S ;  P = P * X / N  ;  S = S + P
      end N
   end                              /* X() and L() used for */
   return S / 1                     /* R() root when needed */

L: procedure expose (INIT)          /* ln( X ) log natural. */
   numeric fuzz 0 ;  arg X
   select
      when X < 0 then S = abs()     /* force error by abs() */
      when X < 1 then S = -L( 1 / X )
      when X > E then do            /* ln(x*(e**N))=N+ln(x) */
         N = 0
         do while X > E             /* split x*(e**(2**P)), */
            R = E ;  S = 0          /* 2**P <= 999,999,999: */
            do P = 1 while R < X & datatype( '1E' || S, 'n' )
               R = E ** ( 2 ** P )  ;  parse var R . 'E' S
               if S = '' then S = 0 ;  else S = S * S + 1
            end P                   /* ln(ln(X))/ln(2) <= P */
            P = 2 ** (P-1) ;  S = E ** P
            do until X < S ;  N = N + P   ;  X = X / S   ;  end
         end
         S = N + L( X )             /* 0 <= ln(x) < 1 added */
      end
      otherwise                     /* ln(x) for 1 <= X < e */
         P = (X-1) / (X+1) ;  X = P*P  ;  S = 0
         do N = 1 by 2 until R = S
            R = S ;  S = S + 2 * P / N ;  P = P * X
         end N
   end
   return S / 1

G: procedure expose (INIT)          /* greatest common div. */
   arg X, Y                         /* G() used by root R() */
   do while Y <> 0   ;  Z = Y ;  Y = X // Y  ;  X = Z ;  end
   return abs( X )   /* smallest multiple X * Y % G( X, Y ) */