/*  Copyright:      (c) 1997-2010 James A.D.W. Anderson. All rights reserved.
 >
 >  File:           $poplocal/local/lib/transarith.p        (Implementation.)
 >
 >  Related files:  $poplocal/local/ref/transarith          (Documentation. )
 >
 >  Documentation:  transarith (REF file).
 >                  http://www.bookofparagon.com
 >                      See under the tabs marked "papers" and "books."
 >                      See the companion paper, "Perspex Machine XII:
 >                      Transfloating-point and Transcomplex Arithmetic with
 >                      Applications in Mathematical Physics," by the author
 >                      of this code.
 >                           `
 >  Language:       This source code is in Pop11 which is part of the Poplog
 >                  programming environment. See:
 >
 >                  http://www.cs.bham.ac.uk/research/projects/poplog/
 >                       freepoplog.html
 >
 >  Purpose:        Provides an implementation of transcomplex arithmetic in
 >                  order to demonstrate the theory of this arithmetic.
 >
 >                  Provides transreal and transcomplex numbers, and
 >                  operations on them. Transreal numbers are implemented as
 >                  Pop11's real numbers and the special objects ninf, pinf,
 >                  nuly. Transcomplex numbers are implemented as records of
 >                  a transreal radius, r, transreal cosine, c, and transreal
 >                  sine, s. Various test procedures are included which
 >                  evaluate formulae discussed in the companion paper.
 >
 >                  W A R N I N G.
 >
 >                  This source code follows transcomplex computational paths,
 >                  even where more accurate, real, computational paths exist.
 >                  The code follows simple computational paths, even where
 >                  more sophisticated computational paths, based on numerical
 >                  analysis, exist. These losses of accuracy makes this code
 >                  unsuitable for end-user applications.
 >
 >                  Note that Pop11 does not explicitly support the IEEE
 >                  floating-point standard so the present source code does
 >                  not trap overflow to signed infinities, does not handle
 >                  minus zero, and does not handle NaNs. Nor does the code
 >                  trap overflow of multilength integer and rational numbers
 >                  on a memory fault. Instead, a system supplied exception
 >                  occurs in these events. In other words, this
 >                  implementation is a model of a total arithmetic, but is
 >                  partial because of hardware faults. It would be better
 >                  to provide a total implementation on transreal hardware.
 >
 >                  Y O U  H A V E  B E E N  W A R N E D.
 >
 >                  The transcomplex arithmetic implemented here contains the
 >                  whole of Cartesian-complex arithmetic, polar-complex
 >                  arithmetic in the principal range, and the whole of both
 >                  transreal and real arithmetic. It is isotropic and is,
 >                  therefore, suitable for use in mathematical physics -
 >                  subject to the above warning about end-user applications.
 >
 >                  In future, a system level implementation of
 >                  transarithmetic would be better, but this would have
 >                  massive impact on the Poplog system's source code.
 >
 >  Note:           Pop11 is a dynamically typed language with an open stack.
 >                  Its builtin arithmetical procedures do a great deal of
 >                  type checking so little extra dynamic checking is needed
 >                  here.
 >
 >                  Prefix "tr_" means "trans_real."
 >                  Prefix "tc_" means "trans_complex."
 >                  Prefix "trans_" means polymorphic "trans_real" or
 >                         "transcomplex."
 >                  Prefix "test_" means testing procedure for system
 >                         developers and readers of the companion paper.
*/




;;; Compiler settings for this file.
compile_mode :pop11 +strict;

section $-library  => transarith
    nuly pinf ninf isstrictlytransreal istransreal checkstrictlytransreal
    checktransreal constransreal desttransreal tr_num tr_den
    ##= ##> ##>= ##< ##<= ##<> ##<=> ##/= ##/> ##/>= ##/< ##/<= ##/<> ##/<=>
    tr_abs tr_sgn tr_sqrt tr_cos tr_sin stdrcs transcomplex_key tc_rad tc_cos
    tc_sin istranscomplex desttranscomplex constranscomplex istransnumber
    isstrictlytranscomplex checktranscomplex checkstrictlytranscomplex
    checktransnumber casttransreal casttranscomplex tc_nuly tc_pinf tc_ninf
    tc_zero tc_i tr_exp tr_log trans_exp trans_log #- ##+ ##- ##* ##/ ##**
    tr_signeddistance tr_distance trans_abs trans_sgn trans_sqrt
    test_1 test_2 test_3 test_4 test_5 test_6 test_7 test_8 test_9 test_10
    test_11 test_12 test_13 test_14 test_15 test_16 test_17 test_18 test_19
    test_20 test_21 test_22
;

;;; Declare existance of library to support auto loading.
vars transarith = true;

section transarith =>
    nuly pinf ninf isstrictlytransreal istransreal checkstrictlytransreal
    checktransreal constransreal desttransreal tr_num tr_den
    ##= ##> ##>= ##< ##<= ##<> ##<=> ##/= ##/> ##/>= ##/< ##/<= ##/<> ##/<=>
    tr_abs tr_sgn tr_sqrt tr_cos tr_sin stdrcs transcomplex_key tc_rad tc_cos
    tc_sin istranscomplex desttranscomplex constranscomplex istransnumber
    isstrictlytranscomplex checktranscomplex checkstrictlytranscomplex
    checktransnumber casttransreal casttranscomplex tc_nuly tc_pinf tc_ninf
    tc_zero tc_i tr_exp tr_log trans_exp trans_log #- ##+ ##- ##* ##/ ##**
    tr_signeddistance tr_distance trans_abs trans_sgn trans_sqrt
    test_1 test_2 test_3 test_4 test_5 test_6 test_7 test_8 test_9 test_10
    test_11 test_12 test_13 test_14 test_15 test_16 test_17 test_18 test_19
    test_20 test_21 test_22
;

;;; Set built in, trigonometric functions to operate in radians rather than
;;; degrees.
;;;
;;; Note: this assignment propagates to all sections. A system level
;;; implementation could be more conservative.
true -> popradians;

