CRYPTOSYSTEM AND METHOD WITH EFFICIENT ELLIPTIC CURVE OPERATORS FOR AN EXTRACTION OF EiSi COORDINATE SYSTEM

Abstract

A system, method and computer-readable medium provide secure communication between a first and a second computer system based on supersingular isogeny elliptic curve cryptography. The first computer system and the second computer system each determine kernels KA and KB including computing mP+nQ by accessing a lookup table stored in a memory that contains a range of doubles of an end point of the respective kernels, where P and Q are points on the public elliptic curve and m and n are integers. The first computer system and the second computer system compute secret isogenies by determining a respective kernel KBA and KAB using mixed-base multiplicands with a single inversion, including computing the respective kernel KBA and KAB by converting the multiplicands to base 32, and computing scalar multiplications using the base 32 multiplicands.

Description

STATEMENT OF ACKNOWLEDGEMENT

The inventor would like to acknowledge the support provided by the Deanship of Scientific Research at Umm Al-Qura University for supporting this work by Grant Code #20UQU0026DSR.

BACKGROUND

Field of the Invention

The disclosure relates to Supersingular Isogeny Diffie-Helman cryptosystem for secure electronic communication between two parties, and in particular, a system and method that minimizes the number of multiplications and inversions in computing the kernel of the elliptic curve.

Background

An objective of a cryptosystem is to minimize complexity of the cryptographic algorithm while making it computationally difficult to detect a secure key that results from the cryptographic algorithm. Efficiency of secure communication that is impacted by the complexity of the algorithm is especially important to end users that desire fast communication of messages despite security risks. Complexity of a cryptographic algorithm may be minimized by optimizing factors including reducing the number of computations and reducing the amount of computer resources, such as storage space, needed to perform the algorithm. For example, in the case of key exchange for secure communications, reducing the total data size stored in a communications device and the total computation amount by a communications device are important factors. At the same time, security of communication continues to be an ever increasing challenge. For example, specialized computer systems are being developed that are better able to compromise existing cryptosystem algorithms, i.e., can determine the secure key and/or decipher the encrypted message within the regular time frame of message transmission. Attacks that can compromise cryptosystems are not necessarily limited to independent calculations and might also be based on analyzing hardware behaviors such as the case with the Side Channel Attack (SCA). There is an ongoing need to develop cryptosystems that can withstand attacks by future specialized computer systems.

Secure Internet-based communications generally rely on public-key cryptography, which allows entities to communicate without the need for sharing confidential material in advance. Elliptic Curve Cryptography (ECC), proposed in 1985, is still a predominant type of public-key cryptography. See L. B. Oliveira, F. M. Q. Pereira, R. Misoczki, D. F. Aranha, F. Borges, M. Nogueira, M. Wangham, M. Wu, and J. Liu, “The computer for the 21st century: present security & privacy challenges,” Journal of Internet Services and Applications, vol. 9, no. 1, p. 24, 2018, incorporated herein by reference in its entirety. ECC is one of the existing cryptosystem algorithms that is commonly used for encrypted emails, online banking, secure ecommerce websites, digital signatures, and other data transfer applications where the size of the storage space for public keys is an issue. Breaching these would have significant effects on society. The adoption of ECC has been accelerated by recommendations from an array of standardization entities, including, NIST, IETF, and ANSI (NIST, 2016). Compared to competitors like RSA and Elgamal, elliptic curve cryptography introduced some of the most efficient public key cryptosystems (PKC) for desirable security. More recently, while there are known quantum and classical attacks that can breach cryptographic protocols based on supersingular isogeny graphs (SIGs), the Supersingular isogeny Diffie-Hellman (SIDH) technique has been found to be able handle quantum-based attacks.

Specialized computer systems, such as quantum computers, are being developed and can break elliptic curve-based cryptosystems. Developments in quantum computer systems threaten to break elliptic curve and factoring techniques for public key cryptography.

Although quantum computers are currently in their infancy, the ongoing development of quantum computers and their theoretical ability to compromise modern cryptographic protocols (such as TLS/SSL) has prompted the development of post-quantum cryptography.

A Supersingular Isogeny Diffie-Hellman (SIDH) key exchange has been developed to withstand attacks by techniques that use quantum computers. Subsequently, SIDH is one of the post-quantum cryptographic algorithms that can offer secure key exchanges between communicating entities over insecure communication channels. See S. Arpin, C. Camacho-Navarro, K. Lauter, J. Lim, K. Nelson, T. Scholl, and J. Sotáková, “Adventures in supersingularland,” arXiv preprint arXiv: 1909.07779, 2019; and C. Costello, P. Longa, and M. Naehrig, “Efficient algorithms for supersingular isogeny diffie-hellman,” in Annual International Cryptology Conference. Springer, 2016, pp. 572-601, each incorporate herein by reference in their entirety.

In algebraic geometry, supersingular elliptic curves form a certain class of elliptic curves over a field of characteristic p>0 with unusually large endomorphism rings. The term “supersingular” comes from the phrase “singular values of the j-invariant” used for values of the j-invariant for which a complex elliptic curve has complex multiplication. The complex elliptic curves with complex multiplication are those for which the endomorphism ring has the maximal possible rank 2. In positive characteristic it is possible for the endomorphism ring to be even larger: it can be an order in a quaternion algebra of dimension 4, in which case the elliptic curve is supersingular.

The core operations for SIDH are the computation of the isogeny and of its kernel. Basically, Velu's formula is used to compute the isogeny, and the P+k[Q] formula is used to compute the kernel. In elliptic curve point multiplication, P and Q are initial and final points on the curve, and k is the size of the underlying field. In key exchange, P and Q are points on the curve and k is the secret key that is generated by both parties. See D. Jao and L. De Feo, “Towards quantum-resistant cryptosystems from supersingular elliptic curve isogenies,” in International Workshop on Post-Quantum Cryptography. Springer, 2011, pp. 19-34, incorporated herein by reference in its entirety. The complexity of SIDH relies on the difficulty of finding isogenies besides computing the scalar multiplication in the kernel formula. Point multiplication mainly relies on elliptic curve point addition and point doubling operations. Thus, speeding up elliptic curve (EC) computations will not only benefit the applications that rely on ECC, but also has an effective impact on the post quantum cryptosystem SIDH. Moreover, attacks are not necessarily limited to independent calculations and might also be based on analyzing hardware behaviors such as with the Side Channel Attack (SCA). As attackers analyze electrical power consumption patterns, which differ between performing point addition or point doubling, they can recover the secret key. Therefore, the development criteria for EC algorithms and systems can include aspects other than speed. For instance, the addition in the Montgomery coordinates system is resistant to such attacks. See P. L. Montgomery, “Speeding the pollard and elliptic curve methods of factorization,” Mathematics of computation, vol. 48, no. 177, pp. 243-264,1987, incorporated herein by reference in its entirety. A drawback of the Montgomery coordinates system is that it is slower than other coordinates systems such as Projective and Jacobian. See L. C. Washington, Elliptic curves: number theory and cryptography. Chapman and Hall/CRC, 2008, incorporated herein by reference in its entirety.

There is a need for effective algorithms for EC systems that provide computational speedup and resistance to side-channel attacks.

SUMMARY

An aspect is a system for secure communication. The system can include a first computer system, a communication network, and a second computer system. The first computer system and the second computer system are each provided with a public elliptic curve E. The first computer system and the second computer system each independently determine kernels K_A and K_B, respectively, and an isogeny mapping using the public elliptic curve E. The first computer system and the second computer system each independently determine mapped points based on the respective isogeny mappings. The first computer system and the second computer system each exchange the determined mapped points. The first computer system and the second computer system each independently compute the coefficients of two different secret elliptic curves under two different secret isogenies. The first computer system and the second computer system each independently compute a common secret key using the secret elliptic curves and the secret isogenies. The first computer system and the second computer system exchange messages using the common secret key. The first computer system and the second computer system each determine kernels K_A and K_B including computing mP+nQ by accessing a lookup table stored in a memory that contains a range of doubles of an end point of the respective kernels, where P and Q are points on the public elliptic curve and m and n are integers.

A further aspect is a supersingular isogeny-based cryptography method. The method can include providing a first computer system and a second computer system with a public elliptic curve E, independently determining, by the first computer system and the second computer system, kernels K_A and K_B, respectively, and an isogeny mapping using the public elliptic curve E, independently determining, by the first computer system and the second computer system, mapped points based on the respective isogeny mappings, exchanging, between the first computer system and the second computer system, the determined mapped points, independently computing, by the first computer system and the second computer system, coefficients of two different secret elliptic curves under two different secret isogenies, independently determining, by the first computer system and the second computer system, a common secret key using the secret elliptic curves and the secret isogenies, and exchanging, between the first computer system and the second computer system, messages using the common secret key. The first computer system and the second computer system each determine kernels K_A and K_B including computing mP+nQ by accessing a lookup table stored in a memory that contains a range of doubles of an end point of the respective kernels, where P and Q are points on the public elliptic curve and m and n are integers.

A further aspect is a non-transitory computer-readable medium storing instructions that are operable when executed by processing circuitry to perform operations. The operations can include receiving a public elliptic curve E, independently determining kernels K_A and K_B, respectively, and an isogeny mapping using the public elliptic curve E, independently determining mapped points based on the respective isogeny mappings, exchanging the determined mapped points, independently computing coefficients of two different secret elliptic curves under two different secret isogenies, independently determining a common secret key using the secret elliptic curves and the secret isogenies, and exchanging messages using the common secret key, and determining kernels K_A and K_B by computing mP+nQ by accessing a lookup table stored in a memory that contains a range of doubles of an end point of the respective kernels, where P and Q are points on the public elliptic curve and m and n are integers.

BRIEF DESCRIPTION OF THE DRAWINGS

A more complete appreciation of the invention and many of the attendant advantages thereof will be readily obtained as the same becomes better understood by reference to the following detailed description when considered in connection with the accompanying drawings, wherein:

FIG. 1 is a diagram of a communications system providing secure communication;

FIG. 2 is an exemplary computer system in the communications system of FIG. 1;

FIGS. 3A and 3B is a sequence diagram of Supersingular Isogeny Diffie-Hellman (SIDH) key exchange algorithm;

FIG. 4 is a flowchart of a method of computing isogeny mapping and kernel in the SIDH algorithm of FIGS. 3A and 3B;

FIG. 5 is a flowchart of a method of computing the kernel of FIG. 4 including a lookup table in accordance with an exemplary aspect of the disclosure;

FIG. 6 is a flowchart of a further method of computing the kernel of FIG. 4 using direct doubling in accordance with an exemplary aspect of the disclosure;

FIG. 7 is a flowchart of a further method of computing the kernel of FIG. 4 using direct doubling with a single inverse in accordance with an exemplary aspect of the disclosure;

FIG. 8 illustrates a cyclic group of the elliptic curve E: y²≡x³+2·x+2 mod 17;

FIG. 9 is a flowchart of a further method of computing the kernel of FIG. 4 using intermediate scalar multiplication algorithms in accordance with an exemplary aspect of the disclosure;

FIG. 10 is a flowchart of a further method of computing the kernal of FIG. 4 using extraction of coordinates in accordance with an exemplary aspect of the disclosure;

