Quadratic function interpolator (Mathematica)

quad.wl

quadraticInterpolation[f_, vars__] := Module[{
  variables =  List[vars],
  rules,
  constant,
  quadratic,
  final,
  cached,
  mdcache,
  calculated = <||>,
  linear 
},

  (* не знаю, как лучше сделать то, чтобы считать во вторых производных xy и yx один раз вместо двух*)
  (* ну вот тут сортирует индексы и направляет во вторую функцию *)
  cached[expr_, {i_,j_}] := mdcache[expr, Sort[{i,j}]];
  SetAttributes[cached, HoldFirst];

  (* тут казалось бы можно тупо mdcahe:= mdcahe[] = ..., но я побоялся что expr если большой, то долго будет искать шаблон *)
  (* поэтому поиск по ассоциациям, но по сути он там нафиг не нужен *)
  mdcache[expr_, {i_,j_}] := If[KeyExistsQ[calculated, {i,j}],
    calculated[{i,j}]
  ,
    calculated[{i,j}] = expr
  ];

  SetAttributes[mdcache, HoldFirst];

  (* типа в правила где все переменные присвоены некому начальному веткору {x0, x1, x2...} *)
  rules = (#[[1]] -> #[[2]]) &/@ variables;
  
  (* примеяем правила и получаем первый член разложения *)
  constant = Hold[f] /. rules // ReleaseHold;

  (* линейный член разложения *)
  linear = Table[
    With[{
      r = {i[[1]] -> 1.}
    },
      (Hold[f] /. r /. rules // ReleaseHold) - constant
    ]
    
  , {i, variables}];

  (* квадратичный член разложения - вторые производные *)
  (* спизжено https://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h10/kompendiet/kap14.pdf  стр.358 *)
  quadratic = Table[
    cached[With[{
      f11 = With[{
        r = {i[[1]] -> 1., j[[1]] -> 1.}
      },
        (Hold[f] /. r /. rules // ReleaseHold)
      ],

      fm1m1 = With[{
        r = {i[[1]] -> -1., j[[1]] -> -1.}
      },
        (Hold[f] /. r /. rules // ReleaseHold)
      ],

      f1m1 = With[{
        r = {i[[1]] -> 1., j[[1]] -> -1.}
      },
        (Hold[f] /. r /. rules // ReleaseHold)
      ],

      fm11 = With[{
        r = {i[[1]] -> -1., j[[1]] -> 1.}
      },
        (Hold[f] /. r /. rules // ReleaseHold)
      ]      
      
    },

      Print["Checking... ", i, j];

      (*FB[*)((f11 - f1m1 - fm11 + fm1m1)(*,*)/(*,*)(4.))(*]FB*)
      
    ], {i,j}]
    
  , {i, variables}, {j, variables}];

  (* делаем основу для чистой функции-замены *)
  final = Sum[
  
    With[{l = linear[[i]]},
      (* нахуй тут холд? Я надеюсь, что WL будет думать что это packed array, если вдруг у нас функция возвращает тензоры *)
      (* иначе этот Слот проникнет внутрь массивов и будет жопа *)
      Slot[i] Hold[l]
    ]
    
  , {i, Length[linear]}] + Sum[
  
    With[{
      q = (*FB[*)((1)(*,*)/(*,*)(2))(*]FB*)quadratic[[i, j]]
    },
      Slot[i]Slot[j] Hold[q]
    ]
    
  , {i, Length[linear]}, {j, Length[linear]}] + With[{c = constant}, Hold[c] ];

  With[{
    r = final
  },
    (* превращаем в чистую функцию и сносим холды *)
    (r &) /. {Hold -> Identity}
  ]
]

SetAttributes[quadraticInterpolation, HoldFirst]

(* пример с какой-нибудь функцией нетривиальное в том смысле что возвращает листы *)

f = quadraticInterpolation[{x^2, y^2} + {x y, z} + {1, 3}, {x,0}, {y,0}, {z,0}];
f[1,2,3]