jonwittwer/lambdalibrary Secret

## lambdalibrary
/*  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, L_RESCALE would become VLL.L_RESCALE

Here are a few examples of module names you could use that will not conflict
with parameter names: VLL, LIB, LL, MY,
*/

/*
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
* L_SFROUND(1234500,4,-1) = 1234000
*/
L_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),
        )
    ))
));

/**
* Create an S-Curve from Project Task Dates and Values
*/
L_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))
));


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

/**
* Create a sequence of integers or characters between start and end
*/
L_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
* L_LINSPACE(2,3,5) = {2; 2.25; 2.5; 2.75; 3}
*/
L_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
* L_LOGSPACE(1,2,5) = {10; 17.8; 31.6; 56.2; 100}
*/
L_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]
* L_RESCALE({1,2,3},0,1) = {0, 0.5, 1}
*/
L_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)}
*/
L_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.
*/
L_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
*/
L_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)
        )
    )
)));


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

/**
* Removes all of the characters in chars (individually) from text
*/
L_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
*/
L_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
*/
L_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 single text string to a row (dim=1, default) or column (dim=2) of characters
*/
L_TEXT2ARRAY = LAMBDA(text,[dim],
LET(doc,"https://www.vertex42.com/lambda/text2array.html",
    dim,IF(ISOMITTED(dim),1,dim),
    array,MID(text,SEQUENCE(1,LEN(text)),1),
    IF(dim=2,TRANSPOSE(array),array)
));
/**
* Converts a Roman numeral to its integer value: IV = 4, XII = 12, etc.
*/
L_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
*/
L_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
*/
L_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
*/
L_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
*/
L_ROT90 = 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.
*/
L_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
*/
L_REPARRAY = LAMBDA(array,m_vert,[n_horiz],
LET(doc,"https://www.vertex42.com/lambda/reparray.html",
    n_vert,IF(ISBLANK(m_vert),n_horiz,m_vert),
    m_horiz,IF(ISBLANK(n_horiz),n_vert,n_horiz),
    MAKEARRAY(n_vert*ROWS(array),m_horiz*COLUMNS(array),LAMBDA(i,j,
        INDEX(array,1+MOD(i-1,ROWS(array)),1+MOD(j-1,COLUMNS(array)))
    ))
));
/**
* Repeat the elements of an array m times vertically and n times horizontally
*/
L_REPELEM = LAMBDA(array,m_vert,[n_horiz],
LET(doc,"https://www.vertex42.com/lambda/repelem.html",
    n_vert,IF(ISBLANK(m_vert),n_horiz,m_vert),
    m_horiz,IF(ISBLANK(n_horiz),n_vert,n_horiz),
    MAKEARRAY(n_vert*ROWS(array),m_horiz*COLUMNS(array),LAMBDA(i,j,
        INDEX(array,ROUNDUP(i/n_vert,0),ROUNDUP(j/m_horiz,0))
    ))
));
/**
* Replace a block within an array by specifying the starting (i,j) location
*/
L_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
*/
L_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
*/
L_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 two arrays
*/
L_COMBINATIONS = LAMBDA(array_1,array_2,
LET(doc,"https://www.vertex42.com/lambda/combinations.html",
    r_1,ROWS(array_1),r_2,ROWS(array_2),
    first,L_REPELEM(array_1,r_2,1),
    second,L_REPARRAY(array_2,r_1,1),
    HSTACK(first,second)
));
/**
* Return all the combinations of (N Choose K) from an Nx1 array
*/
L_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)
    ))))
));
/**
* Replace a single element of array at location row x col with value
*/
L_REPLACEBYLOC = LAMBDA(array,value,row,col,
LET(doc,"https://www.vertex42.com/lambda/replacebyloc.html",
    MAKEARRAY(ROWS(array),COLUMNS(array),
        LAMBDA(r,c,IF(AND(r=row,c=col),value,INDEX(array,r,c)))
    )
));


/**
* Returns the vector {num_rows;num_columns}
*/
L_SIZE = LAMBDA(array,
    VSTACK(ROWS(array),COLUMNS(array))
);
/**
* Returns TRUE if one of the dimensions is 1
*/
L_ISVECTOR = LAMBDA(array,
    OR(ROWS(array)=1,COLUMNS(array)=1)
);
/**
* Returns TRUE if columns>1 and rows>1
*/
L_ISARRAY = LAMBDA(array,
    OR(COLUMNS(array)>1,ROWS(array)>1)
);
/**
* Returns TRUE if columns=1 and rows>1
*/
L_ISCOLUMN = LAMBDA(array,
    AND(ROWS(array)>1,COLUMNS(array)=1)
);
/**
* Returns TRUE if columns>1 and rows=1
*/
L_ISROW = LAMBDA(array,
    AND(COLUMNS(array)>1,ROWS(array)=1)
);
/**
* Returns TRUE if the matrix is symmetric
*/
L_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
*/
L_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)
    ))
));

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