FIG. 11 illustrates a left-to-right algorithm;

FIG. 12 is a flowchart of a further method of computing the kernel of FIG. 4 using base 32 multiplicands in accordance with an exemplary aspect of the disclosure;

FIG. 13 illustrates computing the point 24P by using remi_func and remi_point functions; FIG. 14 illustrates computing the point 29P by using remi_func and remi_point functions;

FIG. 15 is a graph illustrating present work vs original algorithms in terms of number of multiplications;

FIG. 16 is a graph illustrating present work vs other coordinates algorithms in terms of number of multiplications;

FIG. 17 is a graph illustrating present work vs other coordinates algorithms in terms of number of maximum levels; and

FIG. 18 is a graph illustrating the maximum levels for different number of multipliers.

DETAILED DESCRIPTION

As mentioned above, secure Internet-based communications generally rely on public-key cryptography, which allows entities to communicate without the need for sharing confidential material in advance. FIG. 1 is a diagram of a communications system providing secure communication. In one or more embodiments, a first computer system 104 may communicate with a second computer system 108 over the Internet 102. The first computer system 104 may communicate over a secure communications channel 114. The second computer system 108 may communicate over a secure communication channel 118. A server computer 106 may be accessed through the Internet 102 by way of a communication channel 116. The server computer 106 may be a single computer system that can be accessed via communications protocol such as the Internet, or may be multiple server computers in interconnected arrangements, such as hierarchical, series or parallel. The first computer system 104 may communicate with the second computer system 108 securely by encrypting a message using a secret key independently computed by each computer system 104 and 108 using the Diffie-Hellman technique.

The first computer system 104 and the second computer system 108 may include mobile computing devices and/or desktop computing devices to be used by end users for communication by way of the Internet. Mobile computing devices are not limited to, but may include smartphones, tablet computers, laptop computers, handheld gaming devices and other portable computing devices. Desktop computers may range from all-in-one display-computing devices to computer workstations. Desktop computers may also include gaming stations, such as PlayStation, Nintendo switch, Microsoft Xbox, to name few.

In one implementation, the functions and processes of the first computer system 104 and the second computer system 108 may be implemented by one or more respective processing circuitry 226. Processing circuitry includes a programmed processor as a processor includes circuitry. Processing circuitry may also include devices such as an application specific integrated circuit (ASIC) and conventional circuit components arranged to perform the recited functions. Note that processing circuitry refers to a circuit or system of circuits. Herein, the processing circuitry may be in one computer system (as illustrated in FIG. 2) or may be distributed throughout a network of computer systems. Hence, the processing circuitry of the server computer system 106 for example may be in only one server or distributed among different servers/computers.

Next, a hardware description of the processing circuitry 226 according to exemplary embodiments is described with reference to FIG. 2. In FIG. 2, the processing circuitry 226 may include a Mobile Processing Unit (MPU) 200 which performs the processes described herein. The process data and instructions may be stored in memory 202. These processes and instructions may also be stored on a portable storage medium or may be stored remotely. The processing circuitry 226 may have a replaceable Subscriber Identity Module (SIM) 201 that contains information that is unique to the network service of the mobile device 130.

Further, the claimed advancements are not limited by the form of the computer-readable media on which the instructions of the inventive process are stored. For example, the instructions may be stored in FLASH memory, Secure Digital Random Access Memory (SDRAM), Random Access Memory (RAM), Read Only Memory (ROM), Programmable Read-Only Memory (PROM), Erasable Programmable Read-Only Memory (EPROM), Electrically Erasable Programmable Read Only Memory (EEPROM), solid-state hard disk or any other information processing device with which the processing circuitry 226 communicates, such as a server or computer.

Further, the claimed advancements may be provided as a utility application, background daemon, or component of an operating system, or combination thereof, executing in conjunction with MPU 200 and a mobile operating system such as Android, Microsoft® Windows® 10 Mobile, Apple iOS® and other systems known to those skilled in the art.

In order to achieve the processing circuitry 226, the hardware elements may be realized by various circuitry elements, known to those skilled in the art. For example, MPU 200 may be a Qualcomm mobile processor, a Nvidia mobile processor, a Atom® processor from Intel Corporation of America, a Samsung mobile processor, or a Apple A7 mobile processor, or may be other processor types that would be recognized by one of ordinary skill in the art. Alternatively, the MPU 200 may be implemented on an Field-Programmable Gate Array (FPGA), Application Specific Integrated Circuit (ASIC), Programmable Logic Device (PLD) or using discrete logic circuits, as one of ordinary skill in the art would recognize. Further, MPU 200 may be implemented as multiple processors cooperatively working in parallel to perform the instructions of the inventive processes described above.

The processing circuitry 226 in FIG. 2 also includes a network controller 206, such as an Intel Ethernet PRO network interface card from Intel Corporation of America, for interfacing with network 102. As can be appreciated, the network 102 can be a public network, such as the Internet, or a private network such as LAN or WAN network, or any combination thereof and can also include PSTN or ISDN sub-networks. The network 102 can also be wired, such as an Ethernet network. The processing circuitry 226 may include various types of communications processors for wireless communications including 3G, 4G and 5G wireless modems, WiFi®, Bluetooth®, GPS, or any other wireless form of communication that is known.

The processing circuitry 226 includes a Universal Serial Bus (USB) controller 225 which may be managed by the MPU 200.

The processing circuitry 226 further includes a display controller 208, such as a NVIDIA® GeForce® GTX or Quadro® graphics adaptor from NVIDIA Corporation of America for interfacing with display 210. An I/O interface 212 interfaces with buttons 214, such as for volume control. In addition to the I/O interface 212 and the display 210, the processing circuitry 226 may further include a microphone 241 and one or more cameras 231. The microphone 241 may have associated circuitry 240 for processing the sound into digital signals. Similarly, the camera 231 may include a camera controller 230 for controlling image capture operation of the camera 231. In an exemplary aspect, the camera 231 may include a Charge Coupled Device (CCD). The processing circuitry 226 may include an audio circuit 242 for generating sound output signals, and may include an optional sound output port.

The power management and touch screen controller 220 manages power used by the processing circuitry 226 and touch control. The communication bus 222, which may be an Industry Standard Architecture (ISA), Extended Industry Standard Architecture (EISA), Video Electronics Standards Association (VESA), Peripheral Component Interface (PCI), or similar, for interconnecting all of the components of the processing circuitry 226. A description of the general features and functionality of the display 210, buttons 214, as well as the display controller 208, power management controller 220, network controller 206, and I/O interface 212 is omitted herein for brevity as these features are known.

An aspect is effective algorithms that provide improvements to Elliptic Curve (EC) cryptography systems in several ways, both in terms of computational speedup and resistance to side-channel attacks.

FIGS. 3A and 3B is a sequence diagram of Supersingular Isogeny Diffie-Hellman (SIDH) key exchange algorithm for secure communications between the first computer system 104 and the second computer system 108. The algorithm is such that the first computer system 104 and the second computer system 108 each compute a common secret key which is used for the secure communications. Both computer systems 104 and 108 start by being provided with a public supersingular elliptic curve and they both end up with the same isogenous curve to it, but by performing different computations on an isogeny graph.

Initially, in S302, a trusted computer system 106 may determine a prime p derived from two different small primes l_A, l_B, respective exponents, integers e_A, e_B, and a small cofactor f, and a public supersingular elliptic curve E having two subgroups E[l_A^eA] and E[l_B^eB]. The prime p and public supersingular elliptic curve E may be transmitted to each of the first computer system 104 and the second computer system 108 over respective communication channels 114 and 118, via the Internet 102.

Each computer system 104 and 108 determines points, referred to as bases, on the public elliptic curve E and exchanges the computed points. In S312, first computer system 104 determines fixed elliptic points P_A, Q_A, of the subgroup E[l_A^eA], while in S322, second computer system 108 determines fixed elliptic points P_B, Q_B of the subgroup E[l_B^eB].

In S314, first computer system 104 determines integers (m_A and n_A), while in S324, second computer system 108 determines integers (m_B and n_B).

In S316, first computer system 104 computes a kernel K_A=m_A(P_A)+n_A(Q_A) and uses the point K_A and the Velu's formulas to create isogeny mapping ϕ_A: E->E_A. In S326, second computer 108 computes a kernel K_B=m_B(P_B)+n_B(Q_B) and uses the point K_B and the Velu's formulas to create isogeny mapping ϕ_B: E->E_B.

In S318, first computer system 104 applies ϕ_A to P_B and Q_B to form two points on E_A: ϕ_A(P_B) and ϕ_A(Q_B). In S328, second computer 108 applies ϕ_B to P_A and Q_A to form two points on E_B: ϕ_B(P_A) and ϕ_B(Q_A).

First computer system 104 communicates with second computer system 108 to publically exchange values (which may be referred to as the key exchange phase). In a key exchange phase, the system uses the two points as the basis for a kernel of a new isogeny, and computes the coefficients of two new elliptic curves under two new isogenies.

In S332, first computer system 104 transmits E_A, ϕ_A(P_B) and ϕ_A(Q_B) to second computer system 108. In S334, second computer system 108 transmits E_B, ϕ_B(P_A) and ϕ_B(Q_A) to first computer system 104.

Each computer system then performs further private computations.

In S342, first computer system 104 uses m_A, n_A, ϕ_B(P_A) and ϕ_B(Q_A) to form K_BA=m_A(ϕ_B(P_A))+n_A(ϕ_B(Q_A)). The first computer system 104 uses K_BA and Velu's formulas to create an isogeny mapping ψ_BA.

In a similar manner, in S352, second computer system 108 uses m_B, n_B, ϕ_A(P_B) and ϕ_A(Q_B) to form K_AB=m_B(ϕ_A(P_B))+n_B(ϕ_A(Q_B)). The second computer system 108 uses K_AB and Velu's formulas to create an isogeny mapping ψ_AB.

In S344, first computer system 104 uses ψ_BA to create elliptic curve E_BA.

In a similar manner, in S354, second computer system 108 uses ψ_AB to create elliptic curve E_AB.

In S346, first computer system 104 computes k=j-invariant of the curve E_BA.

In a similar manner, in S356, second computer system 108 computes k=j-invariant of the curve E_AB.

In S348 and S358, a message may be transmitted wherein a function of k is used as the common secret key.

Further enhancements to the Supersingular Isogeny Diffie-Hellman (SIDH) key exchange algorithm are provided to compute a kernel, in particular of the general form mP+nQ, where m and n are small integers, in a manner that optimizes implementation speed with the lowest cost among currently known algorithms using only one inversion.

An aspect is to efficiently compute the kernel by the first computer system 104 and the second computer system 108 by compensating the abscissa x and ordinate y equations of conventional EC systems in the higher doubling orders formulas and find a common factor between all slope's denominators in order to obtain a single inverse for each denominator. Unlike with other coordinates systems, the original affine operations with the Weierstrass elliptic curve form require computing an inverse each time point doubling or addition is performed, i.e., at every iteration of common fast EC scalar multiplication algorithms. See Washington. In general, finding inverses is much slower than big integer multiplication. Thus, as with Projective approaches, one goal is to eliminate inverses.

