算法流程如下：

# ============================================================================
# Copyright (c) 2021 Zhou Jin, Shandong University. All rights reserved.
# Elgamal.py - The core part of this algorithm
# Last edited time: 2021/10/25 16:45
#
# THIS PROGRAM IS FREE SOFTWARE
# If you have any question, email elford@foxmail.com
# ============================================================================

import random


# Find the greatest common divisor
def gcd(a, b):
    if a < b:
        return gcd(b, a)
    elif a % b == 0:
        return b
    else:
        return gcd(b, a % b)


# Quick power & modulo
def power(a, b, c):
    ans = 1
    while b != 0:
        if b & 1:
            ans = (ans * a) % c
        b >>= 1
        a = (a * a) % c
    return ans


# Lucas-lemmer sex test
def Lucas_Lehmer(num) -> bool:  # Quickly check whether POw (2,m)-1 is prime
    if num == 2:
        return True
    if num % 2 == 0:
        return False
    s = 4
    Mersenne = pow(2, num) - 1  # Pow (2, num)-1 is a Mersenne number
    for x in range(1, (num - 2) + 1):  # Num-2 is the number of cycles, and +1 means the right interval is open
        s = ((s * s) - 2) % Mersenne
    if s == 0:
        return True
    else:
        return False


# Detection of large prime numbers
def Miller_Rabin(n):
    a = random.randint(2,n-2) # Randomly find 'a' which is belong to [2,n-2]
    s = 0                     # S is the power of factor 2 in d
    d = n - 1
    while (d & 1) == 0:       # Let's factor out all the 2's in d.
        s += 1
        d >>= 1

    x = power(a, d, n)
    for i in range(s):        # Perform s secondary probes
        newX = power(x, 2, n)
        if newX == 1 and x != 1 and x != n - 1:
            return False      # Using the inverse of the quadratic theorem, n is determined to be composite.
        x = newX

    if x != 1:                # If x=a^(n-1) (mod n), then n is determined to be composite.
        return False

    return True               # Judge by the converse of Fermat's little theorem. The number that survives this test
                              # is most likely prime.


# Extended Euclidean algorithm, ab=1 (mod m), yields the multiplicative inverse b of A under module m
def Extended_Eulid(a: int, m: int) -> int:
    def extended_eulid(a: int, m: int):
        if a == 0:  # boundary condition
            return 1, 0, m
        else:
            x, y, gcd = extended_eulid(m % a, a)  # recursion
            x, y = y, (x - (m // a) * y)  # Recursively, the left end is the upper layer
            return x, y, gcd              # Returns the result of the first level calculation.
        # The final return y value is the multiplication inverse of b in mode a
        # If y is complex, y plus a is the corresponding positive inverse

    n = extended_eulid(a, m)
    if n[1] < 0:
        return n[1] + m
    else:
        return n[1]


# Generate the field parameter p, approximately 512bits in length
def Generate_p() -> int:
    a = random.randint(10**150, 10**160)
    while gcd(a, 2) != 1:
        a = random.randint(10**150, 10**160)
    return a


# Generate the field parameter alpha
def Generate_alpha(p: int) -> int:
    return random.randint(2, p)


# Generate a prime number less than p, approximately 512bits long, as the private key
def Generate_private_key(p: int) -> int:
    pri = random.randint(2, p - 2)
    while gcd(pri, p) != 1:
        pri = random.randint(2, p - 2)
    return pri


# Bob or Alice uses the field parameter P and the resulting mask to encrypt the message
def encrypt(message, p, Key_mask) -> int:
    ciphertext = (message * Key_mask) % p
    return ciphertext


# Bob or Alice decrypts the ciphertext using the field parameter P and the masks obtained by them
def decrypt(ciphertext, p, Key_mask):
    # Inverse mask
    inverse_element = Extended_Eulid(Key_mask, p)
    # Decode
    plaintext = (ciphertext * inverse_element) % p
    return plaintext


def quick_power(a: int, b: int) -> int:
    ans = 1
    while b != 0:
        if b & 1:
            ans = ans * a
        b >>= 1
        a = a * a
    return ans


def Generate_prime(key_size: int) -> int:
    while True:
        num = random.randrange(quick_power(2, key_size - 1), quick_power(2, key_size))
        if Miller_Rabin(num):
            return num


if __name__ == '__main__':
    message = int(input("Message:       "))
    if type(message) != int:
        raise ValueError("Must be an integer!")

    p = Generate_prime(512)
    alpha = Generate_alpha(p)

    # Alice selects her own private key, calculates her own public key, and passes the public key to Bob
    Private_Alice = Generate_private_key(p)
    Public_Alice = power(alpha, Private_Alice, p)

    # Bob selects his private key, calculates his public key, and places the public key on the network
    Private_Bob = Generate_private_key(p)
    Public_Bob = power(alpha, Private_Bob, p)

    # Bob and Alice calculate the mask key separately, and the results should be the same
    Key_mask_Alice = power(Public_Alice, Private_Bob, p)
    Key_mask_Bob = power(Public_Bob, Private_Alice, p)

    # Bob encrypts the message
    ciphertext = encrypt(message, p, Key_mask_Bob)

    # Alice decrypts the message
    plaintext = decrypt(ciphertext, p, Key_mask_Alice)

    print("Domain parameters: ")
    print("p:            ", p)
    print("alpha:        ", alpha)
    print()

    print("Private Keys: ")
    print("Alice:        ", Private_Alice)
    print("Bob:          ", Private_Bob)
    print()

    print("Public keys: ")
    print("Alice:       ", Public_Alice)
    print("Bob:         ", Public_Bob)
    print()

    print("Masked key:  ", Key_mask_Bob)
    print()

    print("Ciphertext:  ", ciphertext)
    print("Plaintext:   ", plaintext)