%% 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/24}{2.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

You are viewing proxied material from ftp.icm.edu.pl. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.