The most popular forms of public-key cryptography for current applications have increasingly been based on Elliptic Curves (ECs). See V. S. Miller, “Use of elliptic curves in cryptography,” in Conference on the theory and application of cryptographic techniques. Springer, 1985, pp. 417-426; and N. Koblitz, “Elliptic curve cryptosystems,”Mathematics of computation, vol. 48, no. 177, pp. 203-209,1987, each incorporated herein by reference in their entirety. With Elliptic Curve Cryptography (ECC), messages and secrets are mapped to points on an elliptic curve, which involves scalar point multiplication. However, computing a multiple of a point on an elliptic curve is of high cost. Scalar point multiplication may be performed more efficiently using a double-and-add algorithm. Point doubling and point addition operations define transitions between points. Scalar point multiplication can use such a sequence of point doubling and point addition operations to optimize repeated addition:

$Q=\left[k\right]P=\underset{\underset{k}{︸}}{P+P+\dots +P}.$

The double-and-add algorithm includes steps of:

convert the scalar k into a binary expansion;

calculate the points representing the powers of two multiplies of the generator and just sum those which form the number in binary form (i.e., have ones at given positions of the binary expansion). In other words, P's are doubled in a loop, and only those doubled P's at 1's positions are added.

Cryptosystems based on ECs rely on the difficulty of solving the Elliptic Curve Discrete Log (ECDL) problem. Namely, for elliptic curves with points P of large order and large k numbers, given the points Q and P in the previous equation it is hard to determine the scalar multiple k. However, with the expected emergence of quantum computers, in the near future cryptosystems whose security relies on the difficulty of ECDL will no longer be safe, since the scalar multiple may be easily recovered using Shor's algorithm. See L. Chen, L. Chen, S. Jordan, Y.-K. Liu, D. Moody, R. Peralta, R. Periner, and D. Smith-Tone, Report on post-quantum cryptography. US Department of Commerce, National Institute of Standards and Technology, 2016; and P. W. Shor, “Algorithms for quantum computation: Discrete logarithms and factoring,” in Proceedings 35th annual symposium on foundations of computer science. Ieee, 1994, pp. 124-134, each incorporated herein by reference in their entirety. Other quantum resilient schemes have been proposed. Furthermore, post-quantum cryptosystems such as Supersingular Isogeny Diffie-Hellman (SIDH) are slow techniques, and speeding up its elliptic curve computation is frequently mentioned as a significant goal.

As mentioned in the background section, a core operation for ECC is the scalar multiplication [k]P whose computation speed is seen as key to improving ciphers. For instant, Eisentrager et al proposed a method for computing the formula S=(2P+Q). See K. Eisenträger, K. Lauter, and P. L. Montgomery, “Fast elliptic curve arithmetic and improved weil pairing evaluation,” in Cryptographers' Track at the RSA Conference. Springer, 2003, pp. 343-354, incorporated herein by reference in its entirety. Their improved procedure saves a field multiplication, when compared to the original algorithm. Later, Ciet et al introduced a faster method for computing the same formula when a field inversion costs more than six field multiplications. See M. Ciet, M. Joye, K. Lauter, and P. L. Montgomery, “Trading inversions for multiplications in elliptic curve cryptography,” Designs, codes and cryptography, vol. 39, no. 2, pp. 189-206,2006, incorporated herein by reference in its entirety. Furthermore, they introduced an efficient method for computing point tripling. Mixed powers system of point doubling and tripling for computing the scalar multiplication was represented later by Dimitrov et al. See V. Dimitrov, L. Imbert, and P. K. Mishra, “Efficient and secure elliptic curve point multiplication using double-base chains,” in International Conference on the Theory and Application of Cryptology and Information Security. Springer, 2005, pp. 59-78, incorporated herein by reference in its entirety. Mishra et al. presented an efficient quintuple formula (5P) and introduced a mixed base algorithm with doubling and tripling. See P. K. Mishra and V. Dimitrov, “Efficient quintuple formulas for elliptic curves and efficient scalar multiplication using multibase number representation,” in International Conference on Information Security. Springer, 2007, pp. 390-406, incorporated herein by reference in its entirety. Further development was introduced by Longa and Miri by computing an efficient method for tripling and quintupling mixed with differential addition. See P. Longa and A. Miri, “New multibase non-adjacent form scalar multiplication and its application to elliptic curve cryptosystems (extended version).” IACR Cryptology ePrint Archive, vol. 2008, p. 52, 2008, incorporated herein by reference in its entirety. They proposed an efficient multibases nonadjacent representation (mbNAF) to reduce the cost. In Longa and Miri the same authors present further optimization in terms of cost for computing the form dP+Q. They have succeeded in implementing the previous forms of mixed double and add algorithm by using a single inversion when applying a new precomputation scheme. More recently, Purohit and Rawat used a multibase representation to propose an efficient scalar multiplication algorithm of doubling, tripling, and septupling, restricted on a non super-singular elliptic curve defined over the field F₂m. See G. Purohit and A. S. Rawat, “Fast scalar multiplication in ecc using the multi base number system,” International Journal of Computer Science Issues (IJCSI), vol. 8, no. 1, p. 131, 2011, incorporated herein by reference in its entirety. In addition, they have compared their work with other existing algorithms to achieve better representation in terms of cost. Therefore, speeding up the scalar multiplication computation in parallel with reducing the cost is a critical task.

Among all applications based on EC, Supersingular Isogeny Diffie-Helman (SIDH). provides a promising approach as a post-quantum cryptosystem Its main weakness is the slow elliptic curve computation speed. For elliptic curve schemes, the computation speed-up also favors attacks, which can however be compensated by increasing the size of the key. Isogeny-based cryptography also utilizes points on an elliptic curve, but its security is instead based on the difficulty of computing isogenies between elliptic curves. An isogeny can be thought of as a unique algebraic mapping between two elliptic curves that satisfy the group law. An algorithm for computing isogenies on ordinary curves in sub-exponential time was presented by Childs et al., rendering the use of cryptosystems based on isogenies on ordinary curves unsafe in the presence of quantum computers. See A. Childs, D. Jao, and V. Soukharev, “Constructing elliptic curve isogenies in quantum subexponential time,” Journal of Mathematical Cryptology, vol. 8, no. 1, pp. 1-29, 2014, incorporated herein by reference in its entirety. However, there is no known algorithm for computing isogenies on supersingular curves in sub-exponential time.

FIG. 4 is a flowchart of a method of computing isogeny mapping and kernel in the SIDH algorithm of FIGS. 3A and 3B (S316, S326, S342, S352). The core operations for SIDH includes, step S402, computing the isogeny using Velu's formula, and, step S404, computing its kernel using the P+k[Q] formula, where P and Q are points on the curve and k is the secret key (see Jao et al.). This operation must be performed in both phases of SIDH. First, this happens in the key generation phase (S316, S326), where the point is known in advance. FIG. 5 is a flowchart of a method of computing the kernel of FIG. 4 (S404) including a lookup table in accordance with an exemplary aspect of the disclosure. In S502, a lookup table stored in a memory of the first computer system 104 and of the second computer system 108 may be constructed that contains all doubles of point Q. The doubles can be reused when needed without having to compute the double. Second, in the key exchange phase (S342, S352), where the point Q is variable, a mixed-base representation (up to 32) may be applied in order to further speed up the calculations, given that all mixed-base formulas are implemented with a single inversion.

According to Gutub, there are various ways to apply elliptic curves in applications of cryptography. See A. Gutub, “Efficient utilization of scalable multipliers in parallel to compute gf (p) elliptic curve cryptographic operations,” Kuwait Journal of Science & Engineering (KJSE), December 2007, vol. 34, no. 2, pp. 165-182, 2007, incorporated herein by reference in its entirety. Getub noted how the algorithm utilized for calculating nP from P is based on the binary representation of n. This is because this is the efficient and practical way to implement in hardware systems. That is, the binary algorithm scans the bits of n and doubles the point Q k-times. Gutub further explained that an extra operation of point addition (Q+P) is essential and needed to perform in every case that a particular bit of n is found.

a) Weierstrass Elliptic Curve: This section represents the equations of the original work that the present algorithm is compared with. Elliptic curves over _p, where p>3 are considered. Such a curve, in the short Weierstrass form in the affine plan, is the set of all pairs (x,y)∈_p which fulfill:

y ² ≡x ³ +ax+b (mod p)   (1)

For P=(x_P, y_P) and Q=(x_Q, y_Q), one can compute P+Q by using the following equations, where the computation of λ differs based on two disjoint cases. See C. Paar and J. Pelzl, Understanding cryptography: a textbook for students and practitioners. Springer Science & Business Media, 2009, incorporated herein by reference in its entirety.

In case of addition where P≠Q:

$\begin{array}{cc}\lambda =\left(\frac{y_Q-y_P}{x_Q-x_P}\right)\mathrm{mod}p{x}_R={\lambda }^2-x_P-x_Q\mathrm{mod}p& \left(2\right)\end{array}$

In case of computing 2*P (doubling of order one) where P has coordinates (x₁, y₁):

$\begin{array}{cc}\lambda =\left(\frac{3x_1^2+a}{2y_1}\right)\mathrm{mod}p& \left(3\right)\end{array}$ $\begin{array}{cc}{x}_2={\lambda }^2-2x_1\mathrm{mod}p& \left(4\right)\end{array}$ $\begin{array}{cc}{y}_2=\lambda \left(x_1-x_2\right)-y_1\mathrm{mod}p& \left(5\right)\end{array}$

Where λ is the slope of the tangent through P, and x₂ and y₂, the affine coordinates after doubling P one time. While a two dimensional projective space can also be used for computations in Weierstrass form, here the focus is on computations in the affine plan.

b) Projective: Projective coordinates is another way of representing an elliptic curve. The elliptic curve F can be described by another equation, in the projective space P². That is, the polynomial defines a curve in the projective space P² which is also known as a Weierstrass equation:

Γ:Y²Z+a₁XYZ+a₃YZ²=X³+a₂X²Z+a₄XZ²+a₆Z³

See N. P. Smart, Cryptography made simple. Springer, 2016, vol. 481, incorporated herein by reference in its entirety.

According to Smart, a definition of a projective n-dimensional space P² over a field F is:

• • the set of (n+1)−tuples (x₀, . . . , x_n)∈Fⁿ⁺¹
  • where at least one x_i per tuple does not equal 0, and;
  • where an equivalence relation between two tuples (x_0,1, . . . , x_n,1) and (x_0,2, . . . , x_n,2) in Pⁿ=holds if ∃U∈F such that ∀i,x_i,1=Ux_i,2.

It is noted that a more general definition would replace the 3^rd condition with: ∀i,x_i,1=g_i(U)x_i,2 for some bijective function g_i.

