%% broydensolve.sty
%% Copyright 2025 Matthias Floré
%
% This work may be distributed and/or modified under the
% conditions of the LaTeX Project Public License, either version 1.3c
% of this license or (at your option) any later version.
% The latest version of this license is in
%   http://www.latex-project.org/lppl.txt
% and version 1.3c or later is part of all distributions of LaTeX
% version 2005/12/01 or later.
%
% This work has the LPPL maintenance status `maintained'.
% 
% The Current Maintainer of this work is Matthias Floré.
%
% This work consists of the files broydensolve-doc.pdf, broydensolve.sty,
% broydensolve-doc.tex and README.md.
\NeedsTeXFormat{LaTeX2e}
\ProvidesExplPackage{broydensolve}{2025/07/20}{1.0}{Solve a system of equations with Broyden's good method}

%%> \subsection{Variables}

\bool_new:N \l__broydensolve_do_bool

\clist_new:N \g__broydensolve_roots_clist

\fp_new:N \l__broydensolve_aux_fp
\fp_new:N \g__broydensolve_coord_x_fp
\fp_new:N \g__broydensolve_coord_y_fp
\fp_new:N \l__broydensolve_det_fp
\fp_new:N \l__broydensolve_norm_h_fp
\fp_new:N \l__broydensolve_norm_x_fp

\int_new:N \g__broydensolve_iterations_int
\int_new:N \l__broydensolve_N_int

\seq_new:N \l__broydensolve_coordinates_seq
\seq_new:N \l__broydensolve_coords_seq
\seq_new:N \l__broydensolve_func_seq
\seq_new:N \l__broydensolve_init_seq
\seq_new:N \l__broydensolve_var_seq
\seq_new:N \l__broydensolve_var_N_seq

\tl_new:N \l__broydensolve_func_tl
\tl_new:N \l__broydensolve_name_tl

%%> \subsection{Keys}

\keys_define:nn { broydensolve }
  {
    abs-approx-error .fp_set:N = \l__broydensolve_abs_approx_error_fp ,
    abs-approx-error .initial:n = 10^-3 ,
    coordinates .bool_set:N = \l__broydensolve_coordinates_bool ,
    func .tl_set:N = \l__broydensolve_funcs_tl ,
    func-error .fp_set:N = \l__broydensolve_func_error_fp ,
    func-error .initial:n = 10^-3 ,
    init .code:n = \exp_args:NNe \seq_set_from_clist:Nn \l__broydensolve_init_seq {#1} ,
    iterations .int_set:N = \l__broydensolve_iterations_int ,
    iterations .initial:n = 5 ,
    rel-approx-error .fp_set:N = \l__broydensolve_rel_approx_error_fp ,
    rel-approx-error .initial:n = 10^-3 ,
    stop-crit .str_set_e:N = \l__broydensolve_stop_crit_str ,
    stop-crit .initial:n = rel-approx-error ,
    var .code:n = \exp_args:NNe \seq_set_from_clist:Nn \l__broydensolve_var_seq {#1} ,
    var .initial:n = { x , y , z }
  }

%%> \subsection{Functions}

\cs_new:Npn \__broydensolve_add:n #1
  { ( #1 < 0 ? ( #1 + 2 * pi ) : ( #1 ) ) }

\cs_new_protected:Npn \__broydensolve_add_name:N #1
  {
    \tl_if_empty:nF {#1}
      {
        \cs_if_exist:cTF { __fp_parse_word_#1:N }
          { \tl_build_put_right:Nn \l__broydensolve_func_tl {#1} }
          {
            \bool_if:NTF \l__broydensolve_coordinates_bool
              {
                \tl_set:Ne \l__broydensolve_name_tl { \tl_range:nnn {#1} { 1 } { -2 } }
                \seq_if_in:NVTF \l__broydensolve_coordinates_seq \l__broydensolve_name_tl
                  { \tl_build_put_right:Nn \l__broydensolve_func_tl { \cs:w l__broydensolve_var_#1_fp \cs_end: } }
                  {
                    \seq_if_in:NVF \l__broydensolve_coords_seq \l__broydensolve_name_tl
                      {
                        \path let \p { l__broydensolve_coord } = ( \l__broydensolve_name_tl ) in
                          [
                            / utils / exec =
                              {
                                \fp_gset:Nn \g__broydensolve_coord_x_fp
                                  {
                                    ( \pgf@yy * \x { l__broydensolve_coord } - \pgf@yx * \y { l__broydensolve_coord } )
                                    / \l__broydensolve_det_fp
                                  }
                                \fp_gset:Nn \g__broydensolve_coord_y_fp
                                  {
                                    ( \pgf@xx * \y { l__broydensolve_coord } - \pgf@xy * \x { l__broydensolve_coord } )
                                    / \l__broydensolve_det_fp
                                  }
                              }
                          ] ;
                        \fp_zero_new:c { l__broydensolve_coord_\l__broydensolve_name_tl x_fp }
                        \fp_zero_new:c { l__broydensolve_coord_\l__broydensolve_name_tl y_fp }
                        \fp_set_eq:cN { l__broydensolve_coord_\l__broydensolve_name_tl x_fp } \g__broydensolve_coord_x_fp
                        \fp_set_eq:cN { l__broydensolve_coord_\l__broydensolve_name_tl y_fp } \g__broydensolve_coord_y_fp
                        \seq_put_right:NV \l__broydensolve_coords_seq \l__broydensolve_name_tl
                      }
                    \tl_build_put_right:Ne \l__broydensolve_func_tl { \fp_use:c { l__broydensolve_coord_#1_fp } }
                  }
              }
              { \tl_build_put_right:Nn \l__broydensolve_func_tl { \cs:w l__broydensolve_var_#1_fp \cs_end: } }
          }
      }
  }

\cs_new:Npn \__broydensolve_atan:nnn #1#2#3
  { atan ( \clist_item:nn {#1} {#2} y - \clist_item:nn {#1} {#3} y , \clist_item:nn {#1} {#2} x - \clist_item:nn {#1} {#3} x ) }

\cs_new_protected:Npn \__broydensolve_init:nn #1#2
  {
    \fp_gzero_new:c { g__broydensolve_x_0_#2_fp }
    \fp_gset:cn { g__broydensolve_x_0_#2_fp } {#1}
    \fp_zero_new:c { l__broydensolve_var_#2_fp }
    \fp_set:cn { l__broydensolve_var_#2_fp } {#1}
  }

\cs_new_protected:Npn \__broydensolve_stop_crit:
  {
    \fp_set:Nn \l__broydensolve_aux_fp { abs ( \seq_use:Nn \l__broydensolve_func_seq { ) + abs ( } ) }
    \fp_compare:nNnTF { \l__broydensolve_aux_fp } > { 0 }
      {
        \str_case:VnF \l__broydensolve_stop_crit_str
          {
            { abs-approx-error }
              {
                \int_if_zero:nF { \g__broydensolve_iterations_int }
                  {
                    \fp_zero:N \l__broydensolve_norm_h_fp
                    \int_step_inline:nn { \l__broydensolve_N_int }
                      { \fp_add:Nn \l__broydensolve_norm_h_fp { abs ( \cs:w l__broydensolve_h_##1_fp \cs_end: ) } }
                    \bool_set:Nn \l__broydensolve_do_bool
                      { \fp_compare_p:nNn { \l__broydensolve_norm_h_fp } < { \l__broydensolve_abs_approx_error_fp } }
                  }
              }
            { func-error }
              {
                \bool_set:Nn \l__broydensolve_do_bool
                  { \fp_compare_p:nNn { \l__broydensolve_aux_fp } < { \l__broydensolve_func_error_fp } }
              }
            { iterations }
              {
                \bool_set:Nn \l__broydensolve_do_bool
                  { \int_compare_p:nNn { \g__broydensolve_iterations_int } = { \l__broydensolve_iterations_int } }
              }
            { rel-approx-error }
              {
                \int_if_zero:nF { \g__broydensolve_iterations_int }
                  {
                    \fp_zero:N \l__broydensolve_norm_h_fp
                    \fp_zero:N \l__broydensolve_norm_x_fp
                    \seq_map_indexed_inline:Nn \l__broydensolve_var_N_seq
                      {
                        \fp_add:Nn \l__broydensolve_norm_h_fp { abs ( \cs:w l__broydensolve_h_##1_fp \cs_end: ) }
                        \fp_add:Nn \l__broydensolve_norm_x_fp { abs ( \cs:w l__broydensolve_var_##2_fp \cs_end: ) }
                      }
                    \bool_set:Nn \l__broydensolve_do_bool
                      {
                        \fp_compare_p:nNn
                          { \l__broydensolve_norm_h_fp } < { \l__broydensolve_rel_approx_error_fp * \l__broydensolve_norm_x_fp }
                      }
                  }
              }
          }
          { \PackageError { broydensolve } { Wrong~value~for~key~stop-crit } {} }
      }
      { \bool_set_true:N \l__broydensolve_do_bool }
  }

%%> \subsection{Document commands}

\NewExpandableDocumentCommand \BroydenIterations {}
  { \int_use:N \g__broydensolve_iterations_int }

\NewExpandableDocumentCommand \BroydenRoot { O { \g__broydensolve_iterations_int } m }
  { \fp_use:c { g__broydensolve_x_\int_eval:n {#1}_#2_fp } }

\NewExpandableDocumentCommand \BroydenRoots {}
  { \g__broydensolve_roots_clist }

\NewDocumentCommand \BroydenSetup { m }
  { \keys_set:nn { broydensolve } {#1} }

\NewDocumentCommand \BroydenSolve { m }
  {
    \group_begin:
      \keys_set:nn { broydensolve } {#1}
      \DeclareExpandableDocumentCommand \ang { r () }
        {
          \__broydensolve_add:n
            {
              \int_case:nnF { \clist_count:n {##1} }
                {
                  { 1 } { atan ( \clist_item:nn {##1} { 1 } y , \clist_item:nn {##1} { 1 } x ) }
                  { 2 } { \__broydensolve_atan:nnn {##1} { 2 } { 1 } }
                  { 3 } { \__broydensolve_atan:nnn {##1} { 3 } { 2 } - \__broydensolve_atan:nnn {##1} { 1 } { 2 } }
                  { 4 } { \__broydensolve_atan:nnn {##1} { 4 } { 3 } - \__broydensolve_atan:nnn {##1} { 2 } { 1 } }
                }
                { \ang {} }
            }
        }
      \DeclareExpandableDocumentCommand \col { r () }
        {
          \int_compare:nNnTF { \clist_count:n {##1} } = { 3 }
            {
              \clist_item:nn {##1} { 1 } x * ( \clist_item:nn {##1} { 2 } y - \clist_item:nn {##1} { 3 } y )
              + \clist_item:nn {##1} { 2 } x * ( \clist_item:nn {##1} { 3 } y - \clist_item:nn {##1} { 1 } y )
              + \clist_item:nn {##1} { 3 } x * ( \clist_item:nn {##1} { 1 } y - \clist_item:nn {##1} { 2 } y )
            }
            { \col {} }
        }
      \DeclareExpandableDocumentCommand \dis { r () }
        {
          \int_case:nnF { \clist_count:n {##1} }
            {
              { 1 } { sqrt ( \clist_item:nn {##1} { 1 } x ^ 2 + \clist_item:nn {##1} { 1 } y ^ 2 ) }
              { 2 }
                {
                  sqrt
                    (
                      ( \clist_item:nn {##1} { 2 } x - \clist_item:nn {##1} { 1 } x ) ^ 2
                      + ( \clist_item:nn {##1} { 2 } y - \clist_item:nn {##1} { 1 } y ) ^ 2
                    )
                }
            }
            { \dis {} }
        }
      \int_set:Nn \l__broydensolve_N_int { \clist_count:e { \l__broydensolve_funcs_tl } }
      \bool_if:NTF \l__broydensolve_coordinates_bool
        {
          \seq_map_indexed_inline:Nn \l__broydensolve_var_seq
            {
              \seq_put_right:Nn \l__broydensolve_var_N_seq { ##2 x }
              \seq_put_right:Nn \l__broydensolve_var_N_seq { ##2 y }
              \seq_put_right:Nn \l__broydensolve_coordinates_seq {##2}
              \int_compare:nNnT {##1} = { \l__broydensolve_N_int / 2 }
                { \seq_map_break: }
            }
          \fp_set:Nn \l__broydensolve_det_fp { \pgf@yy * \pgf@xx - \pgf@yx * \pgf@xy }
        }
        {
          \seq_map_indexed_inline:Nn \l__broydensolve_var_seq
            {
              \seq_put_right:Nn \l__broydensolve_var_N_seq {##2}
              \int_compare:nNnT {##1} = { \l__broydensolve_N_int }
                { \seq_map_break: }
            }
        }
      \exp_args:Ne \clist_map_inline:nn { \l__broydensolve_funcs_tl }
        {
          \tl_build_begin:N \l__broydensolve_func_tl
          \tl_build_begin:N \l__broydensolve_name_tl
          \tl_map_inline:nn {##1}
            {
              \bool_lazy_and:nnTF
                { \bool_lazy_or_p:nn { \int_compare_p:nNn { `####1 } < { `a } } { \int_compare_p:nNn { `z } < { `####1 } } }
                { \bool_lazy_or_p:nn { \int_compare_p:nNn { `####1 } < { `A } } { \int_compare_p:nNn { `Z } < { `####1 } } }
                {
                  \tl_build_end:N \l__broydensolve_name_tl
                  \exp_args:NV \__broydensolve_add_name:N \l__broydensolve_name_tl
                  \tl_build_begin:N \l__broydensolve_name_tl
                  \tl_build_put_right:Nn \l__broydensolve_func_tl {####1}
                }
                { \tl_build_put_right:Nn \l__broydensolve_name_tl {####1} }
            }
          \tl_build_end:N \l__broydensolve_name_tl
          \exp_args:NV \__broydensolve_add_name:N \l__broydensolve_name_tl
          \tl_build_end:N \l__broydensolve_func_tl
          \seq_put_right:NV \l__broydensolve_func_seq \l__broydensolve_func_tl
        }
      \seq_map_pairwise_function:NNN \l__broydensolve_init_seq \l__broydensolve_var_N_seq \__broydensolve_init:nn
      \int_step_inline:nn { \l__broydensolve_N_int }
        {
          \int_step_inline:nn { \l__broydensolve_N_int }
            { \fp_zero_new:c { l__broydensolve_B_##1_####1_fp } }
          %set B to the identity matrix
          \fp_set:cn { l__broydensolve_B_##1_##1_fp } { 1 }
          \fp_zero_new:c { l__broydensolve_w_##1_fp }
          \fp_zero_new:c { l__broydensolve_h_##1_fp }
          \fp_zero_new:c { l__broydensolve_aux_##1_fp }
          \fp_zero_new:c { l__broydensolve_auxi_##1_fp }
        }
      \int_gzero:N \g__broydensolve_iterations_int
      \__broydensolve_stop_crit:
      \bool_until_do:Nn \l__broydensolve_do_bool
        {
          %set w=-f(x)
          \seq_map_indexed_inline:Nn \l__broydensolve_func_seq
            { \fp_set:cn { l__broydensolve_w_##1_fp } { - (##2) } }
          %set h=-B*f(x)=B*w
          \int_step_inline:nn { \l__broydensolve_N_int }
            {
              \fp_zero:c { l__broydensolve_h_##1_fp }
              \int_step_inline:nn { \l__broydensolve_N_int }
                {
                  \fp_add:cn { l__broydensolve_h_##1_fp }
                    { \cs:w l__broydensolve_B_##1_####1_fp \cs_end: * \cs:w l__broydensolve_w_####1_fp \cs_end: }
                }
            }
          %set the variable(s) to x+h
          \seq_map_indexed_inline:Nn \l__broydensolve_var_N_seq
            { \fp_add:cn { l__broydensolve_var_##2_fp } { \cs:w l__broydensolve_h_##1_fp \cs_end: } }
          %add f(x+h) to w so that w=f(x+h)-f(x)
          \seq_map_indexed_inline:Nn \l__broydensolve_func_seq
            { \fp_add:cn { l__broydensolve_w_##1_fp } {##2} }
          %update B
          \int_step_inline:nn { \l__broydensolve_N_int }
            {
              \fp_zero:c { l__broydensolve_aux_##1_fp }
              \int_step_inline:nn { \l__broydensolve_N_int }
                {
                  \fp_add:cn { l__broydensolve_aux_##1_fp }
                    { \cs:w l__broydensolve_B_##1_####1_fp \cs_end: * \cs:w l__broydensolve_w_####1_fp \cs_end: }
                }
            }
          \fp_zero:N \l__broydensolve_aux_fp
          \int_step_inline:nn { \l__broydensolve_N_int }
            {
              \fp_add:Nn \l__broydensolve_aux_fp { \cs:w l__broydensolve_h_##1_fp \cs_end: * \cs:w l__broydensolve_aux_##1_fp \cs_end: }
              \fp_sub:cn { l__broydensolve_aux_##1_fp } { \cs:w l__broydensolve_h_##1_fp \cs_end: }
            }
          \int_step_inline:nn { \l__broydensolve_N_int }
            {
              \fp_zero:c { l__broydensolve_auxi_##1_fp }
              \int_step_inline:nn { \l__broydensolve_N_int }
                {
                  \fp_add:cn { l__broydensolve_auxi_##1_fp }
                    { \cs:w l__broydensolve_h_####1_fp \cs_end: * \cs:w l__broydensolve_B_####1_##1_fp \cs_end: }
                }
            }
          \int_step_inline:nn { \l__broydensolve_N_int }
            {
              \int_step_inline:nn { \l__broydensolve_N_int }
                {
                  \fp_sub:cn { l__broydensolve_B_##1_####1_fp }
                    {
                      \cs:w l__broydensolve_aux_##1_fp \cs_end: * \cs:w l__broydensolve_auxi_####1_fp \cs_end:
                      / \l__broydensolve_aux_fp
                    }
                }
            }
          \int_gincr:N \g__broydensolve_iterations_int
          \seq_map_inline:Nn \l__broydensolve_var_N_seq
            {
              \fp_gzero_new:c { g__broydensolve_x_\int_use:N \g__broydensolve_iterations_int _##1_fp }
              \fp_gset_eq:cc { g__broydensolve_x_\int_use:N \g__broydensolve_iterations_int _##1_fp } { l__broydensolve_var_##1_fp }
            }
          \__broydensolve_stop_crit:
        }
      \clist_gclear:N \g__broydensolve_roots_clist
      \seq_map_inline:Nn \l__broydensolve_var_N_seq
        { \clist_gput_right:Ne \g__broydensolve_roots_clist { \BroydenRoot {##1} } }
      \bool_if:NT \l__broydensolve_coordinates_bool
        {
          \seq_map_inline:Nn \l__broydensolve_coordinates_seq
            { \coordinate (##1) at ( \BroydenRoot { ##1 x } , \BroydenRoot { ##1 y } ) ; }
        }
    \group_end:
  }

\endinput