/**
* L_ONES(array) :: Returns a matrix of 1s the size of array
* L_ONES(m,[n]) :: Returns a matrix of 1s the size of m x n (default n=1)
*/
L_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)
        )
    )
));
/**
* L_ZEROS(array) :: Returns a matrix of 0s the size of array
* L_ZEROS(m,[n]) :: Returns a matrix of 0s the size of m x n (default n=1)
*/
L_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
*/
L_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
*/
L_TRACE = LAMBDA(matrix,
LET(doc,"https://www.vertex42.com/lambda/diag.html",
    SUM( matrix*MUNIT(ROWS(matrix)))
));
/**
* Returns a column vector containing row sums
*/
L_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
*/
L_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
*/
L_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.
*/
L_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
*/
L_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(L_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.
*/
L_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.
*/
L_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
* L_QR(matrix) returns {Q;R} (stacked vertically) where matrix=MMULT(Q,R).
*/
L_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
*/
L_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,L_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)),
            QR,L_QR(new_A-shift*MUNIT(n)),
            Q,TAKE(QR,n),
            R,TAKE(QR,-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
*/
L_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
*/
L_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
*/
L_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),
            L_POLYVAL(L_POLYFIT(INDEX(tab,0,2)-xo,INDEX(tab,0,3),n),x-xo)
        )
    ))
));
/**
* Cubic Spline - Creates a cubic piecewise polynomial by specifying control points and slopes at each point
*/
L_CSPLINE = LAMBDA(known_xs,known_ys,[x],[ms],[c],
LET(doc,"https://www.vertex42.com/lambda/cspline.html",
    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),
    slope,(1-c)*IF(OR(ms="",ROWS(ms)<=2),
        IF(AND(NOT(ms=""),ROWS(ms)=1),SEQUENCE(n,1,ms,0),
        MAKEARRAY(n,1,LAMBDA(i,j,
        IF(i=1,IF(ROWS(ms)=2,INDEX(ms,1),(INDEX(ys,i+1)-INDEX(ys,i))/(INDEX(xs,i+1)-INDEX(xs,i))),
        IF(i=n,IF(ROWS(ms)=2,INDEX(ms,2),(INDEX(ys,i)-INDEX(ys,i-1))/(INDEX(xs,i)-INDEX(xs,i-1))),
            1/2*( (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),
    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))
    ))),
    pp,MAP(HSTACK(xs,pp_coeffs),LAMBDA(cell,IFERROR(cell,""))),
    IF(ISOMITTED(x),pp,
        L_PPVAL(pp,x)
    )
));
/**
* Natural Cubic Spline - Creates a C2 Interpolating Spline with d2y/dx2=0 at the end points
*/
L_NSPLINE = LAMBDA(known_xs,known_ys,[x],
LET(doc,"https://www.vertex42.com/lambda/nspline.html",
    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),
    n,ROWS(xs),
    alpha,3/DROP(h,1,0)*(DROP(ys,2,0)-DROP(DROP(ys,1,0),-1,0))-3/DROP(h,-1,0)*(DROP(DROP(ys,1,0),-1,0)-DROP(ys,-2,0)),
    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),
    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,""))),
    m_1,L_POLYVAL(L_POLYDER(DROP(CHOOSEROWS(pp,1),0,1)),0),
    b_1,INDEX(known_ys,1,1)-m_1*INDEX(known_xs,1,1),
    m_n,L_POLYVAL(L_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),
    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(L_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.
*/
L_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,
            L_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
*/
L_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.
*/
L_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.
*/
L_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
*/
L_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,L_POLYVAL(row,b)-L_POLYVAL(row,a))),
        qx
    )
));
/**
* Add two polynomials a + b = c
*/
L_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
*/
L_POLYSUB = LAMBDA(a,b,[drop_leading_zeros],
LET(doc,"https://www.vertex42.com/lambda/polyadd.html",
    L_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
*/
L_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
*/
L_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,L_POLYADD(dividend,-q),
        VSTACK(L_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
*/
L_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
*/
L_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)),
            L_POLYVAL(INDEX(coeffs,row,0),x-INDEX(breaks,row))
        )
    )))
));
/**
* Definite Integral of a Piecewise Polynomial from a to b
*/
L_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+L_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
*/
L_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,L_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
*/
L_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),L_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(L_CIRCSHIFT(xo,1-i,1),1)))))
            )),
            fn,L_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
*/
L_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)) ),
                L_POLYMULT(HSTACK(1,IMPRODUCT(-1,INDEX(roots,1))),HSTACK(1,IMPRODUCT(-1,INDEX(roots,2)))),
                L_POLYMULT(HSTACK(1,-INDEX(roots,1)),HSTACK(1,-INDEX(roots,2)))
            ),
            IF(OR(SUM(--ISTEXT(acc))>0,ISTEXT(INDEX(roots,i+1))),
                L_POLYMULT(acc,HSTACK(1,IMPRODUCT(-1,INDEX(roots,i+1)))),
                L_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
*/
L_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
*/
L_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
*/
L_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
*/
L_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)))
        ))
    )
));

/* ******************************
* 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
*/
L_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
*/
L_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
*/
L_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,L_POLYFIT(INDEX(tab,0,2)-xo,INDEX(tab,0,3),n+order-1),
            deriv,REDUCE(poly,SEQUENCE(n),LAMBDA(acc,i,L_POLYDER(acc))),
            L_POLYVAL(deriv,x-xo)
        )
    ))
));
/**
* Numerical Integration using the Trapezoidal Rule
* Defaults: dx=1 if x and dx are blank, cumulative=FALSE
*/
L_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
*/
L_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
));