The equivalence class of {U(x₀, . . . , x_n, U∈F} is denoted by [x₀, . . . , x_n], where these x₀, . . . , x_n are known as the homogeneous coordinates of that point (see Smart). Projective coordinates are useful in cases where there is a need to eradicate the performance of costly inversion operations (see Gutub).

Higuchi and Takagi and Okeya et al. noted how randomized projective coordinates on a Montgomery-form elliptic curve are effective in securing systems against side channel attacks. See A. Higuchi and N. Takagi, “A fast addition algorithm for elliptic curve arithmetic in gf (2n) using projective coordinates,” Information processing letters, vol. 76, no. 3, pp. 101-103,2000; and K. Okeya, K. Miyazaki, and K. Sakurai, “A fast scalar multiplication method with randomized projective coordinates on a montgomery-form elliptic curve secure against side channel attacks,” in International Conference on Information Security and Cryptology. Springer, 2001, pp. 428-439, each incorporated herein by reference in their entirety. For example, Okeya et al. recommended a scalar multiplication method that does not incur a higher computational cost for randomized projective coordinates of the Montgomery form of elliptic curves.

Homogeneous projective coordinates correspond to the 2-dimensional space through the substitution x=X/Z and y=Y/Z, so that the general Weierstrass form equates to:

E:Y ² Z+a ₁ XYZ+a ₃ YZ ² =X ³ +a ₂ X ² Z+a ₄ XZ ² +a ₆ Z ³.

Jacobian projective coordinates are obtained by substituting x=X/Z² and y=Y/Z², so that the general Weierstrass form equates to:

E:Y ² +a ₁ XYZ+a ₃ YZ ³ =X ³ +a ₂ X ² Z ² +a ₄ XZ ⁴ +a ₆ Z ⁶.

See M. Eichler and D. Zagier, The theory of Jacobi forms. Springer, 1985, vol. 55; P.-Y. Liardet and N. P. Smart, “Preventing spa/dpa in ecc systems using the jacobi form,” in International Workshop on Cryptographic Hardware and Embedded Systems. Springer, 2001, pp. 391-401; and O. Billet and M. Joye, “The jacobi model of an elliptic curve and side-channel analysis,” in International Symposium on Applied Algebra, Algebraic Algorithms, and Error-Correcting Codes. Springer, 2003, pp. 34-42, each incorporated herein by reference in their entirety.

With the use of a projective coordinates approach, the attacker is unable to predict the appearance of a specific value when the projective coordinates are randomized (see Okeya et al.; Higuchi et al.).

Specifically, Higuchi and Takagi proposed a fast addition algorithm on an elliptic curve over GF(2ⁿ) using projective coordinates:

x=X/Z

y=Y/Z²

According to Higuchi and Takagi, the above projective coordinates have less multiplications than the previously known fastest algorithm. See J. López and R. Dahab, “Improved algorithms for elliptic curve arithmetic in gf (2 n),” in International Workshop on Selected Areas in Cryptography. Springer, 1998, pp. 201-212, incorporated herein by reference in its entirety.

An aspect is to apply a term re-grouping step (referred to as labeling) to minimize the number of multiplications. By using these methods, higher orders doubling can be computed in an efficient way with an algorithm involving a single inverse.

Affine Re-Computation of Multi-Stage Doubling

When computing scalar multiplication of elliptic curve points P using fast algorithms inspired from Horner's rule, it is common to include operations of the type kP, that are referred here as the k^th order double of P. For large numbers, it may be necessary to compute doubles, and in particular high order doubles. Certain high order doubles may be directly computed without performing all steps.

FIG. 6 is a flowchart of a further method of computing the kernel of FIG. 4 (S404) using direct doubling in accordance with an exemplary aspect of the disclosure. In S602, higher order doubling is performed independently by the first computer system 104 and the second computer system 108 without going throw all steps that are required for the original affine coordinates algorithm.

2^kP is referred to as the k^th double of P, and to kP as the k^th order double of P. N_xk and N_yk denote the numerators of the x and y coordinates of the k^th order double (kP), which are denoted x_k and y_k respectively. Namely, rewrite

$x_k=\frac{N_{xk}}{U_k^2}\mathrm{and}{y}_k=\frac{N_{\mathrm{yk}}}{U_k^3},$

where k is the order of the desired double, and U_k denotes the corresponding added projective parameter.

In order to compute 4P, first find N_x2 and N_y2 where they are the numerators of x₂ and y₂, the x and y coordinates of the first double (2P) respectively.

For this, substitute the value of λ in Equation (3) in both equations of x and y coordinates, then multiply the transformed Equations (4) and (5) with (2_y1)², for denominators with U₂=2_y1. The obtained N_x2 and N_y2 expressions are,

N _{x 2} =(3x₁²+a)²−2x₁(2y₁)² mod p   (6)

N _y2=(3x₁²+a) (x₁(2y₁)²−N_x2)−2y₁²(2y₁)² mod p   (7)

Replace the variables x₂ and y₂ in the second double slope, getting:

${\lambda }_4=\frac{3x_2^2+a}{2y_2}\mathrm{mod}p$ ${\lambda }_4=\frac{3{\left(\frac{N_{x_2}}{U_2^2}\right)}^2+a}{2\frac{N_{y_2}}{U_2^3}}\mathrm{mod}p$

Note that U₂ is the denominator of the (2P) slope λ₂=λ. Now, eliminate the inverses by amplifying the fraction of λ₄ with

$\frac{U_2^4}{U_2^4},$

$\begin{array}{cc}{\lambda }_4=\frac{3N_{x_2}^2+{\mathrm{aU}}_2^4}{2N_{y_2}{U}_2}\mathrm{mod}p& \left(8\right)\end{array}$

For simplicity, consider,

W ₄=3N_x2²+aU₂⁴ mod p   (9)

q₄=2N_y2 mod p

The new denominator of the obtained slope λ₄ is:

U₄=q₄U₂ mod p   (10)

Then substitute the new slope equation in the x₄ and y₄ equations,

$x_4={\lambda }_4^2-2_{x_2}\mathrm{mod}p$ $x_4=\left(\frac{W_4}{U_4}\right)-2\frac{N_{x_2}}{U_2^2}\mathrm{mod}p$

Eliminating the inverses in x₄ equation by bringing to common denominator and amplifying the obtained fraction with the value of U₄² where from Equation 10,

$U_4=\left(2y_1\right)q_4\mathrm{get},$ $\begin{array}{cc}{U}_4^2{x}_4=W_4^2-2N_{x_2}{q}_4^2\mathrm{mod}p& \left(11\right)\end{array}$ $x_4=\frac{W_4^2-2N_{x_2}{q}_4^2}{U_4^2}\mathrm{mod}p$

Where to match,

$\begin{array}{cc}{x}_4=\frac{N_{x_4}}{U_4^2}\mathrm{mod}p& \left(12\right)\end{array}$

obtain:

N _{x 4} =W ₄ ²−2N_x2q₄²

Same steps will be applied in order to find and simplify y₄

$y_4={\lambda }_4\left(x_2-x_4\right)-y_2\mathrm{mod}p$ $y_4=\frac{W_4}{U_4}\left(\frac{N_{x_2}}{U_2^2}-\frac{N_{x_4}}{U_4^2}\right)-\frac{N_{y_2}}{U_2^3}\mathrm{mod}p$

Then amplify y₄ by U₄³

$\begin{array}{cc}{U}_4^3{y}_4=W_4\left(N_{x_2}{q}_4^2-N_{x_4}\right)-N_{y_2}{q}_4^2\mathrm{mod}p& \left(13\right)\end{array}$ $y_4=\frac{W_4\left(N_{x_2}{q}_4^2-N_{x_4}\right)-N_{y_2}{q}_4^3}{U_4^3}\mathrm{mod}p$

where to match

$\begin{array}{cc}{y}_4=\frac{N_{y_4}}{U_4^3}\mathrm{mod}p\mathrm{obtain}& \left(14\right)\end{array}$ $N_{y_4}=W_4\left(N_{x_2}{q}_4^2-N_{x_4}\right)-N_{y_2}{q}_4^3$

Furthermore, the equations for higher order double can be generalized for any doubling order. By using this form, in S602, compute N_xn and N_yn and then replace all the variables in the equation that are related to the order of the desired double in order to perform any advance double directly (Direct Doubling). Computing the previous W_n's, U_n's, N_xn's and N_yn's is required but with having N_xn and N_yn formulas, the computations can be done smoothly.

Here is the general form that performs any double:

$\begin{array}{cc}{W}_n=3N_{x_{n/2}}^2+{\mathrm{aU}}_{n/2}^4\mathrm{mod}p& \left(15\right)\end{array}$ $q_n=2N_{y_{n/2}}\mathrm{mod}p$ $\begin{array}{cc}{U}_n=q_n{U}_{n/2}\mathrm{mod}p& \left(16\right)\end{array}$ $\begin{array}{cc}{x}_n=\frac{W_n^2-2N_{x_{n/2}}{q}_n^2}{U_n^2}\mathrm{mod}p& \left(17\right)\end{array}$ $\begin{array}{cc}{x}_n=\frac{N_{x_n}}{U_n^2}\mathrm{mod}p& \left(18\right)\end{array}$ $\begin{array}{cc}{N}_{x_n}=W_n^2-2N_{x_{n/2}}{q}_n^2& \left(19\right)\end{array}$ $\begin{array}{cc}{y}_n=\frac{W_n\left(N_{x_{n/2}}{q}_n^2-N_{x_n}\right)-N_{y_{n/2}}{q}_n^2}{U_n^3}\mathrm{mod}p& \left(20\right)\end{array}$ $\begin{array}{cc}{y}_n=\frac{N_{y_n}}{U_n^3}\mathrm{mod}p& \left(21\right)\end{array}$ $\begin{array}{cc}{N}_{y_n}=W_n\left(N_{x_{n/2}}{q}_n^2-N_{x_n}\right)-N_{y_{n/2}}{q}_n^3& \left(22\right)\end{array}$

Where n is the order of the double and n/2 assigned to the previous power of 2 double.

a) Numerical Examples: In this section, the cyclic group of points on the elliptic curve E in FIG. 7 is used, where the order of E is 19 (see Paar et al.):

E:y ² ≡x ³+2·x+2 mod 17   (23)

As seen in FIG. 7, it starts from the primitive element P=(5,1) to 19P that represents the identity element, then flips to P again as it is the characteristic of any cyclic group. This curve is used in all numerical example sections at the end of each algorithm.

Let P=(5,1) to exemplify Direct Doubling algorithm of S602. First, compute N_x1 and N_y1 that are related to the point 2P=(6,3), then apply another four iterations in order to compute the point 32P mod 17 that is equivalent to the point 13P=(16,4).

N _{x 2} =(3x₁²+a)²−2x₁(2y₁)² mod p

N _{x 2} =(3(5)²+2)²−2(5)(2(1))² mod 17

N _{x 2} =13−6 mod 17

N_x2=7

N _{y 2} =(3x₁²+a)(x₁(2y₁)²−N_x2)−2y₁²(2y₁)² mod p

N _{y 2} =(3(5)²+2)((5)(2(1))²−7)−2(1)²(2(1))² mod 17

N _{y 2} =9(3−7)−8 mod 17

N_y2=7

where U₂=2_y1=2.

Now start the first iteration to find the variables N_x4, N_y4, W₄, q₄ and U₄ that are related to the point 4P.

W ₄=3N_x2²+aU₂⁴ mod P

W ₄=3(7)²+2(2)⁴ mod 17

W₄=9

q₄=2N_y2 mod p

q₄=2(7) mod 17

q₄=14

U₄=q₄U₂ mod p

U₄=14(2) mod 17

U₄=11

Then substitute these values in x₄ and y₄ equations, and get,

