/* Copyright:      (c) 1997-2003 James A.D.W. Anderson. All rights reserved.
 > File:           $useperspex/lib/trigonometry.p
 > Release:        Persepx version 1, Pop11.
 > Purpose:        Provides transrational, trigonometric functions as in the following paper.
 >                 Aditionally provides tanq2 and arctanq2.
 > Documentation:  http://www.bookofparagon.btinternet.co.uk/Mathematics/SPIE.2002.Exact.pdf
 >                 Comments "EQN N" refer to equation number N in the above paper.
 > Note:           The paper ommits the sign convention needed for EQN 18.
 >                 Firstly, the integer square root is signed so that sign(sqrt(x)) = sign(x).
 >                 Secondly, r is always positive, so sign(r/q) = sign(q) and sign(r/p) = sign(p).
 > Related Files:  $useperspex/lib/transrational.p
*/


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

section $-library  => trigonometry
    intsqrt
    principal
    pqr arcpqr
       cosq    sinq    tanq    secq    cscq    cotq
    arccosq arcsinq arctanq arcsecq arccscq arccotq;

;;; Load related file.
uses transrational;

;;; Declare existance of library for 'uses'.
vars trigonometry = true;

section trigonometry =>
    intsqrt
    principal
    pqr arcpqr
       cosq    sinq    tanq    secq    cscq    cotq
    arccosq arcsinq arctanq arcsecq arccscq arccotq;


/***********************************************************************
*                            Base Procedures                           *
***********************************************************************/


;;; Signed integer square root.
;;;
;;; Note: the paper ommits to say that the integer square root is
;;; signed.
;;;
;;; The input argument, num, is an integer or else a strictly
;;; transrational number.
;;;
;;; The output parameter, root, is an integral or strictly transrational
;;; number such that sign(root) = sign(num) and root has the largest
;;; magnitude such that abs(root*root) <= num.
;;;
;;; The output parameter, exact, is true if abs(root*root) = num
;;; and is false otherwise.
;;;
;;; Note: more efficient implementation is possible.
define intsqrt(num) -> root -> exact;
    lvars sgn, lo, hi, tmp, num, root, exact = true;

    ;;; Early-out on transrational number.
    if   num == nullity or num == infinity
    then num  -> root;
         return;
    endif;

    ;;; Check argument type.
    unless isintegral(num)
    then   mishap('INTEGRAL NUMBER NEEDED', [^num])
    endunless;

    ;;; Force number positive and record sign.
    sign(num) -> sgn;
    abs(num)  -> num;

    ;;; Compute bounds on root.
    ;;; Note: bounds are powers of 2, so subsequent operations should be
    ;;; fast.
    if   (integer_length(num) div 2 ->> tmp) = 0
    then 0 -> lo;
         1 -> hi;
    else 1  << (tmp - 1) -> lo;
         lo << 2         -> hi
    endif;

    ;;; Early-out if either bound is exact.
    if lo * lo = num then sgn * lo -> root; return; endif;
    if hi * hi = num then sgn * hi -> root; return; endif;

    ;;; Binary search, early-out on success.
    until hi - lo <= 1 do

          ;;; Compute new estimate of root.
          (lo + hi) div 2 -> root;
          root * root -> tmp;

          ;;; Early-out on success or adjust bounds.
          if     num = tmp
          then   sgn * root -> root;
                 return;
          elseif num < tmp
          then   root -> hi
          else   root -> lo
          endif;
    enduntil;

    ;;; No exact integral root.
    sgn*lo -> root;
    false  -> exact;
enddefine;

;;; EQN 10: map all rational numbers, including infinity, onto the range
;;; (-2,2] and map nullity onto itself.
define mapq(ratio);
    lvars ratio, n, d;
    desttransrational(ratio) -> (n,d);
    if   abs(n) > d
    then 2*sign(n) - d/n
    else n ||/ d
    endif;
enddefine;

;;; EQN 11: map all rational numbers in the range (-2,2] onto the
;;; rational number line augmented with infinity and map nullity onto
;;; itself.
define arcmapq(ratio);
    lvars ratio, n, d;
    desttransrational(ratio) -> (n,d);
    if   abs(n) > d
    then |- ( |/ (ratio ||- 2*sign(n)) )
    else n ||/ d
    endif;
enddefine;

