Fully documented lambda functions for powering Excel. Includes functions for Arrays, Matrices, Polynomials, Interpolation, Calculus, and Linear Algebra.

Vertex42 Lambda Library

/*  VERTEX42 LAMBDA LIBRARY
The functions included here are fully documented at the following site:
https://www.vertex42.com/lambda/

If importing this Gist into a separate module via the Excel Labs add-in, 
the name of the module (such as "VLL") is added to the beginning of the function
names in the Name Manager. For example, RESCALE would become VLL.RESCALE

Here are a few examples of module names you could use: L, VLL, LIB, LL 
*/

/*
MIT License

Copyright (c) Vertex42 LLC, https://www.vertex42.com/

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/

/* *******************
* GENERAL FUNCTIONS
*********************/

/**
* Round a value to a number of significant figures
* round_opt: 0=ROUND, -1=ROUNDDOWN, 1=ROUNDUP
* SFROUND(1234500,4,-1) = 1234000
*/
SFROUND = LAMBDA(value,sig_figs,[round_opt],
LET(doc,"https://www.vertex42.com/lambda/sfround.html",
    round_opt,IF(ISOMITTED(round_opt),0,round_opt),
    exponent,INT(LOG10(ABS(value))),
    roundto,sig_figs-(1+exponent),
    IF(value=0,0,IF(sig_figs<1,"Error: sig_figs<1",
        SWITCH(round_opt,
            0,ROUND(value,roundto),
            1,ROUNDUP(value,roundto),
            -1,ROUNDDOWN(value,roundto),
        )
    ))
));
/**
* Return a running total (cumulative sum) of values in an array
* for each column separately.
*/
CUMULATIVESUM = LAMBDA(array,
LET(doc,"https://www.vertex42.com/lambda/cumulativesum.html",
    MAKEARRAY(ROWS(array),COLUMNS(array),LAMBDA(r,c,
        SUM(INDEX(array,SEQUENCE(r),c))
    ))
));
/**
* Create an S-Curve from Project Task Dates and Values
*/
SCURVE = LAMBDA(startdates,enddates,values,[is_daily],[weekend],[holidays],[percent],
LET(doc,"https://www.vertex42.com/lambda/scurve.html",
    weekend,IF(ISOMITTED(weekend),"0000000",weekend),
    holidays,IF(ISOMITTED(holidays),{0},holidays),
    is_daily,IF(ISOMITTED(is_daily),FALSE,is_daily),
    percent,IF(ISOMITTED(percent),FALSE,percent),
    table,HSTACK(startdates,enddates,values),
    ftable,FILTER(table,
        NOT(ISTEXT(startdates))*(startdates>0)*NOT(ISTEXT(enddates))*(enddates>0)*(enddates>=startdates)
    ),
    start,CHOOSECOLS(ftable,1),
    end,CHOOSECOLS(ftable,2),
    days,BYROW(HSTACK(start,end),LAMBDA(row,
        NETWORKDAYS.INTL(CHOOSECOLS(row,1),CHOOSECOLS(row,2),weekend,holidays)
    )),
    val,CHOOSECOLS(ftable,3),
    daily,IF(is_daily,val,val/days),
    dates,SEQUENCE(MAX(end)-MIN(start)+1,1,MIN(start),1),
    date_total,BYROW(dates,LAMBDA(date,
        SUMPRODUCT((NETWORKDAYS.INTL(date,date,weekend,holidays)=1)*(start<=date)*(date<=end)*(daily))
    )),
    cumulative,SCAN(0,date_total,LAMBDA(acc,f,acc+f)),
    HSTACK(dates,date_total,IF(percent=TRUE,cumulative/MAX(cumulative),cumulative))
));
/**
* Used within HYPERLINK to create a hyperlink to a cell in the workbook
*/
BOOKMARK = LAMBDA(cell,
LET(doc,"https://www.vertex42.com/lambda/bookmark.html",
    "#"&ADDRESS(ROW(cell),COLUMN(cell),4,1,TEXTAFTER(CELL("filename",cell),"]"))
));

/* ******************* 
* SEQUENCES, NUMBERING, GRIDS
* ********************/

/**
* Create a sequence of integers or characters between start and end
*/
SE = LAMBDA(start,end,
LET(doc,"https://www.vertex42.com/lambda/se.html",
    IF(ISNUMBER(start),
        SEQUENCE(ABS(end-start)+1,1,start,SIGN(end-start)),
        UNICHAR(SEQUENCE(ABS(UNICODE(end)-UNICODE(start))+1,1,UNICODE(start),SIGN(UNICODE(end)-UNICODE(start))))
    )
));
/**
* Create a vector on interval [start,end] with n linearly spaced points
* LINSPACE(2,3,5) = {2; 2.25; 2.5; 2.75; 3}
*/
LINSPACE = LAMBDA(start,end,n,
LET(doc,"https://www.vertex42.com/lambda/linspace.html",
    SEQUENCE(n,1,start,(end-start)/(n-1))
));

/**
* Create a vector on interval [10^start,10^end] with n logarithmically spaced points
* LOGSPACE(1,2,5) = {10; 17.8; 31.6; 56.2; 100}
*/
LOGSPACE =LAMBDA(start,end,n,
LET(doc,"https://www.vertex42.com/lambda/logspace.html",
    10^SEQUENCE(n,1,start,(end-start)/(n-1))
));
/**
* Returns an array with values rescaled to [lower,upper]
* RESCALE({1,2,3},0,1) = {0, 0.5, 1}
*/
RESCALE = LAMBDA(array,lower,upper,
LET(doc,"https://www.vertex42.com/lambda/rescale.html",
    min,MIN(array),max,MAX(array),
    lower+(upper-lower)*(array-min)/(max-min)
));
/**
* Returns {TRUE;step} if the vector is uniformly spaced (all the
* steps are identical to a precision of n decimal places)
* If FALSE, returns {FALSE, MAX(steps)-MIN(steps)}
*/
ISUNIFORM = LAMBDA(vector,[precision_n],
LET(doc,"https://www.vertex42.com/lambda/isuniform.html",
    rows,ROWS(vector),cols,COLUMNS(vector),
    vector,TAKE(IF(AND(cols>1,rows=1),TRANSPOSE(vector),vector),,1),
    steps,DROP(vector,1)-DROP(vector,-1),
    unique_steps,UNIQUE(IF(ISOMITTED(precision_n),steps,ROUND(steps,precision_n))),
    IF(ROWS(unique_steps)=1,
        VSTACK(TRUE,unique_steps),
        VSTACK(FALSE,MAX(unique_steps)-MIN(unique_steps))
    )
));
/**
* Creates a set of 2D grid coordinates with x- and y-coordinates defined by 
* vectors x and y. Returns the coordinates as two columns of (x,y) pairs.
*/
MESHGRID = LAMBDA(xvec,yvec,[return_option],
LET(doc,"https://www.vertex42.com/lambda/meshgrid.html",
    xvec,IF(AND(ROWS(xvec)=1,COLUMNS(xvec)>1),TRANSPOSE(xvec),xvec),
    yvec,IF(AND(ROWS(yvec)=1,COLUMNS(yvec)>1),TRANSPOSE(yvec),yvec),
    xn, ROWS(xvec),
    yn, ROWS(yvec),
    xM, MAKEARRAY(yn,xn,LAMBDA(i,j,INDEX(TRANSPOSE(xvec),1,j))),
    yM, MAKEARRAY(yn,xn,LAMBDA(i,j,INDEX(yvec,i,1))),
    IF(return_option="X",xM,
        IF(return_option="Y",yM,
            HSTACK(TOCOL(xM),TOCOL(yM))
        )
    )
));
/**
* Return the first N values in the Fibonacci Sequence
*/
FIBONACCI = LAMBDA(n,[sequence_tf],
IF(OR(ROWS(n)>1,COLUMNS(n)>1,n<0,n<>INT(n)),"Error: n should be an integer > 0",
LET(doc,"https://www.vertex42.com/lambda/fibonacci.html",
    sequence_tf,IF(ISBLANK(sequence_tf),TRUE,sequence_tf),
    IF(sequence_tf=TRUE,
        IF(n=1,0,IF(n=2,{0;1},
            REDUCE({0;1},SEQUENCE(n-2,1,3,1),LAMBDA(acc,i,
                VSTACK(acc,INDEX(acc,i-2)+INDEX(acc,i-1))
            ))
        )),
        LET(
            phi, (1+SQRT(5))/2,
            psi, (1-SQRT(5))/2,
            (phi^n-psi^n)/(phi-psi)
        )
    )
)));

/**
* Check if a number or array of numbers is prime.
*/
/*
* Inputs:
*   n: A single number, vector, or array of numbers to check for primality.
*
* Outputs:
*   Returns a TRUE/FALSE array of the same dimensions as n, where TRUE 
*   indicates the corresponding number is prime.
*
* Notes:
*   - Handles non-numeric inputs, negative numbers, and non-integer values.
*   - For numbers < 1000, uses a precomputed list of primes for faster evaluation.
*   - For numbers ≥ 1000, uses divisors up to SQRT(n) to determine primality.
*   - n limited by size of arrays in Excel. Returns #NUM! when n is too large.
*/
ISPRIME = LAMBDA(n,
LET(doc, "https://w=ww.vertex42.com/lambda/isprime.html",
    version,"1/16/2025 - Updated to account for max n",
    // Step 1: List of primes below 1000 for fast checks.
    primes,{2;3;5;7;11;13;17;19;23;29;31;37;41;43;47;53;59;61;67;71;73;79;83;89;
    97;101;103;107;109;113;127;131;137;139;149;151;157;163;167;173;179;181;191;
    193;197;199;211;223;227;229;233;239;241;251;257;263;269;271;277;281;283;293;
    307;311;313;317;331;337;347;349;353;359;367;373;379;383;389;397;401;409;419;
    421;431;433;439;443;449;457;461;463;467;479;487;491;499;503;509;521;523;541;
    547;557;563;569;571;577;587;593;599;601;607;613;617;619;631;641;643;647;653;
    659;661;673;677;683;691;701;709;719;727;733;739;743;751;757;761;769;773;787;
    797;809;811;821;823;827;829;839;853;857;859;863;877;881;883;887;907;911;919;
    929;937;941;947;953;967;971;977;983;991;997},
    // Step 2: Evaluate each number in the input array.
    MAKEARRAY(ROWS(n),COLUMNS(n),LAMBDA(i,j,
        LET(v,INDEX(n,i,j),
            // Step 3: Check to exclude obvious non-prime candidates.
            IF( OR( NOT(ISNUMBER(v)), v<2, INT(v)<>v, AND(v<>2,ISEVEN(v)) ),
                FALSE,
                IF( v < 1000,
                    // Step 4: If v < 1000, check against list            
                    ISNUMBER(XMATCH(v,primes,0)),
                    // Step 5: Otherwise, test divisors
                    LET(
                        test_divisors, VSTACK(2,3,5,7,11,13,17,19,23,29,
                            TOCOL(SEQUENCE(INT((SQRT(v))/30))*30 + 
                            {1,7,11,13,17,19,23,29})
                        ),
                        MIN(MOD(v/test_divisors,1))>0
                    )
                )
            )
        )
    ))
));

