Single Formula Sudoku Solver

Replace 'puzzle1' with a 9x9 range containing an unsolved sudoku, and this formula will calculate the solution:

How it works: the puzzle is represented as an array of 81 numbers, where each bit indicates whether a candidate is still considered possible. e.g. if a cell could be 1,2, or 5, the number for that cell would be 10011 in binary. The constrain_idx function takes an index to this 81-number array, eliminates naked singles (and specific types of naked n-tuples) and identifies hidden singles, and returns an updated value for that cell. constrain_g is just a modified version of constrain_idx which can be applied to an entire grid (81-number candidate array). The constrain_g2 function calls this function recursively until no more eliminations can be made, or the grid is invalid. The ssolve function calls constrain_g2, then, if the solution does not result, starts trying candidates recursively, placing guesses in the cell with the least number of candidates > 1, until a solution is found.

=LET( input_grid, puzzle1, solve, LAMBDA(grid9x9, LET( all_opts, BIN2DEC("111111111"), seq9, SEQUENCE(9), seq81, SEQUENCE(81), opt_arr, 2^(seq9-1), row_ids, QUOTIENT(seq81-1,9)+1, col_ids, MOD(seq81-1,9)+1, box_ids, 3*QUOTIENT(row_ids-1, 3) + QUOTIENT(col_ids-1,3) + 1, row_mates, LAMBDA(idx,FILTER(seq81,--(row_ids=INDEX(row_ids,idx))*--(seq81<>idx))), col_mates, LAMBDA(idx,FILTER(seq81,--(col_ids=INDEX(col_ids,idx))*--(seq81<>idx))), box_mates, LAMBDA(idx,FILTER(seq81,--(box_ids=INDEX(box_ids,idx))*--(seq81<>idx))), all_mates, LAMBDA(idx,FILTER(seq81,--(seq81<>idx)*( --(row_ids=INDEX(row_ids,idx))+ --(col_ids=INDEX(col_ids,idx))+ --(box_ids=INDEX(box_ids,idx)) ))), popcount, LAMBDA(x,SUM(--(BITAND(opt_arr,x)<>0))), is_solved, LAMBDA(grid,SUM(--(MAP(grid,popcount)<>1))=0), is_error, LAMBDA(grid,SUM(--(grid=0))<>0), cn_funcs, LAMBDA(grid,MAP(VSTACK(row_mates,col_mates,box_mates),LAMBDA(f,LAMBDA(idx,INDEX(grid,f(idx)))))), constrain_idx, LAMBDA(grid,c_funcs,idx,LET( ival, INDEX(grid,idx), IF(popcount(ival)=1, ival, LET( mrvals, LAMBDA(vs,LET(uvs,UNIQUE(vs),FILTER(uvs,MAP(uvs,LAMBDA(v,AND(popcount(v)<=SUM(--(vs=v)),SUM(--(vs=v))>0))),0))), mmasks, MAP(c_funcs,LAMBDA(f,REDUCE(0,mrvals(f(idx)),LAMBDA(t,rv,BITOR(t,rv))))), mask, BITXOR(REDUCE(0,mmasks,LAMBDA(t,v,BITOR(t,v))),all_opts), bvs, MAP(c_funcs,LAMBDA(f,BITXOR(REDUCE(0,f(idx),LAMBDA(t,mv,BITOR(t,mv))),all_opts))), XLOOKUP(1,MAP(bvs,popcount),bvs,BITAND(INDEX(grid,idx),mask))*--(SUM(--(MAP(bvs,popcount)>1))=0) )) )), constrain_g, LAMBDA(grid,LET(c_funcs,cn_funcs(grid),MAP(seq81,LAMBDA(i,constrain_idx(grid,c_funcs,i))))), constrain_g2, LAMBDA(grid,LET( cgf, LAMBDA(f,gridt,LET( n_grid, constrain_g(gridt), is_same_err, OR(SUM(--(gridt<>n_grid))=0,is_error(n_grid)), IF(is_same_err,n_grid,f(f,n_grid)) )), cgf(cgf,grid) )), ssolve, LAMBDA(f,gridt,LET( n_grid, constrain_g2(gridt), IFS( is_solved(n_grid),n_grid, is_error(n_grid),FALSE, TRUE,LET( popcounts, MAP(n_grid,popcount), min_uc, SMALL(UNIQUE(popcounts),2), l_idx, XMATCH(min_uc,popcounts), l_val, INDEX(n_grid,l_idx), options, FILTER(opt_arr,BITAND(opt_arr,l_val)<>0), REDUCE(FALSE,options,LAMBDA(t,opt, IF(ROWS(t)=81,t, f(f,IF(seq81=l_idx,opt,n_grid)) ) )) ) ) )), print_sq, LAMBDA(v,TEXTJOIN(", ",TRUE,TRANSPOSE(FILTER(seq9,BITAND(2^(seq9-1),v)<>0)))), print_grid, LAMBDA(grid,WRAPROWS(MAP(grid,print_sq),9)), grid_, MAP(TOCOL(grid9x9),LAMBDA(x,IF(ISBLANK(x),all_opts,2^(x-1)))), print_grid(ssolve(ssolve,grid_)) )), solve(input_grid) )