The code had to be arranged (lots of work, please insert correctly the code, so that it can be pasted, and so that it works).

• The implemented algorithms did not came as a function, instead, there was code doing the job for some, but not all cases. - Also, there was no description of the algorithms, so there is also a small shape of mathematics and understanding what the code does involved in between.
• if there is a simpler (structural) way to post the code, please do it.
• The following is not optimal, it was partially typed in train, then bug-fixed (not improved) till it worked. Please rearrange it for the own needs.
• Using globals is not a good solution, please rewrite for the own needs. (Depending on needs, either speed or solid style are relevant.)

Here is the code, that did the job for me.

# # # ==============================
# GLOBALS
p = 71
F = GF(p)

A = F(41)
B = F(18)

E = EllipticCurve( F, [A, B] )
ZERO = E(0)

S.<W> = PolynomialRing( F )
K.<a> = GF( p^2, modulus=W^2+W+1 )

R.<z> = PolynomialRing( K, sparse=True )

x = (a + 2  )/3 * (z^p-a*z)
y = (1 + 2*a)/3 * (z^p-  z)
f = y^2 - ( x^3 + A*x + B )

# The following routines are not optimal, please improve
"""To be sure that we always have two coefficients,
also for polynomials (of degree one)
with the one or the other vanishing coefficient.
"""
coefs = pol.coefficients( sparse=False )
return coefs + ( ( padto-len(coefs) ) * [ F(0), ] )
return coefs

"""
if P == ZERO:
return Q
if Q == ZERO:
return P

if P[0] == Q[0]:
if P[1] == Q[1]:
return double(P)
elif P[1] == - Q[1]:
return ZERO
else:
return "*** INTERNAL ERROR"

g = ( + ( Q[0] - P[0] ) * y
- ( Q[1] - P[1] ) * x
- ( P[1]*Q[0] - Q[1]*P[0] ) )

f1 = ( f % g ).monic()
f2 = ( z - P[0] - P[1]*a ) * ( z - Q[0] - Q[1]*a )
f3 = f1 // f2
c  = f3.constant_coefficient()
M0, M1 = get_coefficients( S(c) )

return E.point( ( -M0, M1 ) )

def double( P ):
"""
# print "doubling", P
if P == ZERO or P[1] == 0:
return ZERO

g = ( 2 * P[1] * y
- ( 3 * P[0]^2 + A ) * x
- ( -P[0]^3 + A * P[0] + 2*B ) )

f1 = ( f % g ).monic()
f2 = ( z - P[0] - P[1]*a ) * ( z - P[0] - P[1]*a )
f3 = f1 // f2
c  = f3.constant_coefficient()
M0, M1 = get_coefficients( S(-c) )

return E.point( ( M0, -M1 ) )

def mult( P, k ):
"""ad-hoc multiplication P -> k*P on the elliptic curve E
"""
# if k not in ZZ:
#     return "*** ERROR"
# print "computing %s * %s" % ( k, str(P) )
if k == 0:
return P.parent(0)
if k<0:
return -mult(P, -k)
rest = k%2
khalf = k//2
if khalf:
Q = double( mult(P, khalf) )
else:
Q = P.parent(0)
if rest:
return Q

# # # ================================================
print "TEST SUITE:"
k   = 21378964192843
print "k =", k
count = 0
for P in E.points():
count += 1
Q1, Q2 = mult(P, k), k*P
print ( "%2d Point P = %13s | mult(P,k) = %13s | k*P = %13s | %s"
% ( count, P, Q1, Q2, bool(Q1==Q2) ) )

if Q1 != Q2:

else:
print "OK"


Results for the tests:

TEST SUITE:
k = 21378964192843
1 Point P =   (0 : 1 : 0) | mult(P,k) =   (0 : 1 : 0) | k*P =   (0 : 1 : 0) | True
2 Point P =  (0 : 35 : 1) | mult(P,k) = (43 : 63 : 1) | k*P = (43 : 63 : 1) | True
3 Point P =  (0 : 36 : 1) | mult(P,k) =  (43 : 8 : 1) | k*P =  (43 : 8 : 1) | True
4 Point P =  (1 : 29 : 1) | mult(P,k) =  (0 : 35 : 1) | k*P =  (0 : 35 : 1) | True
5 Point P =  (1 : 42 : 1) | mult(P,k) =  (0 : 36 : 1) | k*P =  (0 : 36 : 1) | True
6 Point P =  (2 : 26 : 1) | mult(P,k) = (27 : 17 : 1) | k*P = (27 : 17 : 1) | True
7 Point P =  (2 : 45 : 1) | mult(P,k) = (27 : 54 : 1) | k*P = (27 : 54 : 1) | True
8 Point P =   (5 : 8 : 1) | mult(P,k) = (68 : 62 : 1) | k*P = (68 : 62 : 1) | True
9 Point P =  (5 : 63 : 1) | mult(P,k) =  (68 : 9 : 1) | k*P =  (68 : 9 : 1) | True
10 Point P =  (6 : 14 : 1) | mult(P,k) =   (5 : 8 : 1) | k*P =   (5 : 8 : 1) | True
11 Point P =  (6 : 57 : 1) | mult(P,k) =  (5 : 63 : 1) | k*P =  (5 : 63 : 1) | True
12 Point P =   (7 : 3 : 1) | mult(P,k) = (45 : 64 : 1) | k*P = (45 : 64 : 1) | True
13 Point P =  (7 : 68 : 1) | mult(P,k) =  (45 : 7 : 1) | k*P =  (45 : 7 : 1) | True
14 Point P =  (8 : 19 : 1) | mult(P,k) = (50 : 64 : 1) | k*P = (50 : 64 : 1) | True
15 Point P =  (8 : 52 : 1) | mult(P,k) =  (50 : 7 : 1) | k*P =  (50 : 7 : 1) | True
16 Point P = (10 : 24 : 1) | mult(P,k) = (24 : 49 : 1) | k*P = (24 : 49 : 1) | True
17 Point P = (10 : 47 : 1) | mult(P,k) = (24 : 22 : 1) | k*P = (24 : 22 : 1) | True
18 Point P =  (11 : 5 : 1) | mult(P,k) =  (1 : 42 : 1) | k*P =  (1 : 42 : 1) | True
19 Point P = (11 : 66 : 1) | mult(P,k) =  (1 : 29 : 1) | k*P =  (1 : 29 : 1) | True
20 Point P = (12 : 26 : 1) | mult(P,k) = (38 : 43 : 1) | k*P = (38 : 43 : 1) | True
21 Point P = (12 : 45 : 1) | mult(P,k) = (38 : 28 : 1) | k*P = (38 : 28 : 1) | True
22 Point P = (13 : 11 : 1) | mult(P,k) = (15 : 48 : 1) | k*P = (15 : 48 : 1) | True
23 Point P = (13 : 60 : 1) | mult(P,k) = (15 : 23 : 1) | k*P = (15 : 23 : 1) | True
24 Point P = (15 : 23 : 1) | mult(P,k) =  (56 : 2 : 1) | k*P =  (56 : 2 : 1) | True
25 Point P = (15 : 48 : 1) | mult(P,k) = (56 : 69 : 1) | k*P = (56 : 69 : 1) | True
26 Point P = (17 : 27 : 1) | mult(P,k) = (12 : 26 : 1) | k*P = (12 : 26 : 1) | True
27 Point P = (17 : 44 : 1) | mult(P,k) = (12 : 45 : 1) | k*P = (12 : 45 : 1) | True
28 Point P = (21 : 22 : 1) | mult(P,k) = (51 : 59 : 1) | k*P = (51 : 59 : 1) | True
29 Point P = (21 : 49 : 1) | mult(P,k) = (51 : 12 : 1) | k*P = (51 : 12 : 1) | True
30 Point P =  (23 : 8 : 1) | mult(P,k) = (26 : 49 : 1) | k*P = (26 : 49 : 1) | True
31 Point P = (23 : 63 : 1) | mult(P,k) = (26 : 22 : 1) | k*P = (26 : 22 : 1) | True
32 Point P = (24 : 22 : 1) | mult(P,k) = (10 : 47 : 1) | k*P = (10 : 47 : 1) | True
33 Point P = (24 : 49 : 1) | mult(P,k) = (10 : 24 : 1) | k*P = (10 : 24 : 1) | True
34 Point P = (25 : 14 : 1) | mult(P,k) = (48 : 16 : 1) | k*P = (48 : 16 : 1) | True
35 Point P = (25 : 57 : 1) | mult(P,k) = (48 : 55 : 1) | k*P = (48 : 55 : 1) | True
36 Point P = (26 : 22 : 1) | mult(P,k) = (13 : 60 : 1) | k*P = (13 : 60 : 1) | True
37 Point P = (26 : 49 : 1) | mult(P,k) = (13 : 11 : 1) | k*P = (13 : 11 : 1) | True
38 Point P = (27 : 17 : 1) | mult(P,k) = (36 : 15 : 1) | k*P = (36 : 15 : 1) | True
39 Point P = (27 : 54 : 1) | mult(P,k) = (36 : 56 : 1) | k*P = (36 : 56 : 1) | True
40 Point P = (28 : 16 : 1) | mult(P,k) =  (47 : 7 : 1) | k*P =  (47 : 7 : 1) | True
41 Point P = (28 : 55 : 1) | mult(P,k) = (47 : 64 : 1) | k*P = (47 : 64 : 1) | True
42 Point P =  (29 : 6 : 1) | mult(P,k) = (67 : 28 : 1) | k*P = (67 : 28 : 1) | True
43 Point P = (29 : 65 : 1) | mult(P,k) = (67 : 43 : 1) | k*P = (67 : 43 : 1) | True
44 Point P = (32 : 35 : 1) | mult(P,k) = (17 : 27 : 1) | k*P = (17 : 27 : 1) | True
45 Point P = (32 : 36 : 1) | mult(P,k) = (17 : 44 : 1) | k*P = (17 : 44 : 1) | True
46 Point P = (35 : 33 : 1) | mult(P,k) = (29 : 65 : 1) | k*P = (29 : 65 : 1) | True
47 Point P = (35 : 38 : 1) | mult(P,k) =  (29 : 6 : 1) | k*P =  (29 : 6 : 1) | True
48 Point P = (36 : 15 : 1) | mult(P,k) = (37 : 28 : 1) | k*P = (37 : 28 : 1) | True
49 Point P = (36 : 56 : 1) | mult(P,k) = (37 : 43 : 1) | k*P = (37 : 43 : 1) | True
50 Point P = (37 : 28 : 1) | mult(P,k) = (40 : 57 : 1) | k*P = (40 : 57 : 1) | True
51 Point P = (37 : 43 : 1) | mult(P,k) = (40 : 14 : 1) | k*P = (40 : 14 : 1) | True
52 Point P = (38 : 28 : 1) | mult(P,k) = (39 : 36 : 1) | k*P = (39 : 36 : 1) | True
53 Point P = (38 : 43 : 1) | mult(P,k) = (39 : 35 : 1) | k*P = (39 : 35 : 1) | True
54 Point P = (39 : 35 : 1) | mult(P,k) =  (23 : 8 : 1) | k*P =  (23 : 8 : 1) | True
55 Point P = (39 : 36 : 1) | mult(P,k) = (23 : 63 : 1) | k*P = (23 : 63 : 1) | True
56 Point P = (40 : 14 : 1) | mult(P,k) = (63 : 32 : 1) | k*P = (63 : 32 : 1) | True
57 Point P = (40 : 57 : 1) | mult(P,k) = (63 : 39 : 1) | k*P = (63 : 39 : 1) | True
58 Point P =  (42 : 0 : 1) | mult(P,k) =  (42 : 0 : 1) | k*P =  (42 : 0 : 1) | True
59 Point P =  (43 : 8 : 1) | mult(P,k) = (52 : 30 : 1) | k*P = (52 : 30 : 1) | True
60 Point P = (43 : 63 : 1) | mult(P,k) = (52 : 41 : 1) | k*P = (52 : 41 : 1) | True
61 Point P =  (45 : 7 : 1) | mult(P,k) = (35 : 33 : 1) | k*P = (35 : 33 : 1) | True
62 Point P = (45 : 64 : 1) | mult(P,k) = (35 : 38 : 1) | k*P = (35 : 38 : 1) | True
63 Point P =  (47 : 7 : 1) | mult(P,k) =  (8 : 19 : 1) | k*P =  (8 : 19 : 1) | True
64 Point P = (47 : 64 : 1) | mult(P,k) =  (8 : 52 : 1) | k*P =  (8 : 52 : 1) | True
65 Point P = (48 : 16 : 1) | mult(P,k) = (58 : 46 : 1) | k*P = (58 : 46 : 1) | True
66 Point P = (48 : 55 : 1) | mult(P,k) = (58 : 25 : 1) | k*P = (58 : 25 : 1) | True
67 Point P =  (50 : 7 : 1) | mult(P,k) = (57 : 26 : 1) | k*P = (57 : 26 : 1) | True
68 Point P = (50 : 64 : 1) | mult(P,k) = (57 : 45 : 1) | k*P = (57 : 45 : 1) | True
69 Point P = (51 : 12 : 1) | mult(P,k) = (25 : 57 : 1) | k*P = (25 : 57 : 1) | True
70 Point P = (51 : 59 : 1) | mult(P,k) = (25 : 14 : 1) | k*P = (25 : 14 : 1) | True
71 Point P = (52 : 30 : 1) | mult(P,k) = (28 : 55 : 1) | k*P = (28 : 55 : 1) | True
72 Point P = (52 : 41 : 1) | mult(P,k) = (28 : 16 : 1) | k*P = (28 : 16 : 1) | True
73 Point P =  (56 : 2 : 1) | mult(P,k) = (32 : 36 : 1) | k*P = (32 : 36 : 1) | True
74 Point P = (56 : 69 : 1) | mult(P,k) = (32 : 35 : 1) | k*P = (32 : 35 : 1) | True
75 Point P = (57 : 26 : 1) | mult(P,k) =  (11 : 5 : 1) | k*P =  (11 : 5 : 1) | True
76 Point P = (57 : 45 : 1) | mult(P,k) = (11 : 66 : 1) | k*P = (11 : 66 : 1) | True
77 Point P = (58 : 25 : 1) | mult(P,k) = (21 : 49 : 1) | k*P = (21 : 49 : 1) | True
78 Point P = (58 : 46 : 1) | mult(P,k) = (21 : 22 : 1) | k*P = (21 : 22 : 1) | True
79 Point P = (63 : 32 : 1) | mult(P,k) = (64 : 13 : 1) | k*P = (64 : 13 : 1) | True
80 Point P = (63 : 39 : 1) | mult(P,k) = (64 : 58 : 1) | k*P = (64 : 58 : 1) | True
81 Point P = (64 : 13 : 1) | mult(P,k) =  (6 : 14 : 1) | k*P =  (6 : 14 : 1) | True
82 Point P = (64 : 58 : 1) | mult(P,k) =  (6 : 57 : 1) | k*P =  (6 : 57 : 1) | True
83 Point P = (66 : 16 : 1) | mult(P,k) = (66 : 55 : 1) | k*P = (66 : 55 : 1) | True
84 Point P = (66 : 55 : 1) | mult(P,k) = (66 : 16 : 1) | k*P = (66 : 16 : 1) | True
85 Point P = (67 : 28 : 1) | mult(P,k) =   (7 : 3 : 1) | k*P =   (7 : 3 : 1) | True
86 Point P = (67 : 43 : 1) | mult(P,k) =  (7 : 68 : 1) | k*P =  (7 : 68 : 1) | True
87 Point P =  (68 : 9 : 1) | mult(P,k) =  (2 : 26 : 1) | k*P =  (2 : 26 : 1) | True
88 Point P = (68 : 62 : 1) | mult(P,k) =  (2 : 45 : 1) | k*P =  (2 : 45 : 1) | True
OK

more