/**
* Primes(n) returns the first n prime numbers less than 10000.
* Primes(,100,200) returns all prime numbers within [100,200].
* Primes(n,10,50) returns the first n primes within [10,50] sorted in ascending order.
* Primes(-n,10,50) returns the last n primes within [10,50] sorted in descending order.
*/
/*
* Inputs:
*   [n]: Optional. Number of primes to return. Defaults to all primes within the range [gte, lte].
*   [gte]: Optional. Lower bound of the range. Defaults to 0.
*   [lte]: Optional. Upper bound of the range. Defaults to 10000.
*
* Outputs:
*   Returns an array of prime numbers based on the input parameters:
*     - If n is omitted, returns all primes in [gte, lte].
*     - If n > 0, returns the first n primes in ascending order.
*     - If n < 0, returns the last n primes in descending order.
*
* Notes:
*   - Relies on the variable `the_array_of_primes_lt_10000`, which should contain all primes < 10000.
*   - Uses the FILTER function to narrow results to the specified range.
*/
Primes = LAMBDA([n], [gte], [lte],
LET(doc, "https://www.vertex42.com/lambda/isprime.html",
    version,"1.0.0 1/13/2025 - Original version",
    // Step 1: Handle default values for optional parameters.
    lte, IF(ISOMITTED(lte), 10000, lte),
    gte, IF(ISOMITTED(gte), 0, gte),
    // Step 2: Define the array of primes and handle n defaults.
    array,{2;3;5;7;11;13;17;19;23;29;31;37;41;43;47;53;59;61;67;71;73;79;83;89;97;101;103;107;109;113;127;131;137;139;149;151;157;163;167;173;179;181;191;193;197;199;211;223;227;229;233;239;241;251;257;263;269;271;277;281;283;293;307;311;313;317;331;337;347;349;353;359;367;373;379;383;389;397;401;409;419;421;431;433;439;443;449;457;461;463;467;479;487;491;499;503;509;521;523;541;547;557;563;569;571;577;587;593;599;601;607;613;617;619;631;641;643;647;653;659;661;673;677;683;691;701;709;719;727;733;739;743;751;757;761;769;773;787;797;809;811;821;823;827;829;839;853;857;859;863;877;881;883;887;907;911;919;929;937;941;947;953;967;971;977;983;991;997;1009;1013;1019;1021;1031;1033;1039;1049;1051;1061;1063;1069;1087;1091;1093;1097;1103;1109;1117;1123;1129;1151;1153;1163;1171;1181;1187;1193;1201;1213;1217;1223;1229;1231;1237;1249;1259;1277;1279;1283;1289;1291;1297;1301;1303;1307;1319;1321;1327;1361;1367;1373;1381;1399;1409;1423;1427;1429;1433;1439;1447;1451;1453;1459;1471;1481;1483;1487;1489;1493;1499;1511;1523;1531;1543;1549;1553;1559;1567;1571;1579;1583;1597;1601;1607;1609;1613;1619;1621;1627;1637;1657;1663;1667;1669;1693;1697;1699;1709;1721;1723;1733;1741;1747;1753;1759;1777;1783;1787;1789;1801;1811;1823;1831;1847;1861;1867;1871;1873;1877;1879;1889;1901;1907;1913;1931;1933;1949;1951;1973;1979;1987;1993;1997;1999;2003;2011;2017;2027;2029;2039;2053;2063;2069;2081;2083;2087;2089;2099;2111;2113;2129;2131;2137;2141;2143;2153;2161;2179;2203;2207;2213;2221;2237;2239;2243;2251;2267;2269;2273;2281;2287;2293;2297;2309;2311;2333;2339;2341;2347;2351;2357;2371;2377;2381;2383;2389;2393;2399;2411;2417;2423;2437;2441;2447;2459;2467;2473;2477;2503;2521;2531;2539;2543;2549;2551;2557;2579;2591;2593;2609;2617;2621;2633;2647;2657;2659;2663;2671;2677;2683;2687;2689;2693;2699;2707;2711;2713;2719;2729;2731;2741;2749;2753;2767;2777;2789;2791;2797;2801;2803;2819;2833;2837;2843;2851;2857;2861;2879;2887;2897;2903;2909;2917;2927;2939;2953;2957;2963;2969;2971;2999;3001;3011;3019;3023;3037;3041;3049;3061;3067;3079;3083;3089;3109;3119;3121;3137;3163;3167;3169;3181;3187;3191;3203;3209;3217;3221;3229;3251;3253;3257;3259;3271;3299;3301;3307;3313;3319;3323;3329;3331;3343;3347;3359;3361;3371;3373;3389;3391;3407;3413;3433;3449;3457;3461;3463;3467;3469;3491;3499;3511;3517;3527;3529;3533;3539;3541;3547;3557;3559;3571;3581;3583;3593;3607;3613;3617;3623;3631;3637;3643;3659;3671;3673;3677;3691;3697;3701;3709;3719;3727;3733;3739;3761;3767;3769;3779;3793;3797;3803;3821;3823;3833;3847;3851;3853;3863;3877;3881;3889;3907;3911;3917;3919;3923;3929;3931;3943;3947;3967;3989;4001;4003;4007;4013;4019;4021;4027;4049;4051;4057;4073;4079;4091;4093;4099;4111;4127;4129;4133;4139;4153;4157;4159;4177;4201;4211;4217;4219;4229;4231;4241;4243;4253;4259;4261;4271;4273;4283;4289;4297;4327;4337;4339;4349;4357;4363;4373;4391;4397;4409;4421;4423;4441;4447;4451;4457;4463;4481;4483;4493;4507;4513;4517;4519;4523;4547;4549;4561;4567;4583;4591;4597;4603;4621;4637;4639;4643;4649;4651;4657;4663;4673;4679;4691;4703;4721;4723;4729;4733;4751;4759;4783;4787;4789;4793;4799;4801;4813;4817;4831;4861;4871;4877;4889;4903;4909;4919;4931;4933;4937;4943;4951;4957;4967;4969;4973;4987;4993;4999;5003;5009;5011;5021;5023;5039;5051;5059;5077;5081;5087;5099;5101;5107;5113;5119;5147;5153;5167;5171;5179;5189;5197;5209;5227;5231;5233;5237;5261;5273;5279;5281;5297;5303;5309;5323;5333;5347;5351;5381;5387;5393;5399;5407;5413;5417;5419;5431;5437;5441;5443;5449;5471;5477;5479;5483;5501;5503;5507;5519;5521;5527;5531;5557;5563;5569;5573;5581;5591;5623;5639;5641;5647;5651;5653;5657;5659;5669;5683;5689;5693;5701;5711;5717;5737;5741;5743;5749;5779;5783;5791;5801;5807;5813;5821;5827;5839;5843;5849;5851;5857;5861;5867;5869;5879;5881;5897;5903;5923;5927;5939;5953;5981;5987;6007;6011;6029;6037;6043;6047;6053;6067;6073;6079;6089;6091;6101;6113;6121;6131;6133;6143;6151;6163;6173;6197;6199;6203;6211;6217;6221;6229;6247;6257;6263;6269;6271;6277;6287;6299;6301;6311;6317;6323;6329;6337;6343;6353;6359;6361;6367;6373;6379;6389;6397;6421;6427;6449;6451;6469;6473;6481;6491;6521;6529;6547;6551;6553;6563;6569;6571;6577;6581;6599;6607;6619;6637;6653;6659;6661;6673;6679;6689;6691;6701;6703;6709;6719;6733;6737;6761;6763;6779;6781;6791;6793;6803;6823;6827;6829;6833;6841;6857;6863;6869;6871;6883;6899;6907;6911;6917;6947;6949;6959;6961;6967;6971;6977;6983;6991;6997;7001;7013;7019;7027;7039;7043;7057;7069;7079;7103;7109;7121;7127;7129;7151;7159;7177;7187;7193;7207;7211;7213;7219;7229;7237;7243;7247;7253;7283;7297;7307;7309;7321;7331;7333;7349;7351;7369;7393;7411;7417;7433;7451;7457;7459;7477;7481;7487;7489;7499;7507;7517;7523;7529;7537;7541;7547;7549;7559;7561;7573;7577;7583;7589;7591;7603;7607;7621;7639;7643;7649;7669;7673;7681;7687;7691;7699;7703;7717;7723;7727;7741;7753;7757;7759;7789;7793;7817;7823;7829;7841;7853;7867;7873;7877;7879;7883;7901;7907;7919;7927;7933;7937;7949;7951;7963;7993;8009;8011;8017;8039;8053;8059;8069;8081;8087;8089;8093;8101;8111;8117;8123;8147;8161;8167;8171;8179;8191;8209;8219;8221;8231;8233;8237;8243;8263;8269;8273;8287;8291;8293;8297;8311;8317;8329;8353;8363;8369;8377;8387;8389;8419;8423;8429;8431;8443;8447;8461;8467;8501;8513;8521;8527;8537;8539;8543;8563;8573;8581;8597;8599;8609;8623;8627;8629;8641;8647;8663;8669;8677;8681;8689;8693;8699;8707;8713;8719;8731;8737;8741;8747;8753;8761;8779;8783;8803;8807;8819;8821;8831;8837;8839;8849;8861;8863;8867;8887;8893;8923;8929;8933;8941;8951;8963;8969;8971;8999;9001;9007;9011;9013;9029;9041;9043;9049;9059;9067;9091;9103;9109;9127;9133;9137;9151;9157;9161;9173;9181;9187;9199;9203;9209;9221;9227;9239;9241;9257;9277;9281;9283;9293;9311;9319;9323;9337;9341;9343;9349;9371;9377;9391;9397;9403;9413;9419;9421;9431;9433;9437;9439;9461;9463;9467;9473;9479;9491;9497;9511;9521;9533;9539;9547;9551;9587;9601;9613;9619;9623;9629;9631;9643;9649;9661;9677;9679;9689;9697;9719;9721;9733;9739;9743;9749;9767;9769;9781;9787;9791;9803;9811;9817;9829;9833;9839;9851;9857;9859;9871;9883;9887;9901;9907;9923;9929;9931;9941;9949;9967;9973},
    n, IF(ISOMITTED(n), ROWS(array), n),
    // Step 3: Filter primes within the range [gte, lte].
    filtered, FILTER( array, 
        1*(array <= lte) * (array >= gte)
    ),
    // Step 4: Return the first n or last n primes based on the sign of n.
    IF( n < ROWS(filtered),
        SORT( TAKE(filtered, n, 1), 1, SIGN(n)*1 ),
        filtered
    )
));

/* ******************* 
* TEXT FORMULAS
* ********************/

/**
* Removes all of the characters in chars (individually) from text
*/
REMOVECHARS = LAMBDA(text,remove_chars,
LET(doc,"https://www.vertex42.com/lambda/removechars.html",
    TRIM(REDUCE(text, MID(remove_chars,SEQUENCE(LEN(remove_chars)),1),
        LAMBDA(i,a,SUBSTITUTE(i,a,""))
    ))
));

/**
* Adds a single space before each capital letter A through Z
*/
ADDSPACEBEFORECAP = LAMBDA(text,
LET(doc,"https://www.vertex42.com/lambda/addspacebeforecap.html",
    TRIM(REDUCE(text, CHAR(SEQUENCE(1,26,65)),
        LAMBDA(i,a,SUBSTITUTE(i,a," "&a))
    ))
));

/**
* Count the number of instances of a character or string within text
*/
COUNTCHAR = LAMBDA(char,within_text,
LET(doc,"https://www.vertex42.com/lambda/countchar.html",
    (LEN(within_text)-LEN(SUBSTITUTE(within_text,char,""))) / LEN(char)
));
/**
* Converts a text string to a row (dim=1, default) or column (dim=2) of characters
*/
TEXT2ARRAY = LAMBDA(text,[dim],
LET(doc,"https://www.vertex42.com/lambda/text2array.html",
    dim,IF(ISOMITTED(dim),1,dim),
    text_array,IF(COLUMNS(text)>1,TRANSPOSE(text),text),
    result,REDUCE("",SEQUENCE(ROWS(text_array)),LAMBDA(acc,i,
        LET(cell,INDEX(text_array,i,1),
            //splitArray,IF(cell="","",MID(cell,SEQUENCE(1,LEN(cell)),1)),
            splitArray,IF(cell="","",REGEXEXTRACT(cell,".",1)),
            IF(i=1,splitArray,VSTACK(acc,splitArray))
        )
    )),
    cleaned,MAP(result,LAMBDA(_v,IF(ISNA(_v),"",_v))),
    IF(dim=2,TRANSPOSE(cleaned),cleaned)
));
/**
* Converts a Roman numeral to its integer value: IV = 4, XII = 12, etc.
*/
ROMAN2INT = LAMBDA(text,
LET(doc,"https://www.vertex42.com/lambda/roman2int.html",
    MAP(text,LAMBDA(cell,LET(
        strarray,TRANSPOSE(MID(UPPER(cell),SEQUENCE(LEN(cell),1),1)),
        roman_char,{"I","V","X","L","C","D","M"},
        values,{1,5,10,50,100,500,1000},
        strvalues,XLOOKUP(strarray,roman_char,values,"#Error: Invalid Roman Numeral"),
        IF(LEN(cell)=1,strvalues,
            SUM(strvalues*HSTACK(2*(DROP(strvalues,,1)<=DROP(strvalues,,-1))-1,1)
        )
    ))))
));


/* ****************
* ARRAY FORMULAS
******************* */

