1	# The FP32 sum of the library: value i goes to lane i % 8, each lane keeps its
2	# rounded sum and the error of every addition (TwoSum), and the lanes merge
3	# in a fixed tree, (0,1),(2,3),(4,5),(6,7), then pairs of those.
4	import Base
5	import ./kernel_tile.bend as Tile
6	
7	# A rounded sum and the error it left out.
8	type Parts is Data:
9	  Parts{hi: F32,lo: F32}
10	
11	# hi = a + b rounded, and lo with hi + lo = a + b exactly when nothing overflows.
12	def two_sum(+a: F32,+b: F32) -> Parts:
13	  +hi = (a + b : F32)
14	  +z = (hi - a : F32)
15	  Parts{hi,((a - (hi - z : F32) : F32) + (b - z : F32) : F32)}
16	
17	# The eight lanes, lane 0 next. After n steps the lane of value n is in front:
18	# slot k holds lane (n + k) % 8, and phase counts the steps.
19	type Lanes is Data:
20	  Lanes{hi: Tile.Tile8,lo: Tile.Tile8,phase: U32}
21	
22	def zero() -> Lanes:
23	  Lanes{Tile.filled(0.0),Tile.filled(0.0),0}
24	
25	def turned(+h1: F32,+h2: F32,+h3: F32,+h4: F32,+h5: F32,+h6: F32,+h7: F32,l1: F32,l2: F32,l3: F32,l4: F32,l5: F32,l6: F32,l7: F32,l0: F32,
26	  +phase: U32,parts: Parts) -> Lanes:
27	  Parts{hi,error} = parts
28	  Lanes{Tile.Tile{h1,h2,h3,h4,h5,h6,h7,hi},Tile.Tile{l1,l2,l3,l4,l5,l6,l7,(l0 + error : F32)},(phase + 1 : U32)}
29	
30	def stepped(+value: F32,lanes: Lanes) -> Lanes:
31	  Lanes{Tile.Tile{+h0,+h1,+h2,+h3,+h4,+h5,+h6,+h7},Tile.Tile{l0,l1,l2,l3,l4,l5,l6,l7},+phase} = lanes
32	  turned(h1,h2,h3,h4,h5,h6,h7,l1,l2,l3,l4,l5,l6,l7,l0,phase,two_sum(h0,value))
33	
34	# The front lane takes the value, then the lanes turn by one.
35	def step(lanes: Lanes,value: F32) -> Lanes:
36	  stepped(value,lanes)
37	
38	def merged(+lo: F32,parts: Parts) -> Parts:
39	  Parts{+hi,error} = parts
40	  two_sum(hi,(lo + error : F32))
41	
42	def merge(left: Parts,right: Parts) -> Parts:
43	  Parts{+ah,al} = left
44	  Parts{+bh,bl} = right
45	  merged((al + bl : F32),two_sum(ah,bh))
46	
47	def tree(p0: Parts,p1: Parts,p2: Parts,p3: Parts,p4: Parts,p5: Parts,p6: Parts,p7: Parts) -> Parts:
48	  merge(merge(merge(p0,p1),merge(p2,p3)),merge(merge(p4,p5),merge(p6,p7)))
49	
50	# Lane j sits in slot (j - turns) % 8.
51	def in_order(turns: Nat,+p0: Parts,+p1: Parts,+p2: Parts,+p3: Parts,+p4: Parts,+p5: Parts,+p6: Parts,+p7: Parts) -> Parts:
52	  match turns:
53	    case 1n: tree(p7,p0,p1,p2,p3,p4,p5,p6)
54	    case 2n: tree(p6,p7,p0,p1,p2,p3,p4,p5)
55	    case 3n: tree(p5,p6,p7,p0,p1,p2,p3,p4)
56	    case 4n: tree(p4,p5,p6,p7,p0,p1,p2,p3)
57	    case 5n: tree(p3,p4,p5,p6,p7,p0,p1,p2)
58	    case 6n: tree(p2,p3,p4,p5,p6,p7,p0,p1)
59	    case 7n: tree(p1,p2,p3,p4,p5,p6,p7,p0)
60	    case turns: tree(p0,p1,p2,p3,p4,p5,p6,p7)
61	
62	def rounded(parts: Parts) -> F32:
63	  Parts{hi,lo} = parts
64	  (0.0 + (hi + lo : F32) : F32)
65	
66	# The sum of the values stepped into the lanes.
67	def total(lanes: Lanes) -> F32:
68	  Lanes{Tile.Tile{h0,h1,h2,h3,h4,h5,h6,h7},Tile.Tile{l0,l1,l2,l3,l4,l5,l6,l7},phase} = lanes
69	  rounded(in_order(U32.to_nat(U32.and(phase,7)),Parts{h0,l0},Parts{h1,l1},Parts{h2,l2},Parts{h3,l3},Parts{h4,l4},Parts{h5,l5},Parts{h6,l6},Parts{h7,l7}))