$x_4=\frac{W_4^2-2N_{x_2}{q}_4^2}{U_4^2}\mathrm{mod}p$ $x_4=\frac{9^2-2\left(7\right){\left(14\right)}^2}{2}\mathrm{mod}17$ $x_4=\frac{6}{2}\mathrm{mod}17$ $x_4=3$

where the inverse of 2 is 9 and N_x4=6.

$y_4=\frac{W_4\left(N_{x_2}{q}_4^2-N_{x_4}\right)-N_{y_2}{q}_4^3}{U_4^3}\mathrm{mod}p$ $y_4=\frac{9\left(7{\left(14\right)}^2-6\right)-7{\left(14\right)}^3}{{11}^3}\mathrm{mod}17$ $y_4=\frac{5}{5}\mathrm{mod}17$ $y_4=1$

wherein the inverse of 5 is 7 and N_y4=5.

Now the inputs for the next iteration are ready in order to compute the point 8P.

W ₈=3N_x4²+aU₄⁴ mod p

W ₈=3(6)²+2(11)⁴ mod 17

W₈=14

q₈2N_y4 mod p

q₈2(5) mod 17

q₈=10

U₈=q₈U₄ mod p

U ₈=10(11) mod 17

U₈=8

Then substitute these values in x₈ and y₈ equations, and get,

$x_8=\frac{W_8^2-2N_{x_4}{q}_8^2}{U_8^2}\mathrm{mod}p$ $x_8=\frac{{14}^2-2\left(6\right){\left(10\right)}^2}{8^2}\mathrm{mod}17$ $x_8=\frac{16}{13}\mathrm{mod}17$ $x_8=13$

where the inverse of 13 is 4 and N_x8=16.

$y_8=\frac{W_8\left(N_{x_4}{q}_8^2-N_{x_8}\right)-N_{y_4}{q}_8^3}{U_8^3}\mathrm{mod}p$ $y_8=\frac{14\left(6{\left(10\right)}^2-16\right)-5{\left(10\right)}^3}{8^3}\mathrm{mod}17$ $y_8=\frac{14}{2}\mathrm{mod}17$ $y_8=7$

where the inverse of 2 is 9 and N_y8=14.

Now substitute with the new values of N_x, N_y and U in the next iteration equations in order to compute the point 16P.

W ₁₆=3N_x8²+aU₈⁴ mod p

W ₁₆=3(16)²+2(8)⁴ mod 17

W₁₆=1

q₁₆=2N_y8 mod p

q ₁₆=2(14) mod 17

q₁₆=11

U₁₆=q₁₆U₈ mod p

U ₁₆=11(8) mod 17

U₁₆=3

Then substitute these values in x₁₆ and y₁₆ equations, and get,

$x_{16}=\frac{W_{16}^2-2N_{x_8}{q}_{16}^2}{U_{16}^2}\mathrm{mod}p$ $x_{16}=\frac{1^2-2\left(16\right){\left(11\right)}^2}{3^2}\mathrm{mod}17$ $x_{16}=\frac{5}{9}\mathrm{mod}17$ $x_{16}=10$

The inverse of 9 is 2 and N_x16=5.

$y_{16}=\frac{W_{16}\left(N_{x_8}{q}_{16}^2-N_{x_{16}}\right)-N_{y_8}{q}_{16}^3}{U_{16}^3}\mathrm{mod}p$ $y_{16}=\frac{1\left(16{\left(11\right)}^2-5\right)-14{\left(11\right)}^3}{3^3}\mathrm{mod}17$ $y_{16}=\frac{8}{10}\mathrm{mod}17$ $y_{16}=11$

The inverse of 10 is 12 and N_y16=8.

Now substitute with the new values of N_x, N_y and U in the last iteration equations in order to compute the desired point 32P.

W ₃₂=3N_x16²+aU₁₆⁴ mod p

W ₃₂=3(5)²+2(3)⁴ mod 17

W₃₂=16

q₃₂=2N_y16 mod p

q ₃₂=2(8) mod 17

q₃=16

U₃₂=q₃₂U₁₆ mod p

U ₃₂=16(3) mod 17

U₃₂=14

Then substitute these values in x₃₂ and y₃₂ equations, and get,

$x_{16}=\frac{W_{16}^2-2N_{x_8}{q}_{16}^2}{U_{16}^2}\mathrm{mod}p$ $x_{32}=\frac{{\left(16\right)}^2-2\left(5\right){\left(16\right)}^2}{{\left(14\right)}^2}\mathrm{mod}17$ $x_{32}=\frac{8}{9}\mathrm{mod}17$ $x_{32}=16$

where the inverse of 9 is 2 and N_x32=8. Further,

$y_{32}=\frac{W_{32}\left(N_{x_{16}}{q}_{32}^2-N_{x_{32}}\right)-N_{y_{16}}{q}_{32}^2}{U_{32}^3}\mathrm{mod}p$ $y_{32}=\frac{16\left(5{\left(16\right)}^2-8\right)-8{\left(16\right)}^3}{{\left(14\right)}^3}\mathrm{mod}17$ $y_{32}=\frac{11}{7}\mathrm{mod}17$ $y_{32}=4$

where the inverse of 7 is 5 and N_y32=11.

Intermediate Operations

As it is important to calculate the binary multiplicative 2ⁿ for points Q to compute a large degree isogeny, the algorithm may be enhanced by determining the intermediate steps like 3P, 5P, and 7P etc.

Subramanya Rao have worked on Montgomery curves and found an efficient technique to find point tripling. See S. R. S. Rao, “Three dimensional montgomery ladder, differential point tripling on montgomery curves and point quintupling on weierstrass' and edwards curves,” in International Conference on Cryptology in Africa. Springer, 2016, pp. 84-106, incorporated herein by reference in its entirety. Simply, the present algorithm optimizes an application of a single double to some point P then performs a point addition. This technique could be applied to all intermediate steps. A set of general forms are presented through which represent the interstitial points up to 31P.

A. Fast 2ⁿQ+P

As mentioned earlier in the background section, the complexity of the SIDH cryptosystem relies on the efficiency of computing isogenies between points on the elliptic curve. Thus, a further optimization may be accomplished by reducing the number of computations in performing the kernel Equation P+[k]Q. FIG. 8 is a flowchart of a further method of computing the kernel of FIG. 4 (S404) using direct doubling with a single inverse in accordance with an exemplary aspect of the disclosure. As an advanced exponent may be performed of a point on a curve with a single inverse, it would have been needed to compute an extra inverse for a differential point addition. Therefore, in this section a further optimization is introduced that includes mixing the advanced doubling equations with the addition and performs the addition with a single inverse.

The following equations have some variables like N_x, N_y, and U that are replaced with the variables related to each double.

The optimization includes substituting the value of x and y coordinates of the point 2ⁿP in Equations 18 and 21 respectively in the addition slope equation in 2.

${\lambda }_n=\frac{\frac{N_{y_n}}{U_n^3}-y_1}{\frac{N_{x_n}}{U_n^2}-x_1}\mathrm{mod}p$

Multiplying with U_n³ to eliminate the inverses,

$\begin{array}{cc}{\lambda }_{n+m}=\frac{N_{y_n}-y_1{U}_n^3}{N_{x_n}{U}_n-x_1{U}_n^3}\mathrm{mod}p& \left(24\right)\end{array}$ $\begin{array}{cc}{W}_{n+m}=N_{y_n}-y_1{U}_n^3\mathrm{mod}p& \left(25\right)\end{array}$ $\begin{array}{cc}{q}_{n+m}=N_{x_n}-x_1{U}_n^2\mathrm{mod}p& \left(26\right)\end{array}$ $\begin{array}{cc}{U}_{n+m}=U_n{q}_{n+m}\mathrm{mod}p& \left(27\right)\end{array}$

Substitute λ_n+m in the equations for x_n+m and y_n+m,

$x_{n+m}={\left(\frac{W_{n+m}}{U_{n+m}}\right)}^2-x_1-\frac{N_{x_n}}{U_n^2}\mathrm{mod}p$

Multiplying with U_n+m²,

$\begin{array}{cc}{x}_{n+m}=\frac{W_{n+m}^2-x_1{U}_{n+m}^2-N_{x_n}{q}_{n+m}^2}{U_{n+m}^2}\mathrm{mod}p& \left(28\right)\end{array}$ $x_{n+m}=\frac{N_{x_{n+m}}}{U_{n+m}^2}\mathrm{mod}p$

Now find y_n+m,

$y_{n+m}=\frac{W_{n+m}}{U_{n+m}}\left(x_1-\frac{N_{x_{n+m}}}{U_{n+m}^2}\right)-y_1\mathrm{mod}p$

Multiplying with U_n+m³,

$\begin{array}{cc}{y}_{n+m}=\frac{W_{n+m}\left(x_1{U}_{n+m}^2-N_{x_{n+m}}\right)-y_1{U}_{n+m}^3}{U_{n+m}^3}\mathrm{mod}p& \left(29\right)\end{array}$ $y_{n+m}=\frac{N_{y_{n+m}}}{U_{n+m}^3}\mathrm{mod}p$

a) Numerical Examples: Let P=(5,1), then apply the 2ⁿP+P algorithm to compute the new x and y coordinates. In this example 2²P+P may be applied in order to find the point 5P. Consider the values previously computed in the numerical example of affine re-computation of multi-stage doubling for the point 4P where,

N_x4=6

N_y4=5

U₄=11

Substitute these values in the Equations 28 and 29 then get,

$W_5=N_{y_4}-y_1{U}_4^3\mathrm{mod}p$ $W_5=5-1{\left(11\right)}^3\mathrm{mod}17$ $W_5=0$ $q_5=N_{x_4}-x_1{U}_4^2\mathrm{mod}p$ $q_5=6-5{\left(11\right)}^2\mathrm{mod}17$ $q_5=13$ $U_5=U_4{q}_5\mathrm{mod}p$ $U_5=11\left(13\right)\mathrm{mod}17$ $U_5=7$ $x_5=\frac{W_5^2-x_1{U}_5^2-N_{x_4}{q}_5^2}{U_5^2}\mathrm{mod}p$ $x_5=\frac{{\left(0\right)}^2-5{\left(7\right)}^2-6{\left(13\right)}^2}{{\left(7\right)}^2}\mathrm{mod}17$ $x_5=\frac{16}{15}\mathrm{mod}17$ $x_5=9$

Where the inverse of 15 is 8 and N_x5=16.