/**
* Reverses the order of the rows of an array, unless it is a single row or 
* dimension=2, then reverses the columns
*/
FLIP = LAMBDA(array,[dimension],
LET(doc,"https://www.vertex42.com/lambda/flip.html",
    dimension,IF(ISOMITTED(dimension),1,dimension),
    rows,ROWS(array),
    cols,COLUMNS(array),
    IF(OR(dimension=2,rows=1),
        CHOOSECOLS(array,SEQUENCE(1,cols,cols,-1)),
        CHOOSEROWS(array,SEQUENCE(rows,1,rows,-1))
    )
));
/**
* Reverses the order of the columns of an array
*/
FLIPLR = LAMBDA(array,
LET(doc,"https://www.vertex42.com/lambda/flip.html",
    cols,COLUMNS(array),
    CHOOSECOLS(array,SEQUENCE(1,cols,cols,-1))
));
/**
* Reverses the order of the rows of an array
*/
FLIPUD = LAMBDA(array,
LET(doc,"https://www.vertex42.com/lambda/flip.html",
    rows,ROWS(array),
    CHOOSEROWS(array,SEQUENCE(rows,1,rows,-1))
));
/**
* Rotate an array 90 degrees counterclockwise n times
*/
ROT_90 = LAMBDA(array,[n],
LET(doc,"https://www.vertex42.com/lambda/rot90.html",
    r9Flip,LAMBDA(arr,dim,
        IF(dim=2,
            CHOOSECOLS(arr,SEQUENCE(1,COLUMNS(arr),COLUMNS(arr),-1)),
            CHOOSEROWS(arr,SEQUENCE(ROWS(arr),1,ROWS(arr),-1))
        )
    ),
    IF(ISOMITTED(n),
        r9Flip(TRANSPOSE(array),1),
        CHOOSE(1+MOD(INT(n)-1,4),
            r9Flip(TRANSPOSE(array),1),
            r9Flip(r9Flip(array,2),1),
            r9Flip(TRANSPOSE(array),2),
            array
        )
    )
));
/**
* Shift the rows (dim=1) or columns (dim=2) of an array circularly n times.
*/
CIRCSHIFT = LAMBDA(array,n,[dimension],
LET(doc,"https://www.vertex42.com/lambda/circshift.html",
    rows,ROWS(array),cols,COLUMNS(array),
    by_col,dimension=2,
    dim,IF(by_col,cols,rows),
    m,IF(by_col,1+MOD(INT(n)-1,cols),1+MOD(INT(n)-1,rows)),
    indices,VSTACK(SEQUENCE(m,1,dim-m+1),SEQUENCE(dim-m,1)),
    IF(dim=m,array,
        IF(by_col,
            CHOOSECOLS(array,indices),
            CHOOSEROWS(array,indices)
        )
    )
));
/**
* Repeat an array m times vertically and n times horizontally
*/
/*
* Inputs:
*   array: The input array or range to repeat.
*   [m_vert]: (Optional, Default=1) The number of times to repeat the array vertically.
*   [n_horiz]: (Optional, Default=1) The number of times to repeat the array horizontally.
*
* Outputs:
*   Returns a new array where the input array is repeated vertically and horizontally.
*
* Notes:
*   - m_vert and n_horiz should be whole number scalar values
*/
REPARRAY = LAMBDA(array, [m_vert], [n_horiz],
LET(doc, "https://www.vertex42.com/lambda/reparray.html",
    version, "1/13/2025 - Updated commenting",
    // Handle defaults for m_vert and n_horiz
    m_vert, IF(ISBLANK(m_vert), 1, m_vert),
    n_horiz, IF(ISBLANK(n_horiz), 1, n_horiz),
    // Create the repeated array using MAKEARRAY
    MAKEARRAY(
        m_vert * ROWS(array),      // Total rows in the output array
        n_horiz * COLUMNS(array),  // Total columns in the output array
        LAMBDA(i, j,
            INDEX(
                array,
                1 + MOD(i - 1, ROWS(array)),    // Cyclic row index
                1 + MOD(j - 1, COLUMNS(array))  // Cyclic column index
            )
        )
    )
));
/**
* Repeat the elements of an array m times vertically and n times horizontally
*/
/*
* Inputs:
*   array: The input array or range whose elements are to be repeated.
*   [m_vert]: (Optional, Default=1) A scalar or vector specifying the number of times to repeat each row vertically.
*   [n_horiz]: (Optional, Default=1) A scalar or vector specifying the number of times to repeat each column horizontally.
*
* Outputs:
*   Returns a new array where each element of the input array is repeated as specified by m_vert and n_horiz.
*
* Notes:
*   - If m_vert or n_horiz is a scalar, all rows or columns are repeated equally.
*   - If m_vert or n_horiz is a vector, the corresponding rows or columns are repeated by the values in the vector.
*/
REPELEM = LAMBDA(array, [m_vert], [n_horiz],
LET(doc, "https://www.vertex42.com/lambda/repelem.html",
    version, "1/13/2025 - Updated commenting",
    // Handle defaults and reshape m_vert to column and n_horiz to row
    m_vert, IF(ISOMITTED(m_vert), 1, IF(COLUMNS(m_vert) > 1, TRANSPOSE(m_vert), m_vert)),
    n_horiz, IF(ISOMITTED(n_horiz), 1, IF(ROWS(n_horiz) > 1, TRANSPOSE(n_horiz), n_horiz)),
    // Get the dimensions of the input array
    rows, ROWS(array),
    cols, COLUMNS(array),
    // Create a vector specifying the number of repetitions for each row
    m_array, IF(ROWS(m_vert) > 1, m_vert, SEQUENCE(rows, 1, INDEX(m_vert, 1, 1), 0)),
    // Create a vector specifying the number of repetitions for each column
    n_array, IF(COLUMNS(n_horiz) > 1, n_horiz, SEQUENCE(1, cols, INDEX(n_horiz, 1, 1), 0)),
    // Calculate the total size of the output array
    total_rows, SUM(m_array),
    total_cols, SUM(n_array),
    // Compute cumulative sums for row and column repetitions
    csum_vert, SCAN(0, m_array, LAMBDA(a, b, a + b)),
    csum_horiz, SCAN(0, n_array, LAMBDA(a, b, a + b)),
    // Map each output index to the corresponding input index
    row_indices, XMATCH(SEQUENCE(total_rows), csum_vert, 1),
    col_indices, XMATCH(SEQUENCE(total_cols), csum_horiz, 1),
    // Generate the final array
    MAKEARRAY(total_rows, total_cols, LAMBDA(i, j,
        INDEX(array, INDEX(row_indices, i, 1), INDEX(col_indices, j, 1))
    ))
));
/**
* Replace a block within an array by specifying the starting (i,j) location
*/
REPLACEBLOCK = LAMBDA(array,i,j,new_block,
LET(doc,"https://www.vertex42.com/lambda/replaceblock.html",
    mrows,ROWS(new_block),
    ncols,COLUMNS(new_block),
    MAKEARRAY(ROWS(array),COLUMNS(array),
        LAMBDA(r,c,
            IF(AND((r-i+1)>0,(r-i)<mrows,(c-j+1)>0,(c-j)<ncols),
                INDEX(new_block,r-i+1,c-j+1),
                INDEX(array,r,c)
            )
        )
    )
));

/**
* Return a subset of rows and columns from an array
*/
SLICE = LAMBDA(array,row_start,row_end,[col_start],[col_end],
LET(doc,"https://www.vertex42.com/lambda/slice.html",
    rows,ROWS(array),
    cols,COLUMNS(array),
    row_start,IF(ISBLANK(row_start),1,
        IF(row_start<0,IF(ABS(row_start)>=rows,1,rows+row_start+1),row_start)
    ),
    row_end,IF(OR(ISBLANK(row_end),row_end>rows),rows,
        IF(row_end<0,IF(ABS(row_end)>=rows,1,rows+row_end+1),row_end)
    ),
    col_start,IF(ISBLANK(col_start),1,
        IF(col_start<0,IF(ABS(col_start)>=cols,1,cols+col_start+1),col_start)
    ),
    col_end,IF(OR(ISBLANK(col_end),col_end>cols),cols,
        IF(col_end<0,IF(ABS(col_end)>=cols,1,cols+col_end+1),col_end)
    ),
    new_array,CHOOSEROWS(array,SEQUENCE(row_end-row_start+1,,row_start,1)),
    CHOOSECOLS(new_array,SEQUENCE(col_end-col_start+1,,col_start,1))
));

/**
* Splice an array without optional delete and insert
*/
SPLICE = LAMBDA(array,start_index,delete_count,[insert_array],[by_col],
LET(doc,"https://www.vertex42.com/lambda/splice.html",
    start,IF(start_index<1,1,start_index),
    by_col,IF(ISOMITTED(by_col),FALSE,by_col),
    no_insert,ISOMITTED(insert_array),
    no_first,start<=1,
    first,IF(by_col,TAKE(array,,start-1),TAKE(array,start-1,)),
    no_last,delete_count>IF(by_col,COLUMNS(array),ROWS(array))-start,
    last,IF(by_col,DROP(array,,start+delete_count-1),DROP(array,start+delete_count-1,)),
    tf,CONCATENATE(
        IF(no_first,"T","F"),
        IF(no_insert,"T","F"),
        IF(no_last,"T","F"),
        IF(by_col,"T","F")
    ),
    SWITCH(tf,
        "TTTT",NA(),"TTTF",NA(),
        "TTFT",last,"TTFF",last,
        "TFTT",insert_array,"TFTF",insert_array,
        "TFFT",HSTACK(insert_array,last),
        "TFFF",VSTACK(insert_array,last),
        "FTTT",first,"FTTF",first,
        "FTFT",HSTACK(first,last),
        "FTFF",VSTACK(first,last),
        "FFTT",HSTACK(first,insert_array),
        "FFTF",VSTACK(first,insert_array),
        "FFFT",HSTACK(first,insert_array,last),
        "FFFF",VSTACK(first,insert_array,last),
    )
));

/**
* Return an array of all combinations of rows from multiple arrays.
*/
/* Inputs:
*   array_1, array_2: Required. The first two arrays to combine.
*   [array_3], [array_4], [array_5]: Optional. Additional arrays to include in the combinations.
*
* Outputs:
*   Returns an array where each row is a unique combination of rows from the input arrays.
*
* Requires: REPELEM and REPARRAY
*/
COMBINATIONS = LAMBDA(array_1, array_2, [array_3], [array_4], [array_5],
LET(doc,"https://www.vertex42.com/lambda/combinations.html",
    version, "1.0.1 1/13/2025 - Multiple optional arrays",
    // Combine two arrays into all possible row combinations
    combine_two_arrays, LAMBDA(a, b,
        LET(
            r_1, ROWS(a), r_2, ROWS(b),
            first, REPELEM(a, r_2, 1),
            second, REPARRAY(b, r_1, 1),
            HSTACK(first, second)
        )
    ),
    // Initialize with the first two arrays
    result, combine_two_arrays(array_1, array_2),
    // Extend the combinations by adding optional arrays if provided
    result2, IF(ISOMITTED(array_3), result, combine_two_arrays(result, array_3)),
    result3, IF(ISOMITTED(array_4), result2, combine_two_arrays(result2, array_4)),
    result4, IF(ISOMITTED(array_5), result3, combine_two_arrays(result3, array_5)),
    // Return the final combinations
    result4
));

/**
* Permutations of values in an array chosen n at a time with repetition
*/
PERMUTATIONS = LAMBDA(array,n_chosen,
LET(doc,"https://www.vertex42.com/lambda/permutations.html",
    array,IF(AND(COLUMNS(array)>1,ROWS(array)=1),TRANSPOSE(array),array),
    REDUCE("",SEQUENCE(n_chosen-1),LAMBDA(acc,i,
        IF(i=1,COMBINATIONS(array,array),COMBINATIONS(array,acc))
    ))
));
/**
* Return all the combinations of (N Choose K) from an Nx1 array
*/
COMBINR = LAMBDA(array,k,
LET(doc,"https://www.vertex42.com/lambda/combinr.html",
    array,IF(AND(COLUMNS(array)>1,ROWS(array)=1),TRANSPOSE(array),array),
    cols,COLUMNS(array),
    n,ROWS(array),
    IF(cols>1,"Error: Array size must be n x 1",
    IF(n<3,"Error: n must be >= 3",
    IF(k>n,"Error: k must be < Rows",
    LET(
    m,COMBIN(n,k),
    ms,SEQUENCE(m),
    mx,SEQUENCE(1,k,n-k+1,1),
    combos,REDUCE(SEQUENCE(m,k,0,0),ms,
    LAMBDA(acc,i,
        IF(i=1,SEQUENCE(1,k),
        LET(prev_row,INDEX(acc,i-1,0),
            tf,prev_row=mx,
            col2inc,IFERROR(MATCH(TRUE,tf,0),k+1)-1,
            new_row,IF(col2inc=1,
                SEQUENCE(1,k-col2inc+1,INDEX(prev_row,1,col2inc)+1),
                HSTACK(
                    CHOOSECOLS(prev_row,SEQUENCE(1,col2inc-1)),
                    SEQUENCE(1,k-col2inc+1,INDEX(prev_row,1,col2inc)+1)
                )
            ),
            VSTACK(acc,new_row)
        ))
    )),
    INDEX(array,combos)
    ))))
));

/**
* Returns the vector {num_rows;num_columns}
*/
SIZE = LAMBDA(array,
    VSTACK(ROWS(array),COLUMNS(array))
);
/**
* Returns TRUE if one of the dimensions is 1
*/
ISVECTOR = LAMBDA(array,
    OR(ROWS(array)=1,COLUMNS(array)=1)
);
/**
* Returns TRUE if columns>1 and rows>1
*/
ISARRAY = LAMBDA(array,
    OR(COLUMNS(array)>1,ROWS(array)>1)
);
/**
* Returns TRUE if columns=1 and rows>1
*/
ISCOLUMN = LAMBDA(array,
    AND(ROWS(array)>1,COLUMNS(array)=1)
);
/**
* Returns TRUE if columns>1 and rows=1
*/
ISROW = LAMBDA(array,
    AND(COLUMNS(array)>1,ROWS(array)=1)
);
/**
* Returns TRUE if the matrix is symmetric
*/
ISSYMMETRIC = LAMBDA(matrix,
LET(doc,"https://www.vertex42.com/lambda/",
    IF(ROWS(matrix)<>COLUMNS(matrix),
        "Error: matrix not square",
        SUM(--(matrix<>TRANSPOSE(matrix)))=0
    )
));
/**
* Returns TRUE if the matrix is positive definite
* Checks if all upper-left determinants are greater than zero
*/
ISPOSDEF = LAMBDA(matrix,
LET(doc,"https://www.vertex42.com/lambda/",
    REDUCE(TRUE,SEQUENCE(ROWS(matrix)),LAMBDA(acc,i,
        AND(acc,MDETERM(CHOOSEROWS(CHOOSECOLS(matrix,SEQUENCE(i)),SEQUENCE(i)))>0)
    ))
));

