func pollardRhoParallel(h, g, p, order *big.Int) (*big.Int, error) {
	n := new(big.Int).Set(order)
	iterate := func(x, a, b *big.Int) {
		r := new(big.Int)
		s := r.Mod(x, big.NewInt(3)).Int64()
		switch s {
		case 0:
			x.Exp(x, big.NewInt(2), p)
			a.Mod(r.Mul(a, big.NewInt(2)), n)
			b.Mod(r.Mul(b, big.NewInt(2)), n)
		case 1:
			x.Mod(r.Mul(x, g), p)
			a.Mod(r.Add(a, big.NewInt(1)), n)
		case 2:
			x.Mod(r.Mul(x, h), p)
			b.Mod(r.Add(b, big.NewInt(1)), n)
		}
	}
	c := make(chan triple, 8)
	quit := make(chan struct{})
	triples := make(map[string]triple)
	r, q, t := new(big.Int), new(big.Int), new(big.Int)
	sampler := sample.NewUniformRange(big.NewInt(1), new(big.Int).Add(n, big.NewInt(1)))
	for i := 0; i < 8; i++ {
		a0, err := sampler.Sample()
		if err != nil {
			return nil, err
		}
		b0, err := sampler.Sample()
		if err != nil {
			return nil, err
		}
		x0 := new(big.Int).Mod(new(big.Int).Mul(new(big.Int).Exp(g, a0, p), new(big.Int).Exp(h, b0, p)), p)
		go process(x0, a0, b0, iterate, c, quit)
	}
	for {
		triple := <-c
		str := string(triple.x.Bytes())
		if el, ok := triples[str]; ok {
			a1 := triple.a
			a2 := el.a
			b1 := triple.b
			b2 := el.b
			r.Mod(q.Sub(b2, b1), n)
			t.Mod(q.Sub(a1, a2), n)
			if r.Cmp(big.NewInt(0)) == 0 {
				triples[str] = triple
				continue
			}
			close(quit)
			d := new(big.Int).GCD(nil, nil, r, n)
			if d.Cmp(big.NewInt(1)) == 0 {
				q.ModInverse(r, n)
				return q.Mod(q.Mul(q, t), n), nil
			}
			r.Div(r, d)
			t.Div(t, d)
			nDivD := new(big.Int).Div(n, d)
			q.ModInverse(r, nDivD)
			q.Mod(q.Mul(q, t), nDivD)
			one := big.NewInt(1)
			for j := big.NewInt(0); j.Cmp(d) == -1; j.Add(j, one) {
				if h.Cmp(new(big.Int).Exp(g, q, p)) == 0 {
					return q, nil
				}
				q.Add(q, nDivD)
			}
		} else {
			triples[str] = triple
		}
	}
}
Parallelized Pollard's rho algorithm.
package dlog
import (
	"fmt"
	"math"
	"math/big"
	"github.com/fentec-project/gofe/sample"
)
func pollardRho(h, g, p, order *big.Int) (*big.Int, error) {
	n := new(big.Int).Set(order)
	sampler := sample.NewUniformRange(big.NewInt(1), new(big.Int).Add(n, big.NewInt(1)))
	a1, err := sampler.Sample()
	if err != nil {
		return nil, err
	}
	b1, err := sampler.Sample()
	if err != nil {
		return nil, err
	}
	x1 := new(big.Int).Mod(new(big.Int).Mul(new(big.Int).Exp(g, a1, p), new(big.Int).Exp(h, b1, p)), p)
	x2, a2, b2 := new(big.Int).Set(x1), new(big.Int).Set(a1), new(big.Int).Set(b1)
	r, q, t := new(big.Int), new(big.Int), new(big.Int)
	iterate := func(x, a, b *big.Int) {
		s := r.Mod(x, big.NewInt(3)).Int64()
		switch s {
		case 0:
			x.Exp(x, big.NewInt(2), p)
			a.Mod(r.Mul(a, big.NewInt(2)), n)
			b.Mod(r.Mul(b, big.NewInt(2)), n)
		case 1:
			x.Mod(r.Mul(x, g), p)
			a.Mod(r.Add(a, big.NewInt(1)), n)
		case 2:
			x.Mod(r.Mul(x, h), p)
			b.Mod(r.Add(b, big.NewInt(1)), n)
		}
	}
	iterations := math.MaxInt32
	one := big.NewInt(1)
	for i := 0; i < iterations; i++ {
		iterate(x1, a1, b1)
		iterate(x2, a2, b2)
		iterate(x2, a2, b2)
		if x1.Cmp(x2) == 0 {
			r.Mod(q.Sub(b2, b1), n)
			t.Mod(q.Sub(a1, a2), n)
			if r.Cmp(big.NewInt(0)) == 0 {
				return nil, fmt.Errorf("error in Pollard rho: no modular inverse for factor")
			}
			d := new(big.Int).GCD(nil, nil, r, n)
			if d.Cmp(big.NewInt(1)) == 0 {
				q.ModInverse(r, n)
				return q.Mod(q.Mul(q, t), n), nil
			}
			r.Div(r, d)
			t.Div(t, d)
			nDivD := new(big.Int).Div(n, d)
			q.ModInverse(r, nDivD)
			q.Mod(q.Mul(q, t), nDivD)
			for j := big.NewInt(0); j.Cmp(d) == -1; j.Add(j, one) {
				if h.Cmp(new(big.Int).Exp(g, q, p)) == 0 {
					return q, nil
				}
				q.Add(q, nDivD)
			}
		}
	}
	return nil, fmt.Errorf("failed to find discrete logarithm in %d iterations", iterations)
}
type triple struct {
	x *big.Int
	a *big.Int
	b *big.Int
}
func process(x, a, b *big.Int, iterate func(x, a, b *big.Int), c chan<- triple, quit <-chan struct{}) {
	for {
		select {
		case <-quit:
			return
		default:
		}
		iterate(x, a, b)
		xab := triple{
			x: new(big.Int).Set(x),
			a: new(big.Int).Set(a),
			b: new(big.Int).Set(b),
		}
		c <- xab
	}
}