$y_5=\frac{W_5\left(x_1{U}_5^2-N_{x_4}\right)}{U_5^3}\mathrm{mod}p$ $y_5=\frac{0(5{\left(7\right)}^2-\left(16\right)-1{\left(7\right)}^3}{{\left(7\right)}^3}\mathrm{mod}p$ $y_5=\frac{14}{3}\mathrm{mod}17$ $y_5=16$

Where the inverse of 3 is 6 and N_y5=14.

B. Other General Forms

FIG. 9 is a flowchart of a further method of computing the kernel of FIG. 4 (S404) using intermediate scalar multiplication algorithms in accordance with an exemplary aspect of the disclosure. The same steps as in S602 may be followed in S902 to find intermediate scalar multiplication algorithms, multiplying an EC point P with scalars up to 31P. The above description of intermediate operations illustrate the importance of computing the intermediate equations for the overall speed-up of the present algorithms. Table I lists the general forms illustrated to implement all the points up to 31P.

TABLE I
INTERMEDIATE ALGORITHMS GENERAL FORMS.
Forms     Algorithms
2n P + P  Wn+m = Nyn − y1Un3
          qn+m = Nxn − x1Un2
          Un+m = Unqn+m
          x                                              n                        +                        m                                                              =                                                                                            W                                                      n                            +                            m                                                    2                                                -                                                                              x                            1                                                    ⁢                                                      U                                                          n                              +                              m                                                        2                                                                          -                                                                              N                                                          x                              n                                                                                ⁢                                                      y                            1                                                    ⁢                                                      q                                                          n                              +                              m                                                        2                                                                                                                      U                                                  n                          +                          m                                                2
          y                                              n                        +                        m                                                              =                                                                                                                        W                                                          n                              +                              m                                                                                (                                                                                                                    x                                1                                                            ⁢                                                              U                                                                  n                                  +                                  m                                                                2                                                                                      -                                                          N                                                              x                                                                  n                                  +                                  m                                                                                                                                              )                                                -                                                                              y                            1                                                    ⁢                                                      U                                                          n                              +                              m                                                        3                                                                                                                      U                                                  n                          +                          m                                                3
2n P + 2P Wn+m = Nyn − 8Nym4
          qn+m = Nxn − 4NxmNym2
          Un+m = Unqn+m
          x                                              n                        +                        m                                                              =                                                                                            W                                                      n                            +                            m                                                    2                                                -                                                  4                          ⁢                                                      N                                                          x                              m                                                                                ⁢                                                      N                                                          y                              m                                                        2                                                    ⁢                                                      q                                                          n                              +                              m                                                        2                                                                          -                                                                              N                                                          x                              n                                                                                ⁢                                                      y                            1                                                    ⁢                                                      q                                                          n                              +                              m                                                        2                                                                                                                      U                                                  n                          +                          m                                                2
          y                                              n                        +                        m                                                              =                                                                                                                        W                                                          n                              +                              m                                                                                (                                                                                    4                              ⁢                                                              N                                                                  x                                  m                                                                                            ⁢                                                              N                                                                  y                                  m                                                                2                                                            ⁢                                                              q                                                                  n                                  +                                  m                                                                2                                                                                      -                                                          N                                                              x                                                                  n                                  +                                  m                                                                                                                                              )                                                -                                                  8                          ⁢                          N                                                      ?                                                                                q                                                          n                              +                              m                                                        3                                                                                                                      U                                                  n                          +                          m                                                3
2(nP)     W2n = 3Nxn2 + aUn4
          q2n = 2Nyn
          U2n = Unq2n
          x                                              2                        ⁢                        n                                                              =                                                                                            W                                                      2                            ⁢                            n                                                    2                                                -                                                  2                          ⁢                                                      N                            xn                                                    ⁢                                                      q                                                          2                              ⁢                              n                                                        2                                                                                                                      U                                                  2                          ⁢                          n                                                2
          y                                              2                        ⁢                        n                                                              =                                                                                                                        W                                                          2                              ⁢                              n                                                                                (                                                                                                                    N                                                                  x                                  n                                                                                            ⁢                                                              q                                                                  2                                  ⁢                                  n                                                                2                                                                                      -                                                          N                                                              x                                                                  2                                  ⁢                                  n                                                                                                                                              )                                                -                                                                              N                                                          y                              n                                                                                ⁢                                                      q                                                          2                              ⁢                              n                                                        3                                                                                                                      U                                                  2                          ⁢                          n                                                3
2(nP) + P W2n+1 = Ny2n − y1U2n3
          q2n+1 = Nx2n − x1U2n2
          U2n+1 = U2nq2n+1
          x                                                                        2                          ⁢                          n                                                +                        1                                                              =                                                                                            W                                                                                    2                              ⁢                              n                                                        +                            1                                                    2                                                -                                                                              x                            1                                                    ⁢                                                      U                                                                                          2                                ⁢                                n                                                            +                              1                                                        2                                                                          -                                                                              N                                                          x                                                              2                                ⁢                                n                                                                                                              ⁢                                                      q                                                                                          2                                ⁢                                n                                                            +                              1                                                        2                                                                                                                      U                                                                              2                            ⁢                            n                                                    +                          1                                                2
          y                                                                        2                          ⁢                          n                                                +                        1                                                              =                                                                                                                        W                                                                                          2                                ⁢                                n                                                            +                              1                                                                                (                                                                                                                    x                                1                                                            ⁢                                                              U                                                                                                      2                                    ⁢                                    n                                                                    +                                  1                                                                2                                                                                      -                                                          N                                                              x                                                                                                      2                                    ⁢                                    n                                                                    +                                  1                                                                                                                                              )                                                -                                                                              y                            1                                                    ⁢                                                      U                                                                                          2                                ⁢                                n                                                            +                              1                                                        3                                                                                                                      U                                                                              2                            ⁢                            n                                                    +                          1                                                3
?                                        indicates text missing or illegible when filed

As it is known, the non-adjacent form (NAF) aims to reduce the number of one bit in the binary representation and thus reduce the number of operations, here in Table II the mathematical structure is relied on in representing all points up to 31P with the fastest and most efficient possible form.

TABLE II
STRUCTURE AND COST FOR ALL
INTERMEDIATE ALGORITHMS.
	Algorithms	Structure	Number of Multiplications
	3P	2P + P	19
	5P	4P + P	28
	6P	4P + 2P	26
	7P/9p	8P ± P	38
	10P	2(5P)	38
	11P	2(5P) + P	51
	12P	2(6P)	36
	13P	2(6P) + P	49
	14P/18P	2(7/9P)	48
	15P/17P	16P ± P	48
	19P	2(9P) + P	61
	20P	2(10P)	46
	21P	2(10P) + P	61
	22P	2(11P)	62
	23P/25P	2(12P) ± P	59
	24P	2(12P)	46
	26P	2(13P)	59
	27P	2(14P) − P	71
	28P	2(14P)	58
	29P	2(15P) − P	71
	30P	2(15P)	58
	31P	32P − P	58

Extraction Of Coordinates (EiSi Coordinates) Coordinate System

An aspect is a coordinate system that uses only x and y coordinates and performs a single inverse along the secret key size. FIG. 10 is a flowchart of a further method of computing the kernal of FIG. 4 (S404) using extraction of coordinates in accordance with an exemplary aspect of the disclosure. The EiSi coordinate system has competitive properties when compared with other coordinate systems. The EiSi coordinate system also provides a set of new and efficient operators in the Jacobian coordinate space. Different competing general scalar EC multiplication algorithms are introduced and used to compare between the scalar multiplication algorithms using multiple coordinate systems in order to illustrate the strengths of the EiSi coordinate system.

The EiSi coordinate system can be seen as a modified version of either the affine or Jacobian spaces with different operators. Also, the EiSi space operators offer faster arithmetic. Similarly to Projective techniques, this form of elliptic curve is represented with a single inversion at the last iteration. In S1002, each point P_A=(N_xA: N_yA: U_A) is represented in affine coordinates as(N_xA/U_A²,N_yA/U_A³).

Let P_A and P_B be points on an elliptic curve then, in affine space,

(X_A:Y_A)+(X_B:Y_B)=(X_C:Y_C)

At the first iteration consider U_A=U_B=1 to get,

$\left(N_{x_A}:N_{y_A}\right)+\left(N_{x_B}:N_{y_B}\right)=\left(\frac{N_{x_c}}{U_c^2}:\frac{N_{y_c}}{U_c^3}\right)$

Where, doubling S1004 is performed,

N _{x C} =(3x_A²+a)²−2x_A(2y_A)² mod p   (30)

N _{y C} =(3x_A²+a)(x_A(2y_A)²−N_xA)−2y_A²(2y_A)² mod p   (31)

U_c=2y_A mod p   (32)

Additionally, in S1006, two points are added after the first iteration, where the base point will be changed. A modified version of point addition algorithm is as follows. Let P_A and P_B be a point on the elliptic curve where (N_xA: N_yA: U_A) and (N_xB: N_yB: U_B) are the projective EiSi points representation respectively. Then,

(N_xA:N_yA:U_A)+(N_xB:N_yB:U_B)=(N_xC:N_yC:N_C)

In case P₁≠±P₂ (Addition),

W _C =N _yB U _A ³ −N _yA U _B ³ mod p   (33)

q _C =N _{x B} U _A ² −N _{x A} U _B ² mod p   (34)

U_C=U_AU_Bq_C mod p   (35)

N _{x C} =W _C ² −N _{x A} U _B ² q _C ² −N _{x B} U _A ² q _C ² mod p   (36)

N _{y C} =W _C(N_xAU_B²q_C²−N_xC)−N_yAU_N³q_C³ mod p   (37)

In case P_A=P_B (Higher Order Doubling), let P₁=P_A. Then, as proved in the description o affine re-computation of multi-stage doubling, recursively

W _n=3N_xn/2²+aU_n/2⁴ mod p   (38)

q_n=2N_yn/2 mod p   (39)

U_n=q_nU_n/2 mod p   (40)

N _{x n} =W _n ²−2N_xn/2q_n² mod p   (41)

N _{y n} =W _n(N_xn/2q_n²−N_xn)−N_yn/2q_n³ mod p   (42)

Using the extraction of coordinates, the algorithms may be rewritten to receive N_xn, N_yn and U_n instead of x_n and y_n values. By applying this method, computing the inverse at each iteration may be dispensed with.

Since all algorithms start with finding N_x2 and N_y2 values that are related to the point 2P, some adjustments are made to these algorithms in terms of the inputs, then:

N _{x 2} =(3N_xin+aU_in⁴)²−8N_xinN_yin² mod p   (43)

N _{y 2} =(3N_xin²+aU_in⁴)(4N_xinN_yin²−N_x2)−8N_yin⁴ mod p   (44)

U₂=2N_yinU_in mod p   (45)

Where N_xin and N_yin and U_in are the inputs that represent the point (X₁: Y₁) at the first iteration, where U_m equals one and N_x2, N_x2 and U₂ are the outputs that represent the point 2P.

a) Numerical Example: In this section the same cyclic group that was introduced in the description of affine re-computation of multi-stage doubling is used and some of the values that were previously computed in the previous numerical examples sections are considered in order to illustrate how the new coordinates system finds point doubling and addition correctly with a single inverse along the key size.

Assume that a key size of 4 bits represents the number 10₁₀=(1010)₂. Then the algorithm Left-to-Right is applied in order to compute the new x and y coordinates for the point 10P=(7,11).

First, the algorithm scans from left to right to process the second one-bit. Each 1-bit is represented as doubling and addition while each 0-bit is represented as only doubling. Thus, the algorithm performs doubling in order to find the point 2P.

As in the description of affine re-computation of multi-stage doubling, the flowing values are considered,

U₂=2

N_x2=7

N_y2=7

Then, the next bit is scanned from the left which is 1. Another double and add are applied to get 2(2P)+P=5P.

$\mathrm{KEY}:101\dots$ $\mathrm{Operations}:2P5P\dots$

As in the description of intermediate operations, the values for 5P are considered as well,

U₅=7

N_x5=16

N_y5=14

Note: N_x2, N_y2, N_x5 and N_y5 were computed with no inversion operation.