/**
* Rotate an array 45 degrees counterclockwise, making diagonals become rows
*/
/*
* Inputs:
*   array: The 2D array to rotate.
*   [n]: (Optional) Number of times to rotate by 45 degrees counterclockwise. Default is 1.
*        Accepts values from 1 to 8 or -1 to -8 for multiples of 45° rotations (e.g., 1=45°, 2=90°).
*        The rotation is modular, so values beyond 8 are reduced to their remainder modulo 8.
*   [fill_with]: (Optional) Choose to fill empty spaces with a value. Default is NA().
*
* Outputs: 
*   - Returns the rotated array, depending on the number of rotations.
*   - For odd rotations (e.g., 45° or 135°), diagonals become rows and columns.
*   - For even rotations (e.g., 90° or 180°), uses the ROT_90 function.
*
* Notes:
*   - Uses modular rotation logic to handle rotations beyond 360° (e.g., 9 = 45°).
*   - Empty spaces are filled with #N/A.
*   - Requires supporting functions `ROT_90` and `FLIPLR`.
*   - IMPORTANT: This algorithm does not work in sequence. ROT_45(ROT_45(array)) is not ROT_90(array).
*/
ROT_45 = LAMBDA(array, [n], [fill_with],
LET(doc, "https://www.vertex42.com/lambda/rot_45.html",
    version,"1.0.0 1/13/2025 - Original version",
    // Default values for optional parameters
    n, IF(ISBLANK(n), 1, MOD(n, 8)), // Rotations cycle every 8 steps

    // Array dimensions
    rows, ROWS(array),
    cols, COLUMNS(array),
    new_height, rows + cols - 1,

    // Generate the rotated array based on the rotation step (n)
    rotated, SWITCH(n,
        // 45° counterclockwise rotation
        1, MAKEARRAY(new_height, cols,
                LAMBDA(r, c,
                    IFERROR(INDEX(array, IF(r - cols + c = 0, -1, r - cols + c), c), NA())
                )
            ),
        // 90° rotation
        2, ROT_90(array, 1),
        // 135° counterclockwise rotation. Same as 315° flipped.
        3, FLIPLR(MAKEARRAY(new_height, cols,
                LAMBDA(r, c,
                    IFERROR(INDEX(array, IF(r - c + 1 = 0, -1, r - c + 1), c), NA())
                )
            )),
        // 180° rotation
        4, ROT_90(array, 2),
        // 225° counterclockwise rotation. Same as 45° flipped.
        5, FLIPLR(MAKEARRAY(new_height, cols,
                LAMBDA(r, c,
                    IFERROR(INDEX(array, IF(r - cols + c = 0, -1, r - cols + c), c), NA())
                )
            )),
        // 270° rotation
        6, ROT_90(array, 3),
        // 315° counterclockwise rotation
        7, MAKEARRAY(new_height, cols,
                LAMBDA(r, c,
                    IFERROR(INDEX(array, IF(r - c + 1 = 0, -1, r - c + 1), c), NA())
                )
            ),
        // 0° or 360° just returns the original array
        array
    ),
    result,IF(ISOMITTED(fill_with), rotated, IFNA(rotated, fill_with)),
    result
));


/* ****************
* MATRIX FORMULAS
******************* */

/**
* ONES(array) :: Returns a matrix of 1s the size of array
* ONES(m,[n]) :: Returns a matrix of 1s the size of m x n (default n=1)
*/
ONES = LAMBDA(m_rows,[n_columns],
LET(doc,"https://www.vertex42.com/lambda/ones-and-zeros.html",
    rows,ROWS(m_rows),cols,COLUMNS(m_rows),
    IF(OR(rows>1,cols>1),
        SEQUENCE(rows,cols,1,0),
        IF(ISOMITTED(n_columns),
            SEQUENCE(m_rows,1,1,0),
            SEQUENCE(m_rows,n_columns,1,0)
        )
    )
));
/**
* ZEROS(array) :: Returns a matrix of 0s the size of array
* ZEROS(m,[n]) :: Returns a matrix of 0s the size of m x n (default n=1)
*/
ZEROS = LAMBDA(m_rows,[n_columns],
LET(doc,"https://www.vertex42.com/lambda/ones-and-zeros.html",
    rows,ROWS(m_rows),cols,COLUMNS(m_rows),
    IF(OR(rows>1,cols>1),
        SEQUENCE(rows,cols,0,0),
        IF(ISOMITTED(n_columns),
            SEQUENCE(m_rows,1,0,0),
            SEQUENCE(m_rows,n_columns,0,0)
        )
    )
));
/**
* Convert a vector to a Diagonal matrix or vice versa
*/
DIAG = LAMBDA(array,
LET(doc,"https://www.vertex42.com/lambda/diag.html",
    r, ROWS(array),
    c, COLUMNS(array),
    IF(c=1,
        MAKEARRAY(r,r,LAMBDA(i,j,IF(i=j,INDEX(array,i),0))),
        IF(r <> c, "Error: Not Square",
        MAKEARRAY(r,1,LAMBDA(i,j,INDEX(array,i,i)))
        )
    )
));
/**
* Sum of the Diagonal of a square matrix
*/
TRACE = LAMBDA(matrix,
LET(doc,"https://www.vertex42.com/lambda/diag.html",
    SUM( matrix*MUNIT(ROWS(matrix)))
));
/**
* Returns a column vector containing row sums
*/
ROWSUM = LAMBDA(matrix,
LET(doc,"https://www.vertex42.com/lambda/rowsum-and-colsum.html",
    BYROW(matrix,LAMBDA(row,SUM(row)))
));
/**
* Returns a row vector containing column sums
*/
COLSUM = LAMBDA(matrix,
LET(doc,"https://www.vertex42.com/lambda/rowsum-and-colsum.html",
    BYCOL(matrix,LAMBDA(col,SUM(col)))
));
/**
* Returns the sum of the element-wise multiplication of two vectors
* or the column sum of the element-wise multiplication of two matrices
*/
DOT = LAMBDA(a,b,
LET(doc,"https://www.vertex42.com/lambda/dot.html",
    a,IF(AND(ROWS(a)=1,COLUMNS(a)>1),TRANSPOSE(a),a),
    b,IF(AND(ROWS(b)=1,COLUMNS(b)>1),TRANSPOSE(b),b),
    IF(OR(ROWS(a)<>ROWS(b),COLUMNS(a)<>COLUMNS(b)),
        "Error: a and b must be the same size",
        BYCOL(a*b,LAMBDA(col,SUM(col)))
    )
));
/**
* Returns the Cross Product of two 3x1 vectors a and b.
*/
CROSS = LAMBDA(vector_a,vector_b,
LET(doc,"https://www.vertex42.com/lambda/cross.html",
    toColVec3,LAMBDA(vec,LET(
        cvec,IF(AND(ROWS(vec)=1,COLUMNS(vec)>1),TRANSPOSE(vec),vec),
        IF(ROWS(cvec)=2,VSTACK(cvec,0),cvec)
    )),
    vec_a,toColVec3(vector_a),
    vec_b,toColVec3(vector_b),
    arr,VSTACK(TRANSPOSE(vec_a),TRANSPOSE(vec_b)),
    IF( OR(ROWS(vec_a)<>3,ROWS(vec_b)<>3,COLUMNS(vec_a)>1,COLUMNS(vec_b)>1),
        "Error in vector size",
        VSTACK(
            MDETERM(CHOOSECOLS(arr,2,3)),
            MDETERM(CHOOSECOLS(arr,3,1)),
            MDETERM(CHOOSECOLS(arr,1,2))
        )
    )
));
/**
* Returns the magnitude (2-norm) of a vector or the magnitude of each column of a matrix
*/
MAGNITUDE = LAMBDA(vector,
LET(doc,"https://www.vertex42.com/lambda/magnitude.html",
    vec,IF(AND(ROWS(vector)=1,COLUMNS(vector)>1),TRANSPOSE(vector),vector),
    IF(COLUMNS(vec)>1,
        SQRT(COLSUM(vec*vec)),
        SQRT(SUM(vec*vec))
    )
));
/**
* Returns a Pascal matrix of size nxn. Optionally specify
* type="L" or type="U" for the lower or upper triangle form.
*/
PASCAL = LAMBDA(n,[type],
LET(doc,"https://www.vertex42.com/lambda/pascal.html",
IF(n<2,"Error:n<2",
    SWITCH(type,
    "L",MAKEARRAY(n,n,LAMBDA(r,c,LET(i,r-1,j,c-1,
        IF(j<=i,FACT(i)/(FACT(j)*FACT(i-j)),0)
    ))),
    "U",MAKEARRAY(n,n,LAMBDA(r,c,LET(i,r-1,j,c-1,
        IF(i<=j,FACT(j)/(FACT(i)*FACT(j-i)),0)
    ))),
    REDUCE(SEQUENCE(1,n,1,0),SEQUENCE(n-1),LAMBDA(acc,row,
        VSTACK(acc,HSTACK(1,
            SCAN(1,SEQUENCE(1,n-1),
                LAMBDA(csum,i,csum+INDEX(acc,row,i+1))
            )
        ))
    )))
)));
/**
* Return the Hessenberg of a square matrix using the Householder transformation algorithm.
*/
HESS = LAMBDA(matrix,[output_form],
IF(ROWS(matrix)<>COLUMNS(matrix),"Error: not square",
IF(ROWS(matrix)<3,"Must be > 2x2",
LET(doc,"https://www.vertex42.com/labmda/hess.html",
    output_form,IF(ISOMITTED(output_form),"H",output_form),
    n,ROWS(matrix),
    ks,SEQUENCE(n-2),
    HUk,REDUCE(VSTACK(matrix,MUNIT(n)),ks,LAMBDA(acc,k,
        LET(Amat,DROP(acc,-n),
            x,DROP(INDEX(Amat,0,k),k),
            v,x+SIGN(INDEX(x,1))*VSTACK(SQRT(SUM(x^2)),SEQUENCE(ROWS(x)-1,1,0,0)),
            Pk,MUNIT(ROWS(x))-2*MMULT(v,TRANSPOSE(v))/MMULT(TRANSPOSE(v),v),
            Uk,VSTACK(HSTACK(MUNIT(k),SEQUENCE(k,n-k,0,0)),HSTACK(SEQUENCE(n-k,k,0,0),Pk)),
            new_A,MMULT(MMULT(Uk,Amat),Uk),
            VSTACK(new_A,MMULT(Uk,DROP(acc,n)))
        ))
    ),
    H,DROP(HUk,-n,0),
    IF(output_form="HU",VSTACK(H,TRANSPOSE(DROP(HUk,n,0))),H)
))));
/**
* Returns the QR decomposition of matrix A using the Householder transformation process
* QR(matrix) returns {Q;R} (stacked vertically) where matrix=MMULT(Q,R).
*/
QR = LAMBDA(matrix,
LET(doc,"https://www.vertex42.com/lambda/qr.html",
    m,ROWS(matrix),n,COLUMNS(matrix),
    ks,SEQUENCE(n-1,1),
    Q3Q2Q1,REDUCE(SEQUENCE(n,n,0,0),ks,LAMBDA(acc,k,
        LET(   
            Ahat,IF(k=1,matrix,MMULT(acc,matrix)),
            ek,VSTACK(1,SEQUENCE(n-k,1,0,0)),
            xk,DROP(INDEX(Ahat,0,k),k-1),
            alpha,-SIGN(INDEX(xk,1))*SQRT(SUM(xk^2)),
            uk,xk-alpha*ek,
            vk,uk/SQRT(SUM(uk^2)),
            Qkhat,MUNIT(ROWS(ek))-2*MMULT(vk,TRANSPOSE(vk)),
            Qk,IF(k=1,Qkhat,
                VSTACK(
                HSTACK(MUNIT(k-1), SEQUENCE(k-1,n-k+1,0,0)),
                HSTACK(SEQUENCE(n-k+1,k-1,0,0), Qkhat))
                ),
            IF(k=1,Qk,MMULT(Qk,acc))
        )
    )),
    VSTACK(TRANSPOSE(Q3Q2Q1),MMULT(Q3Q2Q1,matrix))
));
/**
* Attempts to find Eigenvalues of a square matrix with n iterations
* using a QR algorithm with Rayleigh shifts and optional Hessenberg
*/
EIGENVALUE = LAMBDA(matrix,[iterations],[output_form],[use_hess],
LET(doc,"https://www.vertex42.com/lambda/eigenvalue.html",
    n,ROWS(matrix),
    iterations,IF(ISBLANK(iterations),42,iterations),
    output_form,IF(ISBLANK(output_form),"U",output_form),
    initialHU,IF(use_hess=TRUE,HESS(matrix,"HU"),""),
    initialA,IF(use_hess=TRUE,DROP(initialHU,-n,0),matrix),
    UQ,REDUCE(VSTACK(initialA,MUNIT(n)),SEQUENCE(iterations),LAMBDA(acc,k,
        LET(new_A,DROP(acc,-n,0),
            shift,IF(INDEX(new_A,n,n)=0,0.0001,INDEX(new_A,n,n)),
            QRmat,QR(new_A-shift*MUNIT(n)),
            Q,TAKE(QRmat,n),
            R,TAKE(QRmat,-n),
            next_A,MMULT(R,Q)+shift*MUNIT(n),
            VSTACK(next_A,MMULT(DROP(acc,n,0),Q))
        )
    )),
    eigMat,DROP(UQ,-n,0),
    qMat,IF(use_hess=TRUE,
        MMULT(DROP(initialHU,n,0),DROP(UQ,n,0)),
        DROP(UQ,n,0)
    ),
    IF(output_form="UQ",UQ,
    IF(output_form="test",
        LET(
            eigRow,MAKEARRAY(1,n,LAMBDA(i,j,INDEX(eigMat,j,j))),
            label,LAMBDA(label,HSTACK(label,MAKEARRAY(1,n-1,LAMBDA(i,j,"")))),
            VSTACK(
                eigMat,
                label("Q"),qMat,
                label("λ"),eigRow,
                label("det(A-λI)"),BYCOL(eigRow,LAMBDA(λ,MDETERM(matrix-λ*MUNIT(n)))),
                label("Av/v"),MMULT(matrix,qMat)/qMat
                )
        ),
        eigMat
    ))
));
/**
* Returns the Cholesky decomposition for a symmetric positive-definite matrix
*/
CHOLESKY = LAMBDA(matrix,
LET(doc,"https://www.vertex42.com/lambda/cholesky.html",
    n,ROWS(matrix),
    IF(ROWS(matrix)<>COLUMNS(matrix),
        "Error: matrix not square",
    IF(SUM(--(matrix<>TRANSPOSE(matrix)))<>0,
        "Error: not symmetric",
    IF(MDETERM(matrix)<=0,
        "Error: not positive definite",
    REDUCE(0,SEQUENCE(n),LAMBDA(Lmat,j,
        IF(j=1,
            LET(Ljj,SQRT(INDEX(matrix,j,j)),
                Lij,CHOOSEROWS(INDEX(matrix,0,1),SEQUENCE(n-1,1,2,1))/Ljj,
                VSTACK(Ljj,Lij)
            ),
        IF(j=n,
            LET(Ljj,SQRT(INDEX(matrix,j,j)-SUM(CHOOSEROWS(Lmat,j)^2)),
                HSTACK(Lmat,VSTACK(SEQUENCE(n-1,1,0,0),Ljj))
            ),
            LET(Ljj,SQRT(INDEX(matrix,j,j)-SUM(CHOOSEROWS(Lmat,j)^2)),
                Aij,CHOOSEROWS(INDEX(matrix,0,j),SEQUENCE(n-j,1,j+1,1)),
                Lik,CHOOSEROWS(Lmat,SEQUENCE(n-j,1,j+1,1)),
                Ljk,CHOOSEROWS(Lmat,j),
                Lij,1/Ljj*(Aij-MMULT(Lik,TRANSPOSE(Ljk))),
                HSTACK(Lmat,VSTACK(SEQUENCE(j-1,1,0,0),Ljj,Lij))
            )
        ))
    ))
    )))
));

