Hensel's Algorithm in Two Variables

1678 days ago by bardg

f(x,y) = (3*x + 4*x*y + 5*y^2)*( (x-715)^2 + (y-232)^2 ) 
       
g(x,y) = (x - 2*x*y + 7*y^2)*( (x-715)^2 + (y-232)^2 ) 
       
g( 715, 232 ) 
       
0
0
f( 715, 232 ) 
       
0
0
f(x,y).expand() 
       
4*x^3*y + 5*x^2*y^2 + 4*x*y^3 + 3*x^3 + 5*y^4 - 5720*x^2*y - 9003*x*y^2
- 4290*x^2 - 2320*y^3 + 2258804*x*y + 2825245*y^2 + 1695147*x
4*x^3*y + 5*x^2*y^2 + 4*x*y^3 + 3*x^3 + 5*y^4 - 5720*x^2*y - 9003*x*y^2 - 4290*x^2 - 2320*y^3 + 2258804*x*y + 2825245*y^2 + 1695147*x
g(x,y).expand() 
       
-2*x^3*y + 7*x^2*y^2 - 2*x*y^3 + x^3 + 7*y^4 + 2860*x^2*y - 9081*x*y^2 -
1430*x^2 - 3248*y^3 - 1130562*x*y + 3955343*y^2 + 565049*x
-2*x^3*y + 7*x^2*y^2 - 2*x*y^3 + x^3 + 7*y^4 + 2860*x^2*y - 9081*x*y^2 - 1430*x^2 - 3248*y^3 - 1130562*x*y + 3955343*y^2 + 565049*x
f_x( x, y ) = diff( f(x,y), x) 
       
f_x() 
       
(4*y + 3)*((y - 232)^2 + (x - 715)^2) + 2*(x - 715)*(4*x*y + 5*y^2 +
3*x)
(4*y + 3)*((y - 232)^2 + (x - 715)^2) + 2*(x - 715)*(4*x*y + 5*y^2 + 3*x)
g_x( x, y ) = diff( g(x,y), x) 
       
g_x() 
       
-(2*y - 1)*((y - 232)^2 + (x - 715)^2) - 2*(x - 715)*(2*x*y - 7*y^2 - x)
-(2*y - 1)*((y - 232)^2 + (x - 715)^2) - 2*(x - 715)*(2*x*y - 7*y^2 - x)
f_y( x,y ) = diff( f(x,y), y) 
       
f_y() 
       
2*(y - 232)*(4*x*y + 5*y^2 + 3*x) + 2*(2*x + 5*y)*((y - 232)^2 + (x -
715)^2)
2*(y - 232)*(4*x*y + 5*y^2 + 3*x) + 2*(2*x + 5*y)*((y - 232)^2 + (x - 715)^2)
g_y( x,y ) = diff( g(x,y), y ) 
       
g_y() 
       
-2*(y - 232)*(2*x*y - 7*y^2 - x) - 2*(x - 7*y)*((y - 232)^2 + (x -
715)^2)
-2*(y - 232)*(2*x*y - 7*y^2 - x) - 2*(x - 7*y)*((y - 232)^2 + (x - 715)^2)
def my_jacobian( x, y ): nw = f_x( x, y) ne = f_y( x, y) sw = g_x( x, y) se = g_y( x, y) new_matrix = matrix( 2, 2, [ nw, ne, sw, se ] ) return new_matrix 
       
M1 = my_jacobian( 2, 3 ) 
       
M1 % 7 
       
[1 0]
[4 3]
[1 0]
[4 3]
f_x( 2, 3 ) 
       
8305200
8305200
f_y( 2, 3 ) 
       
21276430
21276430
g_x( 2, 3 ) 
       
-2879628
-2879628
g_y( 2, 3 ) 
       
21286506
21286506
det( M1 ) % 7 
       
3
3
def simple_inverse( the_matrix, the_n ): the_det = det( the_matrix ) % the_n if (the_det == 0): print "Sorry, I can't invert a matrix with determinant zero." return det_inverse = inverse_mod( the_det, the_n ) print "debug: det = ", print the_det print "debug: det_inverse = ", print det_inverse old_nw = the_matrix[0][0] old_ne = the_matrix[0][1] old_sw = the_matrix[1][0] old_se = the_matrix[1][1] new_nw = (old_se * det_inverse) % the_n new_ne = (- old_ne * det_inverse) % the_n new_sw = (- old_sw * det_inverse) % the_n new_se = (old_nw * det_inverse) % the_n inverse_matrix = matrix( 2, 2, [ new_nw, new_ne, new_sw, new_se ] ) return inverse_matrix 
       
M2 = simple_inverse( M1, 7 ) 
       
debug: det =  3
debug: det_inverse =  5
debug: det =  3
debug: det_inverse =  5
M2 
       
[1 0]
[1 5]
[1 0]
[1 5]
(M1*M2) % 7 
       
[1 0]
[0 1]
[1 0]
[0 1]
(M2*M1) % 7 
       
[1 0]
[0 1]
[1 0]
[0 1]
f( 2, 3 ) % 7 
       
4
4
g( 2, 3) % 7 
       
6
6
v_vec = matrix( 2, 1, [2, 3] ) 
       
v_vec 
       
[2]
[3]
[2]
[3]
w_vec = matrix( 2, 1, [4, 6]) 
       
w_vec 
       
[4]
[6]
[4]
[6]
y_vec = v_vec - M2*w_vec 
       
y_vec % 7 
       
[5]
[4]
[5]
[4]
f(5, 4) % 7 
       
0
0
g(5, 4) % 7 
       
0
0
f(5, 4) % 49 
       
14
14
g(5, 4) % 49 
       
14
14
f_x(5, 4) % 49 
       
48
48
f_y(5, 4) % 49 
       
30
30
g_x(5, 4) % 49 
       
0
0
g_y(5, 4) % 49 
       
23
23
M3 = my_jacobian( 5, 4) % 49 
       
M3 % 49 
       
[48 30]
[ 0 23]
[48 30]
[ 0 23]
det(M3) % 49 
       
26
26
M4 = simple_inverse( M3, 49 ) 
       
debug: det =  26
debug: det_inverse =  17
debug: det =  26
debug: det_inverse =  17
M4 
       
[48 29]
[ 0 32]
[48 29]
[ 0 32]
(M3 * M4) % 49 
       
[1 0]
[0 1]
[1 0]
[0 1]
(M4 * M3) % 49 
       
[1 0]
[0 1]
[1 0]
[0 1]
(31*19) % 49 
       
1
1
v_vec2 = matrix(2, 1, [5, 4] ) 
       
w_vec2 = matrix(2, 1, [14, 14] ) 
       
y_vec2 = v_vec2 - M4*w_vec2 
       
y_vec2 % 49 
       
[ 5]
[46]
[ 5]
[46]
 
       
y_vec2 = v_vec2 - M4*w_vec2 
       
y_vec2 % 49 
       
[ 5]
[46]
[ 5]
[46]
 
       
f(5 , 46) % 49 
       
0
0
g( 5, 46) % 49 
       
0
0