crypto/fipsmodule/ec/asm/p256_beeu-armv8-asm.pl - boringssl - Git at Google

 # Copyright Amazon.com Inc. or its affiliates. All Rights Reserved.
 #
 # Permission to use, copy, modify, and/or distribute this software for any
 # purpose with or without fee is hereby granted, provided that the above
 # copyright notice and this permission notice appear in all copies.
 #
 # THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 # WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 # MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
 # SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 # WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION
 # OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
 # CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */
 #
 #
 # This code is based on p256_beeu-x86_64-asm.pl (which is based on BN_mod_inverse_odd).
 #

 # The first two arguments should always be the flavour and output file path.
 if ($#ARGV < 1) { die "Not enough arguments provided.
   Two arguments are necessary: the flavour and the output file path."; }

 $flavour = shift;
 $output  = shift;

 $0 =~ m/(.*[\/\\])[^\/\\]+$/; $dir=$1;
 ( $xlate="${dir}arm-xlate.pl" and -f $xlate ) or
 ( $xlate="${dir}../../../perlasm/arm-xlate.pl" and -f $xlate) or
 die "can't locate arm-xlate.pl";

 open OUT,"| \"$^X\" $xlate $flavour $output";
 *STDOUT=*OUT;
 #############################################################################
 # extern int beeu_mod_inverse_vartime(BN_ULONG out[P256_LIMBS],
 #                                     BN_ULONG a[P256_LIMBS],
 #                                     BN_ULONG n[P256_LIMBS]);
 #
 # (Binary Extended GCD (Euclidean) Algorithm.
 #  See A. Menezes, P. vanOorschot, and S. Vanstone's Handbook of Applied Cryptography,
 #  Chapter 14, Algorithm 14.61 and Note 14.64
 #  http://cacr.uwaterloo.ca/hac/about/chap14.pdf)

 # Assumption 1: n is odd for the BEEU
 # Assumption 2: 1 < a < n < 2^256

 # Details
 # The inverse of x modulo y can be calculated using Alg. 14.61, where "a" would be that inverse.
 # In other words,
 # ax == 1 (mod y) (where the symbol “==“ denotes ”congruent“)
 #  a == x^{-1} (mod y)
 #
 # It can be shown that throughout all the iterations of the algorithm, the following holds:
 #    u = Ax + By
 #    v = Cx + Dy
 # The values B and D are not of interest in this case, so they need not be computed by the algorithm.
 # This means the following congruences hold through the iterations of the algorithm.
 #    Ax == u (mod y)
 #    Cx == v (mod y)

 # Now we will modify the notation to match that of BN_mod_inverse_odd()
 # on which beeu_mod_inverse_vartime() in `p256_beeu-x86_64-asm` is based.
 # In those functions:
 #    x, y -> a, n
 #    u, v -> B, A
 #    A, C -> X, Y’, where Y’ = -Y
 # Hence, the following holds throughout the algorithm iterations
 #    Xa == B (mod n)
 #   -Ya == A (mod n)
 #
 # Same algorithm in Python:
 # def beeu(a, n):
 #     X = 1
 #     Y = 0
 #     B = a
 #     A = n
 #     while (B != 0):
 #         while (B % 2) == 0:
 #             B >>= 1
 #             if (X % 2) == 1:
 #                 X = X + n
 #             X >>= 1
 #         while (A % 2) == 0:
 #             A >>= 1
 #             if (Y % 2) == 1:
 #                 Y = Y + n
 #             Y >>= 1
 #         if (B >= A):
 #             B = B - A
 #             X = X + Y
 #         else:
 #             A = A - B
 #             Y = Y + X
 #     if (A != 1):
 #         # error
 #         return 0
 #     else:
 #         while (Y > n):
 #             Y = Y - n
 #         Y = n - Y
 #         return Y


 # For the internal variables,
 # x0-x2, x30 are used to hold the modulus n. The input parameters passed in
 # x1,x2 are copied first before corrupting them. x0 (out) is stored on the stack.
 # x3-x7 are used for parameters, which is not the case in this function, so they are corruptible
 # x8 is corruptible here
 # (the function doesn't return a struct, hence x8 doesn't contain a passed-in address
 #  for that struct).
 # x9-x15 are corruptible registers
 # x19-x28 are callee-saved registers

 # X/Y will hold the inverse parameter
 # Assumption: a,n,X,Y < 2^(256)
 # Initially, X := 1, Y := 0
 #            A := n, B := a

 # Function parameters (as per the Procedure Call Standard)
 my($out, $a_in, $n_in)=map("x$_",(0..2));
 # Internal variables
 my($n0, $n1, $n2, $n3)=map("x$_",(0..2,30));
 my($x0, $x1, $x2, $x3, $x4)=map("x$_",(3..7));
 my($y0, $y1, $y2, $y3, $y4)=map("x$_",(8..12));
 my($shift)=("x13");
 my($t0, $t1, $t2, $t3)=map("x$_",(14,15,19,20));
 my($a0, $a1, $a2, $a3)=map("x$_",(21..24));
 my($b0, $b1, $b2, $b3)=map("x$_",(25..28));

 # if B == 0, jump to end of loop
 sub TEST_B_ZERO {
   return <<___;
     orr     $t0, $b0, $b1
     orr     $t0, $t0, $b2

     // reverse the bit order of $b0. This is needed for clz after this macro
     rbit     $t1, $b0

     orr     $t0, $t0, $b3
     cbz     $t0,.Lbeeu_loop_end
 ___
 }

 # Shift right by 1 bit, adding the modulus first if the variable is odd
 # if least_sig_bit(var0) == 0,
 #     goto shift1_<ctr>
 # else
 #     add n and goto shift1_<ctr>
 # Prerequisite: t0 = 0
 $g_next_label = 0;
 sub SHIFT1 {
   my ($var0, $var1, $var2, $var3, $var4) = @_;
   my $label = ".Lshift1_${g_next_label}";
   $g_next_label++;
   return <<___;
     tbz     $var0, #0, $label
     adds    $var0, $var0, $n0
     adcs    $var1, $var1, $n1
     adcs    $var2, $var2, $n2
     adcs    $var3, $var3, $n3
     adc     $var4, $var4, $t0
 $label:
     // var0 := [var1|var0]<64..1>;
     // i.e. concatenate var1 and var0,
     //      extract bits <64..1> from the resulting 128-bit value
     //      and put them in var0
     extr    $var0, $var1, $var0, #1
     extr    $var1, $var2, $var1, #1
     extr    $var2, $var3, $var2, #1
     extr    $var3, $var4, $var3, #1
     lsr     $var4, $var4, #1
 ___
 }

 # compilation by clang 10.0.0 with -O2/-O3 of
 #      a[0] = (a[0] >> count) | (a[1] << (64-count));
 #      a[1] = (a[1] >> count) | (a[2] << (64-count));
 #      a[2] = (a[2] >> count) | (a[3] << (64-count));
 #      a[3] >>= count;
 # Note: EXTR instruction used in SHIFT1 is similar to x86_64's SHRDQ
 # except that the second source operand of EXTR is only immediate;
 # that's why it cannot be used here where $shift is a variable
 #
 # In the following,
 # t0 := 0 - shift
 #
 # then var0, for example, will be shifted right as follows:
 # var0 := (var0 >> (uint(shift) mod 64)) | (var1 << (uint(t0) mod 64))
 # "uint() mod 64" is from the definition of LSL and LSR instructions.
 #
 # What matters here is the order of instructions relative to certain other
 # instructions, i.e.
 # - lsr and lsl must precede orr of the corresponding registers.
 # - lsl must preced the lsr of the same register afterwards.
 # The chosen order of the instructions overall is to try and maximize
 # the pipeline usage.
 sub SHIFT256 {
   my ($var0, $var1, $var2, $var3) = @_;
   return <<___;
     neg $t0, $shift
     lsr $var0, $var0, $shift
     lsl $t1, $var1, $t0

     lsr $var1, $var1, $shift
     lsl $t2, $var2, $t0

     orr $var0, $var0, $t1

     lsr $var2, $var2, $shift
     lsl $t3, $var3, $t0

     orr $var1, $var1, $t2

     lsr $var3, $var3, $shift

     orr $var2, $var2, $t3
 ___
 }

 $code.=<<___;
 #include "openssl/arm_arch.h"

 .text
 .globl  beeu_mod_inverse_vartime
 .type   beeu_mod_inverse_vartime, %function
 .align  4
 beeu_mod_inverse_vartime:
     // Reserve enough space for 14 8-byte registers on the stack
     // in the first stp call for x29, x30.
     // Then store the remaining callee-saved registers.
     //
     //    | x29 | x30 | x19 | x20 | ... | x27 | x28 |  x0 |  x2 |
     //    ^                                                     ^
     //    sp  <------------------- 112 bytes ----------------> old sp
     //   x29 (FP)
     //
     AARCH64_SIGN_LINK_REGISTER
     stp     x29,x30,[sp,#-112]!
     add     x29,sp,#0
     stp     x19,x20,[sp,#16]
     stp     x21,x22,[sp,#32]
     stp     x23,x24,[sp,#48]
     stp     x25,x26,[sp,#64]
     stp     x27,x28,[sp,#80]
     stp     x0,x2,[sp,#96]

     // B = b3..b0 := a
     ldp     $b0,$b1,[$a_in]
     ldp     $b2,$b3,[$a_in,#16]

     // n3..n0 := n
     // Note: the value of input params are changed in the following.
     ldp     $n0,$n1,[$n_in]
     ldp     $n2,$n3,[$n_in,#16]

     // A = a3..a0 := n
     mov     $a0, $n0
     mov     $a1, $n1
     mov     $a2, $n2
     mov     $a3, $n3

     // X = x4..x0 := 1
     mov     $x0, #1
     eor     $x1, $x1, $x1
     eor     $x2, $x2, $x2
     eor     $x3, $x3, $x3
     eor     $x4, $x4, $x4

     // Y = y4..y0 := 0
     eor     $y0, $y0, $y0
     eor     $y1, $y1, $y1
     eor     $y2, $y2, $y2
     eor     $y3, $y3, $y3
     eor     $y4, $y4, $y4

 .Lbeeu_loop:
     // if B == 0, jump to .Lbeeu_loop_end
     ${\TEST_B_ZERO}

     // 0 < B < |n|,
     // 0 < A <= |n|,
     // (1)      X*a  ==  B   (mod |n|),
     // (2) (-1)*Y*a  ==  A   (mod |n|)

     // Now divide B by the maximum possible power of two in the
     // integers, and divide X by the same value mod |n|.
     // When we're done, (1) still holds.

     // shift := number of trailing 0s in $b0
     // (      = number of leading 0s in $t1; see the "rbit" instruction in TEST_B_ZERO)
     clz     $shift, $t1

     // If there is no shift, goto shift_A_Y
     cbz     $shift, .Lbeeu_shift_A_Y

     // Shift B right by "$shift" bits
     ${\SHIFT256($b0, $b1, $b2, $b3)}

     // Shift X right by "$shift" bits, adding n whenever X becomes odd.
     // $shift--;
     // $t0 := 0; needed in the addition to the most significant word in SHIFT1
     eor     $t0, $t0, $t0
 .Lbeeu_shift_loop_X:
     ${\SHIFT1($x0, $x1, $x2, $x3, $x4)}
     subs    $shift, $shift, #1
     bne     .Lbeeu_shift_loop_X

     // Note: the steps above perform the same sequence as in p256_beeu-x86_64-asm.pl
     // with the following differences:
     // - "$shift" is set directly to the number of trailing 0s in B
     //   (using rbit and clz instructions)
     // - The loop is only used to call SHIFT1(X)
     //   and $shift is decreased while executing the X loop.
     // - SHIFT256(B, $shift) is performed before right-shifting X; they are independent

 .Lbeeu_shift_A_Y:
     // Same for A and Y.
     // Afterwards, (2) still holds.
     // Reverse the bit order of $a0
     // $shift := number of trailing 0s in $a0 (= number of leading 0s in $t1)
     rbit    $t1, $a0
     clz     $shift, $t1

     // If there is no shift, goto |B-A|, X+Y update
     cbz     $shift, .Lbeeu_update_B_X_or_A_Y

     // Shift A right by "$shift" bits
     ${\SHIFT256($a0, $a1, $a2, $a3)}

     // Shift Y right by "$shift" bits, adding n whenever Y becomes odd.
     // $shift--;
     // $t0 := 0; needed in the addition to the most significant word in SHIFT1
     eor     $t0, $t0, $t0
 .Lbeeu_shift_loop_Y:
     ${\SHIFT1($y0, $y1, $y2, $y3, $y4)}
     subs    $shift, $shift, #1
     bne     .Lbeeu_shift_loop_Y

 .Lbeeu_update_B_X_or_A_Y:
     // Try T := B - A; if cs, continue with B > A (cs: carry set = no borrow)
     // Note: this is a case of unsigned arithmetic, where T fits in 4 64-bit words
     //       without taking a sign bit if generated. The lack of a carry would
     //       indicate a negative result. See, for example,
     //       https://community.arm.com/developer/ip-products/processors/b/processors-ip-blog/posts/condition-codes-1-condition-flags-and-codes
     subs    $t0, $b0, $a0
     sbcs    $t1, $b1, $a1
     sbcs    $t2, $b2, $a2
     sbcs    $t3, $b3, $a3
     bcs     .Lbeeu_B_greater_than_A

     // Else A > B =>
     // A := A - B; Y := Y + X; goto beginning of the loop
     subs    $a0, $a0, $b0
     sbcs    $a1, $a1, $b1
     sbcs    $a2, $a2, $b2
     sbcs    $a3, $a3, $b3

     adds    $y0, $y0, $x0
     adcs    $y1, $y1, $x1
     adcs    $y2, $y2, $x2
     adcs    $y3, $y3, $x3
     adc     $y4, $y4, $x4
     b       .Lbeeu_loop

 .Lbeeu_B_greater_than_A:
     // Continue with B > A =>
     // B := B - A; X := X + Y; goto beginning of the loop
     mov     $b0, $t0
     mov     $b1, $t1
     mov     $b2, $t2
     mov     $b3, $t3

     adds    $x0, $x0, $y0
     adcs    $x1, $x1, $y1
     adcs    $x2, $x2, $y2
     adcs    $x3, $x3, $y3
     adc     $x4, $x4, $y4
     b       .Lbeeu_loop

 .Lbeeu_loop_end:
     // The Euclid's algorithm loop ends when A == gcd(a,n);
     // this would be 1, when a and n are co-prime (i.e. do not have a common factor).
     // Since (-1)*Y*a == A (mod |n|), Y>0
     // then out = -Y mod n

     // Verify that A = 1 ==> (-1)*Y*a = A = 1  (mod |n|)
     // Is A-1 == 0?
     // If not, fail.
     sub     $t0, $a0, #1
     orr     $t0, $t0, $a1
     orr     $t0, $t0, $a2
     orr     $t0, $t0, $a3
     cbnz    $t0, .Lbeeu_err

     // If Y>n ==> Y:=Y-n
 .Lbeeu_reduction_loop:
     // x_i := y_i - n_i (X is no longer needed, use it as temp)
     // ($t0 = 0 from above)
     subs    $x0, $y0, $n0
     sbcs    $x1, $y1, $n1
     sbcs    $x2, $y2, $n2
     sbcs    $x3, $y3, $n3
     sbcs    $x4, $y4, $t0

     // If result is non-negative (i.e., cs = carry set = no borrow),
     // y_i := x_i; goto reduce again
     // else
     // y_i := y_i; continue
     csel    $y0, $x0, $y0, cs
     csel    $y1, $x1, $y1, cs
     csel    $y2, $x2, $y2, cs
     csel    $y3, $x3, $y3, cs
     csel    $y4, $x4, $y4, cs
     bcs     .Lbeeu_reduction_loop

     // Now Y < n (Y cannot be equal to n, since the inverse cannot be 0)
     // out = -Y = n-Y
     subs    $y0, $n0, $y0
     sbcs    $y1, $n1, $y1
     sbcs    $y2, $n2, $y2
     sbcs    $y3, $n3, $y3

     // Save Y in output (out (x0) was saved on the stack)
     ldr     x3, [sp,#96]
     stp     $y0, $y1, [x3]
     stp     $y2, $y3, [x3,#16]
     // return 1 (success)
     mov     x0, #1
     b       .Lbeeu_finish

 .Lbeeu_err:
     // return 0 (error)
     eor     x0, x0, x0

 .Lbeeu_finish:
     // Restore callee-saved registers, except x0, x2
     add     sp,x29,#0
     ldp     x19,x20,[sp,#16]
     ldp     x21,x22,[sp,#32]
     ldp     x23,x24,[sp,#48]
     ldp     x25,x26,[sp,#64]
     ldp     x27,x28,[sp,#80]
     ldp     x29,x30,[sp],#112

     AARCH64_VALIDATE_LINK_REGISTER
     ret
 .size beeu_mod_inverse_vartime,.-beeu_mod_inverse_vartime
 ___


 foreach (split("\n",$code)) {
     s/\`([^\`]*)\`/eval $1/ge;

     print $_,"\n";
 }
 close STDOUT or die "error closing STDOUT: $!"; # enforce flush
	# Copyright Amazon.com Inc. or its affiliates. All Rights Reserved.
	#
	# Permission to use, copy, modify, and/or distribute this software for any
	# purpose with or without fee is hereby granted, provided that the above
	# copyright notice and this permission notice appear in all copies.
	#
	# THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
	# WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
	# MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
	# SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
	# WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION
	# OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
	# CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */
	#
	#
	# This code is based on p256_beeu-x86_64-asm.pl (which is based on BN_mod_inverse_odd).
	#

	# The first two arguments should always be the flavour and output file path.
	if ($#ARGV < 1) { die "Not enough arguments provided.
	Two arguments are necessary: the flavour and the output file path."; }

	$flavour = shift;
	$output = shift;

	$0 =~ m/(.*[\/\\])[^\/\\]+$/; $dir=$1;
	( $xlate="${dir}arm-xlate.pl" and -f $xlate ) or
	( $xlate="${dir}../../../perlasm/arm-xlate.pl" and -f $xlate) or
	die "can't locate arm-xlate.pl";

	open OUT,"\| \"$^X\" $xlate $flavour $output";
	STDOUT=OUT;
	#############################################################################
	# extern int beeu_mod_inverse_vartime(BN_ULONG out[P256_LIMBS],
	# BN_ULONG a[P256_LIMBS],
	# BN_ULONG n[P256_LIMBS]);
	#
	# (Binary Extended GCD (Euclidean) Algorithm.
	# See A. Menezes, P. vanOorschot, and S. Vanstone's Handbook of Applied Cryptography,
	# Chapter 14, Algorithm 14.61 and Note 14.64
	# http://cacr.uwaterloo.ca/hac/about/chap14.pdf)

	# Assumption 1: n is odd for the BEEU
	# Assumption 2: 1 < a < n < 2^256

	# Details
	# The inverse of x modulo y can be calculated using Alg. 14.61, where "a" would be that inverse.
	# In other words,
	# ax == 1 (mod y) (where the symbol “==“ denotes ”congruent“)
	# a == x^{-1} (mod y)
	#
	# It can be shown that throughout all the iterations of the algorithm, the following holds:
	# u = Ax + By
	# v = Cx + Dy
	# The values B and D are not of interest in this case, so they need not be computed by the algorithm.
	# This means the following congruences hold through the iterations of the algorithm.
	# Ax == u (mod y)
	# Cx == v (mod y)

	# Now we will modify the notation to match that of BN_mod_inverse_odd()
	# on which beeu_mod_inverse_vartime() in `p256_beeu-x86_64-asm` is based.
	# In those functions:
	# x, y -> a, n
	# u, v -> B, A
	# A, C -> X, Y’, where Y’ = -Y
	# Hence, the following holds throughout the algorithm iterations
	# Xa == B (mod n)
	# -Ya == A (mod n)
	#
	# Same algorithm in Python:
	# def beeu(a, n):
	# X = 1
	# Y = 0
	# B = a
	# A = n
	# while (B != 0):
	# while (B % 2) == 0:
	# B >>= 1
	# if (X % 2) == 1:
	# X = X + n
	# X >>= 1
	# while (A % 2) == 0:
	# A >>= 1
	# if (Y % 2) == 1:
	# Y = Y + n
	# Y >>= 1
	# if (B >= A):
	# B = B - A
	# X = X + Y
	# else:
	# A = A - B
	# Y = Y + X
	# if (A != 1):
	# # error
	# return 0
	# else:
	# while (Y > n):
	# Y = Y - n
	# Y = n - Y
	# return Y


	# For the internal variables,
	# x0-x2, x30 are used to hold the modulus n. The input parameters passed in
	# x1,x2 are copied first before corrupting them. x0 (out) is stored on the stack.
	# x3-x7 are used for parameters, which is not the case in this function, so they are corruptible
	# x8 is corruptible here
	# (the function doesn't return a struct, hence x8 doesn't contain a passed-in address
	# for that struct).
	# x9-x15 are corruptible registers
	# x19-x28 are callee-saved registers

	# X/Y will hold the inverse parameter
	# Assumption: a,n,X,Y < 2^(256)
	# Initially, X := 1, Y := 0
	# A := n, B := a

	# Function parameters (as per the Procedure Call Standard)
	my($out, $a_in, $n_in)=map("x$_",(0..2));
	# Internal variables
	my($n0, $n1, $n2, $n3)=map("x$_",(0..2,30));
	my($x0, $x1, $x2, $x3, $x4)=map("x$_",(3..7));
	my($y0, $y1, $y2, $y3, $y4)=map("x$_",(8..12));
	my($shift)=("x13");
	my($t0, $t1, $t2, $t3)=map("x$_",(14,15,19,20));
	my($a0, $a1, $a2, $a3)=map("x$_",(21..24));
	my($b0, $b1, $b2, $b3)=map("x$_",(25..28));

	# if B == 0, jump to end of loop
	sub TEST_B_ZERO {
	return <<___;
	orr $t0, $b0, $b1
	orr $t0, $t0, $b2

	// reverse the bit order of $b0. This is needed for clz after this macro
	rbit $t1, $b0

	orr $t0, $t0, $b3
	cbz $t0,.Lbeeu_loop_end
	___
	}

	# Shift right by 1 bit, adding the modulus first if the variable is odd
	# if least_sig_bit(var0) == 0,
	# goto shift1_<ctr>
	# else
	# add n and goto shift1_<ctr>
	# Prerequisite: t0 = 0
	$g_next_label = 0;
	sub SHIFT1 {
	my ($var0, $var1, $var2, $var3, $var4) = @_;
	my $label = ".Lshift1_${g_next_label}";
	$g_next_label++;
	return <<___;
	tbz $var0, #0, $label
	adds $var0, $var0, $n0
	adcs $var1, $var1, $n1
	adcs $var2, $var2, $n2
	adcs $var3, $var3, $n3
	adc $var4, $var4, $t0
	$label:
	// var0 := [var1\|var0]<64..1>;
	// i.e. concatenate var1 and var0,
	// extract bits <64..1> from the resulting 128-bit value
	// and put them in var0
	extr $var0, $var1, $var0, #1
	extr $var1, $var2, $var1, #1
	extr $var2, $var3, $var2, #1
	extr $var3, $var4, $var3, #1
	lsr $var4, $var4, #1
	___
	}

	# compilation by clang 10.0.0 with -O2/-O3 of
	# a[0] = (a[0] >> count) \| (a[1] << (64-count));
	# a[1] = (a[1] >> count) \| (a[2] << (64-count));
	# a[2] = (a[2] >> count) \| (a[3] << (64-count));
	# a[3] >>= count;
	# Note: EXTR instruction used in SHIFT1 is similar to x86_64's SHRDQ
	# except that the second source operand of EXTR is only immediate;
	# that's why it cannot be used here where $shift is a variable
	#
	# In the following,
	# t0 := 0 - shift
	#
	# then var0, for example, will be shifted right as follows:
	# var0 := (var0 >> (uint(shift) mod 64)) \| (var1 << (uint(t0) mod 64))
	# "uint() mod 64" is from the definition of LSL and LSR instructions.
	#
	# What matters here is the order of instructions relative to certain other
	# instructions, i.e.
	# - lsr and lsl must precede orr of the corresponding registers.
	# - lsl must preced the lsr of the same register afterwards.
	# The chosen order of the instructions overall is to try and maximize
	# the pipeline usage.
	sub SHIFT256 {
	my ($var0, $var1, $var2, $var3) = @_;
	return <<___;
	neg $t0, $shift
	lsr $var0, $var0, $shift
	lsl $t1, $var1, $t0

	lsr $var1, $var1, $shift
	lsl $t2, $var2, $t0

	orr $var0, $var0, $t1

	lsr $var2, $var2, $shift
	lsl $t3, $var3, $t0

	orr $var1, $var1, $t2

	lsr $var3, $var3, $shift

	orr $var2, $var2, $t3
	___
	}

	$code.=<<___;
	#include "openssl/arm_arch.h"

	.text
	.globl beeu_mod_inverse_vartime
	.type beeu_mod_inverse_vartime, %function
	.align 4
	beeu_mod_inverse_vartime:
	// Reserve enough space for 14 8-byte registers on the stack
	// in the first stp call for x29, x30.
	// Then store the remaining callee-saved registers.
	//
	// \| x29 \| x30 \| x19 \| x20 \| ... \| x27 \| x28 \| x0 \| x2 \|
	// ^ ^
	// sp <------------------- 112 bytes ----------------> old sp
	// x29 (FP)
	//
	AARCH64_SIGN_LINK_REGISTER
	stp x29,x30,[sp,#-112]!
	add x29,sp,#0
	stp x19,x20,[sp,#16]
	stp x21,x22,[sp,#32]
	stp x23,x24,[sp,#48]
	stp x25,x26,[sp,#64]
	stp x27,x28,[sp,#80]
	stp x0,x2,[sp,#96]

	// B = b3..b0 := a
	ldp $b0,$b1,[$a_in]
	ldp $b2,$b3,[$a_in,#16]

	// n3..n0 := n
	// Note: the value of input params are changed in the following.
	ldp $n0,$n1,[$n_in]
	ldp $n2,$n3,[$n_in,#16]

	// A = a3..a0 := n
	mov $a0, $n0
	mov $a1, $n1
	mov $a2, $n2
	mov $a3, $n3

	// X = x4..x0 := 1
	mov $x0, #1
	eor $x1, $x1, $x1
	eor $x2, $x2, $x2
	eor $x3, $x3, $x3
	eor $x4, $x4, $x4

	// Y = y4..y0 := 0
	eor $y0, $y0, $y0
	eor $y1, $y1, $y1
	eor $y2, $y2, $y2
	eor $y3, $y3, $y3
	eor $y4, $y4, $y4

	.Lbeeu_loop:
	// if B == 0, jump to .Lbeeu_loop_end
	${\TEST_B_ZERO}

	// 0 < B < \|n\|,
	// 0 < A <= \|n\|,
	// (1) X*a == B (mod \|n\|),
	// (2) (-1)Ya == A (mod \|n\|)

	// Now divide B by the maximum possible power of two in the
	// integers, and divide X by the same value mod \|n\|.
	// When we're done, (1) still holds.

	// shift := number of trailing 0s in $b0
	// ( = number of leading 0s in $t1; see the "rbit" instruction in TEST_B_ZERO)
	clz $shift, $t1

	// If there is no shift, goto shift_A_Y
	cbz $shift, .Lbeeu_shift_A_Y

	// Shift B right by "$shift" bits
	${\SHIFT256($b0, $b1, $b2, $b3)}

	// Shift X right by "$shift" bits, adding n whenever X becomes odd.
	// $shift--;
	// $t0 := 0; needed in the addition to the most significant word in SHIFT1
	eor $t0, $t0, $t0
	.Lbeeu_shift_loop_X:
	${\SHIFT1($x0, $x1, $x2, $x3, $x4)}
	subs $shift, $shift, #1
	bne .Lbeeu_shift_loop_X

	// Note: the steps above perform the same sequence as in p256_beeu-x86_64-asm.pl
	// with the following differences:
	// - "$shift" is set directly to the number of trailing 0s in B
	// (using rbit and clz instructions)
	// - The loop is only used to call SHIFT1(X)
	// and $shift is decreased while executing the X loop.
	// - SHIFT256(B, $shift) is performed before right-shifting X; they are independent

	.Lbeeu_shift_A_Y:
	// Same for A and Y.
	// Afterwards, (2) still holds.
	// Reverse the bit order of $a0
	// $shift := number of trailing 0s in $a0 (= number of leading 0s in $t1)
	rbit $t1, $a0
	clz $shift, $t1

	// If there is no shift, goto \|B-A\|, X+Y update
	cbz $shift, .Lbeeu_update_B_X_or_A_Y

	// Shift A right by "$shift" bits
	${\SHIFT256($a0, $a1, $a2, $a3)}

	// Shift Y right by "$shift" bits, adding n whenever Y becomes odd.
	// $shift--;
	// $t0 := 0; needed in the addition to the most significant word in SHIFT1
	eor $t0, $t0, $t0
	.Lbeeu_shift_loop_Y:
	${\SHIFT1($y0, $y1, $y2, $y3, $y4)}
	subs $shift, $shift, #1
	bne .Lbeeu_shift_loop_Y

	.Lbeeu_update_B_X_or_A_Y:
	// Try T := B - A; if cs, continue with B > A (cs: carry set = no borrow)
	// Note: this is a case of unsigned arithmetic, where T fits in 4 64-bit words
	// without taking a sign bit if generated. The lack of a carry would
	// indicate a negative result. See, for example,
	// https://community.arm.com/developer/ip-products/processors/b/processors-ip-blog/posts/condition-codes-1-condition-flags-and-codes
	subs $t0, $b0, $a0
	sbcs $t1, $b1, $a1
	sbcs $t2, $b2, $a2
	sbcs $t3, $b3, $a3
	bcs .Lbeeu_B_greater_than_A

	// Else A > B =>
	// A := A - B; Y := Y + X; goto beginning of the loop
	subs $a0, $a0, $b0
	sbcs $a1, $a1, $b1
	sbcs $a2, $a2, $b2
	sbcs $a3, $a3, $b3

	adds $y0, $y0, $x0
	adcs $y1, $y1, $x1
	adcs $y2, $y2, $x2
	adcs $y3, $y3, $x3
	adc $y4, $y4, $x4
	b .Lbeeu_loop

	.Lbeeu_B_greater_than_A:
	// Continue with B > A =>
	// B := B - A; X := X + Y; goto beginning of the loop
	mov $b0, $t0
	mov $b1, $t1
	mov $b2, $t2
	mov $b3, $t3

	adds $x0, $x0, $y0
	adcs $x1, $x1, $y1
	adcs $x2, $x2, $y2
	adcs $x3, $x3, $y3
	adc $x4, $x4, $y4
	b .Lbeeu_loop

	.Lbeeu_loop_end:
	// The Euclid's algorithm loop ends when A == gcd(a,n);
	// this would be 1, when a and n are co-prime (i.e. do not have a common factor).
	// Since (-1)Ya == A (mod \|n\|), Y>0
	// then out = -Y mod n

	// Verify that A = 1 ==> (-1)Ya = A = 1 (mod \|n\|)
	// Is A-1 == 0?
	// If not, fail.
	sub $t0, $a0, #1
	orr $t0, $t0, $a1
	orr $t0, $t0, $a2
	orr $t0, $t0, $a3
	cbnz $t0, .Lbeeu_err

	// If Y>n ==> Y:=Y-n
	.Lbeeu_reduction_loop:
	// x_i := y_i - n_i (X is no longer needed, use it as temp)
	// ($t0 = 0 from above)
	subs $x0, $y0, $n0
	sbcs $x1, $y1, $n1
	sbcs $x2, $y2, $n2
	sbcs $x3, $y3, $n3
	sbcs $x4, $y4, $t0

	// If result is non-negative (i.e., cs = carry set = no borrow),
	// y_i := x_i; goto reduce again
	// else
	// y_i := y_i; continue
	csel $y0, $x0, $y0, cs
	csel $y1, $x1, $y1, cs
	csel $y2, $x2, $y2, cs
	csel $y3, $x3, $y3, cs
	csel $y4, $x4, $y4, cs
	bcs .Lbeeu_reduction_loop

	// Now Y < n (Y cannot be equal to n, since the inverse cannot be 0)
	// out = -Y = n-Y
	subs $y0, $n0, $y0
	sbcs $y1, $n1, $y1
	sbcs $y2, $n2, $y2
	sbcs $y3, $n3, $y3

	// Save Y in output (out (x0) was saved on the stack)
	ldr x3, [sp,#96]
	stp $y0, $y1, [x3]
	stp $y2, $y3, [x3,#16]
	// return 1 (success)
	mov x0, #1
	b .Lbeeu_finish

	.Lbeeu_err:
	// return 0 (error)
	eor x0, x0, x0

	.Lbeeu_finish:
	// Restore callee-saved registers, except x0, x2
	add sp,x29,#0
	ldp x19,x20,[sp,#16]
	ldp x21,x22,[sp,#32]
	ldp x23,x24,[sp,#48]
	ldp x25,x26,[sp,#64]
	ldp x27,x28,[sp,#80]
	ldp x29,x30,[sp],#112

	AARCH64_VALIDATE_LINK_REGISTER
	ret
	.size beeu_mod_inverse_vartime,.-beeu_mod_inverse_vartime
	___


	foreach (split("\n",$code)) {
	s/\`([^\`]*)\`/eval $1/ge;

	print $_,"\n";
	}
	close STDOUT or die "error closing STDOUT: $!"; # enforce flush