Try this for light calcuation?

[curioustorvald]

* * *
* @ *
* * *

main:
	for all blocks in visible_world call calc_main.

subroutine calc_main(@):
	check all nearby: -- might improve performance over air tiles?
		if (their wall) == non-solid && (their air) == non-solid then skip calc_main.
	for all * call spread(@).

subroutine spread(L):
	recursively do:
		a b c
		d @ #
		f g h
		whereas:
			# is a previous cell of the recursion. (can be any of a..h)
			-- basically we skip over #
		break when {a>=K && b>=K && c>=K && d>=K && f>=K && g>=K && h>=K}.
			a>=b (whereas a and b are cvec4) is true iff a.r>=b.r && a.g>=b.g && a.b>=b.b && a.a>=b.a.
		write to {a,b,c,d,f,g,h} with value K whereas:
			let K be (source_light - attenuation).
			attenuation is a value set onto the block.

[curioustorvald]

Regarding the operator >=: if a>b and b-a is smaller than 2/256 in RGB, then a==b.
To put it simply, break the recursion when the difference is tiny.

[curioustorvald]

Another idea: convert RGBA into YOU-A, where the conversion matrix is:

       [ 3  4  1] [R]   [Y]
.125 * [-2 -2  4] [G] = [O]
       [ 4 -3 -1] [B]   [U]

and subsample O and U in the ratio of 4:2:0. You can also subsample Y in the ratio of 4:2:2 and offset to X position by 1. For the A-channel, it can use either of the methods.

When calculating, subsampled O and U channels will have smaller number of tiles to update. However, custom codes need to be written to handle subsampled channels.

[curioustorvald]

The conversion matrix: use YCoCg-R instead. Wikipedia has optimized implementation that cuts number of computations from 9 to 4.

[curioustorvald]

The new light simulator has been deployed (see branch 'NEWLIGHT2'), closing the issue.