#!/usr/local/bin/wish
#This program illustrates the Galois group of a quartic polynomial with
#coefficients in C(z). To vary the parameter z, click on the plane where
#you wish it to go. Do it in small steps to avoid numerical errors. Each
#step uses the previous step as initial data for a Newton iteration to
#compute the roots. As z (represented by a gray point) varies
#in the complex plane, the roots of P vary along with it (represented
#by the red, blue, green and brown points). The pink points represent
#the branch points and by going around them and returning to the origin
#we see the action of Galois. All five possible groups for an equation
#of degree four are illustrated.
#
#This program is based on a program written for a cubic by Philip D. Smolen
#in a class taught by J. Tate. The program was adapted by F. Voloch for
#quartics. Most changes were in the polynomial algebra.

###########################################################################
# Complex number operations.
# Complex numbers are represented as a list of two real numbers, x then y.
###########################################################################

proc c_+ {z w} {
  set x [lindex $z 0]
  set y [lindex $z 1]
  set u [lindex $w 0]
  set v [lindex $w 1]
  list [expr {$x + $u}] [expr {$y + $v}]
}

proc c_- {z w} {
  set x [lindex $z 0]
  set y [lindex $z 1]
  set u [lindex $w 0]
  set v [lindex $w 1]
  list [expr {$x - $u}] [expr {$y - $v}]
}

proc c_* {z w} {
  set x [lindex $z 0]
  set y [lindex $z 1]
  set u [lindex $w 0]
  set v [lindex $w 1]
  list [expr {$x * $u - $y * $v}] [expr {$x * $v + $y * $u}]
}

proc c_/ {z w} {
    set x [lindex $z 0]
    set y [lindex $z 1]
    set u [lindex $w 0]
    set v [lindex $w 1]
    set denom [expr {$u * $u + $v * $v + 0.0}]
    list [expr {($x * $u + $y * $v) / $denom}] [expr {($y * $u - $x * $v) / $denom}]
}



proc c_abs {z} {
    set x [lindex $z 0]
    set y [lindex $z 1]
    expr {sqrt($x * $x + $y * $y)}
}

proc c_as_string {z} {
    set x [lindex $z 0]
    set y [lindex $z 1]
    format {%0.4f%+0.4fi} $x $y
}
    

###########################################################################
# Polynomial operations.
# P(z,w) = w^4-z*w-2 (symmetric) or w^4-z*w^2-1 (klein) or w^4-z-2 (cyclic)
# or w^4-2*w^2-z-3 (dihedral) or 81w^4 - 216w^2 + (192z - 384)w + (192z - 240)
#(alternating)
###########################################################################

#Branch points. Returns a list of complex numbers.

proc bad_points {} {
global g
    if {$g == "s"} {
    set b 2.0867794
    set e [list $b $b]
    set result [list $e [c_* $e {-1 0}] [c_* $e {0 1}] [c_* $e {0 -1}]]
    } elseif {$g == "k"} {
    set result [list {0 2} {0 -2}]
    } elseif {$g == "c"} {
    set result [list {-2 0}]
    } elseif {$g == "d"} {
    set result [list {2 0} {3 0}]
    } elseif {$g == "a"} {
    set result [list {1 0} {2 0}]
    }
return $result
}

#Newton approximation of root with parameter z and initial
#try w, and 8 iterations. Returns a list of complex numbers.

proc Newton {z w} {
set root $w
for {set i 1} {$i < 9} {incr i} {
set root [c_- $root [c_/ [P $z $root] [dP $z $root]]]
}
return $root
}

# This is the function we are finding the roots of.  

proc P {z w} {
global g
     if {$g == "s"} {
     set result [Horner $w [list {1 0} {0 0} {0 0} [c_* $z {-1 0}] {-2 0}]]
     } elseif {$g == "k"} {
     set result [Horner $w [list {1 0} {0 0} [c_* $z {-1 0}] {0 0} {-1 0}]]
     } elseif {$g == "c"} {
     set result [Horner $w [list {1 0} {0 0} {0 0} {0 0} [c_+ [c_* $z {-1 0}] {-2 0}]]]
     } elseif {$g == "d"} {
     set result [Horner $w [list {1 0} {0 0} {-2 0} {0 0} [c_- {3 0} $z]]]
     } elseif {$g == "a"} {
     set result [Horner $w [list {81 0} {0 0} {-216 0} \
     [c_- [c_* $z {192 0}] {384 0}] [c_- [c_* $z {192 0}] {240 0}]]]
     }
return $result
}

#derivative of P wrt w

proc dP {z w} {
global g
     if {$g == "s"} {
      set result [Horner $w [list {4 0} {0 0} {0 0} [c_* $z {-1 0}]]]
     } elseif {$g == "k"} {
     set result [Horner $w [list {4 0} {0 0} [c_* $z {-2 0}] {0 0}]]
     } elseif {$g == "c"} {
     set result [Horner $w [list {4 0} {0 0} {0 0} {0 0}]]
     } elseif {$g == "d"} {
     set result [Horner $w [list {4 0} {0 0} {-4 0} {0 0}]]
     } elseif {$g == "a"} {
     set result [Horner $w [list {324 0} {0 0} {-432 0} [c_- [c_* $z {192 0}] {384 0}]]]
     }
return $result
}

#Horner's rule for evaluation of polynomials.