Now, the last bit is scanned which is 0. A doubling operation is applied to get 2(5P)=10P=(7,11).

$\mathrm{KEY}:1010$ $\mathrm{Operations}:2P5P10P$ $W_{10}=3N_{x_5}^2+{\mathrm{aU}}_5^4\mathrm{mod}p$ $W_{10}=3{\left(16\right)}^2+2{\left(7\right)}^4\mathrm{mod}17$ $W_{10}=11$ $q_{10}=2N_{y_5}\mathrm{mod}p$ $q_{10}=2\left(14\right)\mathrm{mod}17$ $q_{10}=11$ $U_{10}=q_{10}{U}_5\mathrm{mod}p$ $U_{10}=11\left(7\right)\mathrm{mod}17$ $U_{14}=9$ $N_{x_{10}}=W_{10}^2-2N_{x_5}{q}_{10}^2\mathrm{mod}p$ $N_{x_{10}}={\left(11\right)}^2-2\left(16\right){\left(11\right)}^2\mathrm{mod}17$ $N_{x_{10}}=6$ $N_{y_{10}}=W_{10}\left(N_{x_5}{q}_{10}^2-N_{x_{10}}\right)-N_{y_5}{q}_{10}^3\mathrm{mod}p$ $N_{y_{10}}=11\left(\left(16\right){\left(11\right)}^2-6\right)-\left(14\right){\left(11\right)}^3\mathrm{mod}17$ $N_{y_{10}}=12$

At the end of the last iteration, the inverse function is applied in order to find the affine coordinates for the point 10P.

$x_{10}=\frac{N_{x_{10}}}{U_{10}^2}=\frac{6}{13}=7$ $y_{10}=\frac{N_{y_{10}}}{U_{10}^3}=\frac{12}{15}=11$

Where 13⁻¹=4 and 15⁻¹=8.

Fast Multiplication With Mixed Base Multiplicands

The following is a description of a few algorithms that can integrate the fast repeated doubling techniques mentioned so far by applying mixed base multiplicands. With the algorithm mP+nQ one can compute multiplications with scalars up to 31. One can divide m's binary representation into blocks of five bits. In case an obtained block represents one of the unimplemented scalar multiplications, such blocks may be reduced in length.

A. Double and Add Extensions

In the above description of affine re-computation of multi-stage doubling and intermediate operations it is shown how to compute all intermediate exponent and mix doubling with a differential addition with a single inverse. The left-to-right algorithm starts scanning from left the next one-bit considering that the most significant bit is one. Then, it decides whether it applies doubling or doubling and addition depending on the data being read. For instance, if the first two one bits were representing the binary equivalent (101)₂ which is 5₁₀, it will multiply the base by 4 because it was shifted to the left by two bits. Since the last bit scanned is a 1, it also applies a differential addition to the point being doubled with the base point. Thus, the implementation will be 4Q+Q. FIG. 11 shows a practical example for calculating Q⁴⁷. This technique computes Q⁴⁷ with only four inverses, instead of the eight inverses when performing the original equations. However, by applying the EiSi coordinate system, one can compute the whole key with a single inverse similarly to Projective.

Algorithm 1: Double and Add Extensions
	procedure MultiplyL2RKnapsack(k, P) do
	|	if |k| ≤ 0 or k|k|−1 == 0 then return O;
	|	;
	|	D := P;
	|	for (int i := |k| − 2; i ≥ 0; ) do
	|	|	1 := 0;
	|	|	do
	|	|	|	i − −;
	|	|	|	l ++;
	|	|	while (ki == 0 and i ≥ 0);
	|	|	if ki == 1 then
	|	|	|	D := DoubleAndAddKnapsack(D, l, P);
	|	|	else
	|	|	|—	D := DoubleKnapsack(D, l);
	|	|—
	|—	return D;

In Line 1.1 of the pseudocode of the Double-and-Add extensions the DoubleAndAddKnapsack function is applied by taking as a parameter the counter 1 that specifies the current bit location, and the base point P to be added at the end. Otherwise, the DoubleKnapsack function is applied in which the shifting to the left is computed by multiplying the D value with 2^l.

B. Fast Multiplication with Base 32 Multiplicands

A special case algorithm based on base 32 representations of the multiplicands may be used in the scalar multiplication. FIG. 12 is a flowchart of a further method of computing the kernel of FIG. 4 (in S342, s352) using base 32 multiplicands in accordance with an exemplary aspect of the disclosure. In S1202, compute the kernel K by converting multiplicands to base 32. Then, in S1204, for a multiplier of the form qr₃₂ the computation is implemented by as,

32(qP)+rP mod p

For the scalar 10150=27A6₁₆, the obtained algorithm is equivalent to:

(32(9P)+29P)(32)+6P mod p

As noted in the above equation, and similar to Montgomery curve, the key is indistinguishable and can't be recognized by side channel attack since the algorithm applies point doubling and addition each iteration regardless of the key bit value. In addition, by applying direct doubling algorithms the scalar multiplication can be performed with a reduced number of point additions, which costs more than point doubling. Moreover, as noted in the description of EiSi, the EiSi coordinate system operates on two modes. The first is when it receives affine x and y coordinates, while the other deals with N_x, N_y and U as inputs. The base 32 Multiplicands algorithm increases the use of the first mode, which costs fewer multiplications. Thus, the Base 32 Multiplicands algorithm is one of the most efficient scalar multiplication algorithms.

RESULTS AND EXPERIMENTS

Simulation experiments are performed with a Java implementation of the disclosed algorithms. The algorithms have been applied on large parameters defined in the standard curves P-521, P-384, P-256 and P-224 from the National Institute of Standards and Technology (NIST). In addition, 10 different keys were picked that were randomly generated with an appropriate size for the x and y coordinates of each curve. Each algorithm has been executed multiple times and then the average time taken is computed to increase the accuracy of the calculations. Experimentally, the software implementation was tested on BeagleBone Black (BBB) System kit. See G. Coley, “Beaglebone black system reference manual,” Texas Instruments, Dallas, vol. 5, 2013, incorporated herein by reference in its entirety. The BBB has been equipped with a minimum set of features to allow the user to experience the power of the processor.

The system is equipped with one of the ARM Cortex-A8 family, AM3358/9 processor.

A. Functions Description and Properties

In this section, a description is provided for the important functions that are used in a software implementation and their properties. As the EiSi curve receives and returns two different forms of inputs and outputs, (x,y) or (N_x:N_y:U), details of the characteristics of the EiSi curve are clarified.

a) doubling2ⁿN: This special function was designed to receive and return an EiSi point. Basically, it receives the number of doubling of a point and then builds the equation for implementing this doubling. For example, if one wants to compute the point 6P, it requires finding the point 2P then 4P in order to fulfill the constructional equation for 6P, which is 2P+4P.

b) adv_addN2N_N and adv_subN2N_N: These functions receive and return EiSi points. Briefly, they perform point addition and subtraction between two EiSi points.

c) remi_point: This function receives an affine point and return a EiSi point. In addition, it receives the number of doubling of a point and then builds the equation for implementing this doubling. Mainly, it is used in the Base 32 Multiplicands algorithm specifically for computing the remainders, where all of them based on the same base point. It differs from doubling2ⁿN function, where all the doubling algorithms operators and labels are dependent. Basically, the point 4P can't be computed without finding the point 2P. As well as, the point 8P can't be computed without finding the point 2P then 4P. For example, if one wants to compute the third double for the base point 8P=(13,7), that is represented in EiSi coordinates as (16:14:8). The remipoint will compute the N_x, N_y and U values for the point 2(8P) then 4(8P) then return the EiSi point of 8(8P)=7P mod 19 that is represented as (0:14:2).

d) remi_func: This function works as a control for remi_point function. The remi_func has architectures of all the doubling algorithms and how they are implemented. Essentially, it has flags to be checked to avoid repeating any previously computed operations. FIG. 13 and FIG. 14, shows practical examples that were printed from an implementation debug page.

As shown in FIG. 13, the second block 2P parameters were not recomputed as well as at the third block for 8P. Also, in FIG. 14, at the last block it can be seen that 8P and 3P were previously computed and had only performed the point addition.

B. Comparison to Original

In this section, algorithms are compared in terms of number of multiplications, Base 32 Multiplicands and Double and Add, with the original affine algorithm. The original affine equations have been implemented with two different algorithms, Right-to-Left and Left-to-Right. Table III shows the substantial differences in the number of multiplications and inversions between these algorithms.

TABLE III
OUR WORK VS ORIGINAL ALGORITHMS
MEASUREMENTS.
	Number of Operations
	NIST Curve	Algorithm	Mults	Invs
	P-521	RL (Original)	658514	1059
		LR (Original)	478056	778
		DA	9162	1
		Base 32	7921	1
	P-384	RL (Original)	355788	775
		LR (Original)	259177	569
		DA	6714	1
		Base 32	5890	1
	P-256	RL (Original)	162131	519
		LR (Original)	115340	378
		DA	4466	1
		Base 32	4007	1
	P-224	RL (Original)	123461	450
		LR (Original)	88921	332
		DA	3921	1
		Base 32	3529	1

As in Table III, the great difference in number of multiplications and inversions between the present algorithm and the original can be seen as the present algorithm is faster by approximately 35 up to 83 times for the key size 224 bits and 521 bits respectively in case of comparing with RL and 25 up to 60 times in case of LR. This difference is due to the number of inverse operations that the original algorithm requires each point doubling or addition operation. FIG. 15 is a graph that illustrates this difference in chart where the present algorithm appears as a straight line along the x-axis.

C. EiSi Coordinates vs Others

Here, a comparison between the present algorithm and the other coordinates systems, Projective and Jacobian. Table IV shows a comparison of these algorithms in terms of the number of additions, subtractions, multiplications, divisions, modulos, maximum levels of parallelization and elapsed time for implementing them on the NIST standard curves P-521, P-384, P-256 and P-224.

TABLE IV
EISI VS OTHER COORDINATES MEASUREMENTS.
NIST	Number of Operations	Time
Curve	Algorithm	Mults	Divs	ALUs	Mods	MaxIs	ms
P-521	Projective	14900	315	4211	17899	9390	1035
	Jacobian	13312	301	4197	15821	8862	884
	Base 32	7921	296	6518	12290	7123	778
	DA	9162	301	7848	18528	9924	951
P-384	Projective	10901	226	3078	13107	6871	554
	Jacobian	9752	222	3075	11587	6491	480
	Base 32	5890	227	4833	9113	5224	421
	DA	6714	222	5746	13556	7267	508
P-256	Projective	7236	145	2043	8714	4564	286
	Jacobian	6488	147	2046	7708	4316	261
	Base 32	4007	150	3277	6210	3463	227
	DA	4466	147	3820	9028	4833	266
P-224	Projective	6356	127	1796	7656	4010	234
	Jacobian	5697	127	1796	6773	3790	215
	Base 32	3529	132	2886	5472	3041	193
	DA	3921	127	3354	7931	4243	218