/******************** 
* INTERPOLATION
*********************/

/**
* Linearly interpolate between the two nearest points in a table lookup
*/
LINTERP = LAMBDA(xs,known_xs,known_ys,
LET(doc,"https://www.vertex42.com/lambda/linterp.html",
    xs,IF(AND(ROWS(xs)=1,COLUMNS(xs)>1),TRANSPOSE(xs),xs),
    known_xs,IF(AND(ROWS(known_xs)=1,COLUMNS(known_xs)>1),TRANSPOSE(known_xs),known_xs),
    known_ys,IF(AND(ROWS(known_ys)=1,COLUMNS(known_ys)>1),TRANSPOSE(known_ys),known_ys),
    tab,SORT(HSTACK(known_xs,known_ys)),
    BYROW(xs,LAMBDA(x,
        LET(ind,MATCH(x,INDEX(tab,0,1),1),
            pts,
            IF(x<=MIN(known_xs),TAKE(tab,2),
                IF(x>=MAX(known_xs),TAKE(tab,-2),
                    CHOOSEROWS(tab,ind,ind+1)
                )
            ),
            SLOPE(INDEX(pts,,2),INDEX(pts,,1))*x+INTERCEPT(INDEX(pts,,2),INDEX(pts,,1))
        )
    ))
));
/**
* Polynomial interpolation between the n+1 closest points based on distance (x0-x)^2
*/
PINTERP = LAMBDA(xs,known_xs,known_ys,n,
LET(doc,"https://www.vertex42.com/lambda/pinterp.html",
    xs,IF(AND(ROWS(xs)=1,COLUMNS(xs)>1),TRANSPOSE(xs),xs),
    known_xs,IF(AND(ROWS(known_xs)=1,COLUMNS(known_xs)>1),TRANSPOSE(known_xs),known_xs),
    known_ys,IF(AND(ROWS(known_ys)=1,COLUMNS(known_ys)>1),TRANSPOSE(known_ys),known_ys),
    BYROW(xs,LAMBDA(x,
        LET(tab,TAKE(SORT(HSTACK((known_xs-x)^2,known_xs,known_ys)),n+1),
            xo,INDEX(tab,1,2),
            POLYVAL(POLYFIT(INDEX(tab,0,2)-xo,INDEX(tab,0,3),n),x-xo)
        )
    ))
));
/**
* Cubic Spline - Creates a cubic piecewise interpolating polynomial 
* by specifying control points and slopes at each point.
*/
/*
* Inputs:
*   known_xs: A vector of x-coordinates of the control points.
*   known_ys: A vector of y-coordinates of the control points.
*   [x]: (Optional) The x-values at which the spline should be evaluated.
*        If omitted, the function returns the piecewise polynomial data structure.
*   [ms]: (Optional) Slopes at the control points. If omitted, slopes are estimated.
*   [c]: (Optional) Tension parameter for Cardinal splines (default: 0).
*
* Outputs:
*   - Returns either the evaluated spline at x or the piecewise polynomial structure.
*
* Requires: `PPVAL` function to evaluate the piecewise polynomial at x.
*/
CSPLINE = LAMBDA(known_xs,known_ys,[x],[ms],[c],
LET(doc,"https://www.vertex42.com/lambda/cspline.html",
    version, "1.0.1 - Updated commenting",
    // Step 1: Handle defaults and ensure xs, ys and ms are column vectors
    xs, IF(AND(ROWS(known_xs)=1, COLUMNS(known_xs)>1), TRANSPOSE(known_xs), known_xs),
    ys, IF(AND(ROWS(known_ys)=1, COLUMNS(known_ys)>1), TRANSPOSE(known_ys), known_ys),
    ms, IF(ISBLANK(ms), "", IF(AND(ROWS(ms)=1, COLUMNS(ms)>1), TRANSPOSE(ms), ms)),
    c, IF(ISBLANK(c), 0, c),
    n, ROWS(xs),

    // Step 2: Calculate slopes at each control point
    // Slopes are estimated if not provided explicity
    slope, (1-c) * IF(OR( ms="", ROWS(ms)<=2 ),
        IF( AND( NOT(ms=""), ROWS(ms)=1 ),
            SEQUENCE(n, 1, ms, 0), // All slopes are the same if ms is a scalar
            MAKEARRAY(n, 1, LAMBDA(i, j,
                IF( i=1,
                    IF( ROWS(ms)=2,
                        INDEX(ms,1),
                        // Forward difference at the first point
                        (INDEX(ys,i+1)-INDEX(ys,i))/(INDEX(xs,i+1)-INDEX(xs,i))
                    ),
                    IF( i=n,
                        IF( ROWS(ms)=2,
                            INDEX(ms,2),
                            // Backward difference at the last point
                            (INDEX(ys,i)-INDEX(ys,i-1))/(INDEX(xs,i)-INDEX(xs,i-1))
                        ),
                        // Central difference at interior points
                        0.5*(
                            (INDEX(ys,i+1)-INDEX(ys,i))/(INDEX(xs,i+1)-INDEX(xs,i)) +
                            (INDEX(ys,i)-INDEX(ys,i-1))/(INDEX(xs,i)-INDEX(xs,i-1))
                        )
                    )
                )
            ))
        ),
        ms
    ),

    // Step 3: Compute piecewise polynomial coefficients
    // Each interval [x_i, x_i+1] is represented by a cubic polynomial:
    // P(x) = a*(x - x_i)^3 + b*(x - x_i)^2 + c*(x - x_i) + d
    pp_coeffs, REDUCE(0, SEQUENCE(n - 1), LAMBDA(acc, i,
        LET(
            dx, INDEX(xs, i+1) - INDEX(xs, i),
            y_1, INDEX(ys, i),
            y_2, INDEX(ys, i+1),
            m_1, INDEX(slope, i)*dx,
            m_2, INDEX(slope, i+1)*dx,
            coeffs, HSTACK(
                (2*y_1 - 2*y_2 + m_1 + m_2) / dx^3,
                (-3*y_1 + 3*y_2 - 2*m_1 - m_2) / dx^2,
                m_1 / dx,
                y_1
            ),
            IF( i=1, coeffs, VSTACK(acc, coeffs))
        )
    )),

    // Step 4: Assemble the piecewise polynomial data structure
    pp, MAP(HSTACK(xs, pp_coeffs), LAMBDA(cell, IFERROR(cell, ""))),

    // Step 5: Interpolate if x values are provided, or return coefficients
    IF(ISOMITTED(x),
        pp,
        PPVAL(pp,x)
    )
));

/**
* Natural Cubic Spline - Creates a C2 Interpolating Spline with d2y/dx2=0 at the endpoints.
*/
/*
* Inputs:
*   known_xs: A column vector of known x-values (independent variable).
*   known_ys: A column vector of known y-values (dependent variable).
*   [x]: Optional. The value(s) to interpolate. If omitted, returns the piecewise polynomial parameters.
*
* Outputs:
*   Returns the interpolated value(s) for x, or the piecewise polynomial coefficients if x is omitted.
*
* Notes:
*   - Ensures C2 continuity and natural spline conditions (second derivative is zero at endpoints).
*   - Supports multiple x values for interpolation.
*   - Returns 0 for extrapolated values beyond the range of known_xs and known_ys.
*
* Requires: POLYVAL, POLYDER, PPVAL
*/
NSPLINE = LAMBDA(known_xs, known_ys, [x],
LET(doc, "https://www.vertex42.com/lambda/nspline.html", 
    version, "1.0.1 - Updated commenting",
    // Step 1: Ensure xs and ys are column vectors
    xs, IF(AND(ROWS(known_xs) = 1, COLUMNS(known_xs) > 1), TRANSPOSE(known_xs), known_xs),
    ys, IF(AND(ROWS(known_ys) = 1, COLUMNS(known_ys) > 1), TRANSPOSE(known_ys), known_ys),

    // Step 2: Calculate differences between consecutive x-values (h) and other preliminaries
    h, DROP(xs, 1) - DROP(xs, -1),
    n, ROWS(xs),
    alpha, 3 / DROP(h, 1) * (DROP(ys, 2) - DROP(DROP(ys, 1), -1)) 
        - 3 / DROP(h, -1) * (DROP(DROP(ys, 1), -1) - DROP(ys, -2)),

    // Step 3: Solve the tridiagonal system for the spline coefficients
    l_mu_z, REDUCE({1, 0, 0}, SEQUENCE(n - 2, 1, 2, 1), LAMBDA(acc, i,
        LET(
            mu_o, INDEX(acc, i - 1, 2),
            z_o, INDEX(acc, i - 1, 3),
            li, 2 * (INDEX(xs, i + 1) - INDEX(xs, i - 1)) - INDEX(h, i - 1) * mu_o,
            mui, INDEX(h, i) / li,
            zi, (INDEX(alpha, i - 1) - INDEX(h, i - 1) * z_o) / li,
            ret, VSTACK(acc, HSTACK(li, mui, zi)),
            IF(i = n - 1, VSTACK(ret, {1, 0, 0}), ret)
        )
    )),
    li, INDEX(l_mu_z, , 1),
    mui, INDEX(l_mu_z, , 2),
    zi, INDEX(l_mu_z, , 3),

    // Step 4: Back-substitute to compute the cubic coefficients
    d_c_b, REDUCE({0, 0, 0}, SEQUENCE(n - 1, 1, n - 1, -1), LAMBDA(acc, j,
        LET(
            co, INDEX(acc, 1, 2),
            cj, INDEX(zi, j) - INDEX(mui, j) * co,
            bj, (INDEX(ys, j + 1) - INDEX(ys, j)) / INDEX(h, j) - 1 / 3 * INDEX(h, j) * (co + 2 * cj),
            dj, (co - cj) / (3 * INDEX(h, j)),
            ret, VSTACK(HSTACK(dj, cj, bj), acc),
            ret
        )
    )),
    pp, MAP(HSTACK(xs, DROP(HSTACK(d_c_b, ys), -1)), LAMBDA(cell, IFERROR(cell, ""))),

    // Step 5: Handle extrapolation (linear interpolation at the endpoints)
    m_1, POLYVAL(POLYDER(DROP(CHOOSEROWS(pp, 1), 0, 1)), 0),
    b_1, INDEX(known_ys, 1, 1) - m_1 * INDEX(known_xs, 1, 1),
    m_n, POLYVAL(POLYDER(DROP(CHOOSEROWS(pp, n - 1), 0, 1)), INDEX(pp, n, 1) - INDEX(pp, n - 1, 1)),
    b_n, INDEX(known_ys, n, 1) - m_n * INDEX(known_xs, n, 1),

    // Step 6: Interpolate or return coefficients
    IF(ISOMITTED(x), pp,
        (x < INDEX(known_xs, 1, 1)) * (m_1 * x + b_1) +
        (x > INDEX(known_xs, n, 1)) * (m_n * x + b_n) +
        IFERROR(PPVAL(pp, x), 0)
    )
));