proc Horner {z coeff_list} {
set result {0 0}
foreach coeff $coeff_list {
set result [c_+ [c_* $z $result] $coeff]
}
return $result
}




###########################################################################
# GUI
# All I/O is done in complex numbers.
###########################################################################

set border_colors {\#000000 \#ff0000 \#00c000 \#0000ff \#820060}
set fill_colors {\#C0C0C0 \#ff8080 \#80ff80 \#8080ff \#AE0964}

#for now
set degree 4

proc create_shapes {c} {
    global border_colors fill_colors degree
    global path_items start_items end_items
    for {set i 0} {$i <= $degree} {incr i} {
	set path_items($i) [$c create line 0 0 0 0 -fill [lindex $border_colors $i]]
    }
    for {set i 0} {$i <= $degree} {incr i} {
	set start_items($i) [$c create rectangle 0 0 0 0 -outline [lindex $border_colors $i] -fill [lindex $fill_colors $i] -width 2]
    }
    for {set i 0} {$i <= $degree} {incr i} {
	set end_items($i) [$c create oval 0 0 0 0 -outline [lindex $border_colors $i] -fill [lindex $fill_colors $i] -width 2]
    }
}

proc canvas_pos {where x_var y_var} {
    upvar $x_var x
    upvar $y_var y
    set x [expr {[lindex $where 0] * 75.0 + 250}]
    set y [expr {[lindex $where 1] * -75.0 + 250}]
}

proc from_canvas {x y} {
    list [expr {($x - 250) / 75.0}] [expr {($y - 250) / -75.0}]
}

proc reset_pointer {c which where} {
    global path_items start_items end_items path_points
    canvas_pos $where x y
    set path_points($which) [list $x $y]
    $c coords $start_items($which) [expr {$x - 7.5}] [expr {$y - 7.5}] [expr {$x + 7.5}] [expr {$y + 7.5}]
    $c coords $end_items($which) [expr {$x - 5}] [expr {$y - 5}] [expr {$x + 5}] [expr {$y + 5}]
    $c coords $path_items($which) $x $y $x $y
}

proc update_pointer {c which where} {
    global path_items start_items end_items path_points
    canvas_pos $where x y
    lappend path_points($which) $x $y
    $c coords $end_items($which) [expr {$x - 5}] [expr {$y - 5}] [expr {$x + 5}] [expr {$y + 5}]
    eval [concat [list $c coords $path_items($which)] $path_points($which)]
}

proc display_values {} {
    global current_z current_roots degree
    .f.l0 config -text [c_as_string $current_z]
     for {set i 1} {$i <= $degree} {incr i} {
     .f.l$i config -text [c_as_string [lindex $current_roots [expr $i - 1]]]
     }

}

proc reset_z {} {
    global current_z current_roots degree g
    set current_z {0.0 0.0}
    set current_roots ""
    if {$g == "d"} {
    set approximate [list {1 0.5} {-1 0.5} {1 -0.5} {-1 -0.5}]
    } elseif {$g == "a"} {
    set approximate [list {-1 0} {2 0} {-0.5 1} {-0.5 -1}]
    } else {
    set approximate [list {1 0} {-1 0} {0 1} {0 -1}]
    }
    for {set i 0} {$i < $degree} {incr i} {
        lappend current_roots [Newton {0.0 0.0} [lindex $approximate $i]]
        }
    reset_pointer .c 0 $current_z
    for {set a 1} {$a <= $degree} {incr a} {
	reset_pointer .c $a [lindex $current_roots [expr $a - 1]]
    }
    .c delete bad
    foreach z [bad_points] {
	canvas_pos $z x y
	.c create oval [expr {$x - 4}] [expr {$y - 4}] [expr {$x + 4}] [expr {$y + 4}] -fill \#ff00ff -outline {} -tags bad
    }

    display_values
}


proc update2_z {z} {
    global current_z current_roots degree
	for {set i 0} {$i < $degree} {incr i} {
        lappend roots [Newton $z [lindex $current_roots $i]]
	}
    update_pointer .c 0 $z
    for {set a 1} {$a <= $degree} {incr a} {
	update_pointer .c $a [lindex $roots [expr $a - 1]]
    }
    set current_z $z
    set current_roots $roots
    display_values
}




proc init_graphics {} {
    global border_colors degree
    canvas .c -width 500 -height 500 -bg white
    pack .c -side top
    frame .f
    pack .f -side top -fill both
    for {set i 0} {$i <= $degree} {incr i} {
	label .f.l$i -fg [lindex $border_colors $i]
	pack .f.l$i -side left -expand 1
    }
    label .f.err -fg \#ff00ff
    pack .f.err
    create_shapes .c 
    button .b -text Reset -command {reset_z}
    button .bb -text Dismiss -command {destroy .}
    frame .g
    radiobutton .g.s -text S_4 -variable g -value s
    radiobutton .g.a -text A_4 -variable g -value a
    radiobutton .g.d -text D_4 -variable g -value d
    radiobutton .g.k -text K_4 -variable g -value k
    radiobutton .g.c -text C_4 -variable g -value c
    set g "s"
    pack .b -side top -fill both
    pack .bb -side top -fill both
    pack .g -side top
    pack .g.s .g.a .g.d .g.k .g.c -side left
    reset_z
    bind .c <Button-1> {update2_z [from_canvas %x %y]}
}

set g s
init_graphics