The function "sqrtp" computes a square root of x 
modulo a prime p using the Tonelli-Shanks algorithm.
It is assumed that x is a square modulo p.

Library:  function power (x, n: longint): longint;
		(* computes x^n modulo p *)


function sqrtp(var x: longint):longint;
		var
			a, b, c, g, h, hh, q, k, kk, i: longint;

	begin
		p := abs(p);
		a := x mod p;
		q := p - 1;
		k := 0;
		repeat
			k := k + 1;
			q := q div 2
		until odd(q);
		if k = 1 then
			g := -1
		else
			begin
				g := 1;
				repeat
					g := g + 1;
					h := power(g, q);
					hh := power(h, ((p - 1) div q) div 2)
				until hh = (p - 1);
				g := h
			end;
                                                  (* g generates the 2-part *)
		b := power(a, (q - 1) div 2);
		hh := (b*a) mod p;
		c := (b*hh) mod p;
		b := hh;
		while c <> 1 do
			begin
				h := c;
				kk := 0;
				repeat
					kk := kk + 1;
					h := (h*h) mod p
				until h = 1;
				h := 1;
				for i := 1 to (k - kk) do
					h := h * 2;
				g := power(g, (h div 2));
				b := (b*g) mod p;
				g := (g*g) mod p;
				c := (c*g) mod p;
				k := kk
			end;
		sqrtp := b
	end; (*sqrtp*)