;;; Set builtin, floating-point arithmetic to operate in double pecision.
;;; This is optional, but allows the sharpest comparisons with ordinary
;;; implementations. Recall, that Pop11 does not comply with the IEEE
;;; floating-point standards so it is difficult to make a valid comparison
;;; with external code.
;;;
;;; Note: this assignment propagates to all sections. A system level
;;; implementation could be more conservative.
true -> popdprecision;




/*****************************************************************************
*                            Transreal Constants                             *
*****************************************************************************/

;;; Transcomplex constants are given after the implementation of the
;;; transcomplex type.

;;; Strictly transreal constants.
;;;
;;; Note: a system level implementation would be neater, but would have a
;;; massive impact on all of the numerical operations supported by Poplog.
;;; In particular, call-out to conventional languages would be problematical
;;; until those languages support transarithmetic.
constant
    nuly   = consword( '0_/0'),
    pinf   = consword( '1_/0'),
    ninf   = consword('-1_/0'),
;




/*****************************************************************************
*            Transreal Recogniser, Constructor, Destructor, Cast.            *
*****************************************************************************/

;;; Return true if num is strictly transreal and false otherwise.
;;; In other words, return true if num is non-finite and transreal, false
;;; otherwise. It might be better to use "nonfinite," rather than "strictly"
;;; in the name of procedures.
define isstrictlytransreal(num);
    lvars num;
    num == nuly or num == pinf or num == ninf
enddefine;

;;; Return true if num is transreal and false otherwise.
define istransreal(num);
    lvars num;
    isreal(num) or isstrictlytransreal(num)
enddefine;

;;; Typecheck a strictlytransreal number.
define checkstrictlytransreal(num);
    lvars num;
    unless isstrictlytransreal(num)
    then   mishap('STRICTLY TRANSREAL NUMBER NEEDED', [^num])
    endunless;
enddefine;

;;; Typecheck a transreal number.
define checktransreal(num);
    lvars num;
    unless istransreal(num)
    then   mishap('TRANSREAL NUMBER NEEDED', [^num])
    endunless;
enddefine;

;;; desttransreal(num) -> (numerator,denominator)
define desttransreal(num);
    lvars num;

    checktransreal(num);

    ;;; Return numerator and denominator.
    if      num == nuly
    then    0,0
    elseif  num == pinf
    then    1,0
    elseif  num == ninf
    then    -1,0
    elseif  isrational(num)
    then    destratio(num)
    else    num,1
    endif;
enddefine;

;;; Return a transreal number given the numerator and denominator.
;;; Irrational numbers conventionally have denominator d = 1.
;;; constransreal(numerator, denominator) -> num;
define constransreal(n, d);
    lvars n, d, nn, dn, nd, dd;

    ;;; Decompose the given numerator, n.
    desttransreal(n) -> (nn, dn);

    ;;; Decompose the given denominator, d.
    desttransreal(d) -> (nd, dd);

    ;;; Carry the sign of d in dd so that the reciprocal is in canonical form.
    if    nd < 0
    then -nd -> nd;
         -dd -> dd;
    endif;

    ;;; Multiply out by the reciprocal.
    nn*dd -> n;
    dn*nd -> d;

    ;;; Return the transreal number with numerator, n, and denominator, d.
    if   d = 0
    then if     n = 0
         then   nuly
         elseif n > 0
         then   pinf
         else   ninf
         endif
    else n/d
    endif;
enddefine;

;;; Return the numerator, n, of a transreal number, num.
define tr_num(num) -> n;
    lvars num, n, d;
    desttransreal(num) -> (n,d);
enddefine;

;;; Return the denominator, d, of a transreal number, num.
define tr_den(num) -> d;
    lvars n, d, num;
    desttransreal(num) -> (n,d);
enddefine;




/*****************************************************************************
*                      Transreal Arithmetical Functions.                     *
*****************************************************************************/

;;; Return num1 plus num2.
define tr_add(num1, num2);
    lvars n1, n2, d1, d2, num1, num2;

    desttransreal(num1) -> (n1,d1);
    desttransreal(num2) -> (n2,d2);

    if   d1 = 0 and d2 = 0  and abs(n1) = 1 and abs(n2) = 1
    then constransreal(sign(n1 + n2), 0)
    else constransreal(n1*d2 + n2*d1, d1*d2)
    endif;
enddefine;

;;; Negate a transreal number.
define tr_negate(num);
    lvars num, n, d;
    desttransreal(num) -> (n,d);
    negate(n) -> n;
    constransreal(n,d);
enddefine;

;;; Return num1 minus num2.
define tr_sub(num1, num2);
    lvars num1, num2;
    tr_add(num1, tr_negate(num2));
enddefine;

;;; Return the product of two transreal numbers.
define tr_mul(num1, num2);
    lvars n1, n2, d1, d2, num1, num2;
    desttransreal(num1) -> (n1,d1);
    desttransreal(num2) -> (n2,d2);
    constransreal(n1*n2, d1*d2);
enddefine;

;;; Return the reciprocal of a transreal number.
define tr_recip(num);
    lvars n, d, num;
    desttransreal(num) -> (n,d);
    constransreal(d, n);
enddefine;

;;; Return num1 divided by num2.
define tr_div(num1, num2);
    lvars num1, num2;
    tr_mul(num1, tr_recip(num2));
enddefine;




/*****************************************************************************
*                       Transreal Relational Operators.                      *
*****************************************************************************/

;;; Real arithmetic obeys the trichotomy axiom: a number can be greater than,
;;; equal to, or less than zero. But transreal arithmetic obeys the
;;; quadrachotomy axiom: a number can be nullity or else it can be greater
;;; than, equal to, or less than zero. That is, there is one extra case to
;;; consider in transreal arithmetic - a number can be nullity.

;;; Transreal arithmetic provides three, mutually exclusive, Boolean,
;;; ordering relations: less than (<), equal (=), and greater than (>). All
;;; combinations of these operations are implemented here, except the empty
;;; operator, with no occurences of <, =, or >. Thus 2^3 - 1 = 8 - 1 = 7
;;; operators are implemented.

