-
Notifications
You must be signed in to change notification settings - Fork 31
/
SuffixAutomaton.go
2533 lines (2285 loc) · 70 KB
/
SuffixAutomaton.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
// https://www.luogu.com/article/w967d5rp
// https://www.bilibili.com/video/BV1S54y1G7P8
// https://www.cnblogs.com/Linshey/p/14219867.html
// https://maspypy.github.io/library/string/suffix_automaton.hpp
// https://github.com/EndlessCheng/codeforces-go/blob/master/copypasta/sam.go
// https://yutong.site/sam/ 可视化
// Blumber 算法在线构建SAM
//
// !子串是一个前缀的后缀.
//
// -s: aababa
//
// -fail tree
//
// !所有的叶子结点都是一个前缀,所有的前缀都是叶子结点(除开1号结点),所有的叶子结点endPos大小都是1.
// link 指向比当前结点短的最长后缀endPos集.
// 每个节点中的子串,都是以该节点为根的子树中的所有子串的"后缀"
// ```
// 0,1,2,3,4,5 ""
// / \
// / \
// / \
// / \
// 0,1,3,5 "a"(后) 2,4 "b"
// / \ / "ab"
// / \ / \
// / \ / \
// / \ / \
// 1,"aa" 3,5 "ba"(后) 2 "aab" 4 "bab"
// / "aba" "abab"
// / \ "aabab"
// / \
// 3 "aaba" 5 "baba"(后)
// "ababa"
// "aababa"
//
//```
// notes:
// 0. 后缀自动机 (Suffix Automaton, SAM) 是仅接受后缀且状态数最少的 DFA.
// 1. 每一个节点都表示一段子串,所有节点表示的子串们都是唯一的.
// 随着子串长度的减小,它有可能还会出现在其他的地方,于是它的endpos 就会多一些,就会分到其他的状态里。
// 2. len表示的是当前节点的最长长度,当前节点的子串长度范围是 [len-link.len+1, len]
// 3. endPos 集合的大小可以通过topo排序求出来,实际上用桶排实现
// 如果必须要求出 endPos 集合的话,可以用set实现树上自底向上启发式合并,或者用动态开点线段树+线段树合并维护endPos集合
// 4. 节点最短串的最长后缀=父结点最长串
// 5. 两个endPos集合要么包含要么不相交
// 6. 一个子串出现次数就是其对应 endPos 集合的size(注意不是长度范围).
// !由于子串<=>前缀的后缀,
// !可以先通过在 SAM 上找到该子串所处的节点,然后求以该节点为根的子树中,有多少个包含原串前缀的节点
// !另一个含义:从SAM的根到这个结点的转移路径条数。
// 7. 可以把SAM理解为把某个串的所有子串建立AC自动机。
// !8. 设 lcs(i,j) 为前缀i,j的最长公共后缀长度,其等于fail树上 LCA 的 maxLen 值。(反串可以将lcp转化为lcs)
// 9. 一个endpos等价类内的串的长度连续.
// 10.理解
// - 从 SAM 的定义上理解:
// SAM 可以看作一种加强版的 Trie,它可以高度压缩一个字符串的子串信息,
// !一条从根出发到`终止结点`的路径对应了原串的一个后缀,而任意一个从根出发的路径对应了原串一个子串。
// 子串和从根出发的路径一一对应。在这种的理解下,每一个结点的含义并不是固定的,
// 它到底对应哪个子串取决于那条路径是怎么到达它的;而边有着确定的含义。
// - 从 Parent Tree 的角度去理解连边的含义
// 两个不同等价类的Endpos集合要么无交集,要么相包含,因此可以建出一个由 Endpos集合的包含关系连结而成的树——Parent Tree
// 它的连边——后缀链接,若是向下看,是在一个等价类的前面加上一个字符,从而分成若干的其他等价类;
// 向上看,它是指向包含当前集合的最小的集合。
// !而后缀自动机的连边是在一个等价类的后面加上一个字母,看看它会指向谁,显然对于同一个添上的字母,这个指向是唯一确定的。
// - 从结点的含义去理解:
// 每一个结点都对应了一种子串,Parent Tree 的结点与 SAM 的结点一一对应
// 但是, 后缀自动机的边不同于 parent 树上的边
//
// !11. 转移边:parent树往下走代表往前加字符,SAM转移边往后走代表往后加字符
// !12. 子串是什么:
// 从SAM的DAG角度看,子串是后缀的一个前缀;
// !从SAM的Parent Tree角度看,子串是前缀的一个后缀。
// !13. SAM与AC自动机的相似性:
// AC自动机的失配链接和后缀自动机的后缀链接都有性质:
// 指向的两个状态都满足"后者的代表串是前者的代表串的真后缀"。
// 可以把 SAM 理解为把某个串的所有子串建立AC自动机.
// 14. 增量构造中,每次从后面加入一个字符, 有两件事要干:
// 找出能转移到这个状态的状态,建立链接;确定这个状态的min,即找到它在parent树上的父亲。
// 15. 对于SAM任何一个节点u,从根到这个节点的路线有 `maxLen(u)-minLen(u)+1` 条,而这条路线则表示原字符串的一个子串,且各不相同.
// 16. 一般来讲,DAG上可能重复转移,是很难跑计数DP的。
// !但是我们知道后缀自动机的性质 : 任意两个节点的表示集合没有交。
// !所以我们只要统计路径数即可,不需要考虑重复问题。
// !17.可以通过parent树确定SAM的接受状态集合(后缀)。找到MaxLen=n的结点,该结点到根的路径上的所有结点都是接受状态。
//
// applications:
// 1. 查找某个子串位于哪个节点 => s[l:r]是s[:r]的一个后缀,倍增往上跳到maxLen不小于(r-l)的最后一个节点.
// 2. 最长可重叠重复子串 => endPos集合大于等于2的那些节点的最大的范围
// 3. `在线`给出模式串的模式匹配问题(单模式串离线=>KMP,多模式串离线=>AC自动机,多模式串在线=>SAM)
// 一般有固定模式串的字符串处理问题和固定主串的字符串处理问题两大类问题。当固定模式串时,熟知的 AC 自动机算法便可以胜任这类问题。
// 如果主串固定,一般采用对主串构造后缀树、后缀自动机来解决这一类问题。
// 4. 两个字符s和t的最长公共子串 => 对s建立SAM,对t的每个前缀,在SAM中寻找这个前缀的最长后缀,类似AC自动机跳fail.
// 5. 最长不可重叠重复子串 => endPos 集合大于等于2,而且还需要考虑最靠右的那个位置和最靠左的那个位置之间的距离
// if(sz[u] >= 2) res = max(res, min(maxLen[u], r[u] - l[u]));
// 6. 读入字符串时删除首部字符 => 记录已读入的字符串长度,若小于等于当前状态的 parent.MaxLen ,就转移到parent
// 7. 判断子串/后缀 => 建出文本串的SAM,将模式串分别输入SAM,若无法转移到则不是子串,否则是;若转移到接受状态则是后缀,否则不是。
// 9. 线段树合并维护每个点的endPos集合 => DFS 树结构时将父结点的线段树并上儿子的线段树,注意使用不销毁点写法的线段树合并。
// 11. 多次询问区间本质不同子串 => 扫描线+LCT 维护区间子串信息
// 12. 维护endPos等价类最小End、最大End => 后缀链接树自底向上更新End下标信息
// 13. 判断子串s[a:b)是否为某个前缀s[:i)的子串 => 先定位子串到pos结点,然后判断该子串最早结束位置 endPosMinEnd[pos] <= i.
// !14. 求子串s[a:b)在s[c:d)中的出现次数 =>
// 记长度 m=b-a
// - 方法1(waveletMatrix):
// 先定位子串到pos结点,然后根据"子串是前缀的后缀"可知满足其他位置出现的s[a:b)都在pos的子树中。
// 合法的endPos满足:叶子结点的 End 在 [c+m-1,d) 之间.
// - 方法2(线段树合并):
// 先定位子串到pos结点,然后看pos结点的endPos集合在[c+m-1,d)之间的元素个数.
package main
import (
"bufio"
"fmt"
"math/bits"
"os"
"sort"
)
const INF int32 = 1e9 + 10
const SIGMA int32 = 26 // 字符集大小
const OFFSET int32 = 'a' // 字符集的起始字符
type Node struct {
Next [SIGMA]int32 // SAM 转移边
Link int32 // 后缀链接
MaxLen int32 // 当前节点对应的最长子串的长度
End int32 // 最长的字符在原串的下标, 实点下标为非负数, 虚点下标为负数
}
func (n *Node) String() string {
return fmt.Sprintf("{Next:%v, Link:%v, MaxLen:%v, End:%v}", n.Next, n.Link, n.MaxLen, n.End)
}
type SuffixAutomaton struct {
Nodes []*Node
LastPos int32 // 当前插入的字符对应的节点(实点,原串的一个前缀)
n int32 // 当前字符串长度
doubling *DoublingSimple
}
func NewSuffixAutomaton() *SuffixAutomaton {
res := &SuffixAutomaton{}
res.Nodes = append(res.Nodes, res.newNode(-1, 0, -1))
return res
}
// 每次插入会增加一个实点,可能增加一个虚点.
// 返回当前前缀对应的节点编号(lastPos).
func (sam *SuffixAutomaton) Add(char int32) int32 {
c := char - OFFSET
newNode := int32(len(sam.Nodes))
// 新增一个实点以表示当前最长串
sam.Nodes = append(sam.Nodes, sam.newNode(-1, sam.Nodes[sam.LastPos].MaxLen+1, sam.Nodes[sam.LastPos].End+1))
p := sam.LastPos
for p != -1 && sam.Nodes[p].Next[c] == -1 {
sam.Nodes[p].Next[c] = newNode
p = sam.Nodes[p].Link
}
q := int32(0)
if p != -1 {
q = sam.Nodes[p].Next[c]
}
if p == -1 || sam.Nodes[p].MaxLen+1 == sam.Nodes[q].MaxLen {
sam.Nodes[newNode].Link = q
} else {
// 不够用,需要新增一个虚点
newQ := int32(len(sam.Nodes))
sam.Nodes = append(sam.Nodes, sam.newNode(sam.Nodes[q].Link, sam.Nodes[p].MaxLen+1, -abs32(sam.Nodes[q].End)))
sam.Nodes[len(sam.Nodes)-1].Next = sam.Nodes[q].Next
sam.Nodes[q].Link = newQ
sam.Nodes[newNode].Link = newQ
for p != -1 && sam.Nodes[p].Next[c] == q {
sam.Nodes[p].Next[c] = newQ
p = sam.Nodes[p].Link
}
}
sam.n++
sam.LastPos = newNode
return sam.LastPos
}
func (sam *SuffixAutomaton) Size() int32 {
return int32(len(sam.Nodes))
}
// 后缀链接树.也叫 parent tree.
func (sam *SuffixAutomaton) BuildTree() [][]int32 {
n := int32(len(sam.Nodes))
graph := make([][]int32, n)
for v := int32(1); v < n; v++ {
p := sam.Nodes[v].Link
graph[p] = append(graph[p], v)
}
return graph
}
func (sam *SuffixAutomaton) BuildDAG() [][]int32 {
n := int32(len(sam.Nodes))
graph := make([][]int32, n)
for v := int32(0); v < n; v++ {
for _, to := range sam.Nodes[v].Next {
if to != -1 {
graph[v] = append(graph[v], to)
}
}
}
return graph
}
// 将结点按照长度进行计数排序,返回后缀链接树的dfs顺序.
// 注意:后缀链接树上父亲的MaxLen值一定小于儿子,但不能认为编号小的节点MaxLen值也小.
// 常数比建图 + dfs 小.
func (sam *SuffixAutomaton) GetDfsOrder() []int32 {
nodes, size, n := sam.Nodes, sam.Size(), sam.n
counter := make([]int32, n+1)
for i := int32(0); i < size; i++ {
counter[nodes[i].MaxLen]++
}
for i := int32(1); i <= n; i++ {
counter[i] += counter[i-1]
}
order := make([]int32, size)
for i := size - 1; i >= 0; i-- {
v := nodes[i].MaxLen
counter[v]--
order[counter[v]] = i
}
return order
}
// 返回每个节点的endPos集合大小.
// !注意:0号结点(空串)大小为n,有时需要置为0.
func (sam *SuffixAutomaton) GetEndPosSize(dfsOrder []int32) []int32 {
size := sam.Size()
endPosSize := make([]int32, size)
for i := size - 1; i >= 1; i-- {
cur := dfsOrder[i]
if sam.Nodes[cur].End >= 0 { // 实点
endPosSize[cur]++
}
pre := sam.Nodes[cur].Link
endPosSize[pre] += endPosSize[cur]
}
return endPosSize
}
// 每个endPos中(最长串)最小的end.
// 返回每个节点对应子串中在原串中的最小结束位置i(0<=i<n).不存在则为INF.
func (sam *SuffixAutomaton) GetEndPosMinEnd(dfsOrder []int32) []int32 {
size := sam.Size()
res := make([]int32, size)
for i := range res {
res[i] = INF
}
for i := size - 1; i >= 1; i-- {
cur := dfsOrder[i]
res[cur] = abs32(sam.Nodes[cur].End)
pre := sam.Nodes[cur].Link
res[pre] = min32(res[pre], res[cur])
}
return res
}
// 每个endPos中(最长串)最大的end.
// 返回每个节点对应子串中在原串中的最大结束位置i(0<=i<n).不存在则为-1.
func (sam *SuffixAutomaton) GetEndPosMaxEnd(dfsOrder []int32) []int32 {
size := sam.Size()
res := make([]int32, size)
for i := range res {
res[i] = -1
}
for i := size - 1; i >= 1; i-- {
cur := dfsOrder[i]
res[cur] = abs32(sam.Nodes[cur].End)
pre := sam.Nodes[cur].Link
res[pre] = max32(res[pre], res[cur])
}
return res
}
// 线段树合并维护 endPos 集合.
// 将父结点的线段树并上儿子的线段树,使用不销毁点写法的线段树合并。
func (sam *SuffixAutomaton) GetEndPos(dfsOrder []int32) (seg *SegmentTreeOnRange, nodes []*SegNode) {
size := sam.Size()
seg = NewSegmentTreeOnRange(0, sam.n-1)
nodes = make([]*SegNode, size)
for i := int32(0); i < size; i++ {
nodes[i] = seg.Alloc()
if end := sam.Nodes[i].End; end >= 0 {
seg.Set(nodes[i], end, 1)
}
}
for i := size - 1; i >= 1; i-- {
cur := dfsOrder[i]
pre := sam.Nodes[cur].Link
nodes[pre] = seg.Merge(nodes[pre], nodes[cur])
}
return
}
// 给定子串的起始位置和结束位置,返回子串在fail树上的位置.
// 快速定位子串, 可以与其它字符串算法配合使用.
// 倍增往上跳到 MaxLen>=end-start 的最后一个节点.
// start: 子串起始位置, end: 子串结束位置, endPosOfEnd: 子串结束位置在fail树上的位置.
func (sam *SuffixAutomaton) LocateSubstring(start, end int32, endPosOfEnd int32) (pos int32) {
target := end - start
_, pos = sam.Doubling().MaxStep(endPosOfEnd, func(p int32) bool { return sam.Nodes[p].MaxLen >= target })
return
}
// 给定结点编号和子串长度,返回该子串的起始和结束位置.
func (sam *SuffixAutomaton) RecoverSubstring(pos int32, len int32) (start, end int32) {
end = abs32(sam.Nodes[pos].End) + 1
start = end - len
return
}
func (sam *SuffixAutomaton) DistinctSubstringAt(pos int32) int32 {
if pos == 0 {
return 0
}
return sam.Nodes[pos].MaxLen - sam.Nodes[sam.Nodes[pos].Link].MaxLen
}
// 本质不同的子串个数.
func (sam *SuffixAutomaton) DistinctSubstring() int {
res := 0
for i := 1; i < len(sam.Nodes); i++ {
res += int(sam.DistinctSubstringAt(int32(i)))
}
return res
}
// 类似AC自动机转移,输入一个字符,返回(转移后的位置, 转移后匹配的"最长后缀"长度).
// pos: 当前状态, len: 当前匹配的长度, char: 输入字符.
func (sam *SuffixAutomaton) Move(pos, len, char int32) (nextPos, nextLen int32) {
char -= OFFSET
if tmp := sam.Nodes[pos].Next[char]; tmp != -1 {
nextPos = tmp
nextLen = len + 1
} else {
for pos != -1 && sam.Nodes[pos].Next[char] == -1 {
pos = sam.Nodes[pos].Link
}
if pos == -1 {
nextPos = 0
nextLen = 0
} else {
nextPos = sam.Nodes[pos].Next[char]
nextLen = sam.Nodes[pos].MaxLen + 1
}
}
return
}
// 删除当前模式串首部字符,返回(转移后的位置, 转移后匹配的"最长后缀"长度).
// pos: 当前状态, len: 当前匹配的长度, patternLen: 当前模式串长度.
func (sam *SuffixAutomaton) MoveLeft(pos, len, patternLen int32) (nextPos, nextLen int32) {
if len < patternLen { // 没有完全匹配,可以不删字符,匹配到的首字母是模式串某个后缀的首字母
return pos, len
}
if len == 0 {
return 0, 0
}
len--
node := sam.Nodes[pos]
if len == sam.Nodes[node.Link].MaxLen {
pos = node.Link
}
return pos, len
}
// 给定模式串pattern,返回模式串的每个非空前缀s[:i+1]与SAM文本串的最长公共后缀长度.
func (sam *SuffixAutomaton) LongestCommonSuffix(m int32, pattern func(i int32) int32) []int32 {
res := make([]int32, m)
pos, len := int32(0), int32(0)
for i := int32(0); i < m; i++ {
pos, len = sam.Move(pos, len, pattern(i))
res[i] = len
}
return res
}
// 后缀连接树上倍增.
func (sam *SuffixAutomaton) Doubling() *DoublingSimple {
if sam.doubling == nil {
size := sam.Size()
doubling := NewDoubling(size, int(size))
for i := int32(1); i < size; i++ {
doubling.Add(i, sam.Nodes[i].Link)
}
doubling.Build()
sam.doubling = doubling
}
return sam.doubling
}
func (sam *SuffixAutomaton) Print(s string) {
dfsOrder := sam.GetDfsOrder()
tree := sam.BuildTree()
fmt.Println("---")
fmt.Println("Fail Tree")
for i := range tree {
if len(tree[i]) > 0 {
fmt.Println(i, tree[i])
}
}
fmt.Println("---")
fmt.Println("EndPos")
fmt.Println("pos,longest,minLen,maxLen,prefixEnd")
for i := int32(1); i < int32(len(dfsOrder)); i++ {
pos := dfsOrder[i]
link := sam.Nodes[pos].Link
minLen, maxLen := sam.Nodes[link].MaxLen+1, sam.Nodes[pos].MaxLen
start, end := sam.RecoverSubstring(pos, maxLen)
sub := s[start:end]
prefixEnd := sam.Nodes[pos].End
fmt.Println(fmt.Sprintf("%v,%v,%v,%v,%v", pos, sub, minLen, maxLen, prefixEnd))
}
fmt.Println("---")
}
func (sam *SuffixAutomaton) newNode(link, maxLen, end int32) *Node {
res := &Node{Link: link, MaxLen: maxLen, End: end}
for i := int32(0); i < SIGMA; i++ {
res.Next[i] = -1
}
return res
}
type S = int32
// SparseTable 稀疏表: st[j][i] 表示区间 [i, i+2^j) 的贡献值.
type SparseTable struct {
st [][]S
lookup []int
e func() S
op func(S, S) S
n int
}
func NewSparseTable(n int, f func(int) S, e func() S, op func(S, S) S) *SparseTable {
res := &SparseTable{}
b := bits.Len(uint(n))
st := make([][]S, b)
for i := range st {
st[i] = make([]S, n)
}
for i := 0; i < n; i++ {
st[0][i] = f(i)
}
for i := 1; i < b; i++ {
for j := 0; j+(1<<i) <= n; j++ {
st[i][j] = op(st[i-1][j], st[i-1][j+(1<<(i-1))])
}
}
lookup := make([]int, n+1)
for i := 2; i < len(lookup); i++ {
lookup[i] = lookup[i>>1] + 1
}
res.st = st
res.lookup = lookup
res.e = e
res.op = op
res.n = n
return res
}
func NewSparseTableFrom(leaves []S, e func() S, op func(S, S) S) *SparseTable {
return NewSparseTable(len(leaves), func(i int) S { return leaves[i] }, e, op)
}
// 查询区间 [start, end) 的贡献值.
func (st *SparseTable) Query(start, end int) S {
if start >= end {
return st.e()
}
b := st.lookup[end-start]
return st.op(st.st[b][start], st.st[b][end-(1<<b)])
}
// 返回最大的 right 使得 [left,right) 内的值满足 check.
func (st *SparseTable) MaxRight(left int, check func(e S) bool) int {
if left == st.n {
return st.n
}
ok, ng := left, st.n+1
for ok+1 < ng {
mid := (ok + ng) >> 1
if check(st.Query(left, mid)) {
ok = mid
} else {
ng = mid
}
}
return ok
}
// 返回最小的 left 使得 [left,right) 内的值满足 check.
func (st *SparseTable) MinLeft(right int, check func(e S) bool) int {
if right == 0 {
return 0
}
ok, ng := right, -1
for ng+1 < ok {
mid := (ok + ng) >> 1
if check(st.Query(mid, right)) {
ok = mid
} else {
ng = mid
}
}
return ok
}
type E = int32
type SegNode struct {
count E
leftChild, rightChild *SegNode
}
func (n *SegNode) String() string {
return fmt.Sprintf("%v", n.count)
}
type SegmentTreeOnRange struct {
min, max int32
}
// 指定闭区间[min,max]建立权值线段树.
func NewSegmentTreeOnRange(min, max int32) *SegmentTreeOnRange {
return &SegmentTreeOnRange{min: min, max: max}
}
// NewRoot().
func (sm *SegmentTreeOnRange) Alloc() *SegNode {
return &SegNode{}
}
// 权值线段树求第 k 小.
// 调用前需保证 1 <= k <= node.value.
func (sm *SegmentTreeOnRange) Kth(node *SegNode, k int32) (value int32, ok bool) {
if k < 1 || k > sm._eval(node) {
return 0, false
}
return sm._kth(k, node, sm.min, sm.max), true
}
func (sm *SegmentTreeOnRange) Get(node *SegNode, index int32) E {
return sm._get(node, index, sm.min, sm.max)
}
func (sm *SegmentTreeOnRange) Set(node *SegNode, index int32, value E) {
sm._set(node, index, value, sm.min, sm.max)
}
func (sm *SegmentTreeOnRange) Query(node *SegNode, left, right int32) E {
return sm._query(node, left, right, sm.min, sm.max)
}
func (sm *SegmentTreeOnRange) QueryAll(node *SegNode) E {
return sm._eval(node)
}
func (sm *SegmentTreeOnRange) Update(node *SegNode, index int32, count E) {
sm._update(node, index, count, sm.min, sm.max)
}
// 用一个新的节点存合并的结果,会生成重合节点数量的新节点.
func (sm *SegmentTreeOnRange) Merge(a, b *SegNode) *SegNode {
return sm._merge(a, b, sm.min, sm.max)
}
// 把第二棵树直接合并到第一棵树上,比较省空间,缺点是会丢失合并前树的信息.
func (sm *SegmentTreeOnRange) MergeDestructively(a, b *SegNode) *SegNode {
return sm._mergeDestructively(a, b, sm.min, sm.max)
}
func (sm *SegmentTreeOnRange) Enumerate(node *SegNode, f func(i int32, count E)) {
sm._enumerate(node, sm.min, sm.max, f)
}
func (sm *SegmentTreeOnRange) _kth(k int32, node *SegNode, left, right int32) int32 {
if left == right {
return left
}
mid := (left + right) >> 1
if leftCount := sm._eval(node.leftChild); leftCount >= k {
return sm._kth(k, node.leftChild, left, mid)
} else {
return sm._kth(k-leftCount, node.rightChild, mid+1, right)
}
}
func (sm *SegmentTreeOnRange) _get(node *SegNode, index int32, left, right int32) E {
if node == nil {
return 0
}
if left == right {
return node.count
}
mid := (left + right) >> 1
if index <= mid {
return sm._get(node.leftChild, index, left, mid)
} else {
return sm._get(node.rightChild, index, mid+1, right)
}
}
func (sm *SegmentTreeOnRange) _query(node *SegNode, L, R int32, left, right int32) E {
if node == nil {
return 0
}
if L <= left && right <= R {
return node.count
}
mid := (left + right) >> 1
if R <= mid {
return sm._query(node.leftChild, L, R, left, mid)
}
if L > mid {
return sm._query(node.rightChild, L, R, mid+1, right)
}
return sm._query(node.leftChild, L, R, left, mid) + sm._query(node.rightChild, L, R, mid+1, right)
}
func (sm *SegmentTreeOnRange) _set(node *SegNode, index int32, count E, left, right int32) {
if left == right {
node.count = count
return
}
mid := (left + right) >> 1
if index <= mid {
if node.leftChild == nil {
node.leftChild = sm.Alloc()
}
sm._set(node.leftChild, index, count, left, mid)
} else {
if node.rightChild == nil {
node.rightChild = sm.Alloc()
}
sm._set(node.rightChild, index, count, mid+1, right)
}
node.count = sm._eval(node.leftChild) + sm._eval(node.rightChild)
}
func (sm *SegmentTreeOnRange) _update(node *SegNode, index int32, count E, left, right int32) {
if left == right {
node.count += count
return
}
mid := (left + right) >> 1
if index <= mid {
if node.leftChild == nil {
node.leftChild = sm.Alloc()
}
sm._update(node.leftChild, index, count, left, mid)
} else {
if node.rightChild == nil {
node.rightChild = sm.Alloc()
}
sm._update(node.rightChild, index, count, mid+1, right)
}
node.count = sm._eval(node.leftChild) + sm._eval(node.rightChild)
}
func (sm *SegmentTreeOnRange) _enumerate(node *SegNode, left, right int32, f func(i int32, count E)) {
if node == nil {
return
}
if left == right {
f(left, node.count)
return
}
mid := (left + right) >> 1
sm._enumerate(node.leftChild, left, mid, f)
sm._enumerate(node.rightChild, mid+1, right, f)
}
func (sm *SegmentTreeOnRange) _merge(a, b *SegNode, left, right int32) *SegNode {
if a == nil || b == nil {
if a == nil {
return b
}
return a
}
newNode := sm.Alloc()
if left == right {
newNode.count = a.count + b.count
return newNode
}
mid := (left + right) >> 1
newNode.leftChild = sm._merge(a.leftChild, b.leftChild, left, mid)
newNode.rightChild = sm._merge(a.rightChild, b.rightChild, mid+1, right)
newNode.count = sm._eval(newNode.leftChild) + sm._eval(newNode.rightChild)
return newNode
}
func (sm *SegmentTreeOnRange) _mergeDestructively(a, b *SegNode, left, right int32) *SegNode {
if a == nil || b == nil {
if a == nil {
return b
}
return a
}
if left == right {
a.count += b.count
return a
}
mid := (left + right) >> 1
a.leftChild = sm._mergeDestructively(a.leftChild, b.leftChild, left, mid)
a.rightChild = sm._mergeDestructively(a.rightChild, b.rightChild, mid+1, right)
a.count = sm._eval(a.leftChild) + sm._eval(a.rightChild)
return a
}
func (sm *SegmentTreeOnRange) _eval(node *SegNode) E {
if node == nil {
return 0
}
return node.count
}
type DoublingSimple struct {
n int32
log int32
to []int32
}
func NewDoubling(n int32, maxStep int) *DoublingSimple {
res := &DoublingSimple{n: n, log: int32(bits.Len(uint(maxStep)))}
size := n * res.log
res.to = make([]int32, size)
for i := int32(0); i < size; i++ {
res.to[i] = -1
}
return res
}
func (d *DoublingSimple) Add(from, to int32) {
d.to[from] = to
}
func (d *DoublingSimple) Build() {
n := d.n
for k := int32(0); k < d.log-1; k++ {
for v := int32(0); v < n; v++ {
w := d.to[k*n+v]
next := (k+1)*n + v
if w == -1 {
d.to[next] = -1
continue
}
d.to[next] = d.to[k*n+w]
}
}
}
// 求从 `from` 状态开始转移 `step` 次的最终状态的编号。
// 不存在时返回 -1。
func (d *DoublingSimple) Jump(from int32, step int) (to int32) {
to = from
for k := int32(0); k < d.log; k++ {
if to == -1 {
return
}
if step&(1<<k) != 0 {
to = d.to[k*d.n+to]
}
}
return
}
// 求从 `from` 状态开始转移 `step` 次,满足 `check` 为 `true` 的最大的 `step` 以及最终状态的编号。
func (d *DoublingSimple) MaxStep(from int32, check func(next int32) bool) (step int, to int32) {
for k := d.log - 1; k >= 0; k-- {
tmp := d.to[k*d.n+from]
if tmp == -1 {
continue
}
if check(tmp) {
step |= 1 << k
from = tmp
}
}
to = from
return
}
// 注意f(i)>=0.
func NewWaveletMatrix32(n int32, f func(i int32) int32) *WaveletMatrix32 {
dataCopy := make([]int32, n)
max_ := int32(0)
for i := int32(0); i < n; i++ {
v := f(i)
if v > max_ {
max_ = v
}
dataCopy[i] = v
}
maxLog := int32(bits.Len32(uint32(max_)) + 1)
mat := make([]*BitVector32, maxLog)
zs := make([]int32, maxLog)
buff1 := make([]int32, maxLog)
buff2 := make([]int32, maxLog)
ls, rs := make([]int32, n), make([]int32, n)
for dep := int32(0); dep < maxLog; dep++ {
mat[dep] = NewBitVector32(n + 1)
p, q := int32(0), int32(0)
for i := int32(0); i < n; i++ {
k := (dataCopy[i] >> (maxLog - dep - 1)) & 1
if k == 1 {
rs[q] = dataCopy[i]
mat[dep].Set(i)
q++
} else {
ls[p] = dataCopy[i]
p++
}
}
zs[dep] = p
mat[dep].Build()
ls = dataCopy
for i := int32(0); i < q; i++ {
dataCopy[p+i] = rs[i]
}
}
return &WaveletMatrix32{
n: n,
maxLog: maxLog,
mat: mat,
zs: zs,
buff1: buff1,
buff2: buff2,
}
}
type WaveletMatrix32 struct {
n int32
maxLog int32
mat []*BitVector32
zs []int32
buff1, buff2 []int32
}
// [start, end) 内的 value 的個数.
func (w *WaveletMatrix32) Count(start, end, value int32) int32 {
return w.count(value, end) - w.count(value, start)
}
// [start, end) 内 [lower, upper) 的个数.
func (w *WaveletMatrix32) CountRange(start, end, lower, upper int32) int32 {
return w.freqDfs(0, start, end, 0, lower, upper)
}
// 第k(0-indexed)个value的位置.
func (w *WaveletMatrix32) Index(value, k int32) int32 {
w.count(value, w.n)
for dep := w.maxLog - 1; dep >= 0; dep-- {
bit := (value >> uint(w.maxLog-dep-1)) & 1
k = w.mat[dep].IndexWithStart(bit, k, w.buff1[dep])
if k < 0 || k >= w.buff2[dep] {
return -1
}
k -= w.buff1[dep]
}
return k
}
func (w *WaveletMatrix32) IndexWithStart(value, k, start int32) int32 {
return w.Index(value, k+w.count(value, start))
}
// [start, end) 内第k(0-indexed)小的数(kth).
func (w *WaveletMatrix32) Kth(start, end, k int32) int32 {
return w.KthMax(start, end, end-start-k-1)
}
// [start, end) 内第k(0-indexed)大的数.
func (w *WaveletMatrix32) KthMax(start, end, k int32) int32 {
if k < 0 || k >= end-start {
return -1
}
res := int32(0)
for dep := int32(0); dep < w.maxLog; dep++ {
p, q := w.mat[dep].Count(1, start), w.mat[dep].Count(1, end)
if k < q-p {
start = w.zs[dep] + p
end = w.zs[dep] + q
res |= 1 << uint(w.maxLog-dep-1)
} else {
k -= q - p
start -= p
end -= q
}
}
return res
}
// [start, end) 内第k(0-indexed)小的数(kth).
func (w *WaveletMatrix32) KthMin(start, end, k int32) int32 {
return w.KthMax(start, end, end-start-k-1)
}
// [start, end) 中比 value 严格小的数, 不存在返回 -INF.
func (w *WaveletMatrix32) Lower(start, end, value int32) int32 {
k := w.lt(start, end, value)
if k != 0 {
return w.KthMin(start, end, k-1)
}
return -INF
}
// [start, end) 中比 value 严格大的数, 不存在返回 INF.
func (w *WaveletMatrix32) Higher(start, end, value int32) int32 {
k := w.le(start, end, value)
if k == end-start {
return INF
}
return w.KthMin(start, end, k)
}
// [start, end) 中不超过 value 的最大值, 不存在返回 -INF.
func (w *WaveletMatrix32) Floor(start, end, value int32) int32 {
count := w.Count(start, end, value)
if count > 0 {
return value
}
return w.Lower(start, end, value)
}
// [start, end) 中不小于 value 的最小值, 不存在返回 INF.
func (w *WaveletMatrix32) Ceiling(start, end, value int32) int32 {
count := w.Count(start, end, value)
if count > 0 {
return value
}
return w.Higher(start, end, value)
}
func (w *WaveletMatrix32) access(k int32) int32 {
res := int32(0)
for dep := int32(0); dep < w.maxLog; dep++ {
bit := w.mat[dep].Get(k)
res = (res << 1) | bit
k = w.mat[dep].Count(bit, k) + w.zs[dep]*dep
}
return res
}
func (w *WaveletMatrix32) count(value, end int32) int32 {
left, right := int32(0), end
for dep := int32(0); dep < w.maxLog; dep++ {
w.buff1[dep] = left
w.buff2[dep] = right
bit := (value >> (w.maxLog - dep - 1)) & 1
left = w.mat[dep].Count(bit, left) + w.zs[dep]*bit
right = w.mat[dep].Count(bit, right) + w.zs[dep]*bit
}
return right - left
}
func (w *WaveletMatrix32) freqDfs(d, left, right, val, a, b int32) int32 {
if left == right {
return 0
}
if d == w.maxLog {
if a <= val && val < b {
return right - left
}
return 0
}
nv := (1 << (w.maxLog - d - 1)) | val
nnv := ((1 << (w.maxLog - d - 1)) - 1) | nv
if nnv < a || b <= val {
return 0
}
if a <= val && nnv < b {
return right - left
}
lc, rc := w.mat[d].Count(1, left), w.mat[d].Count(1, right)
return w.freqDfs(d+1, left-lc, right-rc, val, a, b) + w.freqDfs(d+1, lc+w.zs[d], rc+w.zs[d], nv, a, b)
}
func (w *WaveletMatrix32) ll(left, right, v int32) (int32, int32) {
res := int32(0)
for dep := int32(0); dep < w.maxLog; dep++ {
w.buff1[dep] = left
w.buff2[dep] = right
bit := v >> uint(w.maxLog-dep-1) & 1
if bit == 1 {
res += right - left + w.mat[dep].Count(1, left) - w.mat[dep].Count(1, right)
}