summaryrefslogtreecommitdiff
path: root/macros/latex/contrib/l3kernel/l3fp-trig.dtx
diff options
context:
space:
mode:
authorNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
committerNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
commite0c6872cf40896c7be36b11dcc744620f10adf1d (patch)
tree60335e10d2f4354b0674ec22d7b53f0f8abee672 /macros/latex/contrib/l3kernel/l3fp-trig.dtx
Initial commit
Diffstat (limited to 'macros/latex/contrib/l3kernel/l3fp-trig.dtx')
-rw-r--r--macros/latex/contrib/l3kernel/l3fp-trig.dtx1698
1 files changed, 1698 insertions, 0 deletions
diff --git a/macros/latex/contrib/l3kernel/l3fp-trig.dtx b/macros/latex/contrib/l3kernel/l3fp-trig.dtx
new file mode 100644
index 0000000000..085641bafd
--- /dev/null
+++ b/macros/latex/contrib/l3kernel/l3fp-trig.dtx
@@ -0,0 +1,1698 @@
+% \iffalse meta-comment
+%
+%% File: l3fp-trig.dtx
+%
+% Copyright (C) 2011-2019 The LaTeX3 Project
+%
+% It may be distributed and/or modified under the conditions of the
+% LaTeX Project Public License (LPPL), either version 1.3c of this
+% license or (at your option) any later version. The latest version
+% of this license is in the file
+%
+% https://www.latex-project.org/lppl.txt
+%
+% This file is part of the "l3kernel bundle" (The Work in LPPL)
+% and all files in that bundle must be distributed together.
+%
+% -----------------------------------------------------------------------
+%
+% The development version of the bundle can be found at
+%
+% https://github.com/latex3/latex3
+%
+% for those people who are interested.
+%
+%<*driver>
+\documentclass[full,kernel]{l3doc}
+\begin{document}
+ \DocInput{\jobname.dtx}
+\end{document}
+%</driver>
+% \fi
+%
+% \title{The \textsf{l3fp-trig} package\\
+% Floating point trigonometric functions}
+% \author{^^A
+% The \LaTeX3 Project\thanks
+% {^^A
+% E-mail:
+% \href{mailto:latex-team@latex-project.org}
+% {latex-team@latex-project.org}^^A
+% }^^A
+% }
+% \date{Released 2019-08-25}
+%
+% \maketitle
+%
+% \begin{documentation}
+%
+% \end{documentation}
+%
+% \begin{implementation}
+%
+% \section{\pkg{l3fp-trig} Implementation}
+%
+% \begin{macrocode}
+%<*initex|package>
+% \end{macrocode}
+%
+% \begin{macrocode}
+%<@@=fp>
+% \end{macrocode}
+%
+% \begin{macro}[EXP]
+% {
+% \@@_parse_word_acos:N ,
+% \@@_parse_word_acosd:N ,
+% \@@_parse_word_acsc:N ,
+% \@@_parse_word_acscd:N ,
+% \@@_parse_word_asec:N ,
+% \@@_parse_word_asecd:N ,
+% \@@_parse_word_asin:N ,
+% \@@_parse_word_asind:N ,
+% \@@_parse_word_cos:N ,
+% \@@_parse_word_cosd:N ,
+% \@@_parse_word_cot:N ,
+% \@@_parse_word_cotd:N ,
+% \@@_parse_word_csc:N ,
+% \@@_parse_word_cscd:N ,
+% \@@_parse_word_sec:N ,
+% \@@_parse_word_secd:N ,
+% \@@_parse_word_sin:N ,
+% \@@_parse_word_sind:N ,
+% \@@_parse_word_tan:N ,
+% \@@_parse_word_tand:N ,
+% }
+% Unary functions.
+% \begin{macrocode}
+\tl_map_inline:nn
+ {
+ {acos} {acsc} {asec} {asin}
+ {cos} {cot} {csc} {sec} {sin} {tan}
+ }
+ {
+ \cs_new:cpx { @@_parse_word_#1:N }
+ {
+ \exp_not:N \@@_parse_unary_function:NNN
+ \exp_not:c { @@_#1_o:w }
+ \exp_not:N \use_i:nn
+ }
+ \cs_new:cpx { @@_parse_word_#1d:N }
+ {
+ \exp_not:N \@@_parse_unary_function:NNN
+ \exp_not:c { @@_#1_o:w }
+ \exp_not:N \use_ii:nn
+ }
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]
+% {
+% \@@_parse_word_acot:N , \@@_parse_word_acotd:N,
+% \@@_parse_word_atan:N , \@@_parse_word_atand:N,
+% }
+% Those functions may receive a variable number of arguments.
+% \begin{macrocode}
+\cs_new:Npn \@@_parse_word_acot:N
+ { \@@_parse_function:NNN \@@_acot_o:Nw \use_i:nn }
+\cs_new:Npn \@@_parse_word_acotd:N
+ { \@@_parse_function:NNN \@@_acot_o:Nw \use_ii:nn }
+\cs_new:Npn \@@_parse_word_atan:N
+ { \@@_parse_function:NNN \@@_atan_o:Nw \use_i:nn }
+\cs_new:Npn \@@_parse_word_atand:N
+ { \@@_parse_function:NNN \@@_atan_o:Nw \use_ii:nn }
+% \end{macrocode}
+% \end{macro}
+%
+% \subsection{Direct trigonometric functions}
+%
+% The approach for all trigonometric functions (sine, cosine, tangent,
+% cotangent, cosecant, and secant), with arguments given in radians or
+% in degrees, is the same.
+% \begin{itemize}
+% \item Filter out special cases ($\pm 0$, $\pm\inf$ and \nan{}).
+% \item Keep the sign for later, and work with the absolute value
+% $\lvert x\rvert$ of the argument.
+% \item Small numbers ($\lvert x\rvert<1$ in radians, $\lvert
+% x\rvert<10$ in degrees) are converted to fixed point numbers (and
+% to radians if $\lvert x\rvert$ is in degrees).
+% \item For larger numbers, we need argument reduction. Subtract a
+% multiple of $\pi/2$ (in degrees,~$90$) to bring the number to the
+% range to $[0, \pi/2)$ (in degrees, $[0,90)$).
+% \item Reduce further to $[0, \pi/4]$ (in degrees, $[0,45]$) using
+% $\sin x = \cos (\pi/2-x)$, and when working in degrees, convert to
+% radians.
+% \item Use the appropriate power series depending on the octant
+% $\lfloor\frac{|x|}{\pi/4}\rfloor \mod 8$ (in degrees, the same
+% formula with $\pi/4\to 45$), the sign, and the function to
+% compute.
+% \end{itemize}
+%
+% \subsubsection{Filtering special cases}
+%
+% \begin{macro}[EXP]{\@@_sin_o:w}
+% This function, and its analogs for \texttt{cos}, \texttt{csc},
+% \texttt{sec}, \texttt{tan}, and \texttt{cot} instead of
+% \texttt{sin}, are followed either by \cs{use_i:nn} and a float in
+% radians or by \cs{use_ii:nn} and a float in degrees. The sine of
+% $\pm 0$ or \nan{} is the same float. The sine of $\pm\infty$ raises
+% an invalid operation exception with the appropriate function name.
+% Otherwise, call the \texttt{trig} function to perform argument
+% reduction and if necessary convert the reduced argument to radians.
+% Then, \cs{@@_sin_series_o:NNwwww} is called to compute the
+% Taylor series: this function receives a sign~|#3|, an initial octant
+% of~$0$, and the function \cs{@@_ep_to_float_o:wwN} which converts the
+% result of the series to a floating point directly rather than taking
+% its inverse, since $\sin(x) = \#3 \sin\lvert x\rvert$.
+% \begin{macrocode}
+\cs_new:Npn \@@_sin_o:w #1 \s_@@ \@@_chk:w #2#3#4; @
+ {
+ \if_case:w #2 \exp_stop_f:
+ \@@_case_return_same_o:w
+ \or: \@@_case_use:nw
+ {
+ \@@_trig:NNNNNwn #1 \@@_sin_series_o:NNwwww
+ \@@_ep_to_float_o:wwN #3 0
+ }
+ \or: \@@_case_use:nw
+ { \@@_invalid_operation_o:fw { #1 { sin } { sind } } }
+ \else: \@@_case_return_same_o:w
+ \fi:
+ \s_@@ \@@_chk:w #2 #3 #4;
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_cos_o:w}
+% The cosine of $\pm 0$ is $1$. The cosine of $\pm\infty$ raises an
+% invalid operation exception. The cosine of \nan{} is itself.
+% Otherwise, the \texttt{trig} function reduces the argument to at
+% most half a right-angle and converts if necessary to radians. We
+% then call the same series as for sine, but using a positive
+% sign~|0| regardless of the sign of~$x$, and with an initial octant
+% of~$2$, because $\cos(x) = + \sin(\pi/2 + \lvert x\rvert)$.
+% \begin{macrocode}
+\cs_new:Npn \@@_cos_o:w #1 \s_@@ \@@_chk:w #2#3; @
+ {
+ \if_case:w #2 \exp_stop_f:
+ \@@_case_return_o:Nw \c_one_fp
+ \or: \@@_case_use:nw
+ {
+ \@@_trig:NNNNNwn #1 \@@_sin_series_o:NNwwww
+ \@@_ep_to_float_o:wwN 0 2
+ }
+ \or: \@@_case_use:nw
+ { \@@_invalid_operation_o:fw { #1 { cos } { cosd } } }
+ \else: \@@_case_return_same_o:w
+ \fi:
+ \s_@@ \@@_chk:w #2 #3;
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_csc_o:w}
+% The cosecant of $\pm 0$ is $\pm \infty$ with the same sign, with a
+% division by zero exception (see \cs{@@_cot_zero_o:Nfw} defined
+% below), which requires the function name. The cosecant of
+% $\pm\infty$ raises an invalid operation exception. The cosecant of
+% \nan{} is itself. Otherwise, the \texttt{trig} function performs
+% the argument reduction, and converts if necessary to radians before
+% calling the same series as for sine, using the sign~|#3|, a starting
+% octant of~$0$, and inverting during the conversion from the fixed
+% point sine to the floating point result, because $\csc(x) = \#3
+% \big( \sin\lvert x\rvert\big)^{-1}$.
+% \begin{macrocode}
+\cs_new:Npn \@@_csc_o:w #1 \s_@@ \@@_chk:w #2#3#4; @
+ {
+ \if_case:w #2 \exp_stop_f:
+ \@@_cot_zero_o:Nfw #3 { #1 { csc } { cscd } }
+ \or: \@@_case_use:nw
+ {
+ \@@_trig:NNNNNwn #1 \@@_sin_series_o:NNwwww
+ \@@_ep_inv_to_float_o:wwN #3 0
+ }
+ \or: \@@_case_use:nw
+ { \@@_invalid_operation_o:fw { #1 { csc } { cscd } } }
+ \else: \@@_case_return_same_o:w
+ \fi:
+ \s_@@ \@@_chk:w #2 #3 #4;
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_sec_o:w}
+% The secant of $\pm 0$ is $1$. The secant of $\pm \infty$ raises an
+% invalid operation exception. The secant of \nan{} is itself.
+% Otherwise, the \texttt{trig} function reduces the argument and turns
+% it to radians before calling the same series as for sine, using a
+% positive sign~$0$, a starting octant of~$2$, and inverting upon
+% conversion, because $\sec(x) = + 1 / \sin(\pi/2 + \lvert x\rvert)$.
+% \begin{macrocode}
+\cs_new:Npn \@@_sec_o:w #1 \s_@@ \@@_chk:w #2#3; @
+ {
+ \if_case:w #2 \exp_stop_f:
+ \@@_case_return_o:Nw \c_one_fp
+ \or: \@@_case_use:nw
+ {
+ \@@_trig:NNNNNwn #1 \@@_sin_series_o:NNwwww
+ \@@_ep_inv_to_float_o:wwN 0 2
+ }
+ \or: \@@_case_use:nw
+ { \@@_invalid_operation_o:fw { #1 { sec } { secd } } }
+ \else: \@@_case_return_same_o:w
+ \fi:
+ \s_@@ \@@_chk:w #2 #3;
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_tan_o:w}
+% The tangent of $\pm 0$ or \nan{} is the same floating point number.
+% The tangent of $\pm\infty$ raises an invalid operation exception.
+% Once more, the \texttt{trig} function does the argument reduction
+% step and conversion to radians before calling
+% \cs{@@_tan_series_o:NNwwww}, with a sign~|#3| and an initial octant
+% of~$1$ (this shift is somewhat arbitrary). See \cs{@@_cot_o:w} for
+% an explanation of the $0$~argument.
+% \begin{macrocode}
+\cs_new:Npn \@@_tan_o:w #1 \s_@@ \@@_chk:w #2#3#4; @
+ {
+ \if_case:w #2 \exp_stop_f:
+ \@@_case_return_same_o:w
+ \or: \@@_case_use:nw
+ {
+ \@@_trig:NNNNNwn #1
+ \@@_tan_series_o:NNwwww 0 #3 1
+ }
+ \or: \@@_case_use:nw
+ { \@@_invalid_operation_o:fw { #1 { tan } { tand } } }
+ \else: \@@_case_return_same_o:w
+ \fi:
+ \s_@@ \@@_chk:w #2 #3 #4;
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_cot_o:w}
+% \begin{macro}[EXP]{\@@_cot_zero_o:Nfw}
+% The cotangent of $\pm 0$ is $\pm \infty$ with the same sign, with a
+% division by zero exception (see \cs{@@_cot_zero_o:Nfw}. The
+% cotangent of $\pm\infty$ raises an invalid operation exception. The
+% cotangent of \nan{} is itself. We use $\cot x = - \tan (\pi/2 +
+% x)$, and the initial octant for the tangent was chosen to be $1$, so
+% the octant here starts at $3$. The change in sign is obtained by
+% feeding \cs{@@_tan_series_o:NNwwww} two signs rather than just the
+% sign of the argument: the first of those indicates whether we
+% compute tangent or cotangent. Those signs are eventually combined.
+% \begin{macrocode}
+\cs_new:Npn \@@_cot_o:w #1 \s_@@ \@@_chk:w #2#3#4; @
+ {
+ \if_case:w #2 \exp_stop_f:
+ \@@_cot_zero_o:Nfw #3 { #1 { cot } { cotd } }
+ \or: \@@_case_use:nw
+ {
+ \@@_trig:NNNNNwn #1
+ \@@_tan_series_o:NNwwww 2 #3 3
+ }
+ \or: \@@_case_use:nw
+ { \@@_invalid_operation_o:fw { #1 { cot } { cotd } } }
+ \else: \@@_case_return_same_o:w
+ \fi:
+ \s_@@ \@@_chk:w #2 #3 #4;
+ }
+\cs_new:Npn \@@_cot_zero_o:Nfw #1#2#3 \fi:
+ {
+ \fi:
+ \token_if_eq_meaning:NNTF 0 #1
+ { \exp_args:NNf \@@_division_by_zero_o:Nnw \c_inf_fp }
+ { \exp_args:NNf \@@_division_by_zero_o:Nnw \c_minus_inf_fp }
+ {#2}
+ }
+% \end{macrocode}
+% \end{macro}
+% \end{macro}
+%
+% \subsubsection{Distinguishing small and large arguments}
+%
+% \begin{macro}[EXP]{\@@_trig:NNNNNwn}
+% The first argument is \cs{use_i:nn} if the operand is in radians and
+% \cs{use_ii:nn} if it is in degrees. Arguments |#2| to~|#5| control
+% what trigonometric function we compute, and |#6| to~|#8| are pieces
+% of a normal floating point number. Call the \texttt{_series}
+% function~|#2|, with arguments |#3|, either a conversion function
+% (\cs{@@_ep_to_float_o:wN} or \cs{@@_ep_inv_to_float_o:wN}) or a sign $0$
+% or~$2$ when computing tangent or cotangent; |#4|, a sign $0$ or~$2$;
+% the octant, computed in an integer expression starting with~|#5| and
+% stopped by a period; and a fixed point number obtained from the
+% floating point number by argument reduction (if necessary) and
+% conversion to radians (if necessary). Any argument reduction
+% adjusts the octant accordingly by leaving a (positive) shift into
+% its integer expression. Let us explain the integer comparison. Two
+% of the four \cs{exp_after:wN} are expanded, the expansion hits the
+% test, which is true if the float is at least~$1$ when working in
+% radians, and at least $10$ when working in degrees. Then one of the
+% remaining \cs{exp_after:wN} hits |#1|, which picks the \texttt{trig}
+% or \texttt{trigd} function in whichever branch of the conditional
+% was taken. The final \cs{exp_after:wN} closes the conditional. At
+% the end of the day, a number is \texttt{large} if it is $\geq 1$ in
+% radians or $\geq 10$ in degrees, and \texttt{small} otherwise. All
+% four \texttt{trig}/\texttt{trigd} auxiliaries receive the operand as
+% an extended-precision number.
+% \begin{macrocode}
+\cs_new:Npn \@@_trig:NNNNNwn #1#2#3#4#5 \s_@@ \@@_chk:w 1#6#7#8;
+ {
+ \exp_after:wN #2
+ \exp_after:wN #3
+ \exp_after:wN #4
+ \int_value:w \@@_int_eval:w #5
+ \exp_after:wN \exp_after:wN \exp_after:wN \exp_after:wN
+ \if_int_compare:w #7 > #1 0 1 \exp_stop_f:
+ #1 \@@_trig_large:ww \@@_trigd_large:ww
+ \else:
+ #1 \@@_trig_small:ww \@@_trigd_small:ww
+ \fi:
+ #7,#8{0000}{0000};
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \subsubsection{Small arguments}
+%
+% \begin{macro}[EXP]{\@@_trig_small:ww}
+% This receives a small extended-precision number in radians and
+% converts it to a fixed point number. Some trailing digits may be
+% lost in the conversion, so we keep the original floating point
+% number around: when computing sine or tangent (or their inverses),
+% the last step is to multiply by the floating point number (as
+% an extended-precision number) rather than the fixed point number.
+% The period serves to end the integer expression for the octant.
+% \begin{macrocode}
+\cs_new:Npn \@@_trig_small:ww #1,#2;
+ { \@@_ep_to_fixed:wwn #1,#2; . #1,#2; }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_trigd_small:ww}
+% Convert the extended-precision number to radians, then call
+% \cs{@@_trig_small:ww} to massage it in the form appropriate for the
+% \texttt{_series} auxiliary.
+% \begin{macrocode}
+\cs_new:Npn \@@_trigd_small:ww #1,#2;
+ {
+ \@@_ep_mul_raw:wwwwN
+ -1,{1745}{3292}{5199}{4329}{5769}{2369}; #1,#2;
+ \@@_trig_small:ww
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \subsubsection{Argument reduction in degrees}
+%
+% \begin{macro}[rEXP]
+% {
+% \@@_trigd_large:ww, \@@_trigd_large_auxi:nnnnwNNNN,
+% \@@_trigd_large_auxii:wNw, \@@_trigd_large_auxiii:www
+% }
+% Note that $25\times 360 = 9000$, so $10^{k+1} \equiv 10^{k}
+% \pmod{360}$ for $k\geq 3$. When the exponent~|#1| is very large, we
+% can thus safely replace it by~$22$ (or even~$19$). We turn the
+% floating point number into a fixed point number with two blocks of
+% $8$~digits followed by five blocks of $4$~digits. The original
+% float is $100\times\meta{block_1}\cdots\meta{block_3}.
+% \meta{block_4}\cdots\meta{block_7}$, or is equal to it modulo~$360$
+% if the exponent~|#1| is very large. The first auxiliary finds
+% $\meta{block_1} + \meta{block_2} \pmod{9}$, a single digit, and
+% prepends it to the $4$~digits of \meta{block_3}. It also unpacks
+% \meta{block_4} and grabs the $4$~digits of \meta{block_7}. The
+% second auxiliary grabs the \meta{block_3} plus any contribution from
+% the first two blocks as~|#1|, the first digit of \meta{block_4}
+% (just after the decimal point in hundreds of degrees) as~|#2|, and
+% the three other digits as~|#3|. It finds the quotient and remainder
+% of |#1#2| modulo~$9$, adds twice the quotient to the integer
+% expression for the octant, and places the remainder (between $0$
+% and~$8$) before |#3| to form a new \meta{block_4}. The resulting
+% fixed point number is $x\in [0, 0.9]$. If $x\geq 0.45$, we add~$1$
+% to the octant and feed $0.9-x$ with an exponent of~$2$ (to
+% compensate the fact that we are working in units of hundreds of
+% degrees rather than degrees) to \cs{@@_trigd_small:ww}. Otherwise,
+% we feed it~$x$ with an exponent of~$2$. The third auxiliary also
+% discards digits which were not packed into the various
+% \meta{blocks}. Since the original exponent~|#1| is at least~$2$,
+% those are all~$0$ and no precision is lost (|#6| and~|#7| are
+% four~$0$ each).
+% \begin{macrocode}
+\cs_new:Npn \@@_trigd_large:ww #1, #2#3#4#5#6#7;
+ {
+ \exp_after:wN \@@_pack_eight:wNNNNNNNN
+ \exp_after:wN \@@_pack_eight:wNNNNNNNN
+ \exp_after:wN \@@_pack_twice_four:wNNNNNNNN
+ \exp_after:wN \@@_pack_twice_four:wNNNNNNNN
+ \exp_after:wN \@@_trigd_large_auxi:nnnnwNNNN
+ \exp_after:wN ;
+ \exp:w \exp_end_continue_f:w
+ \prg_replicate:nn { \int_max:nn { 22 - #1 } { 0 } } { 0 }
+ #2#3#4#5#6#7 0000 0000 0000 !
+ }
+\cs_new:Npn \@@_trigd_large_auxi:nnnnwNNNN #1#2#3#4#5; #6#7#8#9
+ {
+ \exp_after:wN \@@_trigd_large_auxii:wNw
+ \int_value:w \@@_int_eval:w #1 + #2
+ - (#1 + #2 - 4) / 9 * 9 \@@_int_eval_end:
+ #3;
+ #4; #5{#6#7#8#9};
+ }
+\cs_new:Npn \@@_trigd_large_auxii:wNw #1; #2#3;
+ {
+ + (#1#2 - 4) / 9 * 2
+ \exp_after:wN \@@_trigd_large_auxiii:www
+ \int_value:w \@@_int_eval:w #1#2
+ - (#1#2 - 4) / 9 * 9 \@@_int_eval_end: #3 ;
+ }
+\cs_new:Npn \@@_trigd_large_auxiii:www #1; #2; #3!
+ {
+ \if_int_compare:w #1 < 4500 \exp_stop_f:
+ \exp_after:wN \@@_use_i_until_s:nw
+ \exp_after:wN \@@_fixed_continue:wn
+ \else:
+ + 1
+ \fi:
+ \@@_fixed_sub:wwn {9000}{0000}{0000}{0000}{0000}{0000};
+ {#1}#2{0000}{0000};
+ { \@@_trigd_small:ww 2, }
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \subsubsection{Argument reduction in radians}
+%
+% Arguments greater or equal to~$1$ need to be reduced to a range where
+% we only need a few terms of the Taylor series. We reduce to the range
+% $[0,2\pi]$ by subtracting multiples of~$2\pi$, then to the smaller
+% range $[0,\pi/2]$ by subtracting multiples of~$\pi/2$ (keeping track
+% of how many times~$\pi/2$ is subtracted), then to $[0,\pi/4]$ by
+% mapping $x\to \pi/2 - x$ if appropriate. When the argument is very
+% large, say, $10^{100}$, an equally large multiple of~$2\pi$ must be
+% subtracted, hence we must work with a very good approximation
+% of~$2\pi$ in order to get a sensible remainder modulo~$2\pi$.
+%
+% Specifically, we multiply the argument by an approximation
+% of~$1/(2\pi)$ with $\ExplSyntaxOn\int_eval:n { \c__fp_max_exponent_int
+% + 48 }\ExplSyntaxOff$~digits, then discard the integer part of the
+% result, keeping $52$~digits of the fractional part. From the
+% fractional part of $x/(2\pi)$ we deduce the octant (quotient of the
+% first three digits by~$125$). We then multiply by $8$ or~$-8$ (the
+% latter when the octant is odd), ignore any integer part (related to
+% the octant), and convert the fractional part to an extended precision
+% number, before multiplying by~$\pi/4$ to convert back to a value in
+% radians in $[0,\pi/4]$.
+%
+% It is possible to prove that given the precision of floating points
+% and their range of exponents, the $52$~digits may start at most with
+% $24$~zeros. The $5$~last digits are affected by carries from
+% computations which are not done, hence we are left with at least $52 -
+% 24 - 5 = 23$ significant digits, enough to round correctly up to
+% $0.6\cdot\text{ulp}$ in all cases.
+%
+% \begin{variable}[EXP]{\c_@@_trig_intarray}
+% This integer array stores blocks of $8$~decimals of
+% $10^{-16}/(2\pi)$. Each entry is $10^8$ plus an $8$~digit number
+% storing $8$ decimals. In total we store $10112$~decimals of
+% $10^{-16}/(2\pi)$. The number of decimals we really need is the
+% maximum exponent plus the number of digits we later need,~$52$,
+% plus~$12$ ($4-1$~groups of $4$~digits). The memory footprint ($1/2$
+% byte per digit) is the same as an earlier method of storing the data
+% as a control sequence name, but the major advantage is that we can
+% unpack specific subsets of the digits without unpacking the $10112$
+% decimals.
+% \begin{macrocode}
+\intarray_const_from_clist:Nn \c_@@_trig_intarray
+ {
+ 100000000, 100000000, 115915494, 130918953, 135768883, 176337251,
+ 143620344, 159645740, 145644874, 176673440, 158896797, 163422653,
+ 150901138, 102766253, 108595607, 128427267, 157958036, 189291184,
+ 161145786, 152877967, 141073169, 198392292, 139966937, 140907757,
+ 130777463, 196925307, 168871739, 128962173, 197661693, 136239024,
+ 117236290, 111832380, 111422269, 197557159, 140461890, 108690267,
+ 139561204, 189410936, 193784408, 155287230, 199946443, 140024867,
+ 123477394, 159610898, 132309678, 130749061, 166986462, 180469944,
+ 186521878, 181574786, 156696424, 110389958, 174139348, 160998386,
+ 180991999, 162442875, 158517117, 188584311, 117518767, 116054654,
+ 175369880, 109739460, 136475933, 137680593, 102494496, 163530532,
+ 171567755, 103220324, 177781639, 171660229, 146748119, 159816584,
+ 106060168, 103035998, 113391198, 174988327, 186654435, 127975507,
+ 100162406, 177564388, 184957131, 108801221, 199376147, 168137776,
+ 147378906, 133068046, 145797848, 117613124, 127314069, 196077502,
+ 145002977, 159857089, 105690279, 167851315, 125210016, 131774602,
+ 109248116, 106240561, 145620314, 164840892, 148459191, 143521157,
+ 154075562, 100871526, 160680221, 171591407, 157474582, 172259774,
+ 162853998, 175155329, 139081398, 117724093, 158254797, 107332871,
+ 190406999, 175907657, 170784934, 170393589, 182808717, 134256403,
+ 166895116, 162545705, 194332763, 112686500, 126122717, 197115321,
+ 112599504, 138667945, 103762556, 108363171, 116952597, 158128224,
+ 194162333, 143145106, 112353687, 185631136, 136692167, 114206974,
+ 169601292, 150578336, 105311960, 185945098, 139556718, 170995474,
+ 165104316, 123815517, 158083944, 129799709, 199505254, 138756612,
+ 194458833, 106846050, 178529151, 151410404, 189298850, 163881607,
+ 176196993, 107341038, 199957869, 118905980, 193737772, 106187543,
+ 122271893, 101366255, 126123878, 103875388, 181106814, 106765434,
+ 108282785, 126933426, 179955607, 107903860, 160352738, 199624512,
+ 159957492, 176297023, 159409558, 143011648, 129641185, 157771240,
+ 157544494, 157021789, 176979240, 194903272, 194770216, 164960356,
+ 153181535, 144003840, 168987471, 176915887, 163190966, 150696440,
+ 147769706, 187683656, 177810477, 197954503, 153395758, 130188183,
+ 186879377, 166124814, 195305996, 155802190, 183598751, 103512712,
+ 190432315, 180498719, 168687775, 194656634, 162210342, 104440855,
+ 149785037, 192738694, 129353661, 193778292, 187359378, 143470323,
+ 102371458, 137923557, 111863634, 119294601, 183182291, 196416500,
+ 187830793, 131353497, 179099745, 186492902, 167450609, 189368909,
+ 145883050, 133703053, 180547312, 132158094, 131976760, 132283131,
+ 141898097, 149822438, 133517435, 169898475, 101039500, 168388003,
+ 197867235, 199608024, 100273901, 108749548, 154787923, 156826113,
+ 199489032, 168997427, 108349611, 149208289, 103776784, 174303550,
+ 145684560, 183671479, 130845672, 133270354, 185392556, 120208683,
+ 193240995, 162211753, 131839402, 109707935, 170774965, 149880868,
+ 160663609, 168661967, 103747454, 121028312, 119251846, 122483499,
+ 111611495, 166556037, 196967613, 199312829, 196077608, 127799010,
+ 107830360, 102338272, 198790854, 102387615, 157445430, 192601191,
+ 100543379, 198389046, 154921248, 129516070, 172853005, 122721023,
+ 160175233, 113173179, 175931105, 103281551, 109373913, 163964530,
+ 157926071, 180083617, 195487672, 146459804, 173977292, 144810920,
+ 109371257, 186918332, 189588628, 139904358, 168666639, 175673445,
+ 114095036, 137327191, 174311388, 106638307, 125923027, 159734506,
+ 105482127, 178037065, 133778303, 121709877, 134966568, 149080032,
+ 169885067, 141791464, 168350828, 116168533, 114336160, 173099514,
+ 198531198, 119733758, 144420984, 116559541, 152250643, 139431286,
+ 144403838, 183561508, 179771645, 101706470, 167518774, 156059160,
+ 187168578, 157939226, 123475633, 117111329, 198655941, 159689071,
+ 198506887, 144230057, 151919770, 156900382, 118392562, 120338742,
+ 135362568, 108354156, 151729710, 188117217, 195936832, 156488518,
+ 174997487, 108553116, 159830610, 113921445, 144601614, 188452770,
+ 125114110, 170248521, 173974510, 138667364, 103872860, 109967489,
+ 131735618, 112071174, 104788993, 168886556, 192307848, 150230570,
+ 157144063, 163863202, 136852010, 174100574, 185922811, 115721968,
+ 100397824, 175953001, 166958522, 112303464, 118773650, 143546764,
+ 164565659, 171901123, 108476709, 193097085, 191283646, 166919177,
+ 169387914, 133315566, 150669813, 121641521, 100895711, 172862384,
+ 126070678, 145176011, 113450800, 169947684, 122356989, 162488051,
+ 157759809, 153397080, 185475059, 175362656, 149034394, 145420581,
+ 178864356, 183042000, 131509559, 147434392, 152544850, 167491429,
+ 108647514, 142303321, 133245695, 111634945, 167753939, 142403609,
+ 105438335, 152829243, 142203494, 184366151, 146632286, 102477666,
+ 166049531, 140657343, 157553014, 109082798, 180914786, 169343492,
+ 127376026, 134997829, 195701816, 119643212, 133140475, 176289748,
+ 140828911, 174097478, 126378991, 181699939, 148749771, 151989818,
+ 172666294, 160183053, 195832752, 109236350, 168538892, 128468247,
+ 125997252, 183007668, 156937583, 165972291, 198244297, 147406163,
+ 181831139, 158306744, 134851692, 185973832, 137392662, 140243450,
+ 119978099, 140402189, 161348342, 173613676, 144991382, 171541660,
+ 163424829, 136374185, 106122610, 186132119, 198633462, 184709941,
+ 183994274, 129559156, 128333990, 148038211, 175011612, 111667205,
+ 119125793, 103552929, 124113440, 131161341, 112495318, 138592695,
+ 184904438, 146807849, 109739828, 108855297, 104515305, 139914009,
+ 188698840, 188365483, 166522246, 168624087, 125401404, 100911787,
+ 142122045, 123075334, 173972538, 114940388, 141905868, 142311594,
+ 163227443, 139066125, 116239310, 162831953, 123883392, 113153455,
+ 163815117, 152035108, 174595582, 101123754, 135976815, 153401874,
+ 107394340, 136339780, 138817210, 104531691, 182951948, 179591767,
+ 139541778, 179243527, 161740724, 160593916, 102732282, 187946819,
+ 136491289, 149714953, 143255272, 135916592, 198072479, 198580612,
+ 169007332, 118844526, 179433504, 155801952, 149256630, 162048766,
+ 116134365, 133992028, 175452085, 155344144, 109905129, 182727454,
+ 165911813, 122232840, 151166615, 165070983, 175574337, 129548631,
+ 120411217, 116380915, 160616116, 157320000, 183306114, 160618128,
+ 103262586, 195951602, 146321661, 138576614, 180471993, 127077713,
+ 116441201, 159496011, 106328305, 120759583, 148503050, 179095584,
+ 198298218, 167402898, 138551383, 123957020, 180763975, 150429225,
+ 198476470, 171016426, 197438450, 143091658, 164528360, 132493360,
+ 143546572, 137557916, 113663241, 120457809, 196971566, 134022158,
+ 180545794, 131328278, 100552461, 132088901, 187421210, 192448910,
+ 141005215, 149680971, 113720754, 100571096, 134066431, 135745439,
+ 191597694, 135788920, 179342561, 177830222, 137011486, 142492523,
+ 192487287, 113132021, 176673607, 156645598, 127260957, 141566023,
+ 143787436, 129132109, 174858971, 150713073, 191040726, 143541417,
+ 197057222, 165479803, 181512759, 157912400, 125344680, 148220261,
+ 173422990, 101020483, 106246303, 137964746, 178190501, 181183037,
+ 151538028, 179523433, 141955021, 135689770, 191290561, 143178787,
+ 192086205, 174499925, 178975690, 118492103, 124206471, 138519113,
+ 188147564, 102097605, 154895793, 178514140, 141453051, 151583964,
+ 128232654, 106020603, 131189158, 165702720, 186250269, 191639375,
+ 115278873, 160608114, 155694842, 110322407, 177272742, 116513642,
+ 134366992, 171634030, 194053074, 180652685, 109301658, 192136921,
+ 141431293, 171341061, 157153714, 106203978, 147618426, 150297807,
+ 186062669, 169960809, 118422347, 163350477, 146719017, 145045144,
+ 161663828, 146208240, 186735951, 102371302, 190444377, 194085350,
+ 134454426, 133413062, 163074595, 113830310, 122931469, 134466832,
+ 185176632, 182415152, 110179422, 164439571, 181217170, 121756492,
+ 119644493, 196532222, 118765848, 182445119, 109401340, 150443213,
+ 198586286, 121083179, 139396084, 143898019, 114787389, 177233102,
+ 186310131, 148695521, 126205182, 178063494, 157118662, 177825659,
+ 188310053, 151552316, 165984394, 109022180, 163144545, 121212978,
+ 197344714, 188741258, 126822386, 102360271, 109981191, 152056882,
+ 134723983, 158013366, 106837863, 128867928, 161973236, 172536066,
+ 185216856, 132011948, 197807339, 158419190, 166595838, 167852941,
+ 124187182, 117279875, 106103946, 106481958, 157456200, 160892122,
+ 184163943, 173846549, 158993202, 184812364, 133466119, 170732430,
+ 195458590, 173361878, 162906318, 150165106, 126757685, 112163575,
+ 188696307, 145199922, 100107766, 176830946, 198149756, 122682434,
+ 179367131, 108412102, 119520899, 148191244, 140487511, 171059184,
+ 141399078, 189455775, 118462161, 190415309, 134543802, 180893862,
+ 180732375, 178615267, 179711433, 123241969, 185780563, 176301808,
+ 184386640, 160717536, 183213626, 129671224, 126094285, 140110963,
+ 121826276, 151201170, 122552929, 128965559, 146082049, 138409069,
+ 107606920, 103954646, 119164002, 115673360, 117909631, 187289199,
+ 186343410, 186903200, 157966371, 103128612, 135698881, 176403642,
+ 152540837, 109810814, 183519031, 121318624, 172281810, 150845123,
+ 169019064, 166322359, 138872454, 163073727, 128087898, 130041018,
+ 194859136, 173742589, 141812405, 167291912, 138003306, 134499821,
+ 196315803, 186381054, 124578934, 150084553, 128031351, 118843410,
+ 107373060, 159565443, 173624887, 171292628, 198074235, 139074061,
+ 178690578, 144431052, 174262641, 176783005, 182214864, 162289361,
+ 192966929, 192033046, 169332843, 181580535, 164864073, 118444059,
+ 195496893, 153773183, 167266131, 130108623, 158802128, 180432893,
+ 144562140, 147978945, 142337360, 158506327, 104399819, 132635916,
+ 168734194, 136567839, 101281912, 120281622, 195003330, 112236091,
+ 185875592, 101959081, 122415367, 194990954, 148881099, 175891989,
+ 108115811, 163538891, 163394029, 123722049, 184837522, 142362091,
+ 100834097, 156679171, 100841679, 157022331, 178971071, 102928884,
+ 189701309, 195339954, 124415335, 106062584, 139214524, 133864640,
+ 134324406, 157317477, 155340540, 144810061, 177612569, 108474646,
+ 114329765, 143900008, 138265211, 145210162, 136643111, 197987319,
+ 102751191, 144121361, 169620456, 193602633, 161023559, 162140467,
+ 102901215, 167964187, 135746835, 187317233, 110047459, 163339773,
+ 124770449, 118885134, 141536376, 100915375, 164267438, 145016622,
+ 113937193, 106748706, 128815954, 164819775, 119220771, 102367432,
+ 189062690, 170911791, 194127762, 112245117, 123546771, 115640433,
+ 135772061, 166615646, 174474627, 130562291, 133320309, 153340551,
+ 138417181, 194605321, 150142632, 180008795, 151813296, 175497284,
+ 167018836, 157425342, 150169942, 131069156, 134310662, 160434122,
+ 105213831, 158797111, 150754540, 163290657, 102484886, 148697402,
+ 187203725, 198692811, 149360627, 140384233, 128749423, 132178578,
+ 177507355, 171857043, 178737969, 134023369, 102911446, 196144864,
+ 197697194, 134527467, 144296030, 189437192, 154052665, 188907106,
+ 162062575, 150993037, 199766583, 167936112, 181374511, 104971506,
+ 115378374, 135795558, 167972129, 135876446, 130937572, 103221320,
+ 124605656, 161129971, 131027586, 191128460, 143251843, 143269155,
+ 129284585, 173495971, 150425653, 199302112, 118494723, 121323805,
+ 116549802, 190991967, 168151180, 122483192, 151273721, 199792134,
+ 133106764, 121874844, 126215985, 112167639, 167793529, 182985195,
+ 185453921, 106957880, 158685312, 132775454, 133229161, 198905318,
+ 190537253, 191582222, 192325972, 178133427, 181825606, 148823337,
+ 160719681, 101448145, 131983362, 137910767, 112550175, 128826351,
+ 183649210, 135725874, 110356573, 189469487, 154446940, 118175923,
+ 106093708, 128146501, 185742532, 149692127, 164624247, 183221076,
+ 154737505, 168198834, 156410354, 158027261, 125228550, 131543250,
+ 139591848, 191898263, 104987591, 115406321, 103542638, 190012837,
+ 142615518, 178773183, 175862355, 117537850, 169565995, 170028011,
+ 158412588, 170150030, 117025916, 174630208, 142412449, 112839238,
+ 105257725, 114737141, 123102301, 172563968, 130555358, 132628403,
+ 183638157, 168682846, 143304568, 105994018, 170010719, 152092970,
+ 117799058, 132164175, 179868116, 158654714, 177489647, 116547948,
+ 183121404, 131836079, 184431405, 157311793, 149677763, 173989893,
+ 102277656, 107058530, 140837477, 152640947, 143507039, 152145247,
+ 101683884, 107090870, 161471944, 137225650, 128231458, 172995869,
+ 173831689, 171268519, 139042297, 111072135, 107569780, 137262545,
+ 181410950, 138270388, 198736451, 162848201, 180468288, 120582913,
+ 153390138, 135649144, 130040157, 106509887, 192671541, 174507066,
+ 186888783, 143805558, 135011967, 145862340, 180595327, 124727843,
+ 182925939, 157715840, 136885940, 198993925, 152416883, 178793572,
+ 179679516, 154076673, 192703125, 164187609, 162190243, 104699348,
+ 159891990, 160012977, 174692145, 132970421, 167781726, 115178506,
+ 153008552, 155999794, 102099694, 155431545, 127458567, 104403686,
+ 168042864, 184045128, 181182309, 179349696, 127218364, 192935516,
+ 120298724, 169583299, 148193297, 183358034, 159023227, 105261254,
+ 121144370, 184359584, 194433836, 138388317, 175184116, 108817112,
+ 151279233, 137457721, 193398208, 119005406, 132929377, 175306906,
+ 160741530, 149976826, 147124407, 176881724, 186734216, 185881509,
+ 191334220, 175930947, 117385515, 193408089, 157124410, 163472089,
+ 131949128, 180783576, 131158294, 100549708, 191802336, 165960770,
+ 170927599, 101052702, 181508688, 197828549, 143403726, 142729262,
+ 110348701, 139928688, 153550062, 106151434, 130786653, 196085995,
+ 100587149, 139141652, 106530207, 100852656, 124074703, 166073660,
+ 153338052, 163766757, 120188394, 197277047, 122215363, 138511354,
+ 183463624, 161985542, 159938719, 133367482, 104220974, 149956672,
+ 170250544, 164232439, 157506869, 159133019, 137469191, 142980999,
+ 134242305, 150172665, 121209241, 145596259, 160554427, 159095199,
+ 168243130, 184279693, 171132070, 121049823, 123819574, 171759855,
+ 119501864, 163094029, 175943631, 194450091, 191506160, 149228764,
+ 132319212, 197034460, 193584259, 126727638, 168143633, 109856853,
+ 127860243, 132141052, 133076065, 188414958, 158718197, 107124299,
+ 159592267, 181172796, 144388537, 196763139, 127431422, 179531145,
+ 100064922, 112650013, 132686230, 121550837,
+ }
+% \end{macrocode}
+% \end{variable}
+%
+% \begin{macro}[rEXP]
+% {
+% \@@_trig_large:ww,
+% \@@_trig_large_auxi:w,
+% \@@_trig_large_auxii:w,
+% \@@_trig_large_auxiii:w,
+% }
+% The exponent~|#1| is between $1$ and~$\ExplSyntaxOn \int_use:N
+% \c__fp_max_exponent_int$. We wish to look up decimals
+% $10^{\text{\texttt{\#1}}-16}/(2\pi)$ starting from the digit
+% $|#1|+1$. Since they are stored in batches of~$8$, compute
+% $\lfloor|#1|/8\rfloor$ and fetch blocks of $8$ digits starting
+% there. The numbering of items in \cs{c_@@_trig_intarray} starts
+% at~$1$, so the block $\lfloor|#1|/8\rfloor+1$ contains the digit we
+% want, at one of the eight positions. Each call to \cs{int_value:w}
+% \cs{__kernel_intarray_item:Nn} expands the next, until being stopped
+% by \cs{@@_trig_large_auxiii:w} using \cs{exp_stop_f:}. Once all
+% these blocks are unpacked, the \cs{exp_stop_f:} and $0$ to $7$
+% digits are removed by \cs[no-index]{use_none:n\ldots{}n}.
+% Finally, \cs{@@_trig_large_auxii:w} packs $64$ digits (there are
+% between $65$ and $72$ at this point) into groups of~$4$ and the
+% \texttt{auxv} auxiliary is called.
+% \begin{macrocode}
+\cs_new:Npn \@@_trig_large:ww #1, #2#3#4#5#6;
+ {
+ \exp_after:wN \@@_trig_large_auxi:w
+ \int_value:w \@@_int_eval:w (#1 - 4) / 8 \exp_after:wN ,
+ \int_value:w #1 , ;
+ {#2}{#3}{#4}{#5} ;
+ }
+\cs_new:Npn \@@_trig_large_auxi:w #1, #2,
+ {
+ \exp_after:wN \exp_after:wN
+ \exp_after:wN \@@_trig_large_auxii:w
+ \cs:w
+ use_none:n \prg_replicate:nn { #2 - #1 * 8 } { n }
+ \exp_after:wN
+ \cs_end:
+ \int_value:w
+ \__kernel_intarray_item:Nn \c_@@_trig_intarray
+ { \@@_int_eval:w #1 + 1 \scan_stop: }
+ \exp_after:wN \@@_trig_large_auxiii:w \int_value:w
+ \__kernel_intarray_item:Nn \c_@@_trig_intarray
+ { \@@_int_eval:w #1 + 2 \scan_stop: }
+ \exp_after:wN \@@_trig_large_auxiii:w \int_value:w
+ \__kernel_intarray_item:Nn \c_@@_trig_intarray
+ { \@@_int_eval:w #1 + 3 \scan_stop: }
+ \exp_after:wN \@@_trig_large_auxiii:w \int_value:w
+ \__kernel_intarray_item:Nn \c_@@_trig_intarray
+ { \@@_int_eval:w #1 + 4 \scan_stop: }
+ \exp_after:wN \@@_trig_large_auxiii:w \int_value:w
+ \__kernel_intarray_item:Nn \c_@@_trig_intarray
+ { \@@_int_eval:w #1 + 5 \scan_stop: }
+ \exp_after:wN \@@_trig_large_auxiii:w \int_value:w
+ \__kernel_intarray_item:Nn \c_@@_trig_intarray
+ { \@@_int_eval:w #1 + 6 \scan_stop: }
+ \exp_after:wN \@@_trig_large_auxiii:w \int_value:w
+ \__kernel_intarray_item:Nn \c_@@_trig_intarray
+ { \@@_int_eval:w #1 + 7 \scan_stop: }
+ \exp_after:wN \@@_trig_large_auxiii:w \int_value:w
+ \__kernel_intarray_item:Nn \c_@@_trig_intarray
+ { \@@_int_eval:w #1 + 8 \scan_stop: }
+ \exp_after:wN \@@_trig_large_auxiii:w \int_value:w
+ \__kernel_intarray_item:Nn \c_@@_trig_intarray
+ { \@@_int_eval:w #1 + 9 \scan_stop: }
+ \exp_stop_f:
+ }
+\cs_new:Npn \@@_trig_large_auxii:w
+ {
+ \@@_pack_twice_four:wNNNNNNNN \@@_pack_twice_four:wNNNNNNNN
+ \@@_pack_twice_four:wNNNNNNNN \@@_pack_twice_four:wNNNNNNNN
+ \@@_pack_twice_four:wNNNNNNNN \@@_pack_twice_four:wNNNNNNNN
+ \@@_pack_twice_four:wNNNNNNNN \@@_pack_twice_four:wNNNNNNNN
+ \@@_trig_large_auxv:www ;
+ }
+\cs_new:Npn \@@_trig_large_auxiii:w 1 { \exp_stop_f: }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[rEXP]
+% {
+% \@@_trig_large_auxv:www,
+% \@@_trig_large_auxvi:wnnnnnnnn,
+% \@@_trig_large_pack:NNNNNw
+% }
+% First come the first $64$~digits of the fractional part of
+% $10^{\text{\texttt{\#1}}-16}/(2\pi)$, arranged in $16$~blocks
+% of~$4$, and ending with a semicolon. Then a few more digits of the
+% same fractional part, ending with a semicolon, then $4$~blocks of
+% $4$~digits holding the significand of the original argument.
+% Multiply the $16$-digit significand with the $64$-digit fractional
+% part: the \texttt{auxvi} auxiliary receives the significand
+% as~|#2#3#4#5| and $16$~digits of the fractional part as~|#6#7#8#9|,
+% and computes one step of the usual ladder of \texttt{pack} functions
+% we use for multiplication (see \emph{e.g.,} \cs{@@_fixed_mul:wwn}),
+% then discards one block of the fractional part to set things up for
+% the next step of the ladder. We perform $13$~such steps, replacing
+% the last \texttt{middle} shift by the appropriate \texttt{trailing}
+% shift, then discard the significand and remaining $3$~blocks from
+% the fractional part, as there are not enough digits to compute any
+% more step in the ladder. The last semicolon closes the ladder, and
+% we return control to the \texttt{auxvii} auxiliary.
+% \begin{macrocode}
+\cs_new:Npn \@@_trig_large_auxv:www #1; #2; #3;
+ {
+ \exp_after:wN \@@_use_i_until_s:nw
+ \exp_after:wN \@@_trig_large_auxvii:w
+ \int_value:w \@@_int_eval:w \c_@@_leading_shift_int
+ \prg_replicate:nn { 13 }
+ { \@@_trig_large_auxvi:wnnnnnnnn }
+ + \c_@@_trailing_shift_int - \c_@@_middle_shift_int
+ \@@_use_i_until_s:nw
+ ; #3 #1 ; ;
+ }
+\cs_new:Npn \@@_trig_large_auxvi:wnnnnnnnn #1; #2#3#4#5#6#7#8#9
+ {
+ \exp_after:wN \@@_trig_large_pack:NNNNNw
+ \int_value:w \@@_int_eval:w \c_@@_middle_shift_int
+ + #2*#9 + #3*#8 + #4*#7 + #5*#6
+ #1; {#2}{#3}{#4}{#5} {#7}{#8}{#9}
+ }
+\cs_new:Npn \@@_trig_large_pack:NNNNNw #1#2#3#4#5#6;
+ { + #1#2#3#4#5 ; #6 }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[rEXP]
+% {
+% \@@_trig_large_auxvii:w,
+% \@@_trig_large_auxviii:w,
+% }
+% \begin{macro}[EXP]
+% {
+% \@@_trig_large_auxix:Nw,
+% \@@_trig_large_auxx:wNNNNN,
+% \@@_trig_large_auxxi:w
+% }
+% The \texttt{auxvii} auxiliary is followed by $52$~digits and a
+% semicolon. We find the octant as the integer part of $8$~times what
+% follows, or equivalently as the integer part of $|#1#2#3|/125$, and
+% add it to the surrounding integer expression for the octant. We
+% then compute $8$~times the $52$-digit number, with a minus sign if
+% the octant is odd. Again, the last \texttt{middle} shift is
+% converted to a \texttt{trailing} shift. Any integer part (including
+% negative values which come up when the octant is odd) is discarded
+% by \cs{@@_use_i_until_s:nw}. The resulting fractional part should
+% then be converted to radians by multiplying by~$2\pi/8$, but first,
+% build an extended precision number by abusing
+% \cs{@@_ep_to_ep_loop:N} with the appropriate trailing markers.
+% Finally, \cs{@@_trig_small:ww} sets up the argument for the
+% functions which compute the Taylor series.
+% \begin{macrocode}
+\cs_new:Npn \@@_trig_large_auxvii:w #1#2#3
+ {
+ \exp_after:wN \@@_trig_large_auxviii:ww
+ \int_value:w \@@_int_eval:w (#1#2#3 - 62) / 125 ;
+ #1#2#3
+ }
+\cs_new:Npn \@@_trig_large_auxviii:ww #1;
+ {
+ + #1
+ \if_int_odd:w #1 \exp_stop_f:
+ \exp_after:wN \@@_trig_large_auxix:Nw
+ \exp_after:wN -
+ \else:
+ \exp_after:wN \@@_trig_large_auxix:Nw
+ \exp_after:wN +
+ \fi:
+ }
+\cs_new:Npn \@@_trig_large_auxix:Nw
+ {
+ \exp_after:wN \@@_use_i_until_s:nw
+ \exp_after:wN \@@_trig_large_auxxi:w
+ \int_value:w \@@_int_eval:w \c_@@_leading_shift_int
+ \prg_replicate:nn { 13 }
+ { \@@_trig_large_auxx:wNNNNN }
+ + \c_@@_trailing_shift_int - \c_@@_middle_shift_int
+ ;
+ }
+\cs_new:Npn \@@_trig_large_auxx:wNNNNN #1; #2 #3#4#5#6
+ {
+ \exp_after:wN \@@_trig_large_pack:NNNNNw
+ \int_value:w \@@_int_eval:w \c_@@_middle_shift_int
+ #2 8 * #3#4#5#6
+ #1; #2
+ }
+\cs_new:Npn \@@_trig_large_auxxi:w #1;
+ {
+ \exp_after:wN \@@_ep_mul_raw:wwwwN
+ \int_value:w \@@_int_eval:w 0 \@@_ep_to_ep_loop:N #1 ; ; !
+ 0,{7853}{9816}{3397}{4483}{0961}{5661};
+ \@@_trig_small:ww
+ }
+% \end{macrocode}
+% \end{macro}
+% \end{macro}
+%
+% \subsubsection{Computing the power series}
+%
+% \begin{macro}[EXP]
+% {\@@_sin_series_o:NNwwww, \@@_sin_series_aux_o:NNnwww}
+% Here we receive a conversion function \cs{@@_ep_to_float_o:wwN} or
+% \cs{@@_ep_inv_to_float_o:wwN}, a \meta{sign} ($0$ or~$2$), a
+% (non-negative) \meta{octant} delimited by a dot, a \meta{fixed
+% point} number delimited by a semicolon, and an extended-precision
+% number. The auxiliary receives:
+% \begin{itemize}
+% \item the conversion function~|#1|;
+% \item the final sign, which depends on the octant~|#3| and the
+% sign~|#2|;
+% \item the octant~|#3|, which controls the series we use;
+% \item the square |#4 * #4| of the argument as a fixed point number,
+% computed with \cs{@@_fixed_mul:wwn};
+% \item the number itself as an extended-precision number.
+% \end{itemize}
+% If the octant is in $\{1,2,5,6,\ldots{}\}$, we are near an extremum
+% of the function and we use the series
+% \[
+% \cos(x) = 1 - x^2 \bigg( \frac{1}{2!} - x^2 \bigg( \frac{1}{4!}
+% - x^2 \bigg( \cdots \bigg) \bigg) \bigg) .
+% \]
+% Otherwise, the series
+% \[
+% \sin(x) = x \bigg( 1 - x^2 \bigg( \frac{1}{3!} - x^2 \bigg(
+% \frac{1}{5!} - x^2 \bigg( \cdots \bigg) \bigg) \bigg) \bigg)
+% \]
+% is used. Finally, the extended-precision number is converted to a
+% floating point number with the given sign, and \cs{@@_sanitize:Nw}
+% checks for overflow and underflow.
+% \begin{macrocode}
+\cs_new:Npn \@@_sin_series_o:NNwwww #1#2#3. #4;
+ {
+ \@@_fixed_mul:wwn #4; #4;
+ {
+ \exp_after:wN \@@_sin_series_aux_o:NNnwww
+ \exp_after:wN #1
+ \int_value:w
+ \if_int_odd:w \@@_int_eval:w (#3 + 2) / 4 \@@_int_eval_end:
+ #2
+ \else:
+ \if_meaning:w #2 0 2 \else: 0 \fi:
+ \fi:
+ {#3}
+ }
+ }
+\cs_new:Npn \@@_sin_series_aux_o:NNnwww #1#2#3 #4; #5,#6;
+ {
+ \if_int_odd:w \@@_int_eval:w #3 / 2 \@@_int_eval_end:
+ \exp_after:wN \use_i:nn
+ \else:
+ \exp_after:wN \use_ii:nn
+ \fi:
+ { % 1/18!
+ \@@_fixed_mul_sub_back:wwwn {0000}{0000}{0000}{0001}{5619}{2070};
+ #4;{0000}{0000}{0000}{0477}{9477}{3324};
+ \@@_fixed_mul_sub_back:wwwn #4;{0000}{0000}{0011}{4707}{4559}{7730};
+ \@@_fixed_mul_sub_back:wwwn #4;{0000}{0000}{2087}{6756}{9878}{6810};
+ \@@_fixed_mul_sub_back:wwwn #4;{0000}{0027}{5573}{1922}{3985}{8907};
+ \@@_fixed_mul_sub_back:wwwn #4;{0000}{2480}{1587}{3015}{8730}{1587};
+ \@@_fixed_mul_sub_back:wwwn #4;{0013}{8888}{8888}{8888}{8888}{8889};
+ \@@_fixed_mul_sub_back:wwwn #4;{0416}{6666}{6666}{6666}{6666}{6667};
+ \@@_fixed_mul_sub_back:wwwn #4;{5000}{0000}{0000}{0000}{0000}{0000};
+ \@@_fixed_mul_sub_back:wwwn#4;{10000}{0000}{0000}{0000}{0000}{0000};
+ { \@@_fixed_continue:wn 0, }
+ }
+ { % 1/17!
+ \@@_fixed_mul_sub_back:wwwn {0000}{0000}{0000}{0028}{1145}{7254};
+ #4;{0000}{0000}{0000}{7647}{1637}{3182};
+ \@@_fixed_mul_sub_back:wwwn #4;{0000}{0000}{0160}{5904}{3836}{8216};
+ \@@_fixed_mul_sub_back:wwwn #4;{0000}{0002}{5052}{1083}{8544}{1719};
+ \@@_fixed_mul_sub_back:wwwn #4;{0000}{0275}{5731}{9223}{9858}{9065};
+ \@@_fixed_mul_sub_back:wwwn #4;{0001}{9841}{2698}{4126}{9841}{2698};
+ \@@_fixed_mul_sub_back:wwwn #4;{0083}{3333}{3333}{3333}{3333}{3333};
+ \@@_fixed_mul_sub_back:wwwn #4;{1666}{6666}{6666}{6666}{6666}{6667};
+ \@@_fixed_mul_sub_back:wwwn#4;{10000}{0000}{0000}{0000}{0000}{0000};
+ { \@@_ep_mul:wwwwn 0, } #5,#6;
+ }
+ {
+ \exp_after:wN \@@_sanitize:Nw
+ \exp_after:wN #2
+ \int_value:w \@@_int_eval:w #1
+ }
+ #2
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]
+% {\@@_tan_series_o:NNwwww, \@@_tan_series_aux_o:Nnwww}
+% Contrarily to \cs{@@_sin_series_o:NNwwww} which received a
+% conversion auxiliary as~|#1|, here, |#1| is $0$ for tangent
+% and $2$ for
+% cotangent. Consider first the case of the tangent. The octant |#3|
+% starts at $1$, which means that it is $1$ or $2$ for $\lvert
+% x\rvert\in[0,\pi/2]$, it is $3$ or $4$ for $\lvert
+% x\rvert\in[\pi/2,\pi]$, and so on: the intervals on which
+% $\tan\lvert x\rvert\geq 0$ coincide with those for which $\lfloor
+% (|#3| + 1) / 2\rfloor$ is odd. We also have to take into account
+% the original sign of $x$ to get the sign of the final result; it is
+% straightforward to check that the first \cs{int_value:w} expansion
+% produces $0$ for a positive final result, and $2$ otherwise. A
+% similar story holds for $\cot(x)$.
+%
+% The auxiliary receives the sign, the octant, the square of the
+% (reduced) input, and the (reduced) input (an extended-precision
+% number) as arguments. It then
+% computes the numerator and denominator of
+% \[
+% \tan(x) \simeq
+% \frac{x (1 - x^2 (a_1 - x^2 (a_2 - x^2 (a_3 - x^2 (a_4 - x^2 a_5)))))}
+% {1 - x^2 (b_1 - x^2 (b_2 - x^2 (b_3 - x^2 (b_4 - x^2 b_5))))} .
+% \]
+% The ratio is computed by \cs{@@_ep_div:wwwwn}, then converted to a
+% floating point number. For octants~|#3| (really, quadrants) next to
+% a pole of the
+% functions, the fixed point numerator and denominator are exchanged
+% before computing the ratio. Note that this \cs{if_int_odd:w} test
+% relies on the fact that the octant is at least~$1$.
+% \begin{macrocode}
+\cs_new:Npn \@@_tan_series_o:NNwwww #1#2#3. #4;
+ {
+ \@@_fixed_mul:wwn #4; #4;
+ {
+ \exp_after:wN \@@_tan_series_aux_o:Nnwww
+ \int_value:w
+ \if_int_odd:w \@@_int_eval:w #3 / 2 \@@_int_eval_end:
+ \exp_after:wN \reverse_if:N
+ \fi:
+ \if_meaning:w #1#2 2 \else: 0 \fi:
+ {#3}
+ }
+ }
+\cs_new:Npn \@@_tan_series_aux_o:Nnwww #1 #2 #3; #4,#5;
+ {
+ \@@_fixed_mul_sub_back:wwwn {0000}{0000}{1527}{3493}{0856}{7059};
+ #3; {0000}{0159}{6080}{0274}{5257}{6472};
+ \@@_fixed_mul_sub_back:wwwn #3; {0002}{4571}{2320}{0157}{2558}{8481};
+ \@@_fixed_mul_sub_back:wwwn #3; {0115}{5830}{7533}{5397}{3168}{2147};
+ \@@_fixed_mul_sub_back:wwwn #3; {1929}{8245}{6140}{3508}{7719}{2982};
+ \@@_fixed_mul_sub_back:wwwn #3;{10000}{0000}{0000}{0000}{0000}{0000};
+ { \@@_ep_mul:wwwwn 0, } #4,#5;
+ {
+ \@@_fixed_mul_sub_back:wwwn {0000}{0007}{0258}{0681}{9408}{4706};
+ #3;{0000}{2343}{7175}{1399}{6151}{7670};
+ \@@_fixed_mul_sub_back:wwwn #3;{0019}{2638}{4588}{9232}{8861}{3691};
+ \@@_fixed_mul_sub_back:wwwn #3;{0536}{6357}{0691}{4344}{6852}{4252};
+ \@@_fixed_mul_sub_back:wwwn #3;{5263}{1578}{9473}{6842}{1052}{6315};
+ \@@_fixed_mul_sub_back:wwwn#3;{10000}{0000}{0000}{0000}{0000}{0000};
+ {
+ \reverse_if:N \if_int_odd:w
+ \@@_int_eval:w (#2 - 1) / 2 \@@_int_eval_end:
+ \exp_after:wN \@@_reverse_args:Nww
+ \fi:
+ \@@_ep_div:wwwwn 0,
+ }
+ }
+ {
+ \exp_after:wN \@@_sanitize:Nw
+ \exp_after:wN #1
+ \int_value:w \@@_int_eval:w \@@_ep_to_float_o:wwN
+ }
+ #1
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \subsection{Inverse trigonometric functions}
+%
+% All inverse trigonometric functions (arcsine, arccosine, arctangent,
+% arccotangent, arccosecant, and arcsecant) are based on a function
+% often denoted \texttt{atan2}. This function is accessed directly by
+% feeding two arguments to arctangent, and is defined by \(\operatorname{atan}(y, x) =
+% \operatorname{atan}(y/x)\) for generic \(y\) and~\(x\). Its advantages over the
+% conventional arctangent is that it takes values in $[-\pi,\pi]$ rather
+% than $[-\pi/2,\pi/2]$, and that it is better behaved in boundary
+% cases. Other inverse trigonometric functions are expressed in terms
+% of \(\operatorname{atan}\) as
+% \begin{align}
+% \operatorname{acos} x & = \operatorname{atan}(\sqrt{1-x^2}, x) \\
+% \operatorname{asin} x & = \operatorname{atan}(x, \sqrt{1-x^2}) \\
+% \operatorname{asec} x & = \operatorname{atan}(\sqrt{x^2-1}, 1) \\
+% \operatorname{acsc} x & = \operatorname{atan}(1, \sqrt{x^2-1}) \\
+% \operatorname{atan} x & = \operatorname{atan}(x, 1) \\
+% \operatorname{acot} x & = \operatorname{atan}(1, x) .
+% \end{align}
+% Rather than introducing a new function, \texttt{atan2}, the arctangent
+% function \texttt{atan} is overloaded: it can take one or two
+% arguments. In the comments below, following many texts, we call the
+% first argument~$y$ and the second~$x$, because $\operatorname{atan}(y, x) = \operatorname{atan}(y
+% / x)$ is the angular coordinate of the point $(x, y)$.
+%
+% As for direct trigonometric functions, the first step in computing
+% $\operatorname{atan}(y, x)$ is argument reduction. The sign of~$y$ gives that
+% of the result. We distinguish eight regions where the point $(x,
+% \lvert y\rvert)$ can lie, of angular size roughly $\pi/8$,
+% characterized by their \enquote{octant}, between $0$ and~$7$ included. In
+% each region, we compute an arctangent as a Taylor series, then shift
+% this arctangent by the appropriate multiple of $\pi/4$ and sign to get
+% the result. Here is a list of octants, and how we compute the
+% arctangent (we assume $y>0$: otherwise replace $y$ by~$-y$ below):
+% \begin{itemize}
+% \item[0] $0 < \lvert y\rvert < 0.41421 x$, then
+% $\operatorname{atan}\frac{\lvert y\rvert}{x}$
+% is given by a nicely convergent Taylor series;
+% \item[1] $0 < 0.41421 x < \lvert y\rvert < x$, then
+% $\operatorname{atan}\frac{\lvert y\rvert}{x}
+% = \frac{\pi}{4}-\operatorname{atan}\frac{x-\lvert y\rvert}{x+\lvert y\rvert}$;
+% \item[2] $0 < 0.41421 \lvert y\rvert < x < \lvert y\rvert$, then
+% $\operatorname{atan}\frac{\lvert y\rvert}{x}
+% = \frac{\pi}{4}+\operatorname{atan}\frac{-x+\lvert y\rvert}{x+\lvert y\rvert}$;
+% \item[3] $0 < x < 0.41421 \lvert y\rvert$, then
+% $\operatorname{atan}\frac{\lvert y\rvert}{x}
+% = \frac{\pi}{2}-\operatorname{atan}\frac{x}{\lvert y\rvert}$;
+% \item[4] $0 < -x < 0.41421 \lvert y\rvert$, then
+% $\operatorname{atan}\frac{\lvert y\rvert}{x}
+% = \frac{\pi}{2}+\operatorname{atan}\frac{-x}{\lvert y\rvert}$;
+% \item[5] $0 < 0.41421 \lvert y\rvert < -x < \lvert y\rvert$, then
+% $\operatorname{atan}\frac{\lvert y\rvert}{x}
+% =\frac{3\pi}{4}-\operatorname{atan}\frac{x+\lvert y\rvert}{-x+\lvert y\rvert}$;
+% \item[6] $0 < -0.41421 x < \lvert y\rvert < -x$, then
+% $\operatorname{atan}\frac{\lvert y\rvert}{x}
+% =\frac{3\pi}{4}+\operatorname{atan}\frac{-x-\lvert y\rvert}{-x+\lvert y\rvert}$;
+% \item[7] $0 < \lvert y\rvert < -0.41421 x$, then
+% $\operatorname{atan}\frac{\lvert y\rvert}{x}
+% = \pi-\operatorname{atan}\frac{\lvert y\rvert}{-x}$.
+% \end{itemize}
+% In the following, we denote by~$z$ the ratio among
+% $\lvert\frac{y}{x}\rvert$, $\lvert\frac{x}{y}\rvert$,
+% $\lvert\frac{x+y}{x-y}\rvert$, $\lvert\frac{x-y}{x+y}\rvert$ which
+% appears in the right-hand side above.
+%
+% \subsubsection{Arctangent and arccotangent}
+%
+% \begin{macro}[EXP]{\@@_atan_o:Nw, \@@_acot_o:Nw, \@@_atan_default:w}
+% The parsing step manipulates \texttt{atan} and \texttt{acot} like
+% \texttt{min} and \texttt{max}, reading in an array of operands, but
+% also leaves \cs{use_i:nn} or \cs{use_ii:nn} depending on whether the
+% result should be given in radians or in degrees. The helper
+% \cs{@@_parse_function_one_two:nnw} checks that the operand is one or
+% two floating point numbers (not tuples) and leaves its second
+% argument or its tail accordingly (its first argument is used for
+% error messages). More precisely if we are given a single floating
+% point number \cs{@@_atan_default:w} places \cs{c_one_fp} (expanded)
+% after it; otherwise \cs{@@_atan_default:w} is omitted by
+% \cs{@@_parse_function_one_two:nnw}.
+% \begin{macrocode}
+\cs_new:Npn \@@_atan_o:Nw #1
+ {
+ \@@_parse_function_one_two:nnw
+ { #1 { atan } { atand } }
+ { \@@_atan_default:w \@@_atanii_o:Nww #1 }
+ }
+\cs_new:Npn \@@_acot_o:Nw #1
+ {
+ \@@_parse_function_one_two:nnw
+ { #1 { acot } { acotd } }
+ { \@@_atan_default:w \@@_acotii_o:Nww #1 }
+ }
+\cs_new:Npx \@@_atan_default:w #1#2#3 @ { #1 #2 #3 \c_one_fp @ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_atanii_o:Nww, \@@_acotii_o:Nww}
+% If either operand is \texttt{nan}, we return it. If both are
+% normal, we call \cs{@@_atan_normal_o:NNnwNnw}. If both are zero or
+% both infinity, we call \cs{@@_atan_inf_o:NNNw} with argument~$2$,
+% leading to a result among $\{\pm\pi/4, \pm 3\pi/4\}$ (in degrees,
+% $\{\pm 45, \pm 135\}$). Otherwise, one is much bigger than the
+% other, and we call \cs{@@_atan_inf_o:NNNw} with either an argument
+% of~$4$, leading to the values $\pm\pi/2$ (in degrees,~$\pm 90$),
+% or~$0$, leading to $\{\pm 0, \pm\pi\}$ (in degrees, $\{\pm 0,\pm
+% 180\}$). Since $\operatorname{acot}(x, y) = \operatorname{atan}(y, x)$,
+% \cs{@@_acotii_o:ww} simply reverses its two arguments.
+% \begin{macrocode}
+\cs_new:Npn \@@_atanii_o:Nww
+ #1 \s_@@ \@@_chk:w #2#3#4; \s_@@ \@@_chk:w #5 #6 @
+ {
+ \if_meaning:w 3 #2 \@@_case_return_i_o:ww \fi:
+ \if_meaning:w 3 #5 \@@_case_return_ii_o:ww \fi:
+ \if_case:w
+ \if_meaning:w #2 #5
+ \if_meaning:w 1 #2 10 \else: 0 \fi:
+ \else:
+ \if_int_compare:w #2 > #5 \exp_stop_f: 1 \else: 2 \fi:
+ \fi:
+ \exp_stop_f:
+ \@@_case_return:nw { \@@_atan_inf_o:NNNw #1 #3 2 }
+ \or: \@@_case_return:nw { \@@_atan_inf_o:NNNw #1 #3 4 }
+ \or: \@@_case_return:nw { \@@_atan_inf_o:NNNw #1 #3 0 }
+ \fi:
+ \@@_atan_normal_o:NNnwNnw #1
+ \s_@@ \@@_chk:w #2#3#4;
+ \s_@@ \@@_chk:w #5 #6
+ }
+\cs_new:Npn \@@_acotii_o:Nww #1#2; #3;
+ { \@@_atanii_o:Nww #1#3; #2; }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_atan_inf_o:NNNw}
+% This auxiliary is called whenever one number is $\pm 0$ or
+% $\pm\infty$ (and neither is \nan{}). Then the result only depends
+% on the signs, and its value is a multiple of $\pi/4$. We use the
+% same auxiliary as for normal numbers,
+% \cs{@@_atan_combine_o:NwwwwwN}, with arguments the final sign~|#2|;
+% the octant~|#3|; $\operatorname{atan} z/z=1$ as a fixed point number; $z=0$~as a
+% fixed point number; and $z=0$~as an extended-precision number.
+% Given the values we provide, $\operatorname{atan} z$ is computed to be~$0$,
+% and the result is $[|#3|/2]\cdot\pi/4$ if the sign~|#5| of~$x$
+% is positive, and $[(7-|#3|)/2]\cdot\pi/4$ for negative~$x$, where
+% the divisions are rounded up.
+% \begin{macrocode}
+\cs_new:Npn \@@_atan_inf_o:NNNw #1#2#3 \s_@@ \@@_chk:w #4#5#6;
+ {
+ \exp_after:wN \@@_atan_combine_o:NwwwwwN
+ \exp_after:wN #2
+ \int_value:w \@@_int_eval:w
+ \if_meaning:w 2 #5 7 - \fi: #3 \exp_after:wN ;
+ \c_@@_one_fixed_tl
+ {0000}{0000}{0000}{0000}{0000}{0000};
+ 0,{0000}{0000}{0000}{0000}{0000}{0000}; #1
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_atan_normal_o:NNnwNnw}
+% Here we simply reorder the floating point data into a pair of signed
+% extended-precision numbers, that is, a sign, an exponent ending with
+% a comma, and a six-block mantissa ending with a semi-colon. This
+% extended precision is required by other inverse trigonometric
+% functions, to compute things like $\operatorname{atan}(x,\sqrt{1-x^2})$ without
+% intermediate rounding errors.
+% \begin{macrocode}
+\cs_new_protected:Npn \@@_atan_normal_o:NNnwNnw
+ #1 \s_@@ \@@_chk:w 1#2#3#4; \s_@@ \@@_chk:w 1#5#6#7;
+ {
+ \@@_atan_test_o:NwwNwwN
+ #2 #3, #4{0000}{0000};
+ #5 #6, #7{0000}{0000}; #1
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_atan_test_o:NwwNwwN}
+% This receives: the sign~|#1| of~$y$, its exponent~|#2|, its $24$
+% digits~|#3| in groups of~$4$, and similarly for~$x$. We prepare to
+% call \cs{@@_atan_combine_o:NwwwwwN} which expects the sign~|#1|, the
+% octant, the ratio $(\operatorname{atan} z)/z = 1 - \cdots$, and the value of~$z$,
+% both as a fixed point number and as an extended-precision floating
+% point number with a mantissa in $[0.01,1)$. For now, we place |#1|
+% as a first argument, and start an integer expression for the octant.
+% The sign of $x$ does not affect~$z$, so we simply leave
+% a contribution to the octant: $\meta{octant} \to 7 - \meta{octant}$
+% for negative~$x$. Then we order $\lvert y\rvert$ and $\lvert
+% x\rvert$ in a non-decreasing order: if $\lvert y\rvert > \lvert
+% x\rvert$, insert $3-$ in the expression for the octant, and swap the
+% two numbers. The finer test with $0.41421$ is done by
+% \cs{@@_atan_div:wnwwnw} after the operands have been ordered.
+% \begin{macrocode}
+\cs_new:Npn \@@_atan_test_o:NwwNwwN #1#2,#3; #4#5,#6;
+ {
+ \exp_after:wN \@@_atan_combine_o:NwwwwwN
+ \exp_after:wN #1
+ \int_value:w \@@_int_eval:w
+ \if_meaning:w 2 #4
+ 7 - \@@_int_eval:w
+ \fi:
+ \if_int_compare:w
+ \@@_ep_compare:wwww #2,#3; #5,#6; > 0 \exp_stop_f:
+ 3 -
+ \exp_after:wN \@@_reverse_args:Nww
+ \fi:
+ \@@_atan_div:wnwwnw #2,#3; #5,#6;
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[rEXP]{\@@_atan_div:wnwwnw, \@@_atan_near:wwwn}
+% \begin{macro}[EXP]{\@@_atan_near_aux:wwn}
+% This receives two positive numbers $a$ and~$b$ (equal to $\lvert
+% x\rvert$ and~$\lvert y\rvert$ in some order), each as an exponent
+% and $6$~blocks of $4$~digits, such that $0<a<b$. If $0.41421b<a$,
+% the two numbers are \enquote{near}, hence the point $(y,x)$ that we
+% started with is closer to the diagonals $\{\lvert y\rvert = \lvert
+% x\rvert\}$ than to the axes $\{xy = 0\}$. In that case, the octant
+% is~$1$ (possibly combined with the $7-$ and $3-$ inserted earlier)
+% and we wish to compute $\operatorname{atan}\frac{b-a}{a+b}$. Otherwise, the
+% octant is~$0$ (again, combined with earlier terms) and we wish to
+% compute $\operatorname{atan}\frac{a}{b}$. In any case, call \cs{@@_atan_auxi:ww}
+% followed by~$z$, as a comma-delimited exponent and a fixed point
+% number.
+% \begin{macrocode}
+\cs_new:Npn \@@_atan_div:wnwwnw #1,#2#3; #4,#5#6;
+ {
+ \if_int_compare:w
+ \@@_int_eval:w 41421 * #5 < #2 000
+ \if_case:w \@@_int_eval:w #4 - #1 \@@_int_eval_end:
+ 00 \or: 0 \fi:
+ \exp_stop_f:
+ \exp_after:wN \@@_atan_near:wwwn
+ \fi:
+ 0
+ \@@_ep_div:wwwwn #1,{#2}#3; #4,{#5}#6;
+ \@@_atan_auxi:ww
+ }
+\cs_new:Npn \@@_atan_near:wwwn
+ 0 \@@_ep_div:wwwwn #1,#2; #3,
+ {
+ 1
+ \@@_ep_to_fixed:wwn #1 - #3, #2;
+ \@@_atan_near_aux:wwn
+ }
+\cs_new:Npn \@@_atan_near_aux:wwn #1; #2;
+ {
+ \@@_fixed_add:wwn #1; #2;
+ { \@@_fixed_sub:wwn #2; #1; { \@@_ep_div:wwwwn 0, } 0, }
+ }
+% \end{macrocode}
+% \end{macro}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_atan_auxi:ww, \@@_atan_auxii:w}
+% Convert~$z$ from a representation as an exponent and a fixed point
+% number in $[0.01,1)$ to a fixed point number only, then set up the
+% call to \cs{@@_atan_Taylor_loop:www}, followed by the fixed point
+% representation of~$z$ and the old representation.
+% \begin{macrocode}
+\cs_new:Npn \@@_atan_auxi:ww #1,#2;
+ { \@@_ep_to_fixed:wwn #1,#2; \@@_atan_auxii:w #1,#2; }
+\cs_new:Npn \@@_atan_auxii:w #1;
+ {
+ \@@_fixed_mul:wwn #1; #1;
+ {
+ \@@_atan_Taylor_loop:www 39 ;
+ {0000}{0000}{0000}{0000}{0000}{0000} ;
+ }
+ ! #1;
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]
+% {\@@_atan_Taylor_loop:www, \@@_atan_Taylor_break:w}
+% We compute the series of $(\operatorname{atan} z)/z$. A typical intermediate
+% stage has $|#1|=2k-1$, $|#2| =
+% \frac{1}{2k+1}-z^2(\frac{1}{2k+3}-z^2(\cdots-z^2\frac{1}{39}))$, and
+% $|#3|=z^2$. To go to the next step $k\to k-1$, we compute
+% $\frac{1}{2k-1}$, then subtract from it $z^2$ times |#2|. The loop
+% stops when $k=0$: then |#2| is $(\operatorname{atan} z)/z$, and there is a need to
+% clean up all the unnecessary data, end the integer expression
+% computing the octant with a semicolon, and leave the result~|#2|
+% afterwards.
+% \begin{macrocode}
+\cs_new:Npn \@@_atan_Taylor_loop:www #1; #2; #3;
+ {
+ \if_int_compare:w #1 = -1 \exp_stop_f:
+ \@@_atan_Taylor_break:w
+ \fi:
+ \exp_after:wN \@@_fixed_div_int:wwN \c_@@_one_fixed_tl #1;
+ \@@_rrot:www \@@_fixed_mul_sub_back:wwwn #2; #3;
+ {
+ \exp_after:wN \@@_atan_Taylor_loop:www
+ \int_value:w \@@_int_eval:w #1 - 2 ;
+ }
+ #3;
+ }
+\cs_new:Npn \@@_atan_Taylor_break:w
+ \fi: #1 \@@_fixed_mul_sub_back:wwwn #2; #3 !
+ { \fi: ; #2 ; }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]
+% {\@@_atan_combine_o:NwwwwwN, \@@_atan_combine_aux:ww}
+% This receives a \meta{sign}, an \meta{octant}, a fixed point value
+% of $(\operatorname{atan} z)/z$, a fixed point number~$z$, and another
+% representation of~$z$, as an \meta{exponent} and the fixed point
+% number $10^{-\meta{exponent}} z$, followed by either \cs{use_i:nn}
+% (when working in radians) or \cs{use_ii:nn} (when working in
+% degrees). The function computes the floating point result
+% \begin{equation}
+% \meta{sign} \left(
+% \left\lceil\frac{\meta{octant}}{2}\right\rceil
+% \frac{\pi}{4}
+% + (-1)^{\meta{octant}} \frac{\operatorname{atan} z}{z} \cdot z\right) \,,
+% \end{equation}
+% multiplied by $180/\pi$ if working in degrees, and using in any case
+% the most appropriate representation of~$z$. The floating point
+% result is passed to \cs{@@_sanitize:Nw}, which checks for overflow
+% or underflow. If the octant is~$0$, leave the exponent~|#5| for
+% \cs{@@_sanitize:Nw}, and multiply $|#3|=\frac{\operatorname{atan} z}{z}$
+% with~|#6|, the adjusted~$z$. Otherwise, multiply $|#3|=\frac{\operatorname{atan}
+% z}{z}$ with $|#4|=z$, then compute the appropriate multiple of
+% $\frac{\pi}{4}$ and add or subtract the product $|#3|\cdot|#4|$. In
+% both cases, convert to a floating point with
+% \cs{@@_fixed_to_float_o:wN}.
+% \begin{macrocode}
+\cs_new:Npn \@@_atan_combine_o:NwwwwwN #1 #2; #3; #4; #5,#6; #7
+ {
+ \exp_after:wN \@@_sanitize:Nw
+ \exp_after:wN #1
+ \int_value:w \@@_int_eval:w
+ \if_meaning:w 0 #2
+ \exp_after:wN \use_i:nn
+ \else:
+ \exp_after:wN \use_ii:nn
+ \fi:
+ { #5 \@@_fixed_mul:wwn #3; #6; }
+ {
+ \@@_fixed_mul:wwn #3; #4;
+ {
+ \exp_after:wN \@@_atan_combine_aux:ww
+ \int_value:w \@@_int_eval:w #2 / 2 ; #2;
+ }
+ }
+ { #7 \@@_fixed_to_float_o:wN \@@_fixed_to_float_rad_o:wN }
+ #1
+ }
+\cs_new:Npn \@@_atan_combine_aux:ww #1; #2;
+ {
+ \@@_fixed_mul_short:wwn
+ {7853}{9816}{3397}{4483}{0961}{5661};
+ {#1}{0000}{0000};
+ {
+ \if_int_odd:w #2 \exp_stop_f:
+ \exp_after:wN \@@_fixed_sub:wwn
+ \else:
+ \exp_after:wN \@@_fixed_add:wwn
+ \fi:
+ }
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \subsubsection{Arcsine and arccosine}
+%
+% \begin{macro}[EXP]{\@@_asin_o:w}
+% Again, the first argument provided by \pkg{l3fp-parse} is
+% \cs{use_i:nn} if we are to work in radians and \cs{use_ii:nn} for
+% degrees. Then comes a floating point number. The arcsine of $\pm
+% 0$ or \nan{} is the same floating point number. The arcsine of
+% $\pm\infty$ raises an invalid operation exception. Otherwise, call
+% an auxiliary common with \cs{@@_acos_o:w}, feeding it information
+% about what function is being performed (for \enquote{invalid operation}
+% exceptions).
+% \begin{macrocode}
+\cs_new:Npn \@@_asin_o:w #1 \s_@@ \@@_chk:w #2#3; @
+ {
+ \if_case:w #2 \exp_stop_f:
+ \@@_case_return_same_o:w
+ \or:
+ \@@_case_use:nw
+ { \@@_asin_normal_o:NfwNnnnnw #1 { #1 { asin } { asind } } }
+ \or:
+ \@@_case_use:nw
+ { \@@_invalid_operation_o:fw { #1 { asin } { asind } } }
+ \else:
+ \@@_case_return_same_o:w
+ \fi:
+ \s_@@ \@@_chk:w #2 #3;
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_acos_o:w}
+% The arccosine of $\pm 0$ is $\pi / 2$ (in degrees,~$90$). The
+% arccosine of $\pm\infty$ raises an invalid operation exception. The
+% arccosine of \nan{} is itself. Otherwise, call an auxiliary common
+% with \cs{@@_sin_o:w}, informing it that it was called by
+% \texttt{acos} or \texttt{acosd}, and preparing to swap some
+% arguments down the line.
+% \begin{macrocode}
+\cs_new:Npn \@@_acos_o:w #1 \s_@@ \@@_chk:w #2#3; @
+ {
+ \if_case:w #2 \exp_stop_f:
+ \@@_case_use:nw { \@@_atan_inf_o:NNNw #1 0 4 }
+ \or:
+ \@@_case_use:nw
+ {
+ \@@_asin_normal_o:NfwNnnnnw #1 { #1 { acos } { acosd } }
+ \@@_reverse_args:Nww
+ }
+ \or:
+ \@@_case_use:nw
+ { \@@_invalid_operation_o:fw { #1 { acos } { acosd } } }
+ \else:
+ \@@_case_return_same_o:w
+ \fi:
+ \s_@@ \@@_chk:w #2 #3;
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_asin_normal_o:NfwNnnnnw}
+% If the exponent~|#5| is at most $0$, the operand lies
+% within $(-1,1)$ and the operation is permitted: call
+% \cs{@@_asin_auxi_o:NnNww} with the appropriate arguments. If the
+% number is exactly~$\pm 1$ (the test works because we know that
+% $|#5|\geq 1$, $|#6#7|\geq 10000000$, $|#8#9|\geq 0$, with equality
+% only for $\pm 1$), we also call \cs{@@_asin_auxi_o:NnNww}.
+% Otherwise, \cs{@@_use_i:ww} gets rid of the \texttt{asin} auxiliary,
+% and raises instead an invalid operation, because the operand is
+% outside the domain of arcsine or arccosine.
+% \begin{macrocode}
+\cs_new:Npn \@@_asin_normal_o:NfwNnnnnw
+ #1#2#3 \s_@@ \@@_chk:w 1#4#5#6#7#8#9;
+ {
+ \if_int_compare:w #5 < 1 \exp_stop_f:
+ \exp_after:wN \@@_use_none_until_s:w
+ \fi:
+ \if_int_compare:w \@@_int_eval:w #5 + #6#7 + #8#9 = 1000 0001 ~
+ \exp_after:wN \@@_use_none_until_s:w
+ \fi:
+ \@@_use_i:ww
+ \@@_invalid_operation_o:fw {#2}
+ \s_@@ \@@_chk:w 1#4{#5}{#6}{#7}{#8}{#9};
+ \@@_asin_auxi_o:NnNww
+ #1 {#3} #4 #5,{#6}{#7}{#8}{#9}{0000}{0000};
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_asin_auxi_o:NnNww, \@@_asin_isqrt:wn}
+% We compute $x/\sqrt{1-x^2}$. This function is used by \texttt{asin}
+% and \texttt{acos}, but also by \texttt{acsc} and \texttt{asec} after
+% inverting the operand, thus it must manipulate extended-precision
+% numbers. First evaluate $1-x^2$ as $(1+x)(1-x)$: this behaves
+% better near~$x=1$. We do the addition/subtraction with fixed point
+% numbers (they are not implemented for extended-precision floats),
+% but go back to extended-precision floats to multiply and compute the
+% inverse square root $1/\sqrt{1-x^2}$. Finally, multiply by the
+% (positive) extended-precision float $\lvert x\rvert$, and feed the
+% (signed) result, and the number~$+1$, as arguments to the arctangent
+% function. When computing the arccosine, the arguments
+% $x/\sqrt{1-x^2}$ and~$+1$ are swapped by~|#2|
+% (\cs{@@_reverse_args:Nww} in that case) before
+% \cs{@@_atan_test_o:NwwNwwN} is evaluated. Note that the arctangent
+% function requires normalized arguments, hence the need for
+% \texttt{ep_to_ep} and \texttt{continue} after \texttt{ep_mul}.
+% \begin{macrocode}
+\cs_new:Npn \@@_asin_auxi_o:NnNww #1#2#3#4,#5;
+ {
+ \@@_ep_to_fixed:wwn #4,#5;
+ \@@_asin_isqrt:wn
+ \@@_ep_mul:wwwwn #4,#5;
+ \@@_ep_to_ep:wwN
+ \@@_fixed_continue:wn
+ { #2 \@@_atan_test_o:NwwNwwN #3 }
+ 0 1,{1000}{0000}{0000}{0000}{0000}{0000}; #1
+ }
+\cs_new:Npn \@@_asin_isqrt:wn #1;
+ {
+ \exp_after:wN \@@_fixed_sub:wwn \c_@@_one_fixed_tl #1;
+ {
+ \@@_fixed_add_one:wN #1;
+ \@@_fixed_continue:wn { \@@_ep_mul:wwwwn 0, } 0,
+ }
+ \@@_ep_isqrt:wwn
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \subsubsection{Arccosecant and arcsecant}
+%
+% \begin{macro}[EXP]{\@@_acsc_o:w}
+% Cases are mostly labelled by~|#2|, except when |#2| is~$2$: then we
+% use |#3#2|, which is $02=2$ when the number is $+\infty$ and
+% $22$~when the number is $-\infty$. The arccosecant of $\pm 0$
+% raises an invalid operation exception. The arccosecant of
+% $\pm\infty$ is $\pm 0$ with the same sign. The arcosecant of \nan{}
+% is itself. Otherwise, \cs{@@_acsc_normal_o:NfwNnw} does some more
+% tests, keeping the function name (\texttt{acsc} or \texttt{acscd})
+% as an argument for invalid operation exceptions.
+% \begin{macrocode}
+\cs_new:Npn \@@_acsc_o:w #1 \s_@@ \@@_chk:w #2#3#4; @
+ {
+ \if_case:w \if_meaning:w 2 #2 #3 \fi: #2 \exp_stop_f:
+ \@@_case_use:nw
+ { \@@_invalid_operation_o:fw { #1 { acsc } { acscd } } }
+ \or: \@@_case_use:nw
+ { \@@_acsc_normal_o:NfwNnw #1 { #1 { acsc } { acscd } } }
+ \or: \@@_case_return_o:Nw \c_zero_fp
+ \or: \@@_case_return_same_o:w
+ \else: \@@_case_return_o:Nw \c_minus_zero_fp
+ \fi:
+ \s_@@ \@@_chk:w #2 #3 #4;
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_asec_o:w}
+% The arcsecant of $\pm 0$ raises an invalid operation exception. The
+% arcsecant of $\pm\infty$ is $\pi / 2$ (in degrees,~$90$). The
+% arcosecant of \nan{} is itself. Otherwise, do some more tests,
+% keeping the function name \texttt{asec} (or \texttt{asecd}) as an
+% argument for invalid operation exceptions, and a
+% \cs{@@_reverse_args:Nww} following precisely that appearing in
+% \cs{@@_acos_o:w}.
+% \begin{macrocode}
+\cs_new:Npn \@@_asec_o:w #1 \s_@@ \@@_chk:w #2#3; @
+ {
+ \if_case:w #2 \exp_stop_f:
+ \@@_case_use:nw
+ { \@@_invalid_operation_o:fw { #1 { asec } { asecd } } }
+ \or:
+ \@@_case_use:nw
+ {
+ \@@_acsc_normal_o:NfwNnw #1 { #1 { asec } { asecd } }
+ \@@_reverse_args:Nww
+ }
+ \or: \@@_case_use:nw { \@@_atan_inf_o:NNNw #1 0 4 }
+ \else: \@@_case_return_same_o:w
+ \fi:
+ \s_@@ \@@_chk:w #2 #3;
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_acsc_normal_o:NfwNnw}
+% If the exponent is non-positive, the operand is less than~$1$ in
+% absolute value, which is always an invalid operation: complain.
+% Otherwise, compute the inverse of the operand, and feed it to
+% \cs{@@_asin_auxi_o:NnNww} (with all the appropriate arguments). This
+% computes what we want thanks to
+% $\operatorname{acsc}(x)=\operatorname{asin}(1/x)$ and
+% $\operatorname{asec}(x)=\operatorname{acos}(1/x)$.
+% \begin{macrocode}
+\cs_new:Npn \@@_acsc_normal_o:NfwNnw #1#2#3 \s_@@ \@@_chk:w 1#4#5#6;
+ {
+ \int_compare:nNnTF {#5} < 1
+ {
+ \@@_invalid_operation_o:fw {#2}
+ \s_@@ \@@_chk:w 1#4{#5}#6;
+ }
+ {
+ \@@_ep_div:wwwwn
+ 1,{1000}{0000}{0000}{0000}{0000}{0000};
+ #5,#6{0000}{0000};
+ { \@@_asin_auxi_o:NnNww #1 {#3} #4 }
+ }
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macrocode}
+%</initex|package>
+% \end{macrocode}
+%
+% \end{implementation}
+%
+% \PrintChanges
+%
+% \PrintIndex