;;; Pop11 uses the symbol slash (/) to negate the relational operators. This
;;; is cloned here. Hence, 7 negated operators are implemented.

;;; The relational operators implemented here clone the priority of the
;;; corresponding Pop11 relational operators, where the Pop11 operators exist.

;;; Return true if num1 > num2, otherwise return false.
define greatertransreal(num1, num2);
    lvars num1, num2;
    tr_num(tr_sub(num1, num2)) > 0;
enddefine;

define 7 ##= (num1, num2);
    lvars num1, num2;
    checktransreal(num1);
    checktransreal(num2);
    num1 = num2
enddefine;

define 6 ##> (num1, num2);
    lvars num1, num2;
    greatertransreal(num1, num2);
enddefine;

define 6 ##>= (num1, num2);
    lvars num1, num2;
    num1 ##> num2 or num1 = num2;
enddefine;

define 6 ##< (num1, num2);
    lvars num1, num2;
    num2 ##> num1;
enddefine;

define 6 ##<= (num1, num2);
    lvars num1, num2;
    num2 ##>= num1;
enddefine;

define 6 ##<> (num1, num2);
    lvars num1, num2;
    num1 ##< num2 or num1 ##> num2;
enddefine;

define 6 ##<=> (num1, num2);
    lvars num1, num2;
    num1 ##< num2 or num1 = num2 or num1 ##> num2
enddefine;

