diff --git a/OBIPRIMERS_NOTES.md b/OBIPRIMERS_NOTES.md
index e335cc3..d664ced 100644
--- a/OBIPRIMERS_NOTES.md
+++ b/OBIPRIMERS_NOTES.md
@@ -301,3 +301,59 @@ Variantes du protocole :
 
 Décision en suspens : écrire la variante 3 maintenant, ou d'abord finir le pipeline (paires et indices)
 et profiler l'ensemble pour savoir si l'appariement domine encore.
+
+## 12. Mode A, étape 4 : paires et indices (2026-11-03)
+
+Code : `taxa.go` (`Taxa`, `Corpus.AssignTaxa`), `pairs.go` (paramètres, amorces conservées, énumération
+des amplicons, contrainte du 3', comptage), `pairs_eval.go` (`FindPairs`, `PairStats`, Bc, Bs, Yule, Tm),
+test opt-in `TestPairsAgainstEcoPrimers`.
+
+Fidèle à ecoPrimers : amorce conservée = trouvée dans >= floor(|E| x 0,9) exemples ; « good » si en plus
+<= floor(|C| x 0,1) contre-exemples ; une paire a au moins une amorce good ; longueur d'amplicon
+strictement entre `-l` et `-L` ; clé de paire (amorce 1, sens, amorce 2, sens) ordonnée par mot
+canonique ; paire retenue si q > 0,9 et qfp < 0,1 (comparaisons en float32 comme le C) et Bs >= 0,6 ;
+Bc = taxons amplifiés / taxons exemples ; **Bs = taxons identifiés / taxons amplifiés par la paire**
+(et non / tous les taxons : c'est ce que fait le code, l'article est ambigu) ; un taxon est identifié si
+au moins un de ses amplicons n'est partagé avec aucun autre taxon ; tri par Yule x Bs.
+
+Choix d'implémentation (mémoire) : deux phases. Phase 1 ne garde que des compteurs par paire (séquences
+exemples/contre-exemples, longueurs) ; phase 2 ne revisite que les paires retenues (407 / 3 448) et
+n'enregistre qu'un hash 64 bits d'amplicon par occurrence au lieu de stocker les amplicons. ecoPrimers
+conserve un amplicon par occurrence de chaque paire candidate (42 682 paires pour Cyprinidae).
+
+**Contrôle contre les tables complètes d'ecoPrimers** (mêmes paramètres, `-e 3 -3 2 -l 30 -L 150 -c`) :
+
+| Jeu | Paires | Résultat |
+|---|---|---|
+| Cyprinidae | 3 448 / 3 448 | **identique sur toutes les colonnes** : amorces, good, séquences, taxons, Bc, Bs, longueurs, GC, Tm, Tm minimale |
+| Teleostei | 407 / 407 | mêmes paires ; **6 comptent 1 séquence de moins** (NC_033548) |
+
+Cause de l'écart Teleostei **établie par expérience** : l'`extractSite` d'ecoPrimers lit un `N` (code 4)
+dans un registre à décalage ; le 4 déborde sur la base voisine et `~4 & 3` vaut 3 (T), donc un `N` en 3'
+peut être accepté par accident. En émulant exactement ce registre dans le test du 3' (modification
+temporaire, annulée), les **deux tables sont reproduites à 100 %** sur toutes les colonnes. Non conservé :
+une base ambiguë ne peut pas être une correspondance 3' vérifiée. Le test fige donc l'écart (6 paires,
+-1 séquence, rien d'autre ne diffère). Une première émulation plus grossière (N lu comme A) était
+fausse : elle créait un écart dans Cyprinidae.
+
+Autres écarts connus, non testés sur ces données : fenêtres à cheval sur la jonction circulaire, que
+le C lit au-delà du tableau (comportement indéfini) alors que nous reprenons au début ; Tm minimale
+ignorant les sites contenant une base ambiguë (le C en tire des valeurs arbitraires).
+
+**Profil du pipeline complet (16 workers, M3 Max)**, le critère de l'option A :
+
+| Étape | Cyprinidae | Teleostei |
+|---|---|---|
+| lecture + sélection | 0,45 s | 0,44 s |
+| taxons | 0,00 s | 0,00 s |
+| mots stricts | 0,14 s | 1,04 s |
+| **recherche approchée** | **4,14 s** | **10,18 s** |
+| paires + indices + Tm | 1,04 s | 1,12 s |
+| total | ~5,8 s | ~12,8 s |
+
+ecoPrimers complet : 134 s et 350 s. **L'appariement approché reste 72 % à 80 % du total : la variante
+SIMD (variante 3) est justifiée par la mesure.** Mémoire : non isolée (le tas de 2,3 à 2,8 Go du test
+contient la taxonomie NCBI complète).
+
+Non fait : option `-U` (pas de multi-match), séquence de référence `-R`, sets d'amorces, sortie CSV,
+CLI. Prochaine étape : variante 3 (bit-parallèle multi-mots, SIMD).
diff --git a/pkg/obiprimer/pairs.go b/pkg/obiprimer/pairs.go
new file mode 100644
index 0000000..2582bf0
--- /dev/null
+++ b/pkg/obiprimer/pairs.go
@@ -0,0 +1,281 @@
+package obiprimer
+
+import (
+	"cmp"
+	"fmt"
+	"slices"
+	"strings"
+	"sync"
+
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obidefault"
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obikmer"
+)
+
+// PairParams holds the parameters of the primer pair search, named after the
+// ecoPrimers options.
+type PairParams struct {
+	PrimerLength   int     // -O, word size
+	MinLength      int     // -l: the amplicon (primers excluded) must be longer than this
+	MaxLength      int     // -L: and shorter than this
+	ThreePrime     int     // -3: number of 3' bases that must match exactly (0 = no constraint)
+	Sensitivity    float64 // -s: a primer must hit more than this fraction of the examples (0.9)
+	FalsePositive  float64 // -x: and less than this fraction of the counterexamples (0.1)
+	MinSpecificity float64 // -T: pairs with a lower Bs are dropped (0.6)
+	Circular       bool    // -c
+}
+
+func (p PairParams) validate() error {
+	switch {
+	case p.PrimerLength < 2 || p.PrimerLength > 31:
+		return fmt.Errorf("obiprimer: primer length must be between 2 and 31, got %d", p.PrimerLength)
+	case p.MaxLength <= p.MinLength+1:
+		return fmt.Errorf("obiprimer: empty amplicon length range ]%d,%d[", p.MinLength, p.MaxLength)
+	case p.ThreePrime < 0 || p.ThreePrime > p.PrimerLength:
+		return fmt.Errorf("obiprimer: 3' strict bases must be in [0,%d], got %d", p.PrimerLength, p.ThreePrime)
+	case p.Sensitivity <= 0 || p.Sensitivity > 1 || p.FalsePositive < 0 || p.FalsePositive > 1:
+		return fmt.Errorf("obiprimer: quorums must be in ]0,1] and [0,1]")
+	}
+	return nil
+}
+
+// Primer is a conserved word: one found, with at most the allowed errors, in
+// enough example sequences to be a candidate primer.
+type Primer struct {
+	Word            uint64 // canonical word
+	Examples        int    // example sequences where it is found (either strand)
+	CounterExamples int
+	// Good means conserved in the examples and rare enough in the
+	// counterexamples; a pair needs at least one good primer.
+	Good bool
+}
+
+// seqMatch is one occurrence of a primer in one sequence.
+type seqMatch struct {
+	primer  int32
+	pos     uint32
+	reverse bool
+}
+
+// primerSet is the conserved words and, for each corpus sequence, their
+// occurrences sorted by position.
+type primerSet struct {
+	primers []Primer
+	perSeq  [][]seqMatch
+}
+
+// quorumCount is floor(n*fraction), robust to the binary form of the fraction.
+func quorumCount(n int, fraction float64) int { return MinSequences(n, fraction) }
+
+// selectPrimers keeps the words found in at least floor(examples*Sensitivity)
+// example sequences and marks as good those found in at most
+// floor(counterexamples*FalsePositive) counterexamples.
+func selectPrimers(c *Corpus, wm []WordMatches, p PairParams) *primerSet {
+	counts := c.Counts()
+	inQ := quorumCount(counts.Examples, p.Sensitivity)
+	outQ := quorumCount(counts.CounterExamples, p.FalsePositive)
+
+	ps := &primerSet{perSeq: make([][]seqMatch, len(c.Seqs))}
+	for _, w := range wm {
+		ex, co := w.SequencesMatched(c)
+		if ex < inQ {
+			continue
+		}
+		id := int32(len(ps.primers))
+		ps.primers = append(ps.primers, Primer{w.Word, ex, co, co <= outQ})
+		for _, m := range w.Matches {
+			ps.perSeq[m.Seq] = append(ps.perSeq[m.Seq], seqMatch{id, m.Pos, m.Reverse})
+		}
+	}
+	for _, ms := range ps.perSeq {
+		slices.SortStableFunc(ms, func(a, b seqMatch) int { return cmp.Compare(a.pos, b.pos) })
+	}
+	return ps
+}
+
+// pairKey identifies a pair: the two primers, ordered by canonical word, and
+// for each whether the primer is the canonical word as is (true) or its reverse
+// complement (false). It is the key of ecoPrimers' pair tree.
+type pairKey struct {
+	p1, p2 int32
+	a1, a2 bool
+}
+
+// occurrence is one amplicon: a pair of matches in a sequence.
+type occurrence struct {
+	key      pairKey
+	swapped  bool // true when the upstream match is the second primer of the key
+	up, down seqMatch
+	length   int // amplicon length, primers excluded
+}
+
+// forEachOccurrence calls fn for every amplicon of sequence s, applying the
+// length, good-primer and 3' constraints of ecoPrimers.
+func (ps *primerSet) forEachOccurrence(c *Corpus, s int, p PairParams, fn func(occurrence)) {
+	ms := ps.perSeq[s]
+	lp := p.PrimerLength
+	var raw []byte
+	if p.ThreePrime > 0 {
+		raw = c.Seqs[s].Sequence()
+	}
+
+	for i := range ms {
+		for j := i + 1; j < len(ms); j++ {
+			dpos := int(ms[j].pos) - int(ms[i].pos)
+			if dpos <= lp {
+				continue
+			}
+			length := dpos - lp
+			if length >= p.MaxLength {
+				break
+			}
+			up, down := ms[i], ms[j]
+			if length <= p.MinLength || !(ps.primers[up.primer].Good || ps.primers[down.primer].Good) {
+				continue
+			}
+
+			// The upstream match is read as is (primer 1 of the pair when its word
+			// is smaller), the downstream one as the reverse primer.
+			key := pairKey{up.primer, down.primer, !up.reverse, down.reverse}
+			swapped := !(ps.primers[key.p2].Word > ps.primers[key.p1].Word)
+			if swapped {
+				key = pairKey{key.p2, key.p1, key.a2, key.a1}
+			}
+
+			if p.ThreePrime > 0 && !ps.threePrimeMatches(raw, up, down, p) {
+				continue
+			}
+			fn(occurrence{key, swapped, up, down, length})
+		}
+	}
+}
+
+// wordBase returns base t (0 = 5') of a 2-bit packed word of size lp.
+func wordBase(word uint64, lp, t int) byte { return byte(word>>(2*uint(lp-1-t))) & 3 }
+
+// threePrimeMatches reports whether the last p.ThreePrime bases of the forward
+// primer (as found at the upstream match) and of the reverse primer (reverse
+// complement of the downstream window) equal the sequence, i.e. have no
+// mismatch at their 3' end.
+func (ps *primerSet) threePrimeMatches(raw []byte, up, down seqMatch, p PairParams) bool {
+	lp, n := p.PrimerLength, len(raw)
+	at := func(pos int) byte {
+		if p.Circular {
+			pos %= n
+		}
+		return baseCode[raw[pos]]
+	}
+
+	fw := ps.primers[up.primer].Word
+	rv := ps.primers[down.primer].Word
+	for t := lp - p.ThreePrime; t < lp; t++ {
+		// forward primer: the word when found direct, its reverse complement otherwise
+		var f byte
+		if !up.reverse {
+			f = wordBase(fw, lp, t)
+		} else {
+			f = 3 - wordBase(fw, lp, lp-1-t)
+		}
+		if s := at(int(up.pos) + t); s > 3 || s != f {
+			return false
+		}
+
+		// reverse primer: the word when found reverse, its reverse complement otherwise
+		var r byte
+		if down.reverse {
+			r = wordBase(rv, lp, t)
+		} else {
+			r = 3 - wordBase(rv, lp, lp-1-t)
+		}
+		if s := at(int(down.pos) + lp - 1 - t); s > 3 || 3-s != r {
+			return false
+		}
+	}
+	return true
+}
+
+// pairAcc accumulates, over the sequences, the counts ecoPrimers keeps while it
+// builds the pairs.
+type pairAcc struct {
+	examples, counters uint32
+	last               int32 // 1 + last sequence counted, to count each sequence once
+	amplicons          uint32
+	sumLen             uint64
+	minLen, maxLen     uint32
+}
+
+// countPairs enumerates the amplicons of every sequence, in parallel, and
+// returns for each pair the sequences it amplifies and the length statistics
+// of its example amplicons.
+func (ps *primerSet) countPairs(c *Corpus, p PairParams) map[pairKey]*pairAcc {
+	workers := obidefault.ParallelWorkers()
+	partial := make([]map[pairKey]*pairAcc, workers)
+
+	jobs := make(chan int)
+	var wg sync.WaitGroup
+	for w := range workers {
+		wg.Add(1)
+		partial[w] = map[pairKey]*pairAcc{}
+		go func() {
+			defer wg.Done()
+			acc := partial[w]
+			for s := range jobs {
+				isExample := c.Roles[s] == Example
+				ps.forEachOccurrence(c, s, p, func(o occurrence) {
+					a := acc[o.key]
+					if a == nil {
+						a = &pairAcc{minLen: ^uint32(0)}
+						acc[o.key] = a
+					}
+					if a.last != int32(s)+1 {
+						a.last = int32(s) + 1
+						if isExample {
+							a.examples++
+						} else {
+							a.counters++
+						}
+					}
+					if isExample {
+						l := uint32(o.length)
+						a.amplicons++
+						a.sumLen += uint64(l)
+						a.minLen = min(a.minLen, l)
+						a.maxLen = max(a.maxLen, l)
+					}
+				})
+			}
+		}()
+	}
+	for s := range c.Seqs {
+		jobs <- s
+	}
+	close(jobs)
+	wg.Wait()
+
+	// A sequence is handled by one worker only, so the sequence counts add up.
+	total := partial[0]
+	for _, m := range partial[1:] {
+		for k, a := range m {
+			t := total[k]
+			if t == nil {
+				total[k] = a
+				continue
+			}
+			t.examples += a.examples
+			t.counters += a.counters
+			t.amplicons += a.amplicons
+			t.sumLen += a.sumLen
+			t.minLen = min(t.minLen, a.minLen)
+			t.maxLen = max(t.maxLen, a.maxLen)
+		}
+	}
+	return total
+}
+
+// orientedWord returns the upper-case primer string 5'->3': the canonical word
+// when direct, its reverse complement otherwise.
+func orientedWord(word uint64, lp int, direct bool) string {
+	if !direct {
+		word = obikmer.ReverseComplement(word, lp)
+	}
+	return strings.ToUpper(string(obikmer.DecodeKmer(word, lp, nil)))
+}
diff --git a/pkg/obiprimer/pairs_eval.go b/pkg/obiprimer/pairs_eval.go
new file mode 100644
index 0000000..ece6e0a
--- /dev/null
+++ b/pkg/obiprimer/pairs_eval.go
@@ -0,0 +1,360 @@
+package obiprimer
+
+import (
+	"cmp"
+	"fmt"
+	"math"
+	"slices"
+	"sync"
+
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obidefault"
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obithermo"
+)
+
+// PairStats describes one primer pair and its quality indices. Primers are
+// written 5'->3' and ordered so that Primer1 < Primer2 (as in ecoPrimers).
+type PairStats struct {
+	Primer1, Primer2 string
+	Good1, Good2     bool // primer conserved in the examples and rare in the counterexamples
+
+	Examples        int // example sequences amplified
+	CounterExamples int
+	Yule            float32 // Examples/|E| - CounterExamples/|C|
+
+	ExampleTaxa    int     // taxa amplified among the examples
+	CounterTaxa    int     // taxa amplified among the counterexamples
+	Bc             float32 // coverage: ExampleTaxa over all the example taxa
+	Identified     int     // example taxa with an amplicon no other taxon shares
+	Bs             float32 // specificity: Identified over ExampleTaxa
+	MinLen, MaxLen int
+	MeanLen        float64
+	Tm1, Tm2       float64 // melting temperature of the primers against their perfect target
+	MinTm1, MinTm2 float64 // lowest melting temperature over the example amplicons
+	GC1, GC2       int
+}
+
+// Score is the ranking criterion of ecoPrimers.
+func (p PairStats) Score() float32 { return p.Yule * p.Bs }
+
+// ampRecord is an amplicon of a retained pair, in one sequence.
+type ampRecord struct {
+	pair     int32
+	seq      int32
+	hash     uint64
+	tm1, tm2 float32 // NaN when unavailable
+}
+
+// ampHash hashes the amplicon between the two primers read from the side of
+// the first primer of the pair (reverse complemented when the upstream match is
+// the second primer), so that identical amplicons give the same value
+// whatever their strand. Non ACGT bases hash as a fifth symbol.
+func ampHash(raw []byte, from, to int, reverse bool) uint64 {
+	const prime = 1099511628211
+	h := uint64(14695981039346656037) ^ uint64(to-from)
+	h *= prime
+	if !reverse {
+		for _, b := range raw[from:to] {
+			h = (h ^ uint64(baseCode[b])) * prime
+		}
+		return h
+	}
+	for i := to - 1; i >= from; i-- {
+		c := baseCode[raw[i]]
+		if c < 4 {
+			c = 3 - c
+		}
+		h = (h ^ uint64(c)) * prime
+	}
+	return h
+}
+
+// FindPairs searches the primer pairs from the approximate occurrences of the
+// strict words (see MatchApat) and evaluates them.
+//
+// A pair is a pair of conserved primers, at least one of them good, that
+// delimit an amplicon of the allowed length. It is retained when it amplifies
+// more than Sensitivity of the examples and less than FalsePositive of the
+// counterexamples, and its specificity Bs reaches MinSpecificity. The result
+// is sorted by decreasing Yule*Bs.
+//
+// model computes the melting temperatures; nil leaves them NaN.
+func FindPairs(c *Corpus, taxa *Taxa, wm []WordMatches, p PairParams, model *obithermo.Model) ([]PairStats, error) {
+	if err := p.validate(); err != nil {
+		return nil, err
+	}
+	if taxa == nil || len(taxa.Of) != len(c.Seqs) {
+		return nil, fmt.Errorf("obiprimer: taxa do not match the corpus")
+	}
+	ps := selectPrimers(c, wm, p)
+	if len(ps.primers) == 0 {
+		return nil, nil
+	}
+
+	// Phase 1: which pairs amplify enough examples and few enough counterexamples.
+	counts := c.Counts()
+	acc := ps.countPairs(c, p)
+	sens, fpq := float32(p.Sensitivity), float32(p.FalsePositive)
+	var kept []pairKey
+	for k, a := range acc {
+		q := float32(a.examples) / float32(counts.Examples)
+		qfp := float32(0)
+		if counts.CounterExamples > 0 {
+			qfp = float32(a.counters) / float32(counts.CounterExamples)
+		}
+		if q > sens && qfp < fpq {
+			kept = append(kept, k)
+		}
+	}
+	slices.SortFunc(kept, func(a, b pairKey) int {
+		return cmp.Or(cmp.Compare(a.p1, b.p1), cmp.Compare(a.p2, b.p2),
+			boolCmp(a.a1, b.a1), boolCmp(a.a2, b.a2))
+	})
+	index := make(map[pairKey]int32, len(kept))
+	for i, k := range kept {
+		index[k] = int32(i)
+	}
+
+	// Phase 2: the amplicons of the retained pairs only.
+	records := ps.collectAmplicons(c, p, index, model)
+	byPair := make([][]ampRecord, len(kept))
+	for _, r := range records {
+		byPair[r.pair] = append(byPair[r.pair], r)
+	}
+
+	var out []PairStats
+	for i, k := range kept {
+		st, ok := ps.evaluate(c, taxa, p, k, acc[k], byPair[i], model)
+		if ok {
+			out = append(out, st)
+		}
+	}
+	slices.SortFunc(out, func(a, b PairStats) int {
+		return cmp.Or(-cmp.Compare(a.Score(), b.Score()), cmp.Compare(a.Primer1, b.Primer1), cmp.Compare(a.Primer2, b.Primer2))
+	})
+	return out, nil
+}
+
+func boolCmp(a, b bool) int {
+	switch {
+	case a == b:
+		return 0
+	case !a:
+		return -1
+	}
+	return 1
+}
+
+// collectAmplicons gathers, in parallel, one record per amplicon of a retained pair.
+func (ps *primerSet) collectAmplicons(c *Corpus, p PairParams, index map[pairKey]int32, model *obithermo.Model) []ampRecord {
+	lp := p.PrimerLength
+	workers := obidefault.ParallelWorkers()
+	parts := make([][]ampRecord, workers)
+
+	// oriented[primer][0] is the reverse complement of the word, [1] the word itself.
+	oriented := make([][2][]byte, len(ps.primers))
+	for i, pr := range ps.primers {
+		oriented[i][0] = []byte(orientedWord(pr.Word, lp, false))
+		oriented[i][1] = []byte(orientedWord(pr.Word, lp, true))
+	}
+	nan := float32(math.NaN())
+
+	jobs := make(chan int)
+	var wg sync.WaitGroup
+	for w := range workers {
+		wg.Add(1)
+		go func() {
+			defer wg.Done()
+			var site1, site2 [32]byte
+			for s := range jobs {
+				raw := c.Seqs[s].Sequence()
+				n := len(raw)
+				isExample := c.Roles[s] == Example
+				ps.forEachOccurrence(c, s, p, func(o occurrence) {
+					pi, ok := index[o.key]
+					if !ok {
+						return
+					}
+					rec := ampRecord{pair: pi, seq: int32(s), tm1: nan, tm2: nan}
+					if isExample {
+						rec.hash = ampHash(raw, int(o.up.pos)+lp, int(o.down.pos), o.swapped)
+						if model != nil {
+							fp := oriented[o.up.primer][b2i(!o.up.reverse)]
+							rp := oriented[o.down.primer][b2i(o.down.reverse)]
+							fsite, rsite := site1[:lp], site2[:lp]
+							okF := windowAsIs(fsite, raw, int(o.up.pos), n, p.Circular)
+							okR := windowRevComp(rsite, raw, int(o.down.pos), n, p.Circular)
+							tmF, tmR := nan, nan
+							if okF {
+								if t, err := model.Tm(fp, fsite); err == nil {
+									tmF = float32(t)
+								}
+							}
+							if okR {
+								if t, err := model.Tm(rp, rsite); err == nil {
+									tmR = float32(t)
+								}
+							}
+							if o.swapped {
+								rec.tm1, rec.tm2 = tmR, tmF
+							} else {
+								rec.tm1, rec.tm2 = tmF, tmR
+							}
+						}
+					}
+					parts[w] = append(parts[w], rec)
+				})
+			}
+		}()
+	}
+	for s := range c.Seqs {
+		jobs <- s
+	}
+	close(jobs)
+	wg.Wait()
+
+	var all []ampRecord
+	for _, part := range parts {
+		all = append(all, part...)
+	}
+	return all
+}
+
+func b2i(b bool) int {
+	if b {
+		return 1
+	}
+	return 0
+}
+
+const upper = "ACGT"
+
+// windowAsIs copies the len(dst) bases at pos as upper-case ACGT; it reports
+// false if one of them is not A, C, G or T.
+func windowAsIs(dst, raw []byte, pos, n int, circular bool) bool {
+	for t := range dst {
+		i := pos + t
+		if circular {
+			i %= n
+		}
+		c := baseCode[raw[i]]
+		if c > 3 {
+			return false
+		}
+		dst[t] = upper[c]
+	}
+	return true
+}
+
+// windowRevComp writes the reverse complement of the len(dst) bases at pos.
+func windowRevComp(dst, raw []byte, pos, n int, circular bool) bool {
+	l := len(dst)
+	for t := range dst {
+		i := pos + l - 1 - t
+		if circular {
+			i %= n
+		}
+		c := baseCode[raw[i]]
+		if c > 3 {
+			return false
+		}
+		dst[t] = upper[3-c]
+	}
+	return true
+}
+
+// evaluate computes the indices of one retained pair. It reports false when the
+// pair must be dropped (no taxon amplified, or specificity under the threshold).
+func (ps *primerSet) evaluate(c *Corpus, taxa *Taxa, p PairParams, k pairKey, a *pairAcc,
+	recs []ampRecord, model *obithermo.Model) (PairStats, bool) {
+
+	exTaxa := map[int32]struct{}{}
+	coTaxa := map[int32]struct{}{}
+	ampTaxon := map[uint64]int32{} // amplicon -> its taxon, or -2 when shared
+	for _, r := range recs {
+		t := taxa.Of[r.seq]
+		if t < 0 {
+			continue
+		}
+		if c.Roles[r.seq] != Example {
+			coTaxa[t] = struct{}{}
+			continue
+		}
+		exTaxa[t] = struct{}{}
+		if prev, seen := ampTaxon[r.hash]; !seen {
+			ampTaxon[r.hash] = t
+		} else if prev != t {
+			ampTaxon[r.hash] = -2
+		}
+	}
+	if len(exTaxa) == 0 {
+		return PairStats{}, false
+	}
+	identified := map[int32]struct{}{}
+	for _, t := range ampTaxon {
+		if t >= 0 {
+			identified[t] = struct{}{}
+		}
+	}
+
+	counts := c.Counts()
+	st := PairStats{
+		Examples:        int(a.examples),
+		CounterExamples: int(a.counters),
+		ExampleTaxa:     len(exTaxa),
+		CounterTaxa:     len(coTaxa),
+		Identified:      len(identified),
+		MinLen:          int(a.minLen),
+		MaxLen:          int(a.maxLen),
+		MeanLen:         float64(a.sumLen) / float64(a.amplicons),
+	}
+	q := float32(a.examples) / float32(counts.Examples)
+	qfp := float32(0)
+	if counts.CounterExamples > 0 {
+		qfp = float32(a.counters) / float32(counts.CounterExamples)
+	}
+	st.Yule = q - qfp
+	st.Bc = float32(st.ExampleTaxa) / float32(taxa.ExampleTaxa)
+	st.Bs = float32(st.Identified) / float32(st.ExampleTaxa)
+	if st.Bs < float32(p.MinSpecificity) {
+		return PairStats{}, false
+	}
+
+	lp := p.PrimerLength
+	st.Primer1 = orientedWord(ps.primers[k.p1].Word, lp, k.a1)
+	st.Primer2 = orientedWord(ps.primers[k.p2].Word, lp, k.a2)
+	st.Good1, st.Good2 = ps.primers[k.p1].Good, ps.primers[k.p2].Good
+	st.GC1, st.GC2 = obithermo.GCCount([]byte(st.Primer1)), obithermo.GCCount([]byte(st.Primer2))
+
+	st.Tm1, st.Tm2, st.MinTm1, st.MinTm2 = math.NaN(), math.NaN(), math.NaN(), math.NaN()
+	if model != nil {
+		st.Tm1, _ = model.SelfTm([]byte(st.Primer1))
+		st.Tm2, _ = model.SelfTm([]byte(st.Primer2))
+		m1, m2 := float32(math.Inf(1)), float32(math.Inf(1))
+		for _, r := range recs {
+			if c.Roles[r.seq] != Example {
+				continue
+			}
+			if r.tm1 == r.tm1 { // not NaN
+				m1 = min(m1, r.tm1)
+			}
+			if r.tm2 == r.tm2 {
+				m2 = min(m2, r.tm2)
+			}
+		}
+		if !math.IsInf(float64(m1), 1) {
+			st.MinTm1 = float64(m1)
+		}
+		if !math.IsInf(float64(m2), 1) {
+			st.MinTm2 = float64(m2)
+		}
+	}
+
+	if st.Primer2 < st.Primer1 {
+		st.Primer1, st.Primer2 = st.Primer2, st.Primer1
+		st.Good1, st.Good2 = st.Good2, st.Good1
+		st.GC1, st.GC2 = st.GC2, st.GC1
+		st.Tm1, st.Tm2 = st.Tm2, st.Tm1
+		st.MinTm1, st.MinTm2 = st.MinTm2, st.MinTm1
+	}
+	return st, true
+}
diff --git a/pkg/obiprimer/pairs_test.go b/pkg/obiprimer/pairs_test.go
new file mode 100644
index 0000000..8d673d1
--- /dev/null
+++ b/pkg/obiprimer/pairs_test.go
@@ -0,0 +1,237 @@
+package obiprimer
+
+import (
+	"bufio"
+	"fmt"
+	"math"
+	"os"
+	"path/filepath"
+	"runtime"
+	"strconv"
+	"strings"
+	"testing"
+	"time"
+
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiformats"
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obikmer"
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obithermo"
+)
+
+// ecoRow is one result line of ecoPrimers.
+type ecoRow struct {
+	p1, p2         string
+	tm1, mtm1      float64
+	tm2, mtm2      float64
+	gc1, gc2       int
+	good           string
+	in, out        int
+	yule           float64
+	intaxa, outtax int
+	bc             float64
+	well           int
+	bs             float64
+	mind, maxd     int
+	avg            float64
+}
+
+func readEcoRows(t *testing.T, path string) []ecoRow {
+	t.Helper()
+	f, err := os.Open(path)
+	if err != nil {
+		t.Fatal(err)
+	}
+	defer f.Close()
+	var rows []ecoRow
+	sc := bufio.NewScanner(f)
+	sc.Buffer(make([]byte, 1<<20), 1<<20)
+	for sc.Scan() {
+		line := sc.Text()
+		if strings.HasPrefix(line, "#") || strings.TrimSpace(line) == "" {
+			continue
+		}
+		f := strings.Split(line, "\t")
+		fl := func(i int) float64 {
+			v, err := strconv.ParseFloat(strings.TrimSpace(f[i]), 64)
+			if err != nil {
+				t.Fatalf("field %d of %q: %v", i, line, err)
+			}
+			return v
+		}
+		in := func(i int) int { return int(fl(i)) }
+		rows = append(rows, ecoRow{
+			p1: f[1], p2: f[2], tm1: fl(3), mtm1: fl(4), tm2: fl(5), mtm2: fl(6),
+			gc1: in(7), gc2: in(8), good: f[9], in: in(10), out: in(11), yule: fl(12),
+			intaxa: in(13), outtax: in(14), bc: fl(15), well: in(16), bs: fl(17),
+			mind: in(18), maxd: in(19), avg: fl(20),
+		})
+	}
+	if err := sc.Err(); err != nil {
+		t.Fatal(err)
+	}
+	return rows
+}
+
+// Opt-in end-to-end comparison with the complete ecoPrimers result tables of
+// the tutorial: same parameters (-e 3 -3 2 -l 30 -L 150 -c, -r taxid).
+//
+// Cyprinidae is reproduced exactly: its 3448 pairs agree on every column
+// (primers, good flags, sequence and taxon counts, Bc, Bs, lengths, GC, Tm and
+// lowest Tm). Teleostei finds the same 407 pairs, and 6 of them count one
+// sequence less than ecoPrimers (NC_033548, whose downstream site is
+// NGGGGAGGCAAGTCGTAA). The 3' check of ecoPrimers reads an N through the
+// shift register of extractSite, where the code 4 of N spills on the next base
+// and ~4&3 reads as T, so an N can be accepted by accident. Emulating that
+// register exactly in the 3' check was verified to reproduce both tables on
+// every column; it is not kept: an ambiguous base cannot be a verified 3' match.
+func TestPairsAgainstEcoPrimers(t *testing.T) {
+	dir := os.Getenv("OBIPRIMERS_TUTORIAL")
+	if dir == "" {
+		t.Skip("OBIPRIMERS_TUTORIAL not set")
+	}
+	taxonomy, err := obiformats.LoadTaxonomy(filepath.Join(dir, "ncbitaxo_20250211.tgz"), false, false)
+	if err != nil {
+		t.Fatal(err)
+	}
+
+	for _, tc := range []struct {
+		name, taxid string
+		species     int
+		divergent   int
+	}{{"Cyprinidae", "7953", 316, 0}, {"Teleostei", "32443", 3876, 6}} {
+		t.Run(tc.name, func(t *testing.T) {
+			var stages []string
+			tick := time.Now()
+			lap := func(name string) {
+				stages = append(stages, fmt.Sprintf("%s %s", name, time.Since(tick).Round(time.Millisecond)))
+				tick = time.Now()
+			}
+			restrict, err := BelongingTo(taxonomy, []string{tc.taxid})
+			if err != nil {
+				t.Fatal(err)
+			}
+			seqs, err := obiformats.ReadSequencesFromFile(filepath.Join(dir, "mito.vert.fasta.gz"))
+			if err != nil {
+				t.Fatal(err)
+			}
+			corpus, err := LoadCorpus(seqs, Selection{Restrict: restrict})
+			if err != nil {
+				t.Fatal(err)
+			}
+			lap("read+select")
+			taxa, err := corpus.AssignTaxa(taxonomy, "species")
+			if err != nil {
+				t.Fatal(err)
+			}
+			lap("taxa")
+			if taxa.ExampleTaxa != tc.species {
+				t.Errorf("%d example species, ecoPrimers reports %d", taxa.ExampleTaxa, tc.species)
+			}
+
+			sw, err := FindStrictWords(corpus, WordSearch{PrimerLength: 18, Quorum: 0.7, Circular: true, TmpDir: t.TempDir()})
+			if err != nil {
+				t.Fatal(err)
+			}
+			lap("strict words")
+			words := make([]uint64, len(sw.Words))
+			for i, w := range sw.Words {
+				words[i] = w.Kmer
+			}
+			wm, err := matchApat(corpus, words, MatchParams{K: 18, MaxErrors: 3, Circular: true})
+			if err != nil {
+				t.Fatal(err)
+			}
+			lap("approximate search")
+			got, err := FindPairs(corpus, taxa, wm, PairParams{
+				PrimerLength: 18, MinLength: 30, MaxLength: 150, ThreePrime: 2,
+				Sensitivity: 0.9, FalsePositive: 0.1, MinSpecificity: 0.6, Circular: true,
+			}, obithermo.DefaultModel())
+			if err != nil {
+				t.Fatal(err)
+			}
+
+			lap("pairs+indices+Tm")
+			var ms runtime.MemStats
+			runtime.ReadMemStats(&ms)
+			t.Logf("stages: %s; Go heap in use at the end %.0f MB (includes the taxonomy)", strings.Join(stages, ", "), float64(ms.HeapAlloc)/1e6)
+
+			want := readEcoRows(t, filepath.Join(dir, tc.name+".ecoprimers"))
+			t.Logf("%d pairs, ecoPrimers %d", len(got), len(want))
+
+			byKey := map[string]PairStats{}
+			for _, g := range got {
+				byKey[g.Primer1+"/"+g.Primer2] = g
+			}
+			bad := map[string]int{}
+			missing, divergent := 0, 0
+			note := func(field string, w ecoRow, format string, args ...any) {
+				if bad[field]++; bad[field] <= 3 {
+					t.Errorf("%s/%s %s: "+format, append([]any{w.p1, w.p2, field}, args...)...)
+				}
+			}
+			near := func(a, b, tol float64) bool { return math.Abs(a-b) <= tol }
+			for _, w := range want {
+				g, ok := byKey[w.p1+"/"+w.p2]
+				if !ok {
+					if missing++; missing <= 5 {
+						t.Errorf("ecoPrimers pair %s/%s not found", w.p1, w.p2)
+					}
+					continue
+				}
+				delete(byKey, w.p1+"/"+w.p2)
+				good := string("bG"[b2i(g.Good1)]) + string("bG"[b2i(g.Good2)])
+				if good != w.good {
+					note("good", w, "got %s want %s", good, w.good)
+				}
+				// Known divergence, see the comment of the test: ecoPrimers accepts one
+				// amplicon (NC_033548) whose 3' primer site ends with an N. We reject it,
+				// so the counts can be exactly one lower and nothing else may differ.
+				lower := func(got, want int) bool { return got == want || got == want-1 }
+				diverges := g.Examples != w.in || g.ExampleTaxa != w.intaxa || g.Identified != w.well
+				if diverges {
+					divergent++
+				}
+				if !lower(g.Examples, w.in) || g.CounterExamples != w.out {
+					note("sequences", w, "got %d/%d want %d/%d", g.Examples, g.CounterExamples, w.in, w.out)
+				}
+				if !lower(g.ExampleTaxa, w.intaxa) || g.CounterTaxa != w.outtax {
+					note("taxa", w, "got %d/%d want %d/%d", g.ExampleTaxa, g.CounterTaxa, w.intaxa, w.outtax)
+				}
+				if !lower(g.Identified, w.well) {
+					note("identified", w, "got %d want %d", g.Identified, w.well)
+				}
+				tol := 0.0006
+				if diverges {
+					tol = 0.0021
+				}
+				if !near(float64(g.Yule), w.yule, tol) || !near(float64(g.Bc), w.bc, tol) || !near(float64(g.Bs), w.bs, tol) {
+					note("indices", w, "yule/bc/bs got %.4f/%.4f/%.4f want %.3f/%.3f/%.3f", g.Yule, g.Bc, g.Bs, w.yule, w.bc, w.bs)
+				}
+				if g.MinLen != w.mind || g.MaxLen != w.maxd || !near(g.MeanLen, w.avg, 0.006) {
+					note("lengths", w, "got %d/%d/%.2f want %d/%d/%.2f", g.MinLen, g.MaxLen, g.MeanLen, w.mind, w.maxd, w.avg)
+				}
+				if g.GC1 != w.gc1 || g.GC2 != w.gc2 {
+					note("gc", w, "got %d/%d want %d/%d", g.GC1, g.GC2, w.gc1, w.gc2)
+				}
+				// ecoPrimers prints Kelvin-273.0; obithermo uses 273.15
+				if !near(g.Tm1+0.15, w.tm1, 0.051) || !near(g.Tm2+0.15, w.tm2, 0.051) {
+					note("tm", w, "got %.2f/%.2f want %.1f/%.1f", g.Tm1+0.15, g.Tm2+0.15, w.tm1, w.tm2)
+				}
+				if !near(g.MinTm1+0.15, w.mtm1, 0.051) || !near(g.MinTm2+0.15, w.mtm2, 0.051) {
+					note("mintm", w, "got %.2f/%.2f want %.1f/%.1f", g.MinTm1+0.15, g.MinTm2+0.15, w.mtm1, w.mtm2)
+				}
+			}
+			t.Logf("missing from ours: %d; extra in ours: %d; mismatching fields: %v; pairs one sequence below ecoPrimers: %d",
+				missing, len(byKey), bad, divergent)
+			if divergent != tc.divergent {
+				t.Errorf("%d pairs differ from ecoPrimers by one sequence, expected %d", divergent, tc.divergent)
+			}
+			for k := range byKey {
+				if len(byKey) > 0 {
+					t.Errorf("extra pair %s", k)
+				}
+				break
+			}
+		})
+	}
+	_ = obikmer.DecodeKmer
+}
diff --git a/pkg/obiprimer/taxa.go b/pkg/obiprimer/taxa.go
new file mode 100644
index 0000000..bd3417c
--- /dev/null
+++ b/pkg/obiprimer/taxa.go
@@ -0,0 +1,69 @@
+package obiprimer
+
+import (
+	"fmt"
+
+	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obitax"
+)
+
+// Taxa gives every corpus sequence its taxon at one rank (the species by
+// default in ecoPrimers), as a dense index, and the number of distinct taxa
+// among the examples and among the counterexamples.
+//
+// A sequence whose taxon has no ancestor at that rank (or no taxon at all) has
+// the index -1; it is never counted as a taxon.
+type Taxa struct {
+	Rank            string
+	Of              []int32 // Of[i] is the taxon index of Corpus.Seqs[i], or -1
+	ExampleTaxa     int
+	CounterExamples int // distinct taxa among the counterexamples
+	nodes           []*obitax.Taxon
+}
+
+// Name returns the taxon with the given index, for reporting.
+func (t *Taxa) Name(i int32) *obitax.Taxon { return t.nodes[i] }
+
+// AssignTaxa determines the taxon of every sequence of the corpus at rank.
+func (c *Corpus) AssignTaxa(taxonomy *obitax.Taxonomy, rank string) (*Taxa, error) {
+	if taxonomy == nil {
+		return nil, fmt.Errorf("obiprimer: a taxonomy is required to compute taxon statistics")
+	}
+	ranks := taxonomy.RankList()
+	known := false
+	for _, r := range ranks {
+		if r == rank {
+			known = true
+			break
+		}
+	}
+	if !known {
+		return nil, fmt.Errorf("obiprimer: unknown taxonomic rank %q", rank)
+	}
+
+	t := &Taxa{Rank: rank, Of: make([]int32, len(c.Seqs))}
+	index := map[*obitax.TaxNode]int32{}
+	exSeen := map[int32]struct{}{}
+	coSeen := map[int32]struct{}{}
+
+	for i, seq := range c.Seqs {
+		t.Of[i] = -1
+		taxon := seq.Taxon(taxonomy).TaxonAtRank(rank)
+		if taxon == nil {
+			continue
+		}
+		id, ok := index[taxon.Node]
+		if !ok {
+			id = int32(len(t.nodes))
+			index[taxon.Node] = id
+			t.nodes = append(t.nodes, taxon)
+		}
+		t.Of[i] = id
+		if c.Roles[i] == Example {
+			exSeen[id] = struct{}{}
+		} else {
+			coSeen[id] = struct{}{}
+		}
+	}
+	t.ExampleTaxa, t.CounterExamples = len(exSeen), len(coSeen)
+	return t, nil
+}