diff --git a/bargraph.go b/bargraph.go
index f53ef42..9b04d68 100644
--- a/bargraph.go
+++ b/bargraph.go
@@ -11,34 +11,77 @@ import (
 	"fyne.io/fyne/v2/widget"
 )
 
-// barGraphCanvas is a custom widget that renders a bar graph
+// barGraphCanvas is a custom widget that renders a bar graph.
 type barGraphCanvas struct {
 	widget.BaseWidget
 	stats *DiceStatistics
 }
 
 func newBarGraphCanvas(stats *DiceStatistics) *barGraphCanvas {
-	graph := &barGraphCanvas{
-		stats: stats,
-	}
+	graph := &barGraphCanvas{stats: stats}
 	graph.ExtendBaseWidget(graph)
 	return graph
 }
 
 func (b *barGraphCanvas) CreateRenderer() fyne.WidgetRenderer {
 	b.ExtendBaseWidget(b)
-	return &barGraphCanvasRenderer{
-		graph: b,
+
+	background := canvas.NewRectangle(color.NRGBA{R: 20, G: 20, B: 20, A: 255})
+	yAxisLine := canvas.NewLine(color.White)
+	yAxisLine.StrokeWidth = 2
+	xAxisLine := canvas.NewLine(color.White)
+	xAxisLine.StrokeWidth = 2
+	title := canvas.NewText("Probability Distribution", color.White)
+	statsLine1 := canvas.NewText("", color.White)
+	statsLine2 := canvas.NewText("", color.White)
+	yLabel := canvas.NewText("Probability (%)", color.White)
+	xLabel := canvas.NewText("Result Value", color.White)
+
+	renderer := &barGraphCanvasRenderer{
+		graph:      b,
+		background: background,
+		yAxisLine:  yAxisLine,
+		xAxisLine:  xAxisLine,
+		title:      title,
+		statsLine1: statsLine1,
+		statsLine2: statsLine2,
+		yLabel:     yLabel,
+		xLabel:     xLabel,
+		objects: []fyne.CanvasObject{
+			background,
+			yAxisLine,
+			xAxisLine,
+			title,
+			statsLine1,
+			statsLine2,
+			yLabel,
+			xLabel,
+		},
 	}
+
+	renderer.Refresh()
+	return renderer
 }
 
 type barGraphCanvasRenderer struct {
-	graph   *barGraphCanvas
-	objects []fyne.CanvasObject
+	graph       *barGraphCanvas
+	background  *canvas.Rectangle
+	yAxisLine   *canvas.Line
+	xAxisLine   *canvas.Line
+	title       *canvas.Text
+	statsLine1  *canvas.Text
+	statsLine2  *canvas.Text
+	yLabel      *canvas.Text
+	xLabel      *canvas.Text
+	yTicks      []*canvas.Line
+	yTickLabels []*canvas.Text
+	bars        []*canvas.Rectangle
+	barLabels   []*canvas.Text
+	objects     []fyne.CanvasObject
 }
 
 func (r *barGraphCanvasRenderer) Layout(size fyne.Size) {
-	r.Refresh()
+	r.layout(size)
 }
 
 func (r *barGraphCanvasRenderer) MinSize() fyne.Size {
@@ -46,24 +89,75 @@ func (r *barGraphCanvasRenderer) MinSize() fyne.Size {
 }
 
 func (r *barGraphCanvasRenderer) Refresh() {
-	r.objects = []fyne.CanvasObject{}
-
-	if r.graph.stats == nil || len(r.graph.stats.Results) == 0 {
+	stats := r.graph.stats
+	if stats == nil || len(stats.Results) == 0 {
+		r.hideAllButBackground()
+		r.background.Show()
+		r.background.Refresh()
 		return
 	}
 
-	stats := r.graph.stats
-	outcomes := stats.GetSortedOutcomes()
-	maxPercentage := stats.GetMaxPercentage()
-
-	axisMaxPercent, tickStep := calculateYAxisScale(maxPercentage)
-
+	r.background.Show()
+	r.syncData(stats)
 	size := r.graph.Size()
 	if size.Width == 0 || size.Height == 0 {
 		size = fyne.NewSize(900, 550)
 	}
+	r.layout(size)
+}
+
+func (r *barGraphCanvasRenderer) syncData(stats *DiceStatistics) {
+	r.title.Show()
+	r.statsLine1.Show()
+	r.statsLine2.Show()
+	r.yLabel.Show()
+	r.xLabel.Show()
+	r.yAxisLine.Show()
+	r.xAxisLine.Show()
+
+	r.statsLine1.Text = fmt.Sprintf("Range: %d to %d  |  Total Outcomes: %s", stats.MinValue, stats.MaxValue, stats.TotalOutcomesText)
+	r.statsLine2.Text = fmt.Sprintf("Average: %.2f  |  Most Common: %d", stats.Average, stats.MostCommon)
+
+	axisMaxPercent, tickStep := calculateYAxisScale(stats.GetMaxPercentage())
+	numYTicks := int(math.Round(axisMaxPercent / tickStep))
+	r.ensureYTicks(numYTicks + 1)
+	for i := 0; i <= numYTicks; i++ {
+		percent := float64(i) * tickStep
+		r.yTicks[i].Show()
+		r.yTickLabels[i].Text = formatPercentLabel(percent)
+		r.yTickLabels[i].Show()
+	}
+	for i := numYTicks + 1; i < len(r.yTicks); i++ {
+		r.yTicks[i].Hide()
+		r.yTickLabels[i].Hide()
+	}
+
+	outcomes := stats.GetSortedOutcomes()
+	r.ensureBars(len(outcomes))
+	for i := range outcomes {
+		r.bars[i].Show()
+		r.barLabels[i].Hide()
+	}
+	for i := len(outcomes); i < len(r.bars); i++ {
+		r.bars[i].Hide()
+		r.barLabels[i].Hide()
+	}
+
+	r.refreshObjects()
+}
+
+func (r *barGraphCanvasRenderer) layout(size fyne.Size) {
+	r.background.Move(fyne.NewPos(0, 0))
+	r.background.Resize(size)
+
+	stats := r.graph.stats
+	if stats == nil || len(stats.Results) == 0 {
+		return
+	}
+
+	axisMaxPercent, tickStep := calculateYAxisScale(stats.GetMaxPercentage())
+	outcomes := stats.GetSortedOutcomes()
 
-	// Padding and typography scale with the available area so the graph fits the window.
 	topPadding := clamp(size.Height*0.14, 48, 75)
 	bottomPadding := clamp(size.Height*0.15, 56, 80)
 	leftPadding := clamp(size.Width*0.12, 60, 100)
@@ -80,79 +174,43 @@ func (r *barGraphCanvasRenderer) Refresh() {
 	labelTextSize := clamp(size.Height*0.022, 10, 12)
 	tickTextSize := clamp(size.Height*0.018, 8, 10)
 
-	// Background
-	background := canvas.NewRectangle(color.NRGBA{R: 20, G: 20, B: 20, A: 255})
-	background.Move(fyne.NewPos(0, 0))
-	background.Resize(size)
-	r.objects = append(r.objects, background)
+	r.yAxisLine.Move(fyne.NewPos(leftPadding, topPadding))
+	r.yAxisLine.Resize(fyne.NewSize(0, graphHeight))
+	r.xAxisLine.Move(fyne.NewPos(leftPadding, topPadding+graphHeight))
+	r.xAxisLine.Resize(fyne.NewSize(graphWidth, 0))
 
-	// Y-axis
-	yAxisLine := canvas.NewLine(color.White)
-	yAxisLine.StrokeWidth = 2
-	yAxisLine.Move(fyne.NewPos(leftPadding, topPadding))
-	yAxisLine.Resize(fyne.NewSize(0, graphHeight))
-	r.objects = append(r.objects, yAxisLine)
+	r.title.TextSize = titleSize
+	r.title.Move(fyne.NewPos(leftPadding, 5))
 
-	// X-axis
-	xAxisLine := canvas.NewLine(color.White)
-	xAxisLine.StrokeWidth = 2
-	xAxisLine.Move(fyne.NewPos(leftPadding, topPadding+graphHeight))
-	xAxisLine.Resize(fyne.NewSize(graphWidth, 0))
-	r.objects = append(r.objects, xAxisLine)
+	r.statsLine1.TextSize = bodyTextSize
+	r.statsLine1.Move(fyne.NewPos(leftPadding, 5+r.title.MinSize().Height))
 
-	// Title
-	title := canvas.NewText("Probability Distribution", color.White)
-	title.TextSize = titleSize
-	title.Move(fyne.NewPos(leftPadding, 5))
-	r.objects = append(r.objects, title)
+	r.statsLine2.TextSize = bodyTextSize
+	r.statsLine2.Move(fyne.NewPos(leftPadding, 5+r.title.MinSize().Height+r.statsLine1.MinSize().Height))
 
-	// Statistics info line 1
-	statsLine1 := canvas.NewText(fmt.Sprintf("Range: %d to %d  |  Total Outcomes: %d", stats.MinValue, stats.MaxValue, stats.Total), color.White)
-	statsLine1.TextSize = bodyTextSize
-	statsLine1.Move(fyne.NewPos(leftPadding, 5+title.MinSize().Height))
-	r.objects = append(r.objects, statsLine1)
+	r.yLabel.TextSize = labelTextSize
+	r.yLabel.Move(fyne.NewPos(clamp(leftPadding*0.15, 8, 15), topPadding+graphHeight/2-r.yLabel.MinSize().Height/2))
 
-	// Statistics info line 2
-	statsLine2 := canvas.NewText(fmt.Sprintf("Average: %.2f  |  Most Common: %d", stats.Average, stats.MostCommon), color.White)
-	statsLine2.TextSize = bodyTextSize
-	statsLine2.Move(fyne.NewPos(leftPadding, 5+title.MinSize().Height+statsLine1.MinSize().Height))
-	r.objects = append(r.objects, statsLine2)
+	r.xLabel.TextSize = labelTextSize
+	r.xLabel.Move(fyne.NewPos(leftPadding+graphWidth/2-r.xLabel.MinSize().Width/2, topPadding+graphHeight+clamp(bottomPadding*0.35, 18, 28)))
 
-	// Y-axis label
-	yLabel := canvas.NewText("Probability (%)", color.White)
-	yLabel.TextSize = labelTextSize
-	yLabel.Move(fyne.NewPos(clamp(leftPadding*0.15, 8, 15), topPadding+graphHeight/2-yLabel.MinSize().Height/2))
-	r.objects = append(r.objects, yLabel)
-
-	// X-axis label
-	xLabel := canvas.NewText("Result Value", color.White)
-	xLabel.TextSize = labelTextSize
-	xLabel.Move(fyne.NewPos(leftPadding+graphWidth/2-xLabel.MinSize().Width/2, topPadding+graphHeight+clamp(bottomPadding*0.35, 18, 28)))
-	r.objects = append(r.objects, xLabel)
-
-	// Y-axis tick marks and labels
 	numYTicks := int(math.Round(axisMaxPercent / tickStep))
-	for i := 0; i <= numYTicks; i++ {
+	for i := 0; i <= numYTicks && i < len(r.yTicks); i++ {
 		percent := float64(i) * tickStep
-
 		yPos := topPadding + graphHeight - (float32(percent/axisMaxPercent) * graphHeight)
 
-		// Tick mark
-		tick := canvas.NewLine(color.White)
-		tick.StrokeWidth = 1
-		tick.Move(fyne.NewPos(leftPadding-5, yPos))
-		tick.Resize(fyne.NewSize(5, 0))
-		r.objects = append(r.objects, tick)
-
-		// Label
-		label := canvas.NewText(formatPercentLabel(percent), color.White)
-		label.TextSize = tickTextSize
-		label.Move(fyne.NewPos(leftPadding-50, yPos-7))
-		r.objects = append(r.objects, label)
+		r.yTicks[i].Move(fyne.NewPos(leftPadding-5, yPos))
+		r.yTicks[i].Resize(fyne.NewSize(5, 0))
+		r.yTickLabels[i].TextSize = tickTextSize
+		r.yTickLabels[i].Move(fyne.NewPos(leftPadding-50, yPos-7))
 	}
 
-	// Draw bars
 	numBars := len(outcomes)
+	if numBars == 0 {
+		r.refreshObjects()
+		return
+	}
+
 	barSpacing := float32(2)
 	totalSpacing := float32(numBars+1) * barSpacing
 	barWidth := (graphWidth - totalSpacing) / float32(numBars)
@@ -160,45 +218,77 @@ func (r *barGraphCanvasRenderer) Refresh() {
 		barWidth = 1
 	}
 
-	// Calculate label step to prevent overlapping
 	labelStep := calculateLabelStep(graphWidth, numBars)
-
 	for i, value := range outcomes {
 		percentage := stats.Percentages[value]
-
-		// Bar height proportional to percentage
 		barHeight := (float32(percentage) / float32(axisMaxPercent)) * graphHeight
-
-		// X position
 		xPos := leftPadding + barSpacing + float32(i)*(barWidth+barSpacing)
 
-		// Draw bar
-		bar := canvas.NewRectangle(color.NRGBA{R: 100, G: 180, B: 255, A: 255})
-		bar.Move(fyne.NewPos(xPos, topPadding+graphHeight-barHeight))
-		bar.Resize(fyne.NewSize(barWidth, barHeight))
-		r.objects = append(r.objects, bar)
+		r.bars[i].Move(fyne.NewPos(xPos, topPadding+graphHeight-barHeight))
+		r.bars[i].Resize(fyne.NewSize(barWidth, barHeight))
 
-		// X-axis label
-		// Always show first and last label
-		isFirst := i == 0
-		isLast := i == numBars-1
-
-		// Determine if we should show this intermediate label
-		// We show it if it matches the step, BUT we also need to make sure it doesn't clash with the last label
-		// So if we are very close to the end, don't show it (unless it IS the end)
-		showIntermediate := i%labelStep == 0 && i < numBars-labelStep
-
-		if isFirst || isLast || showIntermediate {
-			label := canvas.NewText(fmt.Sprintf("%d", value), color.White)
-			label.TextSize = tickTextSize
-
-			// Center label under bar
-			label.Alignment = fyne.TextAlignCenter
-			label.Move(fyne.NewPos(xPos+barWidth/2-label.MinSize().Width/2, topPadding+graphHeight+10))
-
-			r.objects = append(r.objects, label)
+		r.barLabels[i].TextSize = tickTextSize
+		r.barLabels[i].Alignment = fyne.TextAlignCenter
+		if i == 0 || i == numBars-1 || (labelStep > 0 && i%labelStep == 0 && i < numBars-labelStep) {
+			r.barLabels[i].Text = fmt.Sprintf("%d", value)
+			r.barLabels[i].Show()
+			r.barLabels[i].Move(fyne.NewPos(xPos+barWidth/2-r.barLabels[i].MinSize().Width/2, topPadding+graphHeight+10))
+		} else {
+			r.barLabels[i].Hide()
 		}
 	}
+
+	r.refreshObjects()
+}
+
+func (r *barGraphCanvasRenderer) ensureYTicks(count int) {
+	for len(r.yTicks) < count {
+		tick := canvas.NewLine(color.White)
+		tick.StrokeWidth = 1
+		label := canvas.NewText("", color.White)
+		r.yTicks = append(r.yTicks, tick)
+		r.yTickLabels = append(r.yTickLabels, label)
+		r.objects = append(r.objects, tick, label)
+	}
+}
+
+func (r *barGraphCanvasRenderer) ensureBars(count int) {
+	for len(r.bars) < count {
+		bar := canvas.NewRectangle(color.NRGBA{R: 100, G: 180, B: 255, A: 255})
+		label := canvas.NewText("", color.White)
+		r.bars = append(r.bars, bar)
+		r.barLabels = append(r.barLabels, label)
+		r.objects = append(r.objects, bar, label)
+	}
+}
+
+func (r *barGraphCanvasRenderer) availableGraphWidth() float32 {
+	size := r.graph.Size()
+	if size.Width == 0 {
+		size.Width = 900
+	}
+	leftPadding := clamp(size.Width*0.12, 60, 100)
+	rightPadding := clamp(size.Width*0.03, 16, 24)
+	return size.Width - leftPadding - rightPadding
+}
+
+func (r *barGraphCanvasRenderer) hideAllButBackground() {
+	for _, object := range r.objects {
+		object.Hide()
+	}
+}
+
+func (r *barGraphCanvasRenderer) refreshObjects() {
+	for _, object := range r.objects {
+		object.Refresh()
+	}
+}
+
+func (r *barGraphCanvasRenderer) Objects() []fyne.CanvasObject {
+	return r.objects
+}
+
+func (r *barGraphCanvasRenderer) Destroy() {
 }
 
 func clamp(value, minValue, maxValue float32) float32 {
@@ -212,7 +302,7 @@ func clamp(value, minValue, maxValue float32) float32 {
 }
 
 func calculateLabelStep(graphWidth float32, numBars int) int {
-	labelWidthEstimate := float32(35) // Estimate width of a label
+	labelWidthEstimate := float32(35)
 	maxLabels := int(graphWidth / labelWidthEstimate)
 	if maxLabels < 1 {
 		maxLabels = 1
@@ -258,14 +348,7 @@ func formatPercentLabel(percent float64) string {
 	return fmt.Sprintf("%.1f%%", percent)
 }
 
-func (r *barGraphCanvasRenderer) Objects() []fyne.CanvasObject {
-	return r.objects
-}
-
-func (r *barGraphCanvasRenderer) Destroy() {
-}
-
-// ShowStatisticsWindow creates and shows a statistics window for the given expression
+// ShowStatisticsWindow creates and shows a statistics window for the given expression.
 func ShowStatisticsWindow(expression string) {
 	stats, err := CalculateDiceStatistics(expression)
 	if err != nil {
@@ -273,10 +356,8 @@ func ShowStatisticsWindow(expression string) {
 		return
 	}
 
-	// Create the bar graph
 	graph := newBarGraphCanvas(stats)
 
-	// Create and show the window
 	window := fyne.CurrentApp().NewWindow("Statistics: " + expression)
 	window.SetContent(container.NewMax(graph))
 	window.Resize(fyne.NewSize(900, 550))
diff --git a/statistics.go b/statistics.go
index c0df1e9..56a6477 100644
--- a/statistics.go
+++ b/statistics.go
@@ -3,34 +3,55 @@ package main
 import (
 	"fmt"
 	"math"
+	"math/big"
 	"regexp"
+	"runtime"
 	"sort"
 	"strconv"
 	"strings"
+	"sync"
 )
 
-// DiceStatistics holds the theoretical statistics for a dice roll
+// DiceStatistics holds the theoretical statistics for a dice roll.
 type DiceStatistics struct {
-	MinValue    int
-	MaxValue    int
-	Results     map[int]int     // outcome -> count of ways to achieve it
-	Total       int             // total number of possible outcomes
-	Percentages map[int]float64 // outcome -> percentage
-	Average     float64         // average/mean value
-	MostCommon  int             // most common (median) value
+	MinValue          int
+	MaxValue          int
+	Results           Distribution // outcome -> probability
+	TotalOutcomes     *big.Int     // exact number of equally likely underlying outcomes
+	TotalOutcomesText string       // cached string form for UI
+	Percentages       map[int]float64
+	SortedOutcomes    []int
+	MaxPercentage     float64
+	Average           float64
+	MostCommon        int
 }
 
-// Distribution represents the frequency distribution of outcomes
-type Distribution map[int]int
+// Distribution represents the probability distribution of outcomes.
+type Distribution map[int]float64
+
+type distResult struct {
+	dist          Distribution
+	totalOutcomes *big.Int
+}
+
+type distributionEntry struct {
+	value  int
+	weight float64
+}
+
+const (
+	combineParallelThreshold     = 4096
+	convolutionParallelThreshold = 4096
+	probabilityEpsilon           = 1e-12
+)
 
 // Regex patterns for parsing
 var (
-	diceTokenPattern = regexp.MustCompile(`^([HL])?(\d*)d(\d+)([HL])?`)
-	// Updated numberTokenPattern to include optional decimal part
+	diceTokenPattern   = regexp.MustCompile(`^([HL])?(\d*)d(\d+)([HL])?`)
 	numberTokenPattern = regexp.MustCompile(`^(\d+(\.\d+)?)`)
 )
 
-// CalculateDiceStatistics calculates the theoretical distribution of possible outcomes for a dice expression
+// CalculateDiceStatistics calculates the theoretical distribution of possible outcomes for a dice expression.
 func CalculateDiceStatistics(expression string) (*DiceStatistics, error) {
 	expression = strings.TrimSpace(expression)
 	if expression == "" {
@@ -38,7 +59,7 @@ func CalculateDiceStatistics(expression string) (*DiceStatistics, error) {
 	}
 
 	parser := &statParser{expr: expression, pos: 0}
-	outcomes, err := parser.parseExpression()
+	result, err := parser.parseExpression()
 	if err != nil {
 		return nil, err
 	}
@@ -48,49 +69,18 @@ func CalculateDiceStatistics(expression string) (*DiceStatistics, error) {
 		return nil, fmt.Errorf("unexpected character at position %d: '%c'", parser.pos, parser.expr[parser.pos])
 	}
 
-	if len(outcomes) == 0 {
+	if len(result.dist) == 0 {
 		return nil, fmt.Errorf("no valid outcomes for expression")
 	}
 
-	// Find min and max
-	minVal := 0
-	maxVal := 0
-	first := true
-	totalCount := 0
-
-	for value, count := range outcomes {
-		totalCount += count
-		if first {
-			minVal = value
-			maxVal = value
-			first = false
-		} else {
-			if value < minVal {
-				minVal = value
-			}
-			if value > maxVal {
-				maxVal = value
-			}
-		}
-	}
-
-	// Calculate percentages
-	percentages := make(map[int]float64)
-	for value, count := range outcomes {
-		percentages[value] = (float64(count) / float64(totalCount)) * 100
-	}
-
 	stats := &DiceStatistics{
-		MinValue:    minVal,
-		MaxValue:    maxVal,
-		Results:     outcomes,
-		Total:       totalCount,
-		Percentages: percentages,
+		Results:           normalizeDistribution(result.dist),
+		TotalOutcomes:     cloneBigInt(result.totalOutcomes),
+		TotalOutcomesText: cloneBigInt(result.totalOutcomes).String(),
+		Percentages:       make(map[int]float64, len(result.dist)),
 	}
 
-	// Calculate average and most common value
-	stats.calculateAverageAndMedian()
-
+	stats.populateDerivedFields()
 	return stats, nil
 }
 
@@ -106,11 +96,11 @@ func (p *statParser) skipWhitespace() {
 	}
 }
 
-// parseExpression handles addition and subtraction
-func (p *statParser) parseExpression() (Distribution, error) {
+// parseExpression handles addition and subtraction.
+func (p *statParser) parseExpression() (distResult, error) {
 	left, err := p.parseTerm()
 	if err != nil {
-		return nil, err
+		return distResult{}, err
 	}
 
 	for {
@@ -123,16 +113,16 @@ func (p *statParser) parseExpression() (Distribution, error) {
 			p.pos++
 			right, err := p.parseTerm()
 			if err != nil {
-				return nil, err
+				return distResult{}, err
 			}
-			left = addDist(left, right)
+			left = combineIndependent(left, right, func(a, b int) int { return a + b })
 		} else if p.expr[p.pos] == '-' {
 			p.pos++
 			right, err := p.parseTerm()
 			if err != nil {
-				return nil, err
+				return distResult{}, err
 			}
-			left = subDist(left, right)
+			left = combineIndependent(left, right, func(a, b int) int { return a - b })
 		} else {
 			break
 		}
@@ -141,11 +131,11 @@ func (p *statParser) parseExpression() (Distribution, error) {
 	return left, nil
 }
 
-// parseTerm handles multiplication, division and implicit multiplication
-func (p *statParser) parseTerm() (Distribution, error) {
+// parseTerm handles multiplication, division and implicit multiplication.
+func (p *statParser) parseTerm() (distResult, error) {
 	left, err := p.parsePower()
 	if err != nil {
-		return nil, err
+		return distResult{}, err
 	}
 
 	for {
@@ -159,23 +149,22 @@ func (p *statParser) parseTerm() (Distribution, error) {
 			p.pos++
 			right, err := p.parsePower()
 			if err != nil {
-				return nil, err
+				return distResult{}, err
 			}
-			left = multDist(left, right)
+			left = combineIndependent(left, right, func(a, b int) int { return a * b })
 		} else if c == '/' {
 			p.pos++
 			right, err := p.parsePower()
 			if err != nil {
-				return nil, err
+				return distResult{}, err
 			}
-			left = divDist(left, right)
+			left = combineIndependentFiltered(left, right, divideValues)
 		} else if c == '(' || (c >= '0' && c <= '9') || c == 'd' || c == 'H' || c == 'L' {
-			// Implicit multiplication for things that look like factors
 			right, err := p.parsePower()
 			if err != nil {
-				return nil, err
+				return distResult{}, err
 			}
-			left = multDist(left, right)
+			left = combineIndependent(left, right, func(a, b int) int { return a * b })
 		} else {
 			break
 		}
@@ -184,11 +173,11 @@ func (p *statParser) parseTerm() (Distribution, error) {
 	return left, nil
 }
 
-// parsePower handles exponentiation
-func (p *statParser) parsePower() (Distribution, error) {
+// parsePower handles exponentiation.
+func (p *statParser) parsePower() (distResult, error) {
 	left, err := p.parseFactor()
 	if err != nil {
-		return nil, err
+		return distResult{}, err
 	}
 
 	for {
@@ -201,9 +190,9 @@ func (p *statParser) parsePower() (Distribution, error) {
 			p.pos++
 			right, err := p.parseFactor() // Left-associative to match calculator
 			if err != nil {
-				return nil, err
+				return distResult{}, err
 			}
-			left = powDist(left, right)
+			left = combineIndependent(left, right, powerValues)
 		} else {
 			break
 		}
@@ -212,29 +201,27 @@ func (p *statParser) parsePower() (Distribution, error) {
 	return left, nil
 }
 
-// parseFactor handles parentheses, dice, and numbers
-func (p *statParser) parseFactor() (Distribution, error) {
+// parseFactor handles parentheses, dice, and numbers.
+func (p *statParser) parseFactor() (distResult, error) {
 	p.skipWhitespace()
 	if p.pos >= len(p.expr) {
-		return nil, fmt.Errorf("unexpected end of expression")
+		return distResult{}, fmt.Errorf("unexpected end of expression")
 	}
 
-	// Parentheses
 	if p.expr[p.pos] == '(' {
 		p.pos++
 		dist, err := p.parseExpression()
 		if err != nil {
-			return nil, err
+			return distResult{}, err
 		}
 		p.skipWhitespace()
 		if p.pos >= len(p.expr) || p.expr[p.pos] != ')' {
-			return nil, fmt.Errorf("missing closing parenthesis")
+			return distResult{}, fmt.Errorf("missing closing parenthesis")
 		}
 		p.pos++
 		return dist, nil
 	}
 
-	// Try Dice Pattern
 	remaining := p.expr[p.pos:]
 	if loc := diceTokenPattern.FindStringIndex(remaining); loc != nil {
 		token := remaining[loc[0]:loc[1]]
@@ -242,246 +229,404 @@ func (p *statParser) parseFactor() (Distribution, error) {
 		return parseDiceToken(token)
 	}
 
-	// Try Number Pattern
 	if loc := numberTokenPattern.FindStringIndex(remaining); loc != nil {
 		token := remaining[loc[0]:loc[1]]
 		p.pos += loc[1]
-		// Parse as float then cast to int (truncate/floor) to handle buttons like "."
 		valFloat, err := strconv.ParseFloat(token, 64)
 		if err != nil {
-			return nil, fmt.Errorf("invalid number: %s", token)
+			return distResult{}, fmt.Errorf("invalid number: %s", token)
 		}
-		return Distribution{int(valFloat): 1}, nil
+		return newDistResult(Distribution{int(valFloat): 1}), nil
 	}
 
-	return nil, fmt.Errorf("unexpected character: %c", p.expr[p.pos])
+	return distResult{}, fmt.Errorf("unexpected character: %c", p.expr[p.pos])
 }
 
-func parseDiceToken(token string) (Distribution, error) {
+func parseDiceToken(token string) (distResult, error) {
 	matches := diceTokenPattern.FindStringSubmatch(token)
+	if matches == nil {
+		return distResult{}, fmt.Errorf("invalid dice term: %s", token)
+	}
 
-	if matches != nil {
-		// It is a dice expression
-		prefixModifier := matches[1]
-		countStr := matches[2]
-		sidesStr := matches[3]
-		suffixModifier := matches[4]
+	prefixModifier := matches[1]
+	countStr := matches[2]
+	sidesStr := matches[3]
+	suffixModifier := matches[4]
 
-		count := 1
-		if countStr != "" {
-			c, err := strconv.Atoi(countStr)
-			if err != nil {
-				return nil, err
-			}
-			count = c
-		}
-
-		sides, err := strconv.Atoi(sidesStr)
+	count := 1
+	if countStr != "" {
+		c, err := strconv.Atoi(countStr)
 		if err != nil {
-			return nil, err
+			return distResult{}, err
 		}
-
-		modifier := ""
-		if suffixModifier != "" {
-			modifier = suffixModifier
-		} else if prefixModifier != "" {
-			modifier = prefixModifier
-		}
-
-		return getDiceOutcomes(count, sides, modifier), nil
+		count = c
 	}
 
-	return nil, fmt.Errorf("invalid dice term: %s", token)
-}
-
-// Operations on Distributions
-
-func addDist(a, b Distribution) Distribution {
-	res := make(Distribution)
-	for valA, countA := range a {
-		for valB, countB := range b {
-			res[valA+valB] += countA * countB
-		}
+	sides, err := strconv.Atoi(sidesStr)
+	if err != nil {
+		return distResult{}, err
 	}
-	return res
-}
-
-func subDist(a, b Distribution) Distribution {
-	res := make(Distribution)
-	for valA, countA := range a {
-		for valB, countB := range b {
-			res[valA-valB] += countA * countB
-		}
+	if count <= 0 || sides <= 0 {
+		return distResult{}, fmt.Errorf("dice terms must use positive counts and sides")
 	}
-	return res
-}
 
-func multDist(a, b Distribution) Distribution {
-	res := make(Distribution)
-	for valA, countA := range a {
-		for valB, countB := range b {
-			res[valA*valB] += countA * countB
-		}
+	modifier := ""
+	if suffixModifier != "" {
+		modifier = suffixModifier
+	} else if prefixModifier != "" {
+		modifier = prefixModifier
 	}
-	return res
+
+	return getDiceOutcomes(count, sides, modifier), nil
 }
 
-func divDist(a, b Distribution) Distribution {
-	res := make(Distribution)
-	for valA, countA := range a {
-		for valB, countB := range b {
-			if valB == 0 {
-				continue // Division by zero yields no outcome
+func newDistResult(dist Distribution) distResult {
+	return distResult{
+		dist:          normalizeDistribution(dist),
+		totalOutcomes: big.NewInt(1),
+	}
+}
+
+func cloneBigInt(v *big.Int) *big.Int {
+	if v == nil {
+		return big.NewInt(0)
+	}
+	return new(big.Int).Set(v)
+}
+
+func multiplyOutcomeCounts(a, b *big.Int) *big.Int {
+	return new(big.Int).Mul(cloneBigInt(a), cloneBigInt(b))
+}
+
+func normalizeDistribution(dist Distribution) Distribution {
+	total := 0.0
+	for _, weight := range dist {
+		total += weight
+	}
+	if total == 0 {
+		return dist
+	}
+
+	if math.Abs(total-1) <= probabilityEpsilon {
+		pruned := make(Distribution, len(dist))
+		for value, weight := range dist {
+			if weight != 0 {
+				pruned[value] = weight
 			}
-			res[valA/valB] += countA * countB
+		}
+		return pruned
+	}
+
+	normalized := make(Distribution, len(dist))
+	for value, weight := range dist {
+		probability := weight / total
+		if probability != 0 {
+			normalized[value] = probability
 		}
 	}
-	return res
+	return normalized
 }
 
-func powDist(a, b Distribution) Distribution {
-	res := make(Distribution)
-	for valA, countA := range a {
-		for valB, countB := range b {
-			// Integer exponentiation
-			// Standard behavior for non-negative exponents
-			// Negative exponents with int base result in 0 (unless -1, 1).
-			val := 0
-			if valB >= 0 {
-				val = int(math.Pow(float64(valA), float64(valB)))
-			} else {
-				// Integer division for 1/(a^-b) usually 0
-				val = int(math.Pow(float64(valA), float64(valB)))
+func combineIndependent(a, b distResult, op func(int, int) int) distResult {
+	return distResult{
+		dist:          combineDistributions(a.dist, b.dist, op),
+		totalOutcomes: multiplyOutcomeCounts(a.totalOutcomes, b.totalOutcomes),
+	}
+}
+
+func combineIndependentFiltered(a, b distResult, op func(int, int) (int, bool)) distResult {
+	return distResult{
+		dist:          combineDistributionsFiltered(a.dist, b.dist, op),
+		totalOutcomes: multiplyOutcomeCounts(a.totalOutcomes, b.totalOutcomes),
+	}
+}
+
+func combineDistributions(a, b Distribution, op func(int, int) int) Distribution {
+	if len(a) == 0 || len(b) == 0 {
+		return Distribution{}
+	}
+
+	entriesA := make([]distributionEntry, 0, len(a))
+	for value, weight := range a {
+		entriesA = append(entriesA, distributionEntry{value: value, weight: weight})
+	}
+
+	entriesB := make([]distributionEntry, 0, len(b))
+	for value, weight := range b {
+		entriesB = append(entriesB, distributionEntry{value: value, weight: weight})
+	}
+
+	workSize := len(entriesA) * len(entriesB)
+	if workSize < combineParallelThreshold || len(entriesA) < 2 {
+		return combineDistributionsSerial(entriesA, entriesB, op)
+	}
+
+	workerCount := runtime.GOMAXPROCS(0)
+	if workerCount > len(entriesA) {
+		workerCount = len(entriesA)
+	}
+	if workerCount < 2 {
+		return combineDistributionsSerial(entriesA, entriesB, op)
+	}
+
+	chunkSize := (len(entriesA) + workerCount - 1) / workerCount
+	partials := make([]Distribution, workerCount)
+	var wg sync.WaitGroup
+
+	for worker := 0; worker < workerCount; worker++ {
+		start := worker * chunkSize
+		if start >= len(entriesA) {
+			break
+		}
+		end := start + chunkSize
+		if end > len(entriesA) {
+			end = len(entriesA)
+		}
+
+		wg.Add(1)
+		go func(workerIndex, from, to int) {
+			defer wg.Done()
+			local := make(Distribution)
+			for _, left := range entriesA[from:to] {
+				for _, right := range entriesB {
+					local[op(left.value, right.value)] += left.weight * right.weight
+				}
 			}
-			res[val] += countA * countB
-		}
-	}
-	return res
-}
-
-// getDiceOutcomes returns a map of all possible outcomes for a dice roll and their frequencies
-func getDiceOutcomes(count int, sides int, modifier string) map[int]int {
-	outcomes := make(map[int]int)
-
-	if modifier == "H" {
-		// Keep only the highest die
-		generateHighestOutcomes(count, sides, []int{}, outcomes)
-	} else if modifier == "L" {
-		// Keep only the lowest die
-		generateLowestOutcomes(count, sides, []int{}, outcomes)
-	} else {
-		// Sum all dice
-		generateSumOutcomes(count, sides, []int{}, outcomes)
+			partials[workerIndex] = local
+		}(worker, start, end)
 	}
 
-	return outcomes
+	wg.Wait()
+
+	result := make(Distribution)
+	for _, partial := range partials {
+		for value, weight := range partial {
+			result[value] += weight
+		}
+	}
+	return normalizeDistribution(result)
 }
 
-// generateSumOutcomes recursively generates all sums
-func generateSumOutcomes(remaining int, sides int, current []int, outcomes map[int]int) {
-	if remaining == 0 {
-		sum := 0
-		for _, val := range current {
-			sum += val
+func combineDistributionsFiltered(a, b Distribution, op func(int, int) (int, bool)) Distribution {
+	if len(a) == 0 || len(b) == 0 {
+		return Distribution{}
+	}
+
+	result := make(Distribution)
+	for leftValue, leftWeight := range a {
+		for rightValue, rightWeight := range b {
+			value, ok := op(leftValue, rightValue)
+			if !ok {
+				continue
+			}
+			result[value] += leftWeight * rightWeight
 		}
-		outcomes[sum]++
+	}
+	return normalizeDistribution(result)
+}
+
+func combineDistributionsSerial(entriesA, entriesB []distributionEntry, op func(int, int) int) Distribution {
+	result := make(Distribution)
+	for _, left := range entriesA {
+		for _, right := range entriesB {
+			result[op(left.value, right.value)] += left.weight * right.weight
+		}
+	}
+	return normalizeDistribution(result)
+}
+
+func divideValues(a, b int) (int, bool) {
+	if b == 0 {
+		return 0, false
+	}
+	return a / b, true
+}
+
+func powerValues(a, b int) int {
+	return int(math.Pow(float64(a), float64(b)))
+}
+
+// getDiceOutcomes returns a probability distribution for a dice roll and the exact total outcome count.
+func getDiceOutcomes(count int, sides int, modifier string) distResult {
+	totalOutcomes := new(big.Int).Exp(big.NewInt(int64(sides)), big.NewInt(int64(count)), nil)
+
+	var outcomes Distribution
+	switch modifier {
+	case "H":
+		outcomes = highestDieDistribution(count, sides)
+	case "L":
+		outcomes = lowestDieDistribution(count, sides)
+	default:
+		outcomes = sumDiceDistribution(count, sides)
+	}
+
+	return distResult{
+		dist:          outcomes,
+		totalOutcomes: totalOutcomes,
+	}
+}
+
+func sumDiceDistribution(count int, sides int) Distribution {
+	current := []float64{1}
+	singleDieProbability := 1 / float64(sides)
+
+	for die := 0; die < count; die++ {
+		next := make([]float64, len(current)+sides)
+		if len(current)*sides >= convolutionParallelThreshold {
+			convolveStepParallel(current, next, sides, singleDieProbability)
+		} else {
+			convolveStepSerial(current, next, sides, singleDieProbability)
+		}
+		current = next
+	}
+
+	outcomes := make(Distribution, count*(sides-1)+1)
+	for sum, probability := range current {
+		if probability != 0 {
+			outcomes[sum] = probability
+		}
+	}
+	return normalizeDistribution(outcomes)
+}
+
+func convolveStepSerial(current, next []float64, sides int, singleDieProbability float64) {
+	for sum, probability := range current {
+		if probability == 0 {
+			continue
+		}
+		for face := 1; face <= sides; face++ {
+			next[sum+face] += probability * singleDieProbability
+		}
+	}
+}
+
+func convolveStepParallel(current, next []float64, sides int, singleDieProbability float64) {
+	workerCount := runtime.GOMAXPROCS(0)
+	if workerCount > len(current) {
+		workerCount = len(current)
+	}
+	if workerCount < 2 {
+		convolveStepSerial(current, next, sides, singleDieProbability)
 		return
 	}
 
-	for die := 1; die <= sides; die++ {
-		generateSumOutcomes(remaining-1, sides, append(current, die), outcomes)
-	}
-}
+	chunkSize := (len(current) + workerCount - 1) / workerCount
+	partials := make([][]float64, workerCount)
+	var wg sync.WaitGroup
 
-// generateHighestOutcomes recursively generates all highest-die outcomes
-func generateHighestOutcomes(remaining int, sides int, current []int, outcomes map[int]int) {
-	if remaining == 0 {
-		highest := 0
-		for _, val := range current {
-			if val > highest {
-				highest = val
+	for worker := 0; worker < workerCount; worker++ {
+		start := worker * chunkSize
+		if start >= len(current) {
+			break
+		}
+		end := start + chunkSize
+		if end > len(current) {
+			end = len(current)
+		}
+
+		wg.Add(1)
+		go func(workerIndex, from, to int) {
+			defer wg.Done()
+			local := make([]float64, len(next))
+			for sum := from; sum < to; sum++ {
+				probability := current[sum]
+				if probability == 0 {
+					continue
+				}
+				for face := 1; face <= sides; face++ {
+					local[sum+face] += probability * singleDieProbability
+				}
 			}
+			partials[workerIndex] = local
+		}(worker, start, end)
+	}
+
+	wg.Wait()
+
+	for _, partial := range partials {
+		if partial == nil {
+			continue
 		}
-		outcomes[highest]++
-		return
-	}
-
-	for die := 1; die <= sides; die++ {
-		generateHighestOutcomes(remaining-1, sides, append(current, die), outcomes)
-	}
-}
-
-// generateLowestOutcomes recursively generates all lowest-die outcomes
-func generateLowestOutcomes(remaining int, sides int, current []int, outcomes map[int]int) {
-	if remaining == 0 {
-		lowest := sides + 1
-		for _, val := range current {
-			if val < lowest {
-				lowest = val
-			}
-		}
-		outcomes[lowest]++
-		return
-	}
-
-	for die := 1; die <= sides; die++ {
-		generateLowestOutcomes(remaining-1, sides, append(current, die), outcomes)
-	}
-}
-
-// GetSortedOutcomes returns sorted unique outcomes
-func (s *DiceStatistics) GetSortedOutcomes() []int {
-	var outcomes []int
-	for value := range s.Results {
-		outcomes = append(outcomes, value)
-	}
-	sort.Ints(outcomes)
-	return outcomes
-}
-
-// GetMaxPercentage returns the maximum percentage value
-func (s *DiceStatistics) GetMaxPercentage() float64 {
-	maxPercentage := 0.0
-	for _, percentage := range s.Percentages {
-		if percentage > maxPercentage {
-			maxPercentage = percentage
+		for i, probability := range partial {
+			next[i] += probability
 		}
 	}
-	return maxPercentage
 }
 
-// calculateAverageAndMedian calculates the average and most common value
-func (s *DiceStatistics) calculateAverageAndMedian() {
+func highestDieDistribution(count int, sides int) Distribution {
+	outcomes := make(Distribution, sides)
+	denominator := float64(sides)
+	for value := 1; value <= sides; value++ {
+		current := math.Pow(float64(value)/denominator, float64(count))
+		previous := math.Pow(float64(value-1)/denominator, float64(count))
+		probability := current - previous
+		if probability != 0 {
+			outcomes[value] = probability
+		}
+	}
+	return normalizeDistribution(outcomes)
+}
+
+func lowestDieDistribution(count int, sides int) Distribution {
+	outcomes := make(Distribution, sides)
+	denominator := float64(sides)
+	for value := 1; value <= sides; value++ {
+		current := math.Pow(float64(sides-value+1)/denominator, float64(count))
+		next := math.Pow(float64(sides-value)/denominator, float64(count))
+		probability := current - next
+		if probability != 0 {
+			outcomes[value] = probability
+		}
+	}
+	return normalizeDistribution(outcomes)
+}
+
+func (s *DiceStatistics) populateDerivedFields() {
 	if len(s.Results) == 0 {
-		s.Average = 0
-		s.MostCommon = 0
+		s.TotalOutcomesText = cloneBigInt(s.TotalOutcomes).String()
 		return
 	}
 
-	// Calculate average (mean)
-	sum := 0
-	totalCount := 0
-	for value, count := range s.Results {
-		sum += value * count
-		totalCount += count
-	}
-	s.Average = float64(sum) / float64(totalCount)
+	s.SortedOutcomes = make([]int, 0, len(s.Results))
+	s.Percentages = make(map[int]float64, len(s.Results))
 
-	// Find most common (mode) - the value with highest count
-	maxCount := 0
-	for value, count := range s.Results {
-		if count > maxCount {
-			maxCount = count
+	first := true
+	maxProbability := 0.0
+	for value, probability := range s.Results {
+		s.SortedOutcomes = append(s.SortedOutcomes, value)
+		s.Percentages[value] = probability * 100
+
+		if first {
+			s.MinValue = value
+			s.MaxValue = value
+			s.MostCommon = value
+			first = false
+		} else {
+			if value < s.MinValue {
+				s.MinValue = value
+			}
+			if value > s.MaxValue {
+				s.MaxValue = value
+			}
+		}
+
+		s.Average += float64(value) * probability
+		if probability > maxProbability || (math.Abs(probability-maxProbability) <= probabilityEpsilon && value < s.MostCommon) {
+			maxProbability = probability
 			s.MostCommon = value
 		}
 	}
 
-	// If there are tied values, choose the smallest one
-	if maxCount > 0 {
-		for value, count := range s.Results {
-			if count == maxCount && value < s.MostCommon {
-				s.MostCommon = value
-			}
-		}
-	}
+	sort.Ints(s.SortedOutcomes)
+	s.MaxPercentage = maxProbability * 100
+}
+
+// GetSortedOutcomes returns cached sorted unique outcomes.
+func (s *DiceStatistics) GetSortedOutcomes() []int {
+	return s.SortedOutcomes
+}
+
+// GetMaxPercentage returns the cached maximum percentage value.
+func (s *DiceStatistics) GetMaxPercentage() float64 {
+	return s.MaxPercentage
 }