;;; EQN 12: return the ratio in the principal range (-2,2] U nullity.
define principal(ratio);
    lvars ratio, n, d, d_2, d_4;
    desttransrational(ratio) -> (n,d);

    ;;; Early-out on nullity.
    if ratio == nullity then return(nullity) endif;

    ;;; Early-out on infinity.
    if ratio == infinity then return(2) endif;

    ;;; Store two and four times d.
    d << 1   -> d_2;
    d_2 << 1 -> d_4;

    ;;; Return rational ratio in principal range.
    -(((d_2 - n) mod d_4) - d_2) / d;
enddefine;

;;; EQN 15: returns the Pythagorean tripple p, q, r.
;;; Note: more efficient implementation is possible.
define pqr(ratio) -> (p, q, r);
    lvars n, d, k, ratio, p, q, r;

    ;;; EQN 12: put ratio into its principal range.
    principal(ratio) -> ratio;

    ;;; EQN 11: Decode ratio.
    desttransrational(arcmapq(ratio)) -> (n,d);

    ;;; EQN 13: compute q.
    2*d*n -> q;

    ;;; EQN 13: compute p.
    ;;; Square n and d as a side effect.
    n*n -> n;
    d*d -> d;
    d - n -> p;

    ;;; EQN 13: compute r using side effect.
    n + d  -> r;

    ;;; EQN 14: remove common factors, if any.
    gcd_n(p, q, r, 3) -> k;
    p ||/ k -> p;
    q ||/ k -> q;
    r ||/ k -> r;
enddefine;

;;; Returns the inverse of pqr.
;;;
;;; Note: the checks are defensive; they should be ommitted from secure
;;; code, because checking the validity of Pythagorean tripples is very
;;; expensive.
define arcpqr(p, q, r);
    lvars p2, q2, r2,  n, d, p, q, r;

    ;;; Check type of arguments.
    unless isintegral(p)
    then   mishap('INTEGRAL NUMBER EXPECTED', [^p])
    endunless;
    unless isintegral(q)
    then   mishap('INTEGRAL NUMBER EXPECTED', [^q])
    endunless;
    unless isintegral(r)
    then   mishap('INTEGRAL NUMBER EXPECTED', [^r])
    endunless;

    ;;; Check r non-negative.
    if   r < 0
    then mishap('NON-NEGATIVE RADIUS, R, EXPECTED', [% r %])
    endif;

    ;;; Check true Pythagorean tripples given.
    ;;;
    ;;; Note: where the departure of p, q, r, from pythagorean tripples
    ;;; is small the recovered angle will be close to the true one.
    ;;; So could just set a flag to give a warning.
    unless p*p + q*q = r*r
    then   mishap('INVALID PYTHAGOREAN TRIPPLE', [% p, q, r %])
    endunless;

    ;;; EQN 16: compute inverse.
    if   p >= 0
    then q ||/ (r+p)
    else 2*paritytransrational(q) - (q ||/ (r - p))
    endif;
enddefine;


/***********************************************************************
*                        Trigonometric Functions                       *
***********************************************************************/


;;; All trigonometric funtions with argument nullity return nullity.


;;; EQN 17: return the rational cosine.
;;; Principal range of cosq(x) is 0 < x <= 2.
;;; Note: cosq(infinity) = -1 and cosq(nullity) = nullity.
define cosq(ratio);
    lvars p, q, r, ratio;
    pqr(ratio) -> (p,q,r);
    p ||/ r;
enddefine;

;;; EQN 17: return the rational sine.
;;; Principal range of sinq(x) is -1 < x <= 1.
;;; Note: sinq(infinity) = 0 and sinq(nullity) = nullity.
define sinq(ratio);
    lvars p, q, r, ratio;
    pqr(ratio) -> (p,q,r);
    q ||/ r;
enddefine;

;;; EQN 17: return the rational tangent.
;;; Principal range of tanq(x) is -1 < x <= 1.
;;; Note: tanq(infinity) = 0 and tanq(nullity) = nullity.
define tanq(ratio);
    lvars p, q, r, ratio;
    pqr(ratio) -> (p,q,r);
    q ||/ p;
enddefine;

;;; EQN 17: return the rational secant.
;;; Principal range of secq(x) is 0 < x <= 2.
;;; Note: secq(infinity) = -1 and secq(nullity) = nullity.
define secq(ratio);
    lvars p, q, r, ratio;
    pqr(ratio) -> (p,q,r);
    r ||/ p;
enddefine;

;;; EQN 17: return the rational cosecant.
;;; Principal range of cscq(x) is -1 < x <= 1.
;;; Note: cscq(infinity) = infinity and cscq(nullity) = nullity.
define cscq(ratio);
    lvars p, q, r, ratio;
    pqr(ratio) -> (p,q,r);
    r ||/ q;
enddefine;