As can be seen in Table IV, the present algorithm as represented in the last two algorithms are more efficient when it comes to number of multiplications. Clearly, Base 32 Multiplicands is the optimal algorithm in this case. Moreover, when comparing the maximum level of parallelization, the present algorithm outperforms the other coordinates algorithms as well through Base 32 Multiplicands which makes it the optimal algorithm in terms of both factors. Nevertheless, Double and Add algorithm (DA) which represents original EiSi coordinates appears to be the least efficient in terms of maximum levels. However together with the direct doubling technique of the present algorithm DA outperforms all other algorithms in all aspects. FIG. 16 and FIG. 17 contain graphs that show the comparisons between the present algorithm DA and the other coordinates algorithms in terms of number of multiplications and maximum level of parallelization with different key sizes.

As can be seen in FIG. 16 and FIG. 17, all algorithms are graphed in straight lines of varying slopes which allows the opportunity to apply a straight line equation to any algorithm to predict the expected values, whether it is the number of multiplications or the maximum levels of parallelization at the level of a larger key size. Table V lists all the equations related to each algorithm.

TABLE V
LIST OF ALGORITHMS LINEAR EQUATIONS.
	Algorithm	Number of Mults Eq.	Number of MaxLs Eq.
	Base 32	y = 14.777x + 220.37	y = 13.762x − 52.445
	DA	y = 17.658x − 48.287	y = 19.139x − 60.183
	Jacobian	y = 25.656x − 70.992	y = 17.089x − 52.485
	Projective	y = 28.795x − 121.85	y = 18.131x − 69.095

Predictably, the equations in Table Vare applied on two key sizes of the prime numbers of 751 and 1013. Table VI lists the expected number of multiplications and maximum levels.

TABLE VI
EXPECTED NUMBER OF MULTS AND MAXLS
WITH KEY OF SIZES 751 AND 1013.
	Expected Number of Mulls	Expected Number of MaxLs
Algorithm	751	1013	751	1013
Base 32	11317	15189	10282	13888
DA	13212	17839	14313	19327
Jacobian	19196	25918	12781	17258
Projective	21503	29047	13547	18297

As it can be seen in Table VI, the present algorithm which is represented in Base 32 Multiplicands maintains its place as the optimal algorithm in terms of number of multiplications and maximum levels of parallelization. Likewise, the Jacobian algorithm outperforms Projective in terms of the same factors. As such, another comparison can be made between Base 32 Multiplicands and the Jacobian coordinates algorithm to monitor if the difference in performance will shrink with the size of the key or continue to increase. Despite the slope values in the straight-line equations that show the differences, the delta value, A, is computed which is the difference between the y-axis values along the key size. Table VII shows the comparison between these two algorithms in terms of the same two factors, where,

Δ_i=y_n−y_m   (46)

Where i represents the key size and n and m represent the algorithms labels.

TABLE VII
Δ VALUES FOR BASS 32 VS JACOBIAN.
	Number of Mults	Number of MaxLs
Key Size	Base 32 vs Jacobian	Base 32 vs Jacobian
224	2168	749
256	2481	853
384	3862	1267
521	5391	1739
751	7879	2499
1013	10729	3370

It can be seen from the A values from Table VII, that the results in both cases show that the improvement scales with the size of the input.

D. Number of Multipliers Comparison

After tests and comparisons have proven the efficiency of the present algorithms and overcoming other coordinates systems algorithms, in this section the number of multiplications units each algorithm requires are specified to achieve the maximum levels of parallelism. Table VIII shows the number of multipliers per algorithm in the case of a key size 521.

TABLE VIII
THE NUMBER OF MULTIPLIERS APPROPRIATE TO
ACHIEVE THE HIGHEST LEVEL OF PARALLELISM.
	Algorithm	Number of MaxLs	Multipliers
	DA	9920	3
	Base 32	7133	3
	Projective	9370	6
	Jacobian	8862	4

As it can be seen in Table VIII, the appropriate number of multipliers to achieve the highest level of parallelism varies between algorithms. In addition, it is noted that if the number of multipliers is reduced a little, there may be very close result in terms of maximum levels of prallelization. Thus, another close comparison can be made in which the behavior of each algorithm is monitored in comparison with the others in multiple cases where the number of multipliers is uniform. Table IX shows another comparison between the present optimal algorithm that is specified in the previous sections compared to Jacobian, in terms of the MaxLs at specific number of multipliers.

TABLE IX
MAXL AT DIFFERENT NUMBER OF MULTIPLIERS.
	Multipliers	Algorithm	Number of MaxLs
	1	Bass 32	7982
		Jacobian	13011
	2	Base 32	7134
		Jacobian	8864
	3	Base 32	7133
		Jacobian	8863
	4	Base 32	7133
		Jacobian	8862

As it can be seen in Table IX, the present algorithms out perform the Jacobian algorithm in all levels, starting from a single multiplier, where the Base 32 algorithm appears to be 63% more efficient, up to 4 multipliers, where Jacobian reaches its peak with 24% slower than Base 32. At the 2 multipliers case, the performance of both algorithms improves significantly as the difference in efficiency becomes almost 24% with the preference remaining for the present algorithm. It can be noted that the Base 32 and Jacobian algorithms become highly ineffective as they continue to increase by one parallel level by increasing the number of multiplication units each time until they reach their peak. At the end, and in all cases, whether fewer or more multipliers are used, the efficiency of the present algorithm clearly outweighs the work of other coordinates systems algorithms. FIG. 18 shows the graphical representation of the relationship of the number of multipliers with the maximum levels of parallelism.

Claims

1. A system for secure communication, the system comprising: a first computer system;a communication network; anda second computer system,whereinthe first computer system and the second computer system are each provided with a public elliptic curve E,the first computer system and the second computer system each independently determine kernels KA and KB, respectively, and an isogeny mapping using the public elliptic curve E,the first computer system and the second computer system each independently determine mapped points based on the respective isogeny mappings,the first computer system and the second computer system each exchange the determined mapped points,the first computer system and the second computer system each independently compute the coefficients of two different secret elliptic curves under two different secret isogenies,the first computer system and the second computer system each independently compute a common secret key using the secret elliptic curves and the secret isogenies, andthe first computer system and the second computer system exchange messages using the common secret key, andwherein the first computer system and the second computer system each determine kernels KA and KB including computing mP+nQ by accessing a lookup table stored in a memory that contains a range of doubles of an end point of the respective kernels, where P and Q are points on the public elliptic curve and m and n are integers.

2. The system of claim 1, wherein in the computing the secret isogenies, the first computer system and the second computer system each computes a respective kernel KBA and KAB using mixed-base multiplicands with a single inversion.

3. The system of claim 1, wherein in the determining of the kernels KA and KB, the first computer system and the second computer system each computes the kernel by replacing all variables that are related to the order of the desired double such that an advanced double is performed directly, by determining:

4. The system of claim 3, wherein in the determining of the kernels KA and KB, the first computer system and the second computer system each computes the kernel by mixing the advanced doubling with the addition and the computation is performed with a single inversion, by performing:

5. The system of claim 3, wherein in the determining of the kernels KA and KB, the first computer system and the second computer system each computes intermediate scalar multiplications of the form 2nP+P, 2nP+2P, 2(nP), and 2(nP)+P.

6. The system of claim 1, wherein in the determining of the kernels KA and KB, the first computer system and the second computer system each computes the common secret key using an extraction of coordinates system, where each point PA=(NxA: NyA: UA) is represented in affine coordinates as (NxA/U2A, NyA/U3A).

7. The system of claim 6, further comprising: the first computer system and the second computer system eachperform doubling operations using operators for the affine coordinates, andperform adding using the operators for the affine coordinates.

8. The system of claim 2, wherein the first computer system and the second computer system each compute the respective kernel KBA and KAB by converting the multiplicands to base 32, and compute scalar multiplications using the base 32 multiplicands.

9. A supersingular isogeny-based cryptography method, the method comprising: providing a first computer system and a second computer system with a public elliptic curve E,independently determining, by the first computer system and the second computer system, kernels KA and KB, respectively, and an isogeny mapping using the public elliptic curve E,independently determining, by the first computer system and the second computer system, mapped points based on the respective isogeny mappings,exchanging, between the first computer system and the second computer system, the determined mapped points,independently computing, by the first computer system and the second computer system, coefficients of two different secret elliptic curves under two different secret isogenies,independently determining, by the first computer system and the second computer system, a common secret key using the secret elliptic curves and the secret isogenies, andexchanging, between the first computer system and the second computer system, messages using the common secret key, andwherein the first computer system and the second computer system each determine kernels KA and KB including computing mP+nQ by accessing a lookup table stored in a memory that contains a range of doubles of an end point of the respective kernels, where P and Q are points on the public elliptic curve and m and n are integers.

10. The method of claim 9, wherein the computing the secret isogenies, further includes computing, by the first computer system and the second computer system, a respective kernel KBA and KAB using mixed-base multiplicands with a single inversion.

11. The method of claim 9, wherein the determining of the kernels KA and KB, further includes computing, by the first computer system and the second computer system, the kernel by replacing all variables that are related to the order of the desired double such that an advanced double is performed directly, by determining:

12. The method of claim 11, wherein the determining of the kernels KA and KB, further includes, computing, by the first computer system and the second computer system, the kernel by mixing the advanced doubling with the addition and the computation is performed with a single inversion, by performing:

13. The method of claim 11, wherein the determining of the kernels KA and KB, further includes computing, by the first computer system and the second computer system, intermediate scalar multiplications of the form 2nP+P, 2nP+2P, 2(nP), and 2(nP)+P.

14. The method of claim 9, wherein the determining of the kernels KA and KB, further includes, computing, by the first computer system and the second computer system, the common secret key using an extraction of coordinates system, where each point PA=(NxA: NyA: UA) is represented in affine coordinates as (NxA/U2A, NyA/U3A).

15. The method of claim 14, further comprising: performing doubling operations using operators for the affine coordinates, andperforming adding using the operators for the affine coordinates.

16. The method of claim 10, further comprising, computing, by the first computer system and the second computer system, the respective kernel KBA and KAB by converting the multiplicands to base 32, and computing scalar multiplications using the base 32 multiplicands.

17. A non-transitory computer-readable medium storing instructions that are operable when executed by processing circuitry to perform operations comprising: receiving a public elliptic curve E,independently determining kernels KA and KB, respectively, and an isogeny mapping using the public elliptic curve E,independently determining mapped points based on the respective isogeny mappings,exchanging the determined mapped points,independently computing coefficients of two different secret elliptic curves under two different secret isogenies,independently determining a common secret key using the secret elliptic curves and the secret isogenies, andexchanging messages using the common secret key, anddetermining kernels KA and KB by computing mP+nQ by accessing a lookup table stored in a memory that contains a range of doubles of an end point of the respective kernels, where P and Q are points on the public elliptic curve and m and n are integers.

18. The non-transitory computer-readable medium of claim 17, wherein the computing the secret isogenies, further includes computing a respective kernel KBA and KAB using mixed-base multiplicands with a single inversion.

19. The non-transitory computer-readable medium of claim 18, wherein the determining of the kernels KA and KB, further includes, computing the common secret key using an extraction of coordinates system, where each point PA=(NxA: NyA: UA) is represented in affine coordinates as (NxA/U2A, NyA/U3A).

20. The non-transitory computer-readable medium of claim 18, further comprising, computing the respective kernel KBA and KAB by converting the multiplicands to base 32, and computing scalar multiplications using the base 32 multiplicands.