/**
* C2 Interpolating Cubic Spline with specified end conditions: "not", "free", or slope.
* Solves a system of linear constraint equations. Must use at least 3 control points.
*/
SPLINE = LAMBDA(known_xs,known_ys,[x],[cond_start],[cond_end],[debug],
IF(MAX(ROWS(known_ys),COLUMNS(known_ys))<3,NA(),
LET(doc,"https://www.vertex42.com/lambda/spline.html",
    cond_start,IF(ISOMITTED(cond_start),"not",cond_start),
    cond_end,IF(ISOMITTED(cond_end),"not",cond_end),
    xs,IF(AND(ROWS(known_xs)=1,COLUMNS(known_xs)>1),TRANSPOSE(known_xs),known_xs),
    ys,IF(AND(ROWS(known_ys)=1,COLUMNS(known_ys)>1),TRANSPOSE(known_ys),known_ys),
    h,DROP(xs,1,0)-DROP(xs,-1,0),
    si,DROP(h,-1,0)/DROP(h,1,0),
    n,ROWS(xs),
    dmat,{1,1,1;3,2,1;3,1,0},
    constraints,REDUCE(SEQUENCE(3,3,0,0),SEQUENCE(n-2),LAMBDA(acc,i,LET(
        smat,VSTACK({0,0,0},HSTACK(0,0,-INDEX(si,i,1)),HSTACK(0,-(INDEX(si,i,1)^2),0)),
        IF(i=1,IF(n>3,HSTACK(dmat,smat,SEQUENCE(3,3*(n-3),0,0)),HSTACK(dmat,smat)),
        IF(i=n-2,VSTACK(acc,HSTACK(SEQUENCE(3,3*(n-3),0,0),dmat,smat)),
            VSTACK(acc,HSTACK(SEQUENCE(3,3*(i-1),0,0),dmat,smat,SEQUENCE(3,3*(n-2-i),0,0)))
        ))
    ))),
    amatrix,VSTACK(constraints,
        HSTACK(SEQUENCE(1,3*(n-2),0,0),{1,1,1}),
        SWITCH(cond_start,
            "free",HSTACK(0,1,0,SEQUENCE(1,3*(n-2),0,0)),
            "not",HSTACK(1,0,0,-(INDEX(si,1,1)^3),SEQUENCE(1,3*(n-2)-1,0,0)),
            HSTACK(0,0,1,SEQUENCE(1,3*(n-2),0,0))
        ),
        SWITCH(cond_end,
            "free",HSTACK(SEQUENCE(1,3*(n-2),0,0),3,1,0),
            "not",IF(n>3,
                HSTACK(SEQUENCE(1,3*(n-3),0,0),1,0,0,-(INDEX(si,n-2,1)^3),0,0),
                HSTACK(0,0,0,3,1,0)
            ),
            HSTACK(SEQUENCE(1,3*(n-2),0,0),3,2,1)
        )
    ),
    bmatrix,VSTACK(
        REDUCE({0;0;0},SEQUENCE(n-2),LAMBDA(acc,i,IF(i=1,
            VSTACK(INDEX(ys,2,1)-INDEX(ys,1,1),0,0),
            VSTACK(acc,INDEX(ys,i+1,1)-INDEX(ys,i,1),0,0)
        ))),
        INDEX(ys,n,1)-INDEX(ys,n-1,1),
        SWITCH(cond_start,"free",0,"not",0,cond_start*INDEX(h,1,1)),
        SWITCH(cond_end,"free",0,"not",0,cond_end*INDEX(h,n-1,1))
    ),
    coeffs,HSTACK(WRAPROWS(MMULT(MINVERSE(amatrix),bmatrix),3)/(h^{3,2,1}),DROP(ys,-1,0)),
    pp,MAP(HSTACK(xs,coeffs),LAMBDA(cell,IFERROR(cell,""))),
    res,IF(debug=TRUE,HSTACK(amatrix,bmatrix),
        IF(ISOMITTED(x),pp,
            PPVAL(pp,x)
        )
    ),
    IF(AND(n=3,cond_start="not",cond_end="not"),"Error: redundant knot constraint",res)
)));

/******************** 
* POLYNOMIAL FORMULAS
*********************/

/**
* Returns the coefficients for the nth-degree polynomial fit using LINEST
*/
POLYFIT = LAMBDA(known_xs,known_ys,n,
LET(doc,"https://www.vertex42.com/lambda/polyfit.html",
    LINEST(known_ys,known_xs^SEQUENCE(1,n))
));
/**
* Evaluates a polynomial defined by a row vector of constant coefficients for each value of x.
*/
POLYVAL = LAMBDA(coeffs,x,
LET(doc,"https://www.vertex42.com/lambda/polyval.html",
    coeffs,IF(ROWS(coeffs)>1,TRANSPOSE(coeffs),coeffs),
    n,COLUMNS(coeffs)-1,
    IF(OR(ISTEXT(x)),
        LET(
            matrix,MAKEARRAY(ROWS(x),n+1,LAMBDA(i,j,
                SWITCH(n-(j-1),
                0,INDEX(coeffs,1,j),
                1,IMPRODUCT(INDEX(coeffs,1,j),INDEX(x,i)),
                2,IMPRODUCT(INDEX(coeffs,1,j),INDEX(x,i),INDEX(x,i)),
                IMPRODUCT(INDEX(coeffs,1,j),IMPOWER(INDEX(x,i),n-(j-1)))
                )
            )),
            BYROW(matrix,LAMBDA(row,LET(sum,IMSUM(row),IF(IMAGINARY(sum)=0,IMREAL(sum),sum))))
        ),
        LET(
            X_mat,IF(n=0,SEQUENCE(ROWS(x),1,1,0),
                HSTACK(x^SEQUENCE(1,n,n,-1),SEQUENCE(ROWS(x),1,1,0))
            ),
            BYROW(coeffs*X_mat,LAMBDA(row,SUM(row)))
        )
    )
));
/**
* Returns the derivative of the polynomial defined by coefficients in descending order of power.
*/
POLYDER = LAMBDA(coeffs,
LET(doc,"https://www.vertex42.com/lambda/polyder.html",
    n,COLUMNS(coeffs)-1,
    powers,SEQUENCE(1,n+1,n,-1),
    CHOOSECOLS(coeffs*powers,IF(n=0,1,SEQUENCE(n)))
));
/**
* Returns the integral of the polynomial with integration constant k, optionally evaluated from a to b
*/
POLYINT = LAMBDA(coeffs,[k],[a],[b],
LET(doc,"https://www.vertex42.com/lambda/polyint.html",
    k,IF(ISBLANK(k),SEQUENCE(ROWS(coeffs),1,0,0),IF(ROWS(k)=1,SEQUENCE(ROWS(coeffs),1,k,0),k)),
    n,COLUMNS(coeffs),
    powers,SEQUENCE(1,n,n,-1),
    qx,HSTACK(coeffs/powers,k),
    IF(AND(NOT(ISBLANK(a)),NOT(ISBLANK(b))),
        BYROW(qx,LAMBDA(row,POLYVAL(row,b)-POLYVAL(row,a))),
        qx
    )
));
/**
* Add two polynomials a + b = c
*/
POLYADD = LAMBDA(a,b,[drop_leading_zeros],
LET(doc,"https://www.vertex42.com/lambda/polyadd.html",
    na,COLUMNS(a),
    nb,COLUMNS(b),
    c,IF(na=nb,a+b,
    IF(na>nb,a+HSTACK(SEQUENCE(1,na-nb,0,0),b),
        b+HSTACK(SEQUENCE(1,nb-na,0,0),a)
    )),
    IF(drop_leading_zeros=TRUE,DROP(c,0,MATCH(TRUE,c<>0,0)-1),c)
));
/**
* Subtract polynomial b from polynomial a: a - b = c
*/
POLYSUB = LAMBDA(a,b,[drop_leading_zeros],
LET(doc,"https://www.vertex42.com/lambda/polyadd.html",
    POLYADD(a,-b,drop_leading_zeros)
));
/**
* Multiply two polynomials a and b, represented as row vectors of coefficients
* Uses the convolution algorithm for two vectors
*/
POLYMULT = LAMBDA(a,b,
LET(doc,"https://www.vertex42.com/lambda/polymult.html",
    a,IF(AND(ROWS(a)>1,COLUMNS(a)=1),TRANSPOSE(a),a),
    b,IF(AND(ROWS(b)>1,COLUMNS(b)=1),TRANSPOSE(b),b),
    len_a,COLUMNS(a),len_b,COLUMNS(b),len_c,len_a+len_b-1,
    coeffs,MAKEARRAY(1,len_c,LAMBDA(i,k,
        REDUCE(0,SEQUENCE(MIN(k,len_a)-MAX(1,k+1-len_b)+1,1,MAX(1,k+1-len_b)),
            LAMBDA(acc,j,
                IF(OR(ISTEXT(acc),ISTEXT(INDEX(a,j)),ISTEXT(INDEX(b,k-j+1))),
                    IMSUM(acc,IMPRODUCT(INDEX(a,j),INDEX(b,k-j+1))),
                    acc + INDEX(a,j)*INDEX(b,k-j+1)
                )
            )
        )
    )),
    MAP(coeffs,LAMBDA(val,IF(IMAGINARY(val)=0,IMREAL(val),val)))
));
/**
* Divide two polynomials using long division: a/b = q remainder r
* Returns polynomial q stacked vertically on top of polynomial r
*/
POLYDIV = LAMBDA(a,b,
LET(doc,"https://www.vertex42.com/lambda/polydiv.html",
    na,COLUMNS(a),
    nb,COLUMNS(b),
    result,IF(na<nb,VSTACK(0*b,b),
    REDUCE(VSTACK(0*a,a),SEQUENCE(na-nb+1,1,na-nb,-1),LAMBDA(acc,i,LET(
        dividend,CHOOSEROWS(acc,2),
        t,INDEX(dividend,1,na-nb-i+1)/INDEX(b,1,1),
        tv,IF(i>0,HSTACK(t,SEQUENCE(1,i,0,0)),t),
        q,IF(i>0,HSTACK(t*b,SEQUENCE(1,i,0,0)),t*b),
        r,POLYADD(dividend,-q),
        VSTACK(POLYADD(INDEX(acc,1,0),tv),r)
    )))),
    IF(na>=nb,DROP(result,0,nb-1),result)
));
/**
* Returns the companion matrix for the polynomial defined by coeffs
*/
POLYCOMPAN = LAMBDA(coeffs,
LET(doc,"https://www.vertex42.com/lambda/polycompan.html",
    coeffs,IF(AND(ROWS(coeffs)>1,COLUMNS(coeffs)=1),TRANSPOSE(coeffs),coeffs),
    n,COLUMNS(coeffs)-1,
    m,-CHOOSECOLS(coeffs,SEQUENCE(1,n,2,1))/INDEX(coeffs,1),
    VSTACK(m,HSTACK(MUNIT(n-1),SEQUENCE(n-1,1,0,0)))
));
/**
* Evaluate a Piecewise Polynomial at each value of x
*/
PPVAL = LAMBDA(pp_array,xvec,
LET(doc,"https://www.vertex42.com/lambda/ppval.html",
    breaks,INDEX(pp_array,0,1),
    pieces,ROWS(breaks)-1,
    coeffs,DROP(pp_array,-1,1),
    MAP(xvec,LAMBDA(x,IF(OR(x<MIN(breaks),x>MAX(breaks)),NA(),
        LET(
            row,IF(x=INDEX(breaks,pieces+1,1),pieces,MATCH(x,breaks,1)),
            POLYVAL(INDEX(coeffs,row,0),x-INDEX(breaks,row))
        )
    )))
));
/**
* Definite Integral of a Piecewise Polynomial from a to b
*/
PPINT = LAMBDA(pp_array,[a],[b],[cumulative],
LET(doc,"https://www.vertex42.com/lambda/ppint.html",
    breaks,INDEX(pp_array,0,1),
    pieces,ROWS(breaks)-1,
    coeffs,DROP(pp_array,-1,1),
    a,IF(OR(ISBLANK(a),a<INDEX(breaks,1)),INDEX(breaks,1),a),
    b,IF(OR(ISBLANK(b),b>INDEX(breaks,pieces+1)),INDEX(breaks,pieces+1),b),
    areas,SCAN(0,SEQUENCE(pieces+1),LAMBDA(acc,i,
    IF(i=1,0,LET(
        x_1,INDEX(breaks,i-1),
        x_2,INDEX(breaks,i),
        IF(NOT(AND(a<INDEX(breaks,i),b>INDEX(breaks,i-1))),acc+0,
            acc+POLYINT(CHOOSEROWS(coeffs,i-1),,MAX(0,a-x_1),MIN(b,x_2)-x_1)
        )
    )))),
    IF(cumulative=TRUE,areas,INDEX(areas,pieces+1,1))
));
/**
* Derivative of a Piecewise Polynomial
*/
PPDER = LAMBDA(pp_array,
LET(doc,"https://www.vertex42.com/lambda/ppder.html",
    breaks,INDEX(pp_array,0,1),
    pieces,ROWS(breaks)-1,
    coeffs,DROP(pp_array,-1,1),
    newcoeffs,REDUCE(0,SEQUENCE(pieces),LAMBDA(acc,i,LET(
        deriv,POLYDER(CHOOSEROWS(coeffs,i)),
        IF(i=1,deriv,VSTACK(acc,deriv))
    ))),
    MAP(HSTACK(breaks,newcoeffs),LAMBDA(cell,IFERROR(cell,"")))
));
/**
* Find the roots of a polynomial using the Weierstrauss method
*/
POLYROOTS = LAMBDA(coeffs,[epsilon],[iterations],[xstart],[return_abs],
LET(doc,"https://www.vertex42.com/lambda/polyroots.html",
    iterations,IF(ISBLANK(iterations),42,iterations),
    xstart,IF(ISBLANK(xstart),"0.4+0.9i",xstart),
    return_abs,IF(ISBLANK(return_abs),FALSE,return_abs),
    n,COLUMNS(coeffs)-1,
    xo,TRANSPOSE(IMPOWER(xstart,SEQUENCE(1,n,0))),
    result,REDUCE(xo,SEQUENCE(iterations),LAMBDA(acc,k,
    IF(AND(NOT(ISBLANK(epsilon)),k>1,IFERROR(MAX(INDEX(acc,0,3))<epsilon,FALSE)),acc,
        LET(
            fo,IF(k>1,INDEX(acc,0,2),POLYVAL(coeffs,INDEX(acc,0,1))),
            xo,INDEX(acc,0,1),
            xn,MAKEARRAY(n,1,LAMBDA(i,j,
                IMSUB(INDEX(xo,i),IMDIV(INDEX(fo,i),IMPRODUCT(IMSUB(INDEX(xo,i),DROP(CIRCSHIFT(xo,1-i,1),1)))))
            )),
            fn,POLYVAL(coeffs,xn),
            fabs,BYROW(fn,LAMBDA(value,IMABS(value))),
            HSTACK(xn,fn,fabs)
        )
    ))),
    IF(return_abs=TRUE,CHOOSECOLS(result,{1,3}),INDEX(result,0,1))
));
/**
* Creates a monic polynomial from a vector of roots
*/
POLYFROMROOTS = LAMBDA(roots,
LET(doc,"https://www.vertex42.com/lambda/polyfromroots.html",
    roots,IF(AND(ROWS(roots)=1,COLUMNS(roots)>1),TRANSPOSE(roots),roots),
    poly,REDUCE(1,SEQUENCE(ROWS(roots)-1),LAMBDA(acc,i,
        IF(i=1,
            IF(OR( ISTEXT(INDEX(roots,1)), ISTEXT(INDEX(roots,2)) ),
                POLYMULT(HSTACK(1,IMPRODUCT(-1,INDEX(roots,1))),HSTACK(1,IMPRODUCT(-1,INDEX(roots,2)))),
                POLYMULT(HSTACK(1,-INDEX(roots,1)),HSTACK(1,-INDEX(roots,2)))
            ),
            IF(OR(SUM(--ISTEXT(acc))>0,ISTEXT(INDEX(roots,i+1))),
                POLYMULT(acc,HSTACK(1,IMPRODUCT(-1,INDEX(roots,i+1)))),
                POLYMULT(acc,HSTACK(1,-INDEX(roots,i+1)))
            )
        )
    )),
    realpoly,MAP(poly,LAMBDA(inum,IMREAL(inum))),
    realpoly
));