define 7 ##/= (num1, num2);
    lvars num1, num2;
    not(num1 ##= num2)
enddefine;

define 6 ##/> (num1, num2);
    lvars num1, num2;
    not(num1 ##> num2);
enddefine;

define 6 ##/>= (num1, num2);
    lvars num1, num2;
    not(num1 ##>= num2);
enddefine;

define 6 ##/< (num1, num2);
    lvars num1, num2;
    not(num1 ##< num2);
enddefine;

define 6 ##/<= (num1, num2);
    lvars num1, num2;
    not(num1 ##<= num2);
enddefine;

define 6 ##/<> (num1, num2);
    lvars num1, num2;
    not(num1 ##<> num2);
enddefine;

define 6 ##/<=> (num1, num2);
    lvars num1, num2;
    not(num1 ##<=> num2);
enddefine;




/*****************************************************************************
*           Transreal Functions Needed For Transcomplex Arithmetic.          *
*****************************************************************************/

;;; Return the transreal square root of x.
define tr_sqrt(x);
    lvars x;

    ;;; Check transreal and non-negative.
    if   x ##< 0
    then mishap('NON-NEGATIVE TRANSREAL NUMBER NEEDED', [^x])
    endif;

    if     x == nuly
    then   nuly
    elseif x == pinf
    then   pinf
    else   sqrt(x)
    endif;
enddefine;

;;; Return the radius, cosine and sine in standard form given arbitrary
;;; transreal radius, r, cosine, c, and sine, s.
;;;
;;; In production code it would be better to normalise when the user
;;; creates a transcomplex number, but not to renormalise during builtin
;;; arithmetic and functions, except where numerical analysis suggests that
;;; renormalisation is advisable.
define stdrcs(r, c, s) -> (r, c, s);
    lvars r, c, s, r_tmp;

    checktransreal(r);
    checktransreal(c);
    checktransreal(s);

    ;;; Normalise the cosine and sine so that, apart from rounding errors,
    ;;; c^2 + s^2 = 1.
    tr_sqrt(tr_add(tr_mul(c,c), tr_mul(s,s))) -> r_tmp;
    tr_div(c, r_tmp) -> c;
    tr_div(s, r_tmp) -> s;

    ;;; Normalise the radius, cosine, and sine so that the radius is
    ;;; non-negative.
    if   r ##< 0
    then tr_negate(r) -> r;
         tr_negate(c) -> c;
         tr_negate(s) -> s;
    endif;

    ;;; Normalise the cosine and sine so that all non-finite angles are
    ;;; mapped onto the angle nullity.
    if   isstrictlytransreal(c) or isstrictlytransreal(s)
    then nuly ->> c -> s
    endif;

    ;;; Normalise the cosine and sine so that nullity radius has nullity angle.
    if   r == nuly
    then nuly ->> c -> s
    endif;

    ;;; Normalise the cosine and sine so that zero radius with finite angle
    ;;; is mapped to zero radius with zero angle.
    if   r = 0 and isreal(c)
    then 1 -> c; 0 -> s
    endif;
enddefine;




/*****************************************************************************
*           Transcomplex Recogniser, Constructor, Destructor, Cast.          *
*****************************************************************************/

;;; Implement polar-transcomplex numbers, along with their constructor,
;;; destructor, recogniser (and key).
;;;
;;; It might be better to separate type checking from from number set
;;; membership. For example, 2 is a member of the real, transreal, complex,
;;; and transcomplex numbers; but it has builtin type integer and
;;; polymorphically real. A purist would keep the builtin real and complex
;;; types distinct, but many languages promote real to complex. Types are a
;;; language issue and need to be resolved for each individual language. The
;;; membership of mathematical sets is already clear. So it might be better to
;;; develop a hierarchy of builtin types and then have separate predicates to
;;; test type and finiteness.

defclass transcomplex {tc_rad, tc_cos, tc_sin};

;;; Modify the destructor so that it can destroy a real number.
constant basetranscomplexdest = class_dest(transcomplex_key);

;;; Given num real or transcomplex, return (r, c, s).
;;; This mimics the behaviour of Pop11's destcomplex.
;;; Question: why does Pop11's class_dest not have an updater?
;;; Is there a better way of modifying the class_dest?
define desttranscomplex(num);
    lvars num;
    if   istransreal(num)
    then stdrcs(num, 1, 0)
    else basetranscomplexdest(num)
    endif;
enddefine;

;;; Modify the constructor so that it puts radius, consine and sine in
;;; standard form.
;;;
;;; Question: why does Pop11's class_cons not have an updater?
;;; Is there a better way of modifying the class_cons?
constant basetranscomplexcons = class_cons(transcomplex_key);

define constranscomplex(r, c, s);
    lvars r, c, s;
    basetranscomplexcons(stdrcs(r,c,s));
enddefine;

;;; Return true if num is transreal or transcomplex, false otherwise.
define istransnumber(num);
    lvars num;
    istransreal(num) or istranscomplex(num)
enddefine;

;;; Return true if num is strictly transcomplex, false otherwise.
;;; In other words, return true if num is transcomplex and non-finite, false
;;; otherwise.
define isstrictlytranscomplex(num);
    lvars num;
    if   istranscomplex(num)
    then isstrictlytransreal(tc_rad(num)) or
         isstrictlytransreal(tc_cos(num))
    else false
    endif;
enddefine;

;;; Typecheck a transcomplex number.
define checktranscomplex(num);
    lvars num;
    unless istranscomplex(num)
    then   mishap('TRANSCOMPLEX NUMBER NEEDED', [^num])
    endunless;
enddefine;

;;; Typecheck a strictly transcomplex number.
define checkstrictlytranscomplex(num);
    lvars num;
    unless isstrictlytranscomplex(num)
    then   mishap('STRICTLY TRANSCOMPLEX NUMBER NEEDED', [^num])
    endunless;
enddefine;

;;; Typecheck a transreal or transcomplex number.
define checktransnumber(num);
    lvars num;
    unless istransnumber(num)
    then   mishap('TRANSREAL OR TRANSCOMPLEX NUMBER NEEDED', [^num])
    endunless;
enddefine;

;;; Convert a number to transreal form unless it is transcomplex and does not
;;; lie on the real axis.
;;;
;;; casttransreal(num) -> num.
;;;
;;; For display to the user it might be better to have a second version of
;;; the cast that rounds to real if the number is nearly on the real axis.
;;;
;;; This procedure is sometimes used to avoid transcomplex contagion, but it
;;; might be better either to accept contagion or else have separate
;;; transreal and transcomplex functions, the transreal ones being totalised
;;; by returning nullity when no real solution is available on the given
;;; computational path.
define casttransreal(num);
    lvars num, r, c, s;

    desttranscomplex(num) -> (r,c,s);

    if     c = 1
    then   r
    elseif c = -1
    then   tr_negate(r)
    elseif r == nuly
    then   nuly
    else   num
    endif;
enddefine;

;;; Convert a number to transcomplex form.
define casttranscomplex(num) -> num;
    lvars num;
    checktransnumber(num);
    if     istransreal(num)
    then   constranscomplex(num, 1, 0) -> num
    endif;
enddefine;




/*****************************************************************************
*                           Transcomplex Constants.                          *
*****************************************************************************/

constant tc_nuly = casttranscomplex(nuly),
         tc_pinf = casttranscomplex(pinf),
         tc_ninf = casttranscomplex(ninf),
         tc_zero = casttranscomplex(0),
         tc_i    = constranscomplex(1,0,1),
;



/*****************************************************************************
*                     Transcomplex Arithmetical Functions                    *
*****************************************************************************/

;;; Negate a transcomplex number.
define trans_negate(num);
    lvars num, r, c, s;

    desttranscomplex(num) -> (r, c, s);

    constranscomplex(r, tr_negate(c), tr_negate(s)) -> num;

    ;;; Cast to transreal if possible.
    casttransreal(num);
enddefine;

;;; Return the product of transreal or transcomplex numbers.
;;; All cases are computed as transcomplex so as to exercise the most
;;; general form of the source code.
define trans_mul(num1, num2);
    lvars num1, num2, r1, c1, s1, r2, c2, s2, r3, c3, s3;

    desttranscomplex(num1) -> (r1,c1,s1);
    desttranscomplex(num2) -> (r2,c2,s2);

    tr_mul(r1,r2) -> r3;

    if   r3 = 0
    then 1 -> c3;
         0 -> s3;
    else tr_sub(tr_mul(c1,c2), tr_mul(s1,s2)) -> c3;
         tr_add(tr_mul(s1,c2), tr_mul(c1,s2)) -> s3;
    endif;

    casttransreal(
        constranscomplex(r3,c3,s3)
    );
enddefine;

;;; Return the quotient of transreal or transcomplex numbers.
;;; All cases are computed as transcomplex so as to exercise the most
;;; general form of the source code.
define trans_div(num1, num2);
    lvars num1, num2, r1, c1, s1, r2, c2, s2, r3, c3, s3;

    desttranscomplex(num1) -> (r1,c1,s1);
    desttranscomplex(num2) -> (r2,c2,s2);

    tr_div(r1,r2) -> r3;

    if   r3 = 0
    then 1 -> c3;
         0 -> s3;
    else tr_add(tr_mul(c1,c2), tr_mul(s1,s2)) -> c3;
         tr_sub(tr_mul(s1,c2), tr_mul(c1,s2)) -> s3;
    endif;

    casttransreal(
        constranscomplex(r3,c3,s3)
    );
enddefine;

;;; Return num1 plus num2.
;;; All cases are computed as transcomplex so as to exercise the most
;;; general form of the source code; but a sum involving a vector at angle
;;; nullity involves only a transreal computation.
define trans_add(num1, num2);
    lvars num1, num2, r1, c1, s1, r2, c2, s2, r3, c3, s3, k1, k2, k3, x, y;

    desttranscomplex(num1) -> (r1,c1,s1);
    desttranscomplex(num2) -> (r2,c2,s2);

    tr_add(r1,r2) -> r3;

    if   c1 == nuly or c2 == nuly
    then
         ;;; Phase collapses to nullity, but magnitudes sum (as above at r3).
         nuly ->> c3 -> s3
    else

         if r1 = r3 then 1 else tr_div(r1,r3) endif -> k1;
         if r2 = r3 then 1 else tr_div(r2,r3) endif -> k2;

         tr_add(tr_mul(k1,c1), tr_mul(k2,c2)) -> x;
         tr_add(tr_mul(k1,s1), tr_mul(k2,s2)) -> y;

         tr_sqrt(tr_add(tr_mul(x,x), tr_mul(y,y)))  -> k3;

         if   k3 = 0
         then 1 -> c3;
              0 -> s3;
         else tr_div(x,k3) -> c3;
              tr_div(y,k3) -> s3;
         endif;

         tr_mul(r3,k3) -> r3;
    endif;

    ;;; Cast the sum to transreal if possible.
    casttransreal(constranscomplex(r3,c3,s3));
enddefine;

;;; Return num1 minus num2.
;;; All cases are computed as transcomplex so as to exercise the most
;;; general form of the source code.
define trans_sub(num1, num2);
    lvars num1, num2;
    trans_add(num1, trans_negate(num2));
enddefine;




/*****************************************************************************
*                 Transcomplex Exponential, Logarithm, Power.                *
*****************************************************************************/

;;; Return the transcomplex, including transreal, exponential.
define trans_exp(x);
    lvars x, y, r, c, s;

    desttranscomplex(x) -> (r, c, s);

    if     r == nuly
    then
           ;;; Exponential is identically nullity.
           x

    elseif r == pinf
    then
           if   c ##< 0 or s = -1
           then
                ;;; Exponential is zero.
                tc_zero
           else
                ;;; Exponential is identically infinite.
                x
           endif;

    elseif c == nuly
    then
           ;;; Finite exponential of finite radius at angle identically
           ;;; nullity.
           constranscomplex(exp(r), c, s)

    else
           ;;; Compute Cartesian-complex exponential.
           exp(r*c +: r*s) -> y;

           ;;; Compute Cartesian-complex components.
           realpart(y) -> x;
           imagpart(y) -> y;

           sqrt(x*x + y*y) -> r;

           x/r -> c;
           y/r -> s;

           ;;; Complex exponential.
           constranscomplex(r,c,s);
    endif -> y;

    ;;; Cast to transreal if possible.
    casttransreal(y);
enddefine;

;;; Return the transcomplex, including transreal, logarithm.
define trans_log(x);
    lvars x, y, r, c, s;

    desttranscomplex(x) -> (r,c,s);

    if     r == nuly
    then
           ;;; Logarithm is identically nullity.
           x

    elseif r == pinf
    then
           if   c ##< 0 or s = -1
           then
                ;;; Logarithm of a negative infinity.
                ;;; This is a totalising step.
                tc_nuly

           else
                ;;; Logarithm is identically infinite.
                x
           endif

    elseif r = 0
    then
           ;;; Logarithm is negative infinity.
           ;;; Note that at angle nullity, log(r) = log(1/r) >= 0.
           constranscomplex(ninf, c, s)

    elseif c == nuly
    then
           ;;; Finite logarithm of finite, non-zero radius at angle nullity.
           ;;; Note that at angle nullity, log(r) = log(1/r) >= 0.
           constranscomplex(log(r), c, s)

    else
           ;;; Compute Cartesian-complex logarithm.
           log(r*c +: r*s) -> y;

           ;;; Compute Cartesian-complex components.
           realpart(y) -> x;
           imagpart(y) -> y;

           sqrt(x*x + y*y) -> r;

           if   r = 0
           then 1 -> c;
                0 -> s;
           else x/r -> c;
                y/r -> s;
           endif;

           ;;; Complex logarithm.
           constranscomplex(r,c,s);
    endif -> y;

    ;;; Cast to transreal if possible.
    casttransreal(y);
enddefine;

;;; Compute num1 raised to the power of num2.
;;; All cases are computed as transcomplex so as to exercise the most
;;; general form of the source code.
define trans_power(num1, num2);
    lvars num1, num2, num3, r1, c1, s1, sgnnum1, absnum1;

    desttranscomplex(num1) -> (r1,c1,s1);

    if   r1 == pinf
    then

         ;;; Radius of sgnnum must be positive, finite, and not unity.
         ;;; It is convenient to set the radius equal to the radix of the
         ;;; floating-point arithmetic. In binary arithmetic we set radius to
         ;;; two, in denary arithmetic, we set radius to ten.
         constranscomplex( 2,c1,s1) -> sgnnum1;
         constranscomplex(r1, 1, 0) -> absnum1;

         trans_mul(
             trans_exp(trans_mul(num2, trans_log(sgnnum1))),
             trans_exp(trans_mul(num2, trans_log(absnum1)))
         )
    else
             trans_exp(trans_mul(num2, trans_log(  num1)))
    endif -> num3;

    ;;; Cast the result to transreal if possible.
    casttransreal(num3);
enddefine;




/*****************************************************************************
*                         Transrithmetical Operators.                        *
*****************************************************************************/

;;; The operations defined in this section operate on both transreal and
;;; transcomplex numbers.
;;;
;;; The symbol # denotes a number so #- A is the unary negation of A and
;;; A ##- B is the infix subtraction of B from A. It would be possible to
;;; overload the standard Pop11 operators, but this would severly reduce
;;; performance of the general system.
;;;
;;; Note: transarithmetical operators have the same precedence as the
;;; corresponding Pop11 operators.

define 1 #- (num);
    lvars num;
    trans_negate(num);
enddefine;

define 5 ##+ (num1, num2);
    lvars num1, num2;
    trans_add(num1, num2);
enddefine;

define 5 ##- (num1, num2);
    lvars num1, num2;
    trans_sub(num1, num2);
enddefine;

define 4 ##* (num1, num2);
    lvars num1, num2;
    trans_mul(num1, num2);
enddefine;

define 4 ##/ (num1, num2);
    lvars num1, num2;
    trans_div(num1, num2);
enddefine;

;;; In my view it is a bug that Pop11's power operator has precedence 3.
;;; Pop11 has: 2 ** 1/2 = 1 but 2 ** (1/2) = 1.414214 Nonetheless, I have
;;; cloned the precedence 3 here.
define 3 ##** (num1, num2);
    lvars num1, num2;
    trans_power(num1, num2);
enddefine;



/*****************************************************************************
*                      Transreal Arithmetical Functions.                     *
*****************************************************************************/

;;; The functions supplied here are not needed to implement transcomplex
;;; arithmetic, but users might find them helpful.

;;; Return the transreal absolute value of x.
define tr_abs(x);
    lvars x;

    checktransreal(x);

    if     x == nuly or x == pinf
    then   x
    elseif x == ninf
    then   pinf
    else   abs(x)
    endif;
enddefine;

;;; Return the transreal sign of x.
define tr_sgn(x);
    lvars x;

    checktransreal(x);

    if     x == nuly
    then   nuly
    elseif x == pinf
    then   1
    elseif x == ninf
    then   -1
    else   sign(x)
    endif;
enddefine;

;;; Return the transreal cosine of the angle, a.
define tr_cos(a);
    lvars a;
    if   isstrictlytransreal(a)
    then nuly
    else cos(a)
    endif;
enddefine;

;;; Return the transreal sine of the angle, a.
define tr_sin(a);
    lvars a;
    if   isstrictlytransreal(a)
    then nuly
    else sin(a)
    endif;
enddefine;

;;; Return the signed distance from transreal num2 to transreal num1.
;;; With real numbers this corresponds to num1 - num2.
define tr_signeddistance(num1, num2);
    lvars num1, num2;
    if   num1 = num2
    then 0
    else num1 ##- num2
    endif;
enddefine;

;;; Return the absolute distance between transreal num1 and transreal num2.
define tr_distance(num1, num2);
    lvars num1, num2;
    tr_abs(tr_signeddistance(num1, num2))
enddefine;

;;; Return the transreal exponential: tr_exp(x) -> y.
define tr_exp(x);
    lvars x;

    checktransreal(x);

    if     x == nuly or x == pinf
    then   x
    elseif x == ninf
    then   0
    else   exp(x)
    endif;
enddefine;

;;; Return the transreal logarithm: tr_log(x) -> y.
;;; This is the inverse of the transreal exponential with the totalisation
;;; that all x in the input class, but outside the domain, return y = nuly.
define tr_log(x);
    lvars x;
    if     x ##< 0 or x == nuly
    then   nuly
    elseif x == pinf
    then   pinf
    elseif x = 0
    then   ninf
    else   log(x)
    endif;
enddefine;

;;; Return the transabsolute sign of num1.
define trans_sgn(num);
    lvars num, r, c, s;
    desttranscomplex(num) -> (r,c,s);
    if r == pinf then 2 else 1 endif -> r;
    constranscomplex(r,c,s)
enddefine;

;;; Return the transabsolute value of num1.
define trans_abs(num);
    lvars num, r, c, s;
    desttranscomplex(num) -> (r,c,s);
    constranscomplex(r,1,0)
enddefine;

;;; Return the transcomplex square root of num1.
define trans_sqrt(num);
    lvars num;
    trans_power(num, 1/2);
enddefine;




/*****************************************************************************
*                             Testing Procedures.                            *
*****************************************************************************/

;;; These procedures compute the numerical examples given in the companion
;;; journal paper. The page numebrs related to the text of the paper as
;;; submitted for review. These page numbers might be different in the
;;; published paper. For example, the published paper might not start on page
;;; one.

define test_1();
    'Page  6' =>
    [nuly => %nuly%] =>
    [pinf => %pinf%] =>
    [ninf => %ninf%] =>
    [pinf ##/ pinf => %pinf ##/ pinf%] =>
    [pinf ##- pinf => %pinf ##- pinf%] =>
    [0 ##** 0 => %0 ##** 0%] =>
    [trans_exp(ninf) => %trans_exp(ninf)%] =>
    [trans_log(0) => %trans_log(0)%] =>
    [nuly = 0 ##** 0 => %nuly = 0 ##** 0%] =>
    [    nuly = trans_exp(trans_log(0 ##** 0)) =>
        %nuly = trans_exp(trans_log(0 ##** 0))%
    ] =>
    [    nuly = trans_exp(0 ##* trans_log(0)) =>
        %nuly = trans_exp(0 ##* trans_log(0))%
    ] =>
    nl(1);
enddefine;

define test_2();
    'Page 11' =>
    [#-0 = 0 => % #-0 = 0%] =>
    [ninf = ninf => %ninf = ninf%] =>
    [#-nuly = nuly => % #-nuly = nuly%] =>
    [trans_sqrt(0) = 0 => %trans_sqrt(0) = 0%] =>
    [trans_sqrt(#-0) = 0 => %trans_sqrt(#-0) = 0%] =>
    [0 ##/ 1 = 0 => %0 ##/ 1 = 0%] =>
    [0 ##/ #-1 = 0 => %0 ##/ #-1 = 0%] =>
    [1 ##/ 0 = pinf => %1 ##/ 0 = pinf%] =>
    [#-1 ##/ 0 = ninf => %#-1 ##/ 0 = ninf%] =>
    [0 ##/ 0 = nuly => %0 ##/ 0 = nuly%] =>
    [pinf ##/ pinf = nuly => %pinf ##/ pinf = nuly%] =>
    [ninf ##/ pinf = nuly => %ninf ##/ pinf = nuly%] =>
    [pinf ##/ ninf = nuly => %pinf ##/ ninf = nuly%] =>
    [ninf ##/ ninf = nuly => %ninf ##/ ninf = nuly%] =>
    [pinf ##* 0 = nuly => %pinf ##* 0 = nuly%] =>
    [ninf ##* 0 = nuly => %ninf ##* 0 = nuly%] =>
    [0 ##* pinf = nuly => %0 ##* pinf = nuly%] =>
    [0 ##* ninf = nuly => %0 ##* ninf = nuly%] =>
    [pinf ##- pinf = nuly => %pinf ##- pinf = nuly%] =>
    nl(1);
enddefine;

define test_3();
    'Page 12' =>
    [pinf ##+ 1 = pinf => %pinf ##+ 1 = pinf%] =>
    nl(1);
enddefine;

define test_4();
    'Page 13' =>
    [1 ##/ (pinf ##+ 1) = 0 => %1 ##/ (pinf ##+ 1) = 0%] =>
    [nuly = 0 ##* tr_sin(pinf) => %nuly = 0 ##* tr_sin(pinf)%] =>
    nl(1);
enddefine;

define test_5();
    'Page 15' =>
    [1 = 1 ##/ 1 => %1 = 1 ##/ 1%] =>
    [#-1 = #-1 ##/ 1 => %#-1 = #-1 ##/ 1%] =>
    nl(1);
enddefine;

define test_6();
    'Page 16' =>
    [0 = 0 ##/ 1 => %0 = 0 ##/ 1%] =>
    [pinf = 1 ##/ 0 => %pinf = 1 ##/ 0%] =>
    [ninf = #-1 ##/ 0 => %ninf = #-1 ##/ 0%] =>
    [nuly = 0 ##/ 0 => %nuly = 0 ##/ 0%] =>
    [nuly = pinf ##/ pinf => %nuly = pinf ##/ pinf%] =>
    [ninf = pinf ##/ #-3 => %ninf = pinf ##/ #-3%] =>
    nl(1);
enddefine;

define test_7();
    'Page 17' =>
    [nuly = pinf ##- pinf => %nuly = pinf ##- pinf%] =>
    [nuly = nuly ##- pinf => %nuly = nuly ##- pinf%] =>
    nl(1);
enddefine;

define test_8();
    'Page 18' =>
    [nuly = pinf ##* (2 ##- 2) => %nuly = pinf ##* (2 ##- 2)%] =>
    [nuly = pinf ##* 2 ##- pinf ##* 2 => %nuly = pinf ##* 2 ##- pinf ##* 2%] =>
    [pinf = pinf ##* (2 ##- 1) => %pinf = pinf ##* (2 ##- 1)%] =>
    [nuly = pinf ##* 2 ##- pinf ##* 1 => %nuly = pinf ##* 2 ##- pinf ##* 1%] =>
    [nuly = pinf ##* 0 => %nuly = pinf ##* 0%] =>
    [nuly = pinf ##- pinf => %nuly = pinf ##- pinf%] =>
    nl(1);
enddefine;

define test_9();
    'Page 19' =>
    nl(1);

    'Gravitational attraction between P_1 and P_2.' =>
    [ ((6.7 ##* 10 ##** #-11) ##*
       (3.2 ##* 10 ##** #-27) ##*
       (3.2 ##* 10 ##** #-27)
      ) ##/ (0 ##** 2) =>
     %((6.7 ##* 10 ##** #-11) ##*
       (3.2 ##* 10 ##** #-27) ##*
       (3.2 ##* 10 ##** #-27)
      ) ##/ (0 ##** 2)
     %
    ] =>
    nl(1);

    'Gravitational attraction between singularity and P_3.' =>
    [ ((6.7 ##* 10 ##** #-11) ##*
       (2 ##* 3.2 ##* 10 ##** #-27) ##*
       (3.2 ##* 10 ##** #-27)
      ) ##/ ((1.6 ##* 10 ##* #-35) ##** 2) =>
     %((6.7 ##* 10 ##** #-11) ##*
       (2 ##* 3.2 ##* 10 ##** #-27) ##*
       (3.2 ##* 10 ##** #-27)
      ) ##/ ((1.6 ##* 10 ##** #-35) ##** 2)
     %
    ] =>
    nl(1);

    'Electrostatic repulsion between P_1 and P_2.' =>
    [ ((9 ##* 10 ##** 9) ##*
       (1.6 ##* 10 ##** #-19) ##*
       (1.6 ##* 10 ##** #-19)
      ) ##/ (0 ##** 2) =>
     %((9 ##* 10 ##** 9) ##*
       (1.6 ##* 10 ##** #-19) ##*
       (1.6 ##* 10 ##** #-19)
      ) ##/ (0 ##** 2)
     %
    ] =>
    nl(1);

    'Electrostatic repulsion between singularity and P_3.' =>
    [ ((9 ##* 10 ##** 9) ##*
       (2 ##* 1.6 ##* 10 ##** #-19) ##*
       (1.6 ##* 10 ##** #-19)
      ) ##/ ((1.6 #* 10 ##** #-35) ##** 2) =>
     %((9 ##* 10 ##** 9) ##*
       (2 ##* 1.6 ##* 10 ##** #-19) ##*
       (1.6 ##* 10 ##** #-19)
      ) ##/ ((1.6 ##* 10 ##** #-35) ##** 2)
     %
    ] =>
    nl(1);
enddefine;

define test_10();
    'Page 20' =>
    [nuly = pinf ##- pinf => %nuly = pinf ##- pinf%] =>
    [ 5.4 ##* 10 ##** 6 ##- 1.8 ##* 10 ##** 42 =>
     %5.4 ##* 10 ##** 6 ##- 1.8 ##* 10 ##** 42%
    ] =>
    nl(1);
enddefine;

define test_11();
    'Page 21' =>
    [pinf = pinf ##/ 0 => %pinf = pinf ##/ 0%] =>
    [nuly = pinf ##/ pinf => %nuly = pinf ##/ pinf%] =>
    nl(1);
enddefine;

define test_12();
    lvars F, m, a, row, col, k_m = 2, k_a = 3;

    'Page 22' =>
    'Table 3: F = ma with k_m = 2, k_a = 3' =>

    [% 0, k_m, pinf, nuly %] -> row;
    [% ninf, #-k_a, 0, k_a, pinf, nuly %] -> col;

    for m in row do
        for a in col do
            m ##* a -> F;
            [m = ^m a = ^a F = ^F] =>
        endfor;
    endfor;

    nl(1);
enddefine;

define test_13();
    lvars F, m, a, row, col, k_F = 2, k_m = 3;

    'Page 22' =>
    'Table 4: a = F/m with k_F = 2, k_m = 3' =>

    [% ninf, #-k_F, 0, k_F, pinf, nuly %] -> row;
    [% 0, k_m, pinf, nuly %] -> col;

    for F in row do
        for m in col do
            F ##/ m -> a;
            [F = ^F m = ^m a = ^a] =>
        endfor;
    endfor;

    nl(1);
enddefine;

define test_14();
    lvars F, m, a, row, col, k_F = 2, k_a = 3;

    'Page 23' =>
    'Table 5: m = F/a with k_F = 2, k_a = 3' =>

    [% ninf, #-k_F, 0, k_F, pinf, nuly %] -> row;
    [% ninf, #-k_a, 0, k_a, pinf, nuly %] -> col;

    for F in row do
        for a in col do
            F ##/ a -> m;
            [F = ^F a = ^a m = ^m] =>
        endfor;
    endfor;

    nl(1);
enddefine;

define test_15();
    dlocal pop_pr_places = 16;

    'Page 25' =>

    nl(1);
    'Exact integer computation:' =>
    [2 ** ( 10 + 1) - 2 => %2 ** ( 10 + 1) - 2%] =>
    [2 ** ( 23 + 1) - 2 => %2 ** ( 23 + 1) - 2%] =>
    [2 ** ( 52 + 1) - 2 => %2 ** ( 52 + 1) - 2%] =>
    [2 ** (112 + 1) - 2 => %2 ** (112 + 1) - 2%] =>

    nl(1);
    'Floating-point computation following a transcomplex path:' =>
    [2 ##** ( 10 ##+ 1) ##- 2 => %2 ##** ( 10 ##+ 1) ##- 2%] =>
    [2 ##** ( 23 ##+ 1) ##- 2 => %2 ##** ( 23 ##+ 1) ##- 2%] =>
    [2 ##** ( 52 ##+ 1) ##- 2 => %2 ##** ( 52 ##+ 1) ##- 2%] =>
    [2 ##** (112 ##+ 1) ##- 2 => %2 ##** (112 ##+ 1) ##- 2%] =>

    nl(1);
    'No signed zero:' =>
    [0 = #-0 and 1 ##/ 0 = 1 ##/ #- 0 => %0 = #-0 and 1 ##/ 0 = 1 ##/ #- 0%] =>
    nl(1);
enddefine;

define test_16();
    lvars tc_i_inf = constranscomplex(pinf, tr_cos(pi/2), tr_sin(pi/2));
    'Page 34' =>
    [^tc_pinf + ^tc_i_inf => %tc_pinf ##+ tc_i_inf%] =>
    nl(1);
enddefine;

define test_17();
    'Page 41' =>
    [^pinf ##/ 0 = ^pinf => %tc_pinf ##/ 0 = pinf%] =>
    nl(1);
enddefine;

define test_18();
    'Page 43' =>
    [ nuly ##/ tr_sqrt(nuly##**2 ##+ nuly##**2) = nuly =>
     %nuly ##/ tr_sqrt(nuly##**2 ##+ nuly##**2) = nuly%
    ] =>
    nl(1);
enddefine;

define test_19();
    'Page 46' =>

    [ pinf ##- pinf = nuly =>
     %pinf ##- pinf = nuly%] =>

    [ ninf ##+ pinf = nuly =>
     %ninf ##+ pinf = nuly%] =>

    [ nuly ##- nuly = nuly =>
     %nuly ##- nuly = nuly%] =>

    [ #-nuly ##+ nuly = nuly =>
     %#-nuly ##+ nuly = nuly%] =>

    [ tr_cos(ninf) = nuly and tr_sin(ninf) = nuly =>
     %tr_cos(ninf) = nuly and tr_sin(ninf) = nuly%] =>

    [ tr_cos(pinf) = nuly and tr_sin(pinf) = nuly =>
     %tr_cos(pinf) = nuly and tr_sin(pinf) = nuly%] =>

    [ tr_cos(nuly) = nuly and tr_sin(nuly) = nuly =>
     %tr_cos(nuly) = nuly and tr_sin(nuly) = nuly%] =>

    nl(1);
enddefine;

define test_20();
    lvars e, a1, a2, z1, z2;
    dlocal pop_pr_places = 16;

    'Page 51' =>
    'This is not a numerical example in the paper, but it is important' =>

    ;;; Create variables for opposite numbers.
    constranscomplex(1, nuly, nuly) -> z1;
    constranscomplex(1, nuly, nuly) -> z2;

    ;;; Test over two rotations.
    for e from -2 by 1/8 to 2 do

         ;;; Compute exactly opposite angles.
         pi/2 + e -> a1;
        -pi/2 + e -> a2;

        ;;; Compute exactly opposite numbers.
        tr_cos(a1) -> tc_cos(z1);
        tr_sin(a1) -> tc_sin(z1);
        tr_cos(a2) -> tc_cos(z2);
        tr_sin(a2) -> tc_sin(z2);

        nl(1);

        [z1 = ^z1] =>
        [z2 = ^z2] =>
        [z1 ##+ z2 => %z1 ##+ z2%] =>
    endfor;

    nl(1);
enddefine;

define test_21();
    lvars inf;

    'Page 52' =>
    constranscomplex(pinf, nuly, nuly) -> inf;
    [^inf ##* ^inf = ^inf => %inf ##* inf = inf%] =>
    constranscomplex(pinf, sqrt(2), sqrt(2)) -> inf;
    [istranscomplex(^inf) => %istranscomplex(inf)%] =>
    [istransreal(^inf) => %istransreal(inf)%] =>
    constranscomplex(pinf, 1, 0) -> inf;
    [^inf ##* ^inf = ^pinf => %inf ##* inf = pinf%] =>
    nl(1);
enddefine;

define test_22();
    lvars inf;

    'Page 52' =>
    constranscomplex(pinf, 1, 0) -> inf;
    [^inf ##+ ^inf = ^pinf => %inf ##+ inf = pinf%] =>
    [^inf ##- ^inf = ^nuly => %inf ##- inf = nuly%] =>
    [^inf ##/ ^inf = ^nuly => %inf ##/ inf = nuly%] =>
    [   0 ##/    0 = ^nuly => %  0 ##/   0 = nuly%] =>
    [^inf ##*    0 = ^nuly => %inf ##*   0 = nuly%] =>
    nl(1);
enddefine;

section_cancel(current_section);

endsection;
endsection;