;;; EQN 17: return the rational cotangent.
;;; Principal range of cotq(x) is 0 <= x < 2.
;;; Note: cotq(infinity) = infinity and cotq(nullity) = nullity.
define cotq(ratio);
    lvars p, q, r, ratio;
    pqr(ratio) -> (p,q,r);
    p ||/ q;
enddefine;

;;; EQN 18: inverse of cosq.
define arccosq(c);
    lvars p, q, r, c;

    ;;; Check range of argument.
    if   c ||< -1 or c ||> 1
    then mishap('COSINE OUT OF RANGE', [^c])
    endif;

    ;;; Get p, r from the canonical ratio describing the cosine.
    desttransrational(c) -> (p,r);

    ;;; Compute the inverse.
    intsqrt((r + p)*(r - p)) -> q ->;
    arcpqr(p, q, r);
enddefine;

;;; EQN 18: inverse of sinq.
define arcsinq(s);
    lvars p, q, r, s;

    ;;; Check range of argument.
    if   s ||< -1 or s ||> 1
    then mishap('SINE OUT OF RANGE', [^s])
    endif;

    ;;; Get p, r from the canonical ratio describing the sine.
    desttransrational(s) -> (q,r);

    ;;; Compute the inverse.
    intsqrt((r + q)*(r - q)) -> p ->;
    arcpqr(p, q, r);
enddefine;

;;; EQN 18: inverse of tanq.
;;;
;;; Note EQN 19: arctanq(x) = tan(x/2) when x = tanq(y). Hence arctanq
;;; can be used to create half-angle functions.
;;;
;;; EQN 20: double and integral-n functions can be created simply by
;;; composing rotations and then using arctanq2 to obtain their
;;; parameter.
define arctanq(t);
    lvars p, q, r, t;

    ;;; Get p, q from the canonical ratio describing the tangent.
    desttransrational(t) -> (q,p);

    ;;; Compute the inverse.
    intsqrt(p*p + q*q) -> r ->;
    arcpqr(p, q, r);
enddefine;

;;; EQN 18: inverse of secq.
;;; The pepar ommits to say that r is always positive, so
;;; sign(r/q) = sign(q) and sign(r/p) = sign(p).
define arcsecq(c);
    lvars p, q, r, c;

    ;;; Get p, r from the canonical ratio describing the cosine.
    desttransrational(c) -> (r,p);

    ;;; Force the sign of r/p to be carried by p.
    if   r ||< 0
    then |- p -> p;
         |- r -> r;
    endif;

    ;;; Compute the inverse.
    intsqrt((r + p)*(r - p)) -> q ->;
    arcpqr(p, q, r);
enddefine;

;;; EQN 18: inverse of cscq.
;;; The pepar ommits to say that r is always positive, so
;;; sign(r/q) = sign(q) and sign(r/p) = sign(p).
define arccscq(s);
    lvars p, q, r, s;

    ;;; Get q, r from the canonical ratio describing the sine.
    desttransrational(s) -> (r,q);

    ;;; Force the sign of r/q to be carried by q.
    if   r ||< 0
    then |- q -> q;
         |- r -> r;
    endif;

    ;;; Compute the inverse.
    intsqrt((r + q)*(r - q)) -> p ->;
    arcpqr(p, q, r);
enddefine;

;;; EQN 18: inverse of cotq.
define arccotq(t);
    lvars p, q, r, t;

    ;;; Get p, q from the canonical ratio describing the tangent.
    desttransrational(t) -> (p,q);

    ;;; Compute the inverse.
    intsqrt(p*p + q*q) -> r ->;
    arcpqr(p, q, r);
enddefine;


/***********************************************************************
*               AdditionalTrigonometric Functions                      *
***********************************************************************/


;;; Return the cosine and sine.
;;; tanq2(ratio) -> (cosine, sine);
;;; Principal range of tanq2(x) is -2 < x <= 2.
define tanq2(ratio);
    lvars c, s, ratio, p, q, r;
    pqr(ratio) -> (p,q,r);
    p ||/ r;
    q ||/ r;
enddefine;

;;; Inverse of tanq2.
define arctanq2(c, s);
    lvars c, s, a;

    ;;; Compute positive solution.
    arctanq(abstransrational(s ||/ c)) -> a;

    ;;; Force solution to match the given signs of the cosine and sine.
    if   c ||> 0
    then if    s ||< 0
         then  |- a
         else  a
         endif
    else if    s ||< 0
         then -2 ||+ a
         else  2 ||- a
         endif
    endif;
enddefine;

section_cancel(current_section);

endsection;
endsection;