/* **********************
* EQUATION SOLVERS
********************** */

/**
* Newton-Raphson Method equation solver for finding f(x)=y
*/
SOLVE_NR = LAMBDA(function,xstart,[y],[xstep],[epsilon],[iterations],[info],
LET(doc,"https://www.vertex42.com/lambda/newton-raphson.html",
    n,ROWS(xstart),
    y,IF(ISOMITTED(y),0,y),
    xstep,IF(ISOMITTED(xstep),0.1,xstep),
    epsilon,IF(ISOMITTED(epsilon),1E-14,epsilon),
    iterations,IF(ISOMITTED(iterations),20,iterations),
    ret,REDUCE({"x","Fx","Iterations"},SEQUENCE(n),LAMBDA(acc_row,i,
        VSTACK(acc_row,
        REDUCE(
            HSTACK(INDEX(xstart,i),function(INDEX(xstart,i))-y,0),
            SEQUENCE(iterations),
            LAMBDA(acc,k,
            LET(cur_row,TAKE(acc,-1),
                IF(ABS(INDEX(cur_row,1,2))<ABS(epsilon),acc,
                    LET(
                        xn,INDEX(cur_row,1,1),
                        fn,INDEX(cur_row,1,2),
                        fprimen,(function(xn+xstep)-y-fn)/xstep,
                        xnew,xn-fn/fprimen,
                        fnew,function(xnew)-y,
                        new_row,HSTACK(xnew,fnew,k),
                        output,IF(info="verbose",VSTACK(acc,new_row),new_row),
                        output
                    )
                )
            ))
        ))
    )),
    IF(OR(info=TRUE,info="verbose"),ret,DROP(ret,1,-2))
));
/**
* Quadratic Formula returns the real or complex Roots of ax^2+bx+c=0
*/
QUADRATIC = LAMBDA(a,[b],[c],
LET(doc,"https://www.vertex42.com/lambda/quadratic.html",
    b,IF(ISOMITTED(b),INDEX(a,2),b),
    c,IF(ISOMITTED(c),INDEX(a,3),c),
    a,IF(COLUMNS(a)>1,INDEX(a,1),a),
    vertex,-b/(2*a),
    radicand,(vertex)^2-c/a,
    pm,SQRT(ABS(radicand)),
    rmax,vertex+pm,
    rmin,vertex-pm,
    IF(radicand>=0,
        IF(radicand=0,VSTACK(rmax,rmax),
            IF(vertex>=0,
                VSTACK(rmax,(c/a)/(rmax)),
                VSTACK((c/a)/rmin,rmin)
            )
        ),
        VSTACK(COMPLEX(vertex,pm),COMPLEX(vertex,-pm))
    )
));

/************************
*  UTILITY FUNCTIONS
*************************/
/**
* Runs a function n times in sequence using BYROW and returns the time in milliseconds
* Resolution is about 10 milliseconds
*/
TIMER = LAMBDA(n_iterations,function,[ret_values],[j],[k],
LET(doc,"https://www.vertex42.com/lambda/timer.html",
    started,NOW(),
    seq,SEQUENCE(n_iterations),
    result,BYROW(seq,LAMBDA(i,INDEX(function(i),IF(ISBLANK(j),1,j),IF(ISBLANK(k),1,k)))),
    time_in_ms,(NOW()-started)*24*3600000,
    IF(ret_values=TRUE,result,time_in_ms)
));
/**
* Like BYROW, but can return a multi-column array instead of just a single column
*/
BYROW2D = LAMBDA(array, function,
LET(doc,"https://www.vertex42.com/lambda/byrow2d.html",
    REDUCE("",
        SEQUENCE(ROWS(array)),
        LAMBDA(acc,i,IF(i=1,
            function(CHOOSEROWS(array,i)),
            VSTACK(acc,function(CHOOSEROWS(array,i)))
        ))
    )
));
/**
* Like BYCOL, but can return a multi-row array instead of just a single row
*/
BYCOL2D = LAMBDA(array, function,
LET(doc,"https://www.vertex42.com/lambda/byrow2d.html",
    REDUCE("",
        SEQUENCE(COLUMNS(array)),
        LAMBDA(acc,i,IF(i=1,
            function(CHOOSECOLS(array,i)),
            HSTACK(acc,function(CHOOSECOLS(array,i)))
        ))
    )
));
/**
* HSTACK the result of expression for each value in a vector
* e.g. HSTACK_ForEach(vector, LAMBDA(val, SEQUENCE(val) ) ))
*/
HSTACK_ForEach = LAMBDA(vector, lambda_function,
LET(doc,"https://www.vertex42.com/lambda/byrow2d.html",
    DROP(REDUCE("", vector,
        LAMBDA(acc,val, HSTACK(acc, lambda_function(val) ) )
    ),,1)
));
/**
* VSTACK the result of expression for each value in a vector
* e.g. VSTACK_ForEach(vector, LAMBDA(val, SEQUENCE(val) ) ))
*/
VSTACK_ForEach = LAMBDA(vector, lambda_function,
LET(doc,"https://www.vertex42.com/lambda/byrow2d.html",
    DROP(REDUCE("", vector,
        LAMBDA(acc,val, VSTACK(acc, lambda_function(val) ) )
    ),1)
));

/* ******************************
* Functions for Structure Arrays
* ******************************/

/**
* Returns the field defined by the name
*/
STRUCT_GET = LAMBDA(name,structure,
LET(doc,"https://www.vertex42.com/lambda/structures.html",
    n,COLUMNS(name),
    getStruct,LAMBDA(id,struct_i,
        LET(
            row,MATCH(id,INDEX(struct_i,0,1),0),
            rows,INDEX(struct_i,row,2),
            cols,INDEX(struct_i,row,3),
            INDEX(struct_i,row+1,2):INDEX(struct_i,row+rows,1+cols)
        )
    ),
    IF(n=1,
        getStruct(name,structure),
        REDUCE(structure,SEQUENCE(n),LAMBDA(acc,i,
            getStruct(INDEX(name,i),acc)
        ))
    )
));

/**
* Creates a structure with up to 3 fields.
*/
STRUCT_CREATE = LAMBDA(name_1,array_1,[name_2],[array_2],[name_3],[array_3],[name_4],[array_4],
LET(doc,"https://www.vertex42.com/lambda/structures.html",
    one,VSTACK(HSTACK(name_1,ROWS(array_1),COLUMNS(array_1)),HSTACK("",array_1)),
    include_2,NOT(AND(ISOMITTED(name_2),ISOMITTED(array_2))),
    include_3,NOT(AND(ISOMITTED(name_3),ISOMITTED(array_3))),
    include_4,NOT(AND(ISOMITTED(name_4),ISOMITTED(array_4))),
    two,IF(include_2,
        IFERROR(VSTACK(HSTACK(name_2,ROWS(array_2),COLUMNS(array_2)),HSTACK("",array_2)),""),
        "Empty"
    ),
    three,IF(include_3,
        IFERROR(VSTACK(HSTACK(name_3,ROWS(array_3),COLUMNS(array_3)),HSTACK("",array_3)),""),
        "Empty"
    ),
    four,IF(include_4,
        IFERROR(VSTACK(HSTACK(name_4,ROWS(array_4),COLUMNS(array_4)),HSTACK("",array_4)),""),
        "Empty"
    ),
    ret,IF(NOT(include_2),
        one,
        IF(NOT(include_3),
            VSTACK(one,two),
        IF(NOT(include_4),
            VSTACK(one,two,three),
            VSTACK(one,two,three,four)
        ))
    ),
    MAP(ret,LAMBDA(cell,IF(ISNA(cell),"",cell)))
));

/** 
* Appends the new name and array to the existing structure.
*/
STRUCT_APPEND = LAMBDA(structure,name,array,
LET(doc,"https://www.vertex42.com/lambda/structures.html",
    ret,VSTACK(structure,
        HSTACK(name,ROWS(array),COLUMNS(array)),
        HSTACK("",array)
    ),
    MAP(ret,LAMBDA(cell,IF(ISNA(cell),"",cell)))
));

/**
* Removes a field by name from an existing structure.
*/
STRUCT_REMOVE = LAMBDA(name,structure,
LET(doc,"https://www.vertex42.com/lambda/structures.html",
    row,MATCH(name,INDEX(structure,0,1),0),
    rows,INDEX(structure,row,2),
    filter_rows,VSTACK(SEQUENCE(row-1,1,1,0),
        SEQUENCE(rows+1,1,0,0),
        SEQUENCE(ROWS(structure)-(row+rows),1,1,0)
    ),
    FILTER(structure,filter_rows)
));

/************************
*  IMAGINARY NUMBERS
*************************/
/**
* Returns the Conjugate Transpose of a Complex Matrix
* Equal to TRANSPOSE if values are Real
*/
IM_CTRANSPOSE = LAMBDA(matrix,
LET(doc,"https://www.vertex42.com/lambda/complex.html",
    TRANSPOSE(MAP(matrix,LAMBDA(cell,IMCONJUGATE(cell))))
));
/**
* Returns the Complex Dot Product of two vectors
* a={"1+i";"-i"}, b={"3+i";"1-4i"}
* Equal to SUMPRODUCT if values are Real, but not IM_SUMPRODUCT
*/
IM_DOT = LAMBDA(a_vec,b_vec,
LET(doc,"https://www.vertex42.com/lambda/complex.html",
    IMSUM(MAP(a_vec,b_vec,LAMBDA(ai,bi,IMPRODUCT(ai,IMCONJUGATE(bi)))))
));
/**
* Returns the sum of the element-wise product of complex vectors a and b
* Equal to SUMPRODUCT if values are Real, but not IM_DOT
*/
IM_SUMPRODUCT = LAMBDA(a_vec,b_vec,
LET(doc,"https://www.vertex42.com/lambda/complex.html",
    IMSUM(MAP(a_vec,b_vec,LAMBDA(ai,bi,IMPRODUCT(ai,bi))))
));
/**
* Returns the matrix multiplication of two complex matrices
* Equal to MMULT if values are Real
*/
IM_MMULT = LAMBDA(a_mat,b_mat,
LET(doc,"https://www.vertex42.com/lambda/complex.html",
    m,ROWS(a_mat),
    n,COLUMNS(b_mat),
    IF(COLUMNS(a_mat)<>ROWS(b_mat),"Error: invalid sizes",
        MAKEARRAY(m,n,LAMBDA(i,j,
            IM_SUMPRODUCT(TRANSPOSE(INDEX(a_mat,i,0)),INDEX(b_mat,0,j))
        ))
    )
));
/**
* Checks if there is at least one values in an array that is a complex number
* and if all values can be treated as complex numbers
*/
IM_ISCOMPLEX = LAMBDA(matrix,
LET(doc,"https://www.vertex42.com/lambda/complex.html",
    istext,IF(SUM(--ISTEXT(matrix))>0,TRUE,FALSE),
    isinum,IF(SUM(--ISERROR(MAP(matrix,LAMBDA(cell,IMSUM(cell,cell)))))>0,FALSE,TRUE),
    AND(istext,isinum)
));
/**
* Returns TRUE if the matrix is equal to its conjugate transpose
*/
IM_ISHERMITIAN = LAMBDA(matrix,
LET(doc,"https://www.vertex42.com/lambda/complex.html",
    matrix,IM_CTRANSPOSE(IM_CTRANSPOSE(matrix)),
    IF(ROWS(matrix)<>COLUMNS(matrix),FALSE,SUM(--(matrix<>IM_CTRANSPOSE(matrix)))=0)
));
/**
* If |coefficient| is less than epsilon (default=1E-14), make it equal zero
*/
IM_ROUNDTOZERO = LAMBDA(inumber,[epsilon],
LET(doc,"https://www.vertex42.com/lambda/complex.html",
    epsilon,IF(ISOMITTED(epsilon),1E-14,epsilon),
    MAP(inumber,LAMBDA(inum,
    LET(real,IF(ABS(IMREAL(inum))<epsilon,0,IMREAL(inum)),
        imag,IF(ABS(IMAGINARY(inum))<epsilon,0,IMAGINARY(inum)),
        IF(imag=0,real,COMPLEX(real,imag))
    )))
));
/**
* Rounds the Real and Imaginary coefficients separately to num_digits
*/
IM_ROUND = LAMBDA(inumber,num_digits,
LET(doc,"https://www.vertex42.com/lambda/complex.html",
    num_digits,IF(ISOMITTED(num_digits),{"",""},num_digits),
    rdigits,INDEX(num_digits,1),
    idigits,IF(OR(ROWS(num_digits)>1,COLUMNS(num_digits)>1),INDEX(num_digits,2),rdigits),
    MAP(inumber,LAMBDA(inum,
    LET(real,IF(ISBLANK(rdigits),IMREAL(inum),ROUND(IMREAL(inum),rdigits)),
        imag,IF(ISBLANK(idigits),IMAGINARY(inum),ROUND(IMAGINARY(inum),idigits)),
        IF(imag=0,real,COMPLEX(real,imag))
    )))
));
/**
* Converts a Complex Number z=x+yi to an array of coefficients {x,y}
*/
IM_TOARRAY = LAMBDA(inumber,
LET(doc,"https://www.vertex42.com/lambda/complex.html",
    MAKEARRAY(ROWS(inumber),2,LAMBDA(i,j,
        IF(j=1,IMREAL(INDEX(inumber,i)),IMAGINARY(INDEX(inumber,i)))
    ))
));

/**********************************
* FINITE DIFFERENCE and INTEGRATION
***********************************/

/**
* Calculate the difference between adjacent values of a column vector
*/
DIFF = LAMBDA(x,[n],
LET(doc,"https://www.vertex42.com/lambda/diff.html",
    n,IF(ISBLANK(n),1,n),
    REDUCE(1,SEQUENCE(n),LAMBDA(acc,i,
        IF(i=1,
            DROP(x,1)-DROP(x,-1),
            DROP(acc,1)-DROP(acc,-1)
        )
    ))
));
/**
* Calculate the n-th Finite Difference Derivative assuming uniform spacing h
*/
FDIFF = LAMBDA(ys,[h],[n],
LET(doc,"https://www.vertex42.com/lambda/fdiff.html",
    n,IF(ISBLANK(n),1,n),
    h,IF(ISBLANK(h),1,h),
    m,ROWS(ys),
    MAKEARRAY(m,1,LAMBDA(i,j,
        SWITCH(n,
            1,SWITCH(i,
                1,SUMPRODUCT({-3;4;-1},CHOOSEROWS(ys,1,2,3))/(2*h),
                m,SUMPRODUCT({1;-4;3},CHOOSEROWS(ys,m-2,m-1,m))/(2*h),
                SUMPRODUCT({-1;0;1},CHOOSEROWS(ys,i-1,i,i+1))/(2*h)
            ),
            2,SWITCH(i,
                1,SUMPRODUCT({2;-5;4;-1},CHOOSEROWS(ys,1,2,3,4))/(h^2),
                m,SUMPRODUCT({-1;4;-5;2},CHOOSEROWS(ys,m-3,m-2,m-1,m))/(h^2),
                SUMPRODUCT({1;-2;1},CHOOSEROWS(ys,i-1,i,i+1))/(h^2)
            ),
            3,SWITCH(i,
                1,SUMPRODUCT({-5;18;-24;14;-3},CHOOSEROWS(ys,1,2,3,4,5))/(2*h^3),
                2,SUMPRODUCT({-3;10;-12;6;-1},CHOOSEROWS(ys,1,2,3,4,5))/(2*h^3),
                m-1,SUMPRODUCT({1;-6;12;-10;3},CHOOSEROWS(ys,m-4,m-3,m-2,m-1,m))/(2*h^3),
                m,SUMPRODUCT({3;-14;24;-18;5},CHOOSEROWS(ys,m-4,m-3,m-2,m-1,m))/(2*h^3),
                SUMPRODUCT({-1;2;0;-2;1},CHOOSEROWS(ys,i-2,i-1,i,i+1,i+2))/(2*h^3)
            ),
            4,SWITCH(i,
                1,SUMPRODUCT({3;-14;26;-24;11;-2},CHOOSEROWS(ys,1,2,3,4,5,6))/(h^4),
                2,SUMPRODUCT({2;-9;16;-14;6;-1},CHOOSEROWS(ys,1,2,3,4,5,6))/(h^4),
                m-1,SUMPRODUCT({-1;6;-14;16;-9;2},CHOOSEROWS(ys,m-5,m-4,m-3,m-2,m-1,m))/(h^4),
                m,SUMPRODUCT({-2;11;-24;26;-14;3},CHOOSEROWS(ys,m-5,m-4,m-3,m-2,m-1,m))/(h^4),
                SUMPRODUCT({1;-4;6;-4;1},CHOOSEROWS(ys,i-2,i-1,i,i+1,i+2))/(h^4)
            ),
            "Error: Invalid Integer n, 1<=n<=4"
        )
    ))
));
/**
* Calculate the nth-Derivative of y(x) using a polynomial-based Finite Difference approximation
* Default: Finds the 1st derivative at 2nd-order accuracy
*/
PDIFF = LAMBDA(xs,ys,[nth_deriv],[order],
LET(doc,"https://www.vertex42.com/lambda/pdiff.html",
    n,IF(ISBLANK(nth_deriv),1,nth_deriv),
    order,IF(ISBLANK(order),2,order),
    BYROW(xs,LAMBDA(x,
        LET(tab,TAKE(SORT(HSTACK((xs-x)^2,xs,ys)),n+order),
            xo,INDEX(tab,1,2),
            poly,POLYFIT(INDEX(tab,0,2)-xo,INDEX(tab,0,3),n+order-1),
            deriv,REDUCE(poly,SEQUENCE(n),LAMBDA(acc,i,POLYDER(acc))),
            POLYVAL(deriv,x-xo)
        )
    ))
));
/**
* Numerical Integration using the Trapezoidal Rule
* Defaults: dx=1 if x and dx are blank, cumulative=FALSE
*/
TRAPZ = LAMBDA(y,[x],[dx],[cumulative],
LET(doc,"https://www.vertex42.com/lambda/trapz.html",
    dx,IF(ISBLANK(dx),1,dx),
    areas,IF(ISOMITTED(x),
        (DROP(y,-1)+DROP(y,1))/2*dx,
        ((DROP(y,-1)+DROP(y,1))/2)*(DROP(x,1)-DROP(x,-1))
    ),
    REDUCE("",SEQUENCE(COLUMNS(areas)),
        LAMBDA(acc,i,LET(
            col,CHOOSECOLS(areas,i),
            new,IF(cumulative=TRUE,
                SCAN(0,SEQUENCE(ROWS(col)+1),
                    LAMBDA(acc,i,IF(i=1,0,acc+INDEX(col,i-1,1)))
                ),
                SUM(col)
            ),
            IF(i=1,new,HSTACK(acc,new))
        ))
    )
));
/**
* Numerical Integration using Composite Simpson's 1/3 Rule
*/
SIMPSON = LAMBDA(y,[x],[dx],
LET(doc,"https://www.vertex42.com/lambda/simpson.html",
    dx,IF(ISOMITTED(dx),1,dx),
    n,ROWS(y),
    integral,REDUCE(0,SEQUENCE(INT((n-1)/2),1,1,2),LAMBDA(acc,i,
        IF(ISOMITTED(x),
            acc+(1/3*dx*(INDEX(y,i)+4*INDEX(y,i+1)+INDEX(y,i+2))),
            LET(
                h_1,INDEX(x,i+1)-INDEX(x,i),
                h_2,INDEX(x,i+2)-INDEX(x,i+1),
                area,(h_1+h_2)/6*(
                    (2-h_2/h_1)*INDEX(y,i) + 
                    (h_1+h_2)^2/(h_1*h_2)*INDEX(y,i+1) + 
                    (2-h_1/h_2)*INDEX(y,i+2)),
                acc+area
            )
        )
    )),
    lastinterval,IF(ISODD(n),0,
        IF(ISOMITTED(x),
            dx*(5/12*INDEX(y,n)+2/3*INDEX(y,n-1)-1/12*(INDEX(y,n-2))),
            LET(
                h_1,INDEX(x,n)-INDEX(x,n-1),
                h_2,INDEX(x,n-1)-INDEX(x,n-2),
                alpha,(2*h_1^2+3*h_1*h_2)/(6*(h_2+h_1)),
                beta,(h_1^2+3*h_1*h_2)/(6*h_2),
                nu,(h_1^3)/(6*h_2*(h_2+h_1)),
                alpha*INDEX(y,n)+beta*INDEX(y,n-1)-nu*INDEX(y,n-2)
            )
        )
    ),
    integral+lastinterval
));

/* *******************
* ADDRESS FUNCTIONS
*********************/

/**
* Shortcut version of the ADDRESS function. Uses cell_or_range as the input.
*/
/*
* Inputs:
*   cell_or_range: The reference cell or range for which to return the address(es).
*   [abs_num]: Optional. Determines the address style (e.g., 1=absolute, 4=relative). Defaults to 1.
*   [a1_r1c1]: Optional. Determines the reference style (1=A1, 0=R1C1). Defaults to 1.
*   [include_sheetname]: Optional. If TRUE, includes the sheet name. Defaults to FALSE.
*
* Outputs:
*   Returns the address of the specified cell, optionally including the sheet name.
*/
AD = LAMBDA(cell_or_range, [abs_num], [a1_r1c1], [include_sheetname],
LET(doc, "https://www.vertex42.com/lambda/",
    version, "1.0.1 - Updated parameter name and commenting",
    cell, cell_or_range,
    abs_num, IF(ISOMITTED(abs_num), 1, abs_num),
    a1_r1c1, IF(ISOMITTED(a1_r1c1), 1, a1_r1c1),
    include_sheetname, IF(ISOMITTED(include_sheetname), FALSE, include_sheetname),
    IF(include_sheetname,
        ADDRESS(ROW(cell), COLUMN(cell), abs_num, a1_r1c1, TEXTAFTER(CELL("filename", cell), "]")),
        ADDRESS(ROW(cell), COLUMN(cell), abs_num, a1_r1c1)
    )
));