bignum/i386/bignum_i386_asm.s

/*
 * CDDL HEADER START
 *
 * The contents of this file are subject to the terms of the
 * Common Development and Distribution License (the "License").
 * You may not use this file except in compliance with the License.
 *
 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
 * or http://www.opensolaris.org/os/licensing.
 * See the License for the specific language governing permissions
 * and limitations under the License.
 *
 * When distributing Covered Code, include this CDDL HEADER in each
 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
 * If applicable, add the following below this CDDL HEADER, with the
 * fields enclosed by brackets "[]" replaced with your own identifying
 * information: Portions Copyright [yyyy] [name of copyright owner]
 *
 * CDDL HEADER END
 */
/*
 * Copyright (c) 2004, 2010, Oracle and/or its affiliates. All rights reserved.
 */

#include <sys/asm_linkage.h>
#include <sys/x86_archext.h>
#include <sys/controlregs.h>

#if defined(__lint)

#include <sys/types.h>

uint32_t
bignum_use_sse2()
{ return (0); }

/* Not to be called by C code */
/* ARGSUSED */
uint32_t
big_mul_set_vec_sse2_r()
{ return (0); }

/* Not to be called by C code */
/* ARGSUSED */
uint32_t
big_mul_add_vec_sse2_r()
{ return (0); }

/* ARGSUSED */
uint32_t
big_mul_set_vec_sse2(uint32_t *r, uint32_t *a, int len, uint32_t digit)
{ return (0); }

/* ARGSUSED */
uint32_t
big_mul_add_vec_sse2(uint32_t *r, uint32_t *a, int len, uint32_t digit)
{ return (0); }

/* ARGSUSED */
void
big_mul_vec_sse2(uint32_t *r, uint32_t *a, int alen, uint32_t *b, int blen)
{}

/* ARGSUSED */
void
big_sqr_vec_sse2(uint32_t *r, uint32_t *a, int len)
{}

#if defined(MMX_MANAGE)

/* ARGSUSED */
uint32_t
big_mul_set_vec_sse2_nsv(uint32_t *r, uint32_t *a, int len, uint32_t digit)
{ return (0); }

/* ARGSUSED */
uint32_t
big_mul_add_vec_sse2_nsv(uint32_t *r, uint32_t *a, int len, uint32_t digit)
{ return (0); }

/* Not to be called by C code */
/* ARGSUSED */
void
big_sqr_vec_sse2_fc(uint32_t *r, uint32_t *a, int len)
{}

#endif  /* MMX_MANAGE */

/*
 * UMUL
 *
 */

/* ARGSUSED */
uint32_t
big_mul_set_vec_umul(uint32_t *r, uint32_t *a, int len, uint32_t digit)
{ return (0); }

/* ARGSUSED */
uint32_t
big_mul_add_vec_umul(uint32_t *r, uint32_t *a, int len, uint32_t digit)
{ return (0); }

#else   /* __lint */

#if defined(MMX_MANAGE)

#if defined(_KERNEL)

#define KPREEMPT_DISABLE call kpr_disable
#define KPREEMPT_ENABLE call kpr_enable
#define TEST_TS(reg)                    \
    movl    %cr0, reg;              \
    clts;                       \
    testl   $CR0_TS, reg

#else   /* _KERNEL */

#define KPREEMPT_DISABLE
#define KPREEMPT_ENABLE

#define TEST_TS(reg)                    \
    movl    $0, reg;                \
    testl   $CR0_TS, reg

#endif  /* _KERNEL */

#define MMX_SIZE 8
#define MMX_ALIGN 8

#define SAVE_MMX_PROLOG(sreg, nreg)         \
    subl    $_MUL(MMX_SIZE, nreg + MMX_ALIGN), %esp;    \
    movl    %esp, sreg;             \
    addl    $MMX_ALIGN, sreg;           \
    andl    $-1![MMX_ALIGN-1], sreg;

#define RSTOR_MMX_EPILOG(nreg)              \
    addl    $_MUL(MMX_SIZE, nreg + MMX_ALIGN), %esp;

#define SAVE_MMX_0TO4(sreg)         \
    SAVE_MMX_PROLOG(sreg, 5);       \
    movq    %mm0, 0(sreg);          \
    movq    %mm1, 8(sreg);          \
    movq    %mm2, 16(sreg);         \
    movq    %mm3, 24(sreg);         \
    movq    %mm4, 32(sreg)

#define RSTOR_MMX_0TO4(sreg)            \
    movq    0(sreg), %mm0;          \
    movq    8(sreg), %mm1;          \
    movq    16(sreg), %mm2;         \
    movq    24(sreg), %mm3;         \
    movq    32(sreg), %mm4;         \
    RSTOR_MMX_EPILOG(5)

#endif  /* MMX_MANAGE */

/ Note: this file contains implementations for
/   big_mul_set_vec()
/   big_mul_add_vec()
/   big_mul_vec()
/   big_sqr_vec()
/ One set of implementations is for SSE2-capable models.
/ The other uses no MMX, SSE, or SSE2 instructions, only
/ the x86 32 X 32 -> 64 unsigned multiply instruction, MUL.
/
/ The code for the implementations is grouped by SSE2 vs UMUL,
/ rather than grouping pairs of implementations for each function.
/ This is because the bignum implementation gets "imprinted"
/ on the correct implementation, at the time of first use,
/ so none of the code for the other implementations is ever
/ executed.  So, it is a no-brainer to layout the code to minimize
/ the "footprint" of executed code.

/ Can we use SSE2 instructions?  Return value is non-zero
/ if we can.
/
/ Note:
/   Using the cpuid instruction directly would work equally
/   well in userland and in the kernel, but we do not use the
/   cpuid instruction in the kernel, we use x86_featureset,
/   instead.  This means we honor any decisions the kernel
/   startup code may have made in setting this variable,
/   including disabling SSE2.  It might even be a good idea
/   to honor this kind of setting in userland, as well, but
/   the variable, x86_featureset is not readily available to
/   userland processes.
/
/ uint32_t
/ bignum_use_sse2()

    ENTRY(bignum_use_sse2)
#if defined(_KERNEL)
    xor %eax, %eax
    bt  $X86FSET_SSE2, x86_featureset
    adc     %eax, %eax
#else   /* _KERNEL */
    pushl   %ebx
    movl    $1, %eax        / Get feature information
    cpuid
    movl    %edx, %eax      / set return value
    popl    %ebx
    andl    $CPUID_INTC_EDX_SSE2, %eax
#endif  /* _KERNEL */
    ret
    SET_SIZE(bignum_use_sse2)


/ ------------------------------------------------------------------------
/       SSE2 Implementations
/ ------------------------------------------------------------------------

/ r = a * digit, r and a are vectors of length len
/ returns the carry digit
/ Suitable only for x86 models that support SSE2 instruction set extensions
/
/ uint32_t
/ big_mul_set_vec_sse2_r(uint32_t *r, uint32_t *a, int len, uint32_t digit)
/
/ r %edx
/ a %ebx
/ len   %ecx
/ digit %mm3
/
/ Does not touch the following registers: %esi, %edi, %mm4
/
/ N.B.:
/   This is strictly for internal use.
/   The interface is very light-weight.
/   All parameters are passed in registers.
/   It does not conform to the SYSV x86 ABI.
/   So, don't even think about calling this function directly from C code.
/
/ The basic multiply digit loop is unrolled 8 times.
/ Each comment is preceded by an instance number.
/ Instructions that have been moved retain their original, "natural"
/ instance number.  It should be easier this way to follow
/ the step-wise refinement process that went into constructing
/ the final code.

#define UNROLL      8
#define UNROLL32    32

    ENTRY(big_mul_set_vec_sse2_r)
    xorl    %eax, %eax  / if (len == 0) return (0);
    testl   %ecx, %ecx
    jz  .L17

    pxor    %mm0, %mm0  / cy = 0

.L15:
    cmpl    $UNROLL, %ecx
    jl  .L16
    movd    0(%ebx), %mm1   / 1: mm1 = a[i]
    pmuludq %mm3, %mm1  / 1: mm1 = digit * a[i]
    paddq   %mm1, %mm0  / 1: mm0 = digit * a[i] + cy;
    movd    4(%ebx), %mm1   / 2: mm1 = a[i]
    movd    %mm0, 0(%edx)   / 1: r[i] = product[31..0]
    psrlq   $32, %mm0   / 1: cy = product[63..32]

    pmuludq %mm3, %mm1  / 2: mm1 = digit * a[i]
    paddq   %mm1, %mm0  / 2: mm0 = digit * a[i] + cy;
    movd    8(%ebx), %mm1   / 3: mm1 = a[i]
    movd    %mm0, 4(%edx)   / 2: r[i] = product[31..0]
    psrlq   $32, %mm0   / 2: cy = product[63..32]

    pmuludq %mm3, %mm1  / 3: mm1 = digit * a[i]
    paddq   %mm1, %mm0  / 3: mm0 = digit * a[i] + cy;
    movd    12(%ebx), %mm1  / 4: mm1 = a[i]
    movd    %mm0, 8(%edx)   / 3: r[i] = product[31..0]
    psrlq   $32, %mm0   / 3: cy = product[63..32]

    pmuludq %mm3, %mm1  / 4: mm1 = digit * a[i]
    paddq   %mm1, %mm0  / 4: mm0 = digit * a[i] + cy;
    movd    16(%ebx), %mm1  / 5: mm1 = a[i]
    movd    %mm0, 12(%edx)  / 4: r[i] = product[31..0]
    psrlq   $32, %mm0   / 4: cy = product[63..32]

    pmuludq %mm3, %mm1  / 5: mm1 = digit * a[i]
    paddq   %mm1, %mm0  / 5: mm0 = digit * a[i] + cy;
    movd    20(%ebx), %mm1  / 6: mm1 = a[i]
    movd    %mm0, 16(%edx)  / 5: r[i] = product[31..0]
    psrlq   $32, %mm0   / 5: cy = product[63..32]

    pmuludq %mm3, %mm1  / 6: mm1 = digit * a[i]
    paddq   %mm1, %mm0  / 6: mm0 = digit * a[i] + cy;
    movd    24(%ebx), %mm1  / 7: mm1 = a[i]
    movd    %mm0, 20(%edx)  / 6: r[i] = product[31..0]
    psrlq   $32, %mm0   / 6: cy = product[63..32]

    pmuludq %mm3, %mm1  / 7: mm1 = digit * a[i]
    paddq   %mm1, %mm0  / 7: mm0 = digit * a[i] + cy;
    movd    28(%ebx), %mm1  / 8: mm1 = a[i]
    movd    %mm0, 24(%edx)  / 7: r[i] = product[31..0]
    psrlq   $32, %mm0   / 7: cy = product[63..32]

    pmuludq %mm3, %mm1  / 8: mm1 = digit * a[i]
    paddq   %mm1, %mm0  / 8: mm0 = digit * a[i] + cy;
    movd    %mm0, 28(%edx)  / 8: r[i] = product[31..0]
    psrlq   $32, %mm0   / 8: cy = product[63..32]

    leal    UNROLL32(%ebx), %ebx    / a += UNROLL
    leal    UNROLL32(%edx), %edx    / r += UNROLL
    subl    $UNROLL, %ecx       / len -= UNROLL
    jz  .L17
    jmp .L15

.L16:
    movd    0(%ebx), %mm1   / 1: mm1 = a[i]
    pmuludq %mm3, %mm1  / 1: mm1 = digit * a[i]
    paddq   %mm1, %mm0  / 1: mm0 = digit * a[i] + cy;
    movd    %mm0, 0(%edx)   / 1: r[i] = product[31..0]
    psrlq   $32, %mm0   / 1: cy = product[63..32]
    subl    $1, %ecx
    jz  .L17

    movd    4(%ebx), %mm1   / 2: mm1 = a[i]
    pmuludq %mm3, %mm1  / 2: mm1 = digit * a[i]
    paddq   %mm1, %mm0  / 2: mm0 = digit * a[i] + cy;
    movd    %mm0, 4(%edx)   / 2: r[i] = product[31..0]
    psrlq   $32, %mm0   / 2: cy = product[63..32]
    subl    $1, %ecx
    jz  .L17

    movd    8(%ebx), %mm1   / 3: mm1 = a[i]
    pmuludq %mm3, %mm1  / 3: mm1 = digit * a[i]
    paddq   %mm1, %mm0  / 3: mm0 = digit * a[i] + cy;
    movd    %mm0, 8(%edx)   / 3: r[i] = product[31..0]
    psrlq   $32, %mm0   / 3: cy = product[63..32]
    subl    $1, %ecx
    jz  .L17

    movd    12(%ebx), %mm1  / 4: mm1 = a[i]
    pmuludq %mm3, %mm1  / 4: mm1 = digit * a[i]
    paddq   %mm1, %mm0  / 4: mm0 = digit * a[i] + cy;
    movd    %mm0, 12(%edx)  / 4: r[i] = product[31..0]
    psrlq   $32, %mm0   / 4: cy = product[63..32]
    subl    $1, %ecx
    jz  .L17

    movd    16(%ebx), %mm1  / 5: mm1 = a[i]
    pmuludq %mm3, %mm1  / 5: mm1 = digit * a[i]
    paddq   %mm1, %mm0  / 5: mm0 = digit * a[i] + cy;
    movd    %mm0, 16(%edx)  / 5: r[i] = product[31..0]
    psrlq   $32, %mm0   / 5: cy = product[63..32]
    subl    $1, %ecx
    jz  .L17

    movd    20(%ebx), %mm1  / 6: mm1 = a[i]
    pmuludq %mm3, %mm1  / 6: mm1 = digit * a[i]
    paddq   %mm1, %mm0  / 6: mm0 = digit * a[i] + cy;
    movd    %mm0, 20(%edx)  / 6: r[i] = product[31..0]
    psrlq   $32, %mm0   / 6: cy = product[63..32]
    subl    $1, %ecx
    jz  .L17

    movd    24(%ebx), %mm1  / 7: mm1 = a[i]
    pmuludq %mm3, %mm1  / 7: mm1 = digit * a[i]
    paddq   %mm1, %mm0  / 7: mm0 = digit * a[i] + cy;
    movd    %mm0, 24(%edx)  / 7: r[i] = product[31..0]
    psrlq   $32, %mm0   / 7: cy = product[63..32]

.L17:
    movd    %mm0, %eax  / return (cy)
    / no emms.  caller is responsible for emms
    ret
    SET_SIZE(big_mul_set_vec_sse2_r)


/ r = a * digit, r and a are vectors of length len
/ returns the carry digit
/ Suitable only for x86 models that support SSE2 instruction set extensions
/
/ r      8(%ebp)    %edx
/ a     12(%ebp)    %ebx
/ len       16(%ebp)    %ecx
/ digit     20(%ebp)    %mm3
/
/ In userland, there is just the one function, big_mul_set_vec_sse2().
/ But in the kernel, there are two variations:
/    1. big_mul_set_vec_sse2() which does what is necessary to save and
/       restore state, if necessary, and to ensure that preemtion is
/       disabled.
/    2. big_mul_set_vec_sse2_nsv() which just does the work;
/       it is the caller's responsibility to ensure that MMX state
/       does not need to be saved and restored and that preemption
/       is already disabled.

#if defined(MMX_MANAGE)
    ENTRY(big_mul_set_vec_sse2)
    pushl   %ebp
    movl    %esp, %ebp
    pushl   %ebx
    pushl   %esi
    KPREEMPT_DISABLE
    TEST_TS(%ebx)
    pushl   %ebx
    jnz .setvec_no_save
    pushl   %edi
    SAVE_MMX_0TO4(%edi)
    movl    8(%ebp), %edx
    movl    12(%ebp), %ebx
    movl    16(%ebp), %ecx
    movd    20(%ebp), %mm3
    call    big_mul_set_vec_sse2_r
    movl    %eax, %esi
    RSTOR_MMX_0TO4(%edi)
    popl    %edi
    jmp .setvec_rtn

.setvec_no_save:
    movl    8(%ebp), %edx
    movl    12(%ebp), %ebx
    movl    16(%ebp), %ecx
    movd    20(%ebp), %mm3
    call    big_mul_set_vec_sse2_r
    movl    %eax, %esi

.setvec_rtn:
    emms
    popl    %ebx
    movl    %ebx, %cr0
    KPREEMPT_ENABLE
    movl    %esi, %eax
    popl    %esi
    popl    %ebx
    leave
    ret
    SET_SIZE(big_mul_set_vec_sse2)

    ENTRY(big_mul_set_vec_sse2_nsv)
    pushl   %ebp
    movl    %esp, %ebp
    pushl   %ebx
    movl    8(%ebp), %edx
    movl    12(%ebp), %ebx
    movl    16(%ebp), %ecx
    movd    20(%ebp), %mm3
    call    big_mul_set_vec_sse2_r
    popl    %ebx
    leave
    ret
    SET_SIZE(big_mul_set_vec_sse2_nsv)

#else   /* !defined(MMX_MANAGE) */

/ r = a * digit, r and a are vectors of length len
/ returns the carry digit
/ Suitable only for x86 models that support SSE2 instruction set extensions
/
/ r      8(%ebp)    %edx
/ a     12(%ebp)    %ebx
/ len       16(%ebp)    %ecx
/ digit     20(%ebp)    %mm3

    ENTRY(big_mul_set_vec_sse2)
    pushl   %ebp
    movl    %esp, %ebp
    pushl   %ebx
    movl    8(%ebp), %edx
    movl    12(%ebp), %ebx
    movl    16(%ebp), %ecx
    movd    20(%ebp), %mm3
    call    big_mul_set_vec_sse2_r
    popl    %ebx
    emms
    leave
    ret
    SET_SIZE(big_mul_set_vec_sse2)

#endif  /* MMX_MANAGE */


/ r = r + a * digit, r and a are vectors of length len
/ returns the carry digit
/ Suitable only for x86 models that support SSE2 instruction set extensions
/
/ uint32_t
/ big_mul_add_vec_sse2_r(uint32_t *r, uint32_t *a, int len, uint32_t digit)
/
/ r %edx
/ a %ebx
/ len   %ecx
/ digit %mm3
/
/ N.B.:
/   This is strictly for internal use.
/   The interface is very light-weight.
/   All parameters are passed in registers.
/   It does not conform to the SYSV x86 ABI.
/   So, don't even think about calling this function directly from C code.
/
/ The basic multiply digit loop is unrolled 8 times.
/ Each comment is preceded by an instance number.
/ Instructions that have been moved retain their original, "natural"
/ instance number.  It should be easier this way to follow
/ the step-wise refinement process that went into constructing
/ the final code.

    ENTRY(big_mul_add_vec_sse2_r)
    xorl    %eax, %eax
    testl   %ecx, %ecx
    jz  .L27

    pxor    %mm0, %mm0  / cy = 0

.L25:
    cmpl    $UNROLL, %ecx
    jl  .L26
    movd    0(%ebx), %mm1   / 1: mm1 = a[i]
    movd    0(%edx), %mm2   / 1: mm2 = r[i]
    pmuludq %mm3, %mm1  / 1: mm1 = digit * a[i]
    paddq   %mm1, %mm2  / 1: mm2 = digit * a[i] + r[i]
    movd    4(%ebx), %mm1   / 2: mm1 = a[i]
    paddq   %mm2, %mm0  / 1: mm0 = digit * a[i] + r[i] + cy;
    movd    %mm0, 0(%edx)   / 1: r[i] = product[31..0]
    movd    4(%edx), %mm2   / 2: mm2 = r[i]
    psrlq   $32, %mm0   / 1: cy = product[63..32]

    pmuludq %mm3, %mm1  / 2: mm1 = digit * a[i]
    paddq   %mm1, %mm2  / 2: mm2 = digit * a[i] + r[i]
    movd    8(%ebx), %mm1   / 3: mm1 = a[i]
    paddq   %mm2, %mm0  / 2: mm0 = digit * a[i] + r[i] + cy;
    movd    %mm0, 4(%edx)   / 2: r[i] = product[31..0]
    movd    8(%edx), %mm2   / 3: mm2 = r[i]
    psrlq   $32, %mm0   / 2: cy = product[63..32]

    pmuludq %mm3, %mm1  / 3: mm1 = digit * a[i]
    paddq   %mm1, %mm2  / 3: mm2 = digit * a[i] + r[i]
    movd    12(%ebx), %mm1  / 4: mm1 = a[i]
    paddq   %mm2, %mm0  / 3: mm0 = digit * a[i] + r[i] + cy;
    movd    %mm0, 8(%edx)   / 3: r[i] = product[31..0]
    movd    12(%edx), %mm2  / 4: mm2 = r[i]
    psrlq   $32, %mm0   / 3: cy = product[63..32]

    pmuludq %mm3, %mm1  / 4: mm1 = digit * a[i]
    paddq   %mm1, %mm2  / 4: mm2 = digit * a[i] + r[i]
    movd    16(%ebx), %mm1  / 5: mm1 = a[i]
    paddq   %mm2, %mm0  / 4: mm0 = digit * a[i] + r[i] + cy;
    movd    %mm0, 12(%edx)  / 4: r[i] = product[31..0]
    movd    16(%edx), %mm2  / 5: mm2 = r[i]
    psrlq   $32, %mm0   / 4: cy = product[63..32]

    pmuludq %mm3, %mm1  / 5: mm1 = digit * a[i]
    paddq   %mm1, %mm2  / 5: mm2 = digit * a[i] + r[i]
    movd    20(%ebx), %mm1  / 6: mm1 = a[i]
    paddq   %mm2, %mm0  / 5: mm0 = digit * a[i] + r[i] + cy;
    movd    %mm0, 16(%edx)  / 5: r[i] = product[31..0]
    movd    20(%edx), %mm2  / 6: mm2 = r[i]
    psrlq   $32, %mm0   / 5: cy = product[63..32]

    pmuludq %mm3, %mm1  / 6: mm1 = digit * a[i]
    paddq   %mm1, %mm2  / 6: mm2 = digit * a[i] + r[i]
    movd    24(%ebx), %mm1  / 7: mm1 = a[i]
    paddq   %mm2, %mm0  / 6: mm0 = digit * a[i] + r[i] + cy;
    movd    %mm0, 20(%edx)  / 6: r[i] = product[31..0]
    movd    24(%edx), %mm2  / 7: mm2 = r[i]
    psrlq   $32, %mm0   / 6: cy = product[63..32]

    pmuludq %mm3, %mm1  / 7: mm1 = digit * a[i]
    paddq   %mm1, %mm2  / 7: mm2 = digit * a[i] + r[i]
    movd    28(%ebx), %mm1  / 8: mm1 = a[i]
    paddq   %mm2, %mm0  / 7: mm0 = digit * a[i] + r[i] + cy;
    movd    %mm0, 24(%edx)  / 7: r[i] = product[31..0]
    movd    28(%edx), %mm2  / 8: mm2 = r[i]
    psrlq   $32, %mm0   / 7: cy = product[63..32]

    pmuludq %mm3, %mm1  / 8: mm1 = digit * a[i]
    paddq   %mm1, %mm2  / 8: mm2 = digit * a[i] + r[i]
    paddq   %mm2, %mm0  / 8: mm0 = digit * a[i] + r[i] + cy;
    movd    %mm0, 28(%edx)  / 8: r[i] = product[31..0]
    psrlq   $32, %mm0   / 8: cy = product[63..32]

    leal    UNROLL32(%ebx), %ebx    / a += UNROLL
    leal    UNROLL32(%edx), %edx    / r += UNROLL
    subl    $UNROLL, %ecx       / len -= UNROLL
    jz  .L27
    jmp .L25

.L26:
    movd    0(%ebx), %mm1   / 1: mm1 = a[i]
    movd    0(%edx), %mm2   / 1: mm2 = r[i]
    pmuludq %mm3, %mm1  / 1: mm1 = digit * a[i]
    paddq   %mm1, %mm2  / 1: mm2 = digit * a[i] + r[i]
    paddq   %mm2, %mm0  / 1: mm0 = digit * a[i] + r[i] + cy;
    movd    %mm0, 0(%edx)   / 1: r[i] = product[31..0]
    psrlq   $32, %mm0   / 1: cy = product[63..32]
    subl    $1, %ecx
    jz  .L27

    movd    4(%ebx), %mm1   / 2: mm1 = a[i]
    movd    4(%edx), %mm2   / 2: mm2 = r[i]
    pmuludq %mm3, %mm1  / 2: mm1 = digit * a[i]
    paddq   %mm1, %mm2  / 2: mm2 = digit * a[i] + r[i]
    paddq   %mm2, %mm0  / 2: mm0 = digit * a[i] + r[i] + cy;
    movd    %mm0, 4(%edx)   / 2: r[i] = product[31..0]
    psrlq   $32, %mm0   / 2: cy = product[63..32]
    subl    $1, %ecx
    jz  .L27

    movd    8(%ebx), %mm1   / 3: mm1 = a[i]
    movd    8(%edx), %mm2   / 3: mm2 = r[i]
    pmuludq %mm3, %mm1  / 3: mm1 = digit * a[i]
    paddq   %mm1, %mm2  / 3: mm2 = digit * a[i] + r[i]
    paddq   %mm2, %mm0  / 3: mm0 = digit * a[i] + r[i] + cy;
    movd    %mm0, 8(%edx)   / 3: r[i] = product[31..0]
    psrlq   $32, %mm0   / 3: cy = product[63..32]
    subl    $1, %ecx
    jz  .L27

    movd    12(%ebx), %mm1  / 4: mm1 = a[i]
    movd    12(%edx), %mm2  / 4: mm2 = r[i]
    pmuludq %mm3, %mm1  / 4: mm1 = digit * a[i]
    paddq   %mm1, %mm2  / 4: mm2 = digit * a[i] + r[i]
    paddq   %mm2, %mm0  / 4: mm0 = digit * a[i] + r[i] + cy;
    movd    %mm0, 12(%edx)  / 4: r[i] = product[31..0]
    psrlq   $32, %mm0   / 4: cy = product[63..32]
    subl    $1, %ecx
    jz  .L27

    movd    16(%ebx), %mm1  / 5: mm1 = a[i]
    movd    16(%edx), %mm2  / 5: mm2 = r[i]
    pmuludq %mm3, %mm1  / 5: mm1 = digit * a[i]
    paddq   %mm1, %mm2  / 5: mm2 = digit * a[i] + r[i]
    paddq   %mm2, %mm0  / 5: mm0 = digit * a[i] + r[i] + cy;
    movd    %mm0, 16(%edx)  / 5: r[i] = product[31..0]
    psrlq   $32, %mm0   / 5: cy = product[63..32]
    subl    $1, %ecx
    jz  .L27

    movd    20(%ebx), %mm1  / 6: mm1 = a[i]
    movd    20(%edx), %mm2  / 6: mm2 = r[i]
    pmuludq %mm3, %mm1  / 6: mm1 = digit * a[i]
    paddq   %mm1, %mm2  / 6: mm2 = digit * a[i] + r[i]
    paddq   %mm2, %mm0  / 6: mm0 = digit * a[i] + r[i] + cy;
    movd    %mm0, 20(%edx)  / 6: r[i] = product[31..0]
    psrlq   $32, %mm0   / 6: cy = product[63..32]
    subl    $1, %ecx
    jz  .L27

    movd    24(%ebx), %mm1  / 7: mm1 = a[i]
    movd    24(%edx), %mm2  / 7: mm2 = r[i]
    pmuludq %mm3, %mm1  / 7: mm1 = digit * a[i]
    paddq   %mm1, %mm2  / 7: mm2 = digit * a[i] + r[i]
    paddq   %mm2, %mm0  / 7: mm0 = digit * a[i] + r[i] + cy;
    movd    %mm0, 24(%edx)  / 7: r[i] = product[31..0]
    psrlq   $32, %mm0   / 7: cy = product[63..32]

.L27:
    movd    %mm0, %eax
    / no emms.  caller is responsible for emms
    ret
    SET_SIZE(big_mul_add_vec_sse2_r)


/ r = r + a * digit, r and a are vectors of length len
/ returns the carry digit
/ Suitable only for x86 models that support SSE2 instruction set extensions
/
/ r      8(%ebp)    %edx
/ a     12(%ebp)    %ebx
/ len       16(%ebp)    %ecx
/ digit     20(%ebp)    %mm3
/
/ In userland, there is just the one function, big_mul_add_vec_sse2().
/ But in the kernel, there are two variations:
/    1. big_mul_add_vec_sse2() which does what is necessary to save and
/       restore state, if necessary, and to ensure that preemtion is
/       disabled.
/    2. big_mul_add_vec_sse2_nsv() which just does the work;
/       it is the caller's responsibility to ensure that MMX state
/       does not need to be saved and restored and that preemption
/       is already disabled.


#if defined(MMX_MANAGE)

    ENTRY(big_mul_add_vec_sse2)
    pushl   %ebp
    movl    %esp, %ebp
    pushl   %ebx
    pushl   %esi
    KPREEMPT_DISABLE
    TEST_TS(%ebx)
    pushl   %ebx
    jnz .addvec_no_save
    pushl   %edi
    SAVE_MMX_0TO4(%edi)
    movl    8(%ebp), %edx
    movl    12(%ebp), %ebx
    movl    16(%ebp), %ecx
    movd    20(%ebp), %mm3
    call    big_mul_add_vec_sse2_r
    movl    %eax, %esi
    RSTOR_MMX_0TO4(%edi)
    popl    %edi
    jmp .addvec_rtn

.addvec_no_save:
    movl    8(%ebp), %edx
    movl    12(%ebp), %ebx
    movl    16(%ebp), %ecx
    movd    20(%ebp), %mm3
    call    big_mul_add_vec_sse2_r
    movl    %eax, %esi

.addvec_rtn:
    emms
    popl    %ebx
    movl    %ebx, %cr0
    KPREEMPT_ENABLE
    movl    %esi, %eax
    popl    %esi
    popl    %ebx
    leave
    ret
    SET_SIZE(big_mul_add_vec_sse2)

    ENTRY(big_mul_add_vec_sse2_nsv)
    pushl   %ebp
    movl    %esp, %ebp
    pushl   %ebx
    movl    8(%ebp), %edx
    movl    12(%ebp), %ebx
    movl    16(%ebp), %ecx
    movd    20(%ebp), %mm3
    call    big_mul_add_vec_sse2_r
    popl    %ebx
    leave
    ret
    SET_SIZE(big_mul_add_vec_sse2_nsv)


#else   /* !defined(MMX_MANAGE) */

    ENTRY(big_mul_add_vec_sse2)
    pushl   %ebp
    movl    %esp, %ebp
    pushl   %ebx
    movl    8(%ebp), %edx
    movl    12(%ebp), %ebx
    movl    16(%ebp), %ecx
    movd    20(%ebp), %mm3
    call    big_mul_add_vec_sse2_r
    popl    %ebx
    emms
    leave
    ret
    SET_SIZE(big_mul_add_vec_sse2)

#endif  /* MMX_MANAGE */


/ void
/ big_mul_vec_sse2(uint32_t *r, uint32_t *a, int alen, uint32_t *b, int blen)
/ {
/   int i;
/
/   r[alen] = big_mul_set_vec_sse2(r, a, alen, b[0]);
/   for (i = 1; i < blen; ++i)
/       r[alen + i] = big_mul_add_vec_sse2(r+i, a, alen, b[i]);
/ }


#if defined(MMX_MANAGE)
    ENTRY(big_mul_vec_sse2_fc)
#else
    ENTRY(big_mul_vec_sse2)
#endif
    subl    $0x8, %esp
    pushl   %ebx
    pushl   %ebp
    pushl   %esi
    pushl   %edi
    movl    40(%esp), %eax
    movl    %eax, 20(%esp)
    pushl   (%eax)
    movl    40(%esp), %edi
    pushl   %edi
    movl    40(%esp), %esi
    pushl   %esi
    movl    40(%esp), %ebx
    pushl   %ebx
#if defined(MMX_MANAGE)
    call    big_mul_set_vec_sse2_nsv
#else
    call    big_mul_set_vec_sse2
#endif
    addl    $0x10, %esp
    movl    %eax, (%ebx,%edi,4)
    movl    44(%esp), %eax
    movl    %eax, 16(%esp)
    cmpl    $0x1, %eax
    jle .mulvec_rtn
    movl    $0x1, %ebp

    .align 16
.mulvec_add:
    movl    20(%esp), %eax
    pushl   (%eax,%ebp,4)
    pushl   %edi
    pushl   %esi
    leal    (%ebx,%ebp,4), %eax
    pushl   %eax
#if defined(MMX_MANAGE)
    call    big_mul_add_vec_sse2_nsv
#else
    call    big_mul_add_vec_sse2
#endif
    addl    $0x10, %esp
    leal    (%ebp,%edi), %ecx
    movl    %eax, (%ebx,%ecx,4)
    incl    %ebp
    cmpl    16(%esp), %ebp
    jl  .mulvec_add
.mulvec_rtn:
#if defined(MMX_MANAGE)
    emms
#endif
    popl    %edi
    popl    %esi
    popl    %ebp
    popl    %ebx
    addl    $0x8, %esp
    ret
#if defined(MMX_MANAGE)
    SET_SIZE(big_mul_vec_sse2_fc)
#else
    SET_SIZE(big_mul_vec_sse2)
#endif

#if defined(MMX_MANAGE)

    ENTRY(big_mul_vec_sse2)
    pushl   %ebp
    movl    %esp, %ebp
    subl    $8, %esp
    pushl   %edi
    KPREEMPT_DISABLE
    TEST_TS(%eax)
    movl    %eax, -8(%ebp)
    jnz .mulvec_no_save
    SAVE_MMX_0TO4(%edi)
    movl    %edi, -4(%ebp)
.mulvec_no_save:
    movl    24(%ebp), %eax      / blen
    pushl   %eax
    movl    20(%ebp), %eax      / b
    pushl   %eax
    movl    16(%ebp), %eax      / alen
    pushl   %eax
    movl    12(%ebp), %eax      / a
    pushl   %eax
    movl    8(%ebp), %eax       / r
    pushl   %eax
    call    big_mul_vec_sse2_fc
    addl    $20, %esp
    movl    -8(%ebp), %eax
    testl   $CR0_TS, %eax
    jnz .mulvec_no_rstr
    movl    -4(%ebp), %edi
    RSTOR_MMX_0TO4(%edi)
.mulvec_no_rstr:
    movl    %eax, %cr0
    KPREEMPT_ENABLE
    popl    %edi
    leave
    ret
    SET_SIZE(big_mul_vec_sse2)

#endif  /* MMX_MANAGE */


#undef UNROLL
#undef UNROLL32


/ r = a * a, r and a are vectors of length len
/ Suitable only for x86 models that support SSE2 instruction set extensions
/
/ This function is not suitable for a truly general-purpose multiprecision
/ arithmetic library, because it does not work for "small" numbers, that is
/ numbers of 1 or 2 digits.  big_mul() just uses the ordinary big_mul_vec()
/ for any small numbers.

#if defined(MMX_MANAGE)
    ENTRY(big_sqr_vec_sse2_fc)
#else
    ENTRY(big_sqr_vec_sse2)
    pushl   %ebp
    movl    %esp, %ebp
#endif

    pushl   %ebx
    pushl   %edi
    pushl   %esi

    / r[1..alen] = a[0] * a[1..alen-1]

    movl    8(%ebp), %edi       / r = arg(r)
    movl    12(%ebp), %esi      / a = arg(a)
    movl    16(%ebp), %ecx      / cnt = arg(alen)
    movd    %ecx, %mm4      / save_cnt = arg(alen)
    leal    4(%edi), %edx       / dst = &r[1]
    movl    %esi, %ebx      / src = a
    movd    0(%ebx), %mm3       / mm3 = a[0]
    leal    4(%ebx), %ebx       / src = &a[1]
    subl    $1, %ecx        / --cnt
    call    big_mul_set_vec_sse2_r  / r[1..alen-1] = a[0] * a[1..alen-1]
    movl    %edi, %edx      / dst = r
    movl    %esi, %ebx      / src = a
    movd    %mm4, %ecx      / cnt = save_cnt
    movl    %eax, (%edx, %ecx, 4)   / r[cnt] = cy

/   /* High-level vector C pseudocode */
/   for (i = 1; i < alen-1; ++i)
/       r[2*i + 1 ... ] += a[i] * a[i+1 .. alen-1]
/
/   /* Same thing, but slightly lower level C-like pseudocode */
/   i = 1;
/   r = &arg_r[2*i + 1];
/   a = &arg_a[i + 1];
/   digit = arg_a[i];
/   cnt = alen - 3;
/   while (cnt != 0) {
/       r[cnt] = big_mul_add_vec_sse2_r(r, a, cnt, digit);
/       r += 2;
/       ++a;
/       --cnt;
/   }
/
/   /* Same thing, but even lower level
/    * For example, pointers are raw pointers,
/    * with no scaling by object size.
/    */
/   r = arg_r + 12; /* i == 1; 2i + 1 == 3;  4*3 == 12; */
/   a = arg_a + 8;
/   digit = *(arg_a + 4);
/   cnt = alen - 3;
/   while (cnt != 0) {
/       cy = big_mul_add_vec_sse2_r();
/       *(r + 4 * cnt) = cy;
/       r += 8;
/       a += 4;
/       --cnt;
/   }

    leal    4(%edi), %edi       / r += 4; r = &r[1]
    leal    4(%esi), %esi       / a += 4; a = &a[1]
    movd    %mm4, %ecx      / cnt = save
    subl    $2, %ecx        / cnt = alen - 2; i in 1..alen-2
    movd    %ecx, %mm4      / save_cnt
    jecxz   .L32            / while (cnt != 0) {
.L31:
    movd    0(%esi), %mm3       / digit = a[i]
    leal    4(%esi), %esi       / a += 4; a = &a[1]; a = &a[i + 1]
    leal    8(%edi), %edi       / r += 8; r = &r[2]; r = &r[2 * i + 1]
    movl    %edi, %edx      / edx = r
    movl    %esi, %ebx      / ebx = a
    cmp $1, %ecx        / The last triangle term is special
    jz  .L32
    call    big_mul_add_vec_sse2_r
    movd    %mm4, %ecx      / cnt = save_cnt
    movl    %eax, (%edi, %ecx, 4)   / r[cnt] = cy
    subl    $1, %ecx        / --cnt
    movd    %ecx, %mm4      / save_cnt = cnt
    jmp .L31            / }

.L32:
    movd    0(%ebx), %mm1       / mm1 = a[i + 1]
    movd    0(%edx), %mm2       / mm2 = r[2 * i + 1]
    pmuludq %mm3, %mm1      / mm1 = p = digit * a[i + 1]
    paddq   %mm1, %mm2      / mm2 = r[2 * i + 1] + p
    movd    %mm2, 0(%edx)       / r[2 * i + 1] += lo32(p)
    psrlq   $32, %mm2       / mm2 = cy
    movd    %mm2, 4(%edx)       / r[2 * i + 2] = cy
    pxor    %mm2, %mm2
    movd    %mm2, 8(%edx)       / r[2 * i + 3] = 0

    movl    8(%ebp), %edx       / r = arg(r)
    movl    12(%ebp), %ebx      / a = arg(a)
    movl    16(%ebp), %ecx      / cnt = arg(alen)

    / compute low-order corner
    / p = a[0]**2
    / r[0] = lo32(p)
    / cy   = hi32(p)
    movd    0(%ebx), %mm2       / mm2 = a[0]
    pmuludq %mm2, %mm2      / mm2 = p = a[0]**2
    movd    %mm2, 0(%edx)       / r[0] = lo32(p)
    psrlq   $32, %mm2       / mm2 = cy = hi32(p)

    / p = 2 * r[1]
    / t = p + cy
    / r[1] = lo32(t)
    / cy   = hi32(t)
    movd    4(%edx), %mm1       / mm1 = r[1]
    psllq   $1, %mm1        / mm1 = p = 2 * r[1]
    paddq   %mm1, %mm2      / mm2 = t = p + cy
    movd    %mm2, 4(%edx)       / r[1] = low32(t)
    psrlq   $32, %mm2       / mm2 = cy = hi32(t)

    / r[2..$-3] = inner_diagonal[*]**2 + 2 * r[2..$-3]
    subl    $2, %ecx        / cnt = alen - 2
.L34:
    movd    4(%ebx), %mm0       / mm0 = diag = a[i+1]
    pmuludq %mm0, %mm0      / mm0 = p = diag**2
    paddq   %mm0, %mm2      / mm2 = t = p + cy
    movd    %mm2, %eax
    movd    %eax, %mm1      / mm1 = lo32(t)
    psrlq   $32, %mm2       / mm2 = hi32(t)

    movd    8(%edx), %mm3       / mm3 = r[2*i]
    psllq   $1, %mm3        / mm3 = 2*r[2*i]
    paddq   %mm3, %mm1      / mm1 = 2*r[2*i] + lo32(t)
    movd    %mm1, 8(%edx)       / r[2*i] = 2*r[2*i] + lo32(t)
    psrlq   $32, %mm1
    paddq   %mm1, %mm2

    movd    12(%edx), %mm3      / mm3 = r[2*i+1]
    psllq   $1, %mm3        / mm3 = 2*r[2*i+1]
    paddq   %mm3, %mm2      / mm2 = 2*r[2*i+1] + hi32(t)
    movd    %mm2, 12(%edx)      / r[2*i+1] = mm2
    psrlq   $32, %mm2       / mm2 = cy
    leal    8(%edx), %edx       / r += 2
    leal    4(%ebx), %ebx       / ++a
    subl    $1, %ecx        / --cnt
    jnz .L34

    / Carry from last triangle term must participate in doubling,
    / but this step isn't paired up with a squaring the elements
    / of the inner diagonal.
    / r[$-3..$-2] += 2 * r[$-3..$-2] + cy
    movd    8(%edx), %mm3       / mm3 = r[2*i]
    psllq   $1, %mm3        / mm3 = 2*r[2*i]
    paddq   %mm3, %mm2      / mm2 = 2*r[2*i] + cy
    movd    %mm2, 8(%edx)       / r[2*i] = lo32(2*r[2*i] + cy)
    psrlq   $32, %mm2       / mm2 = cy = hi32(2*r[2*i] + cy)

    movd    12(%edx), %mm3      / mm3 = r[2*i+1]
    psllq   $1, %mm3        / mm3 = 2*r[2*i+1]
    paddq   %mm3, %mm2      / mm2 = 2*r[2*i+1] + cy
    movd    %mm2, 12(%edx)      / r[2*i+1] = mm2
    psrlq   $32, %mm2       / mm2 = cy

    / compute high-order corner and add it in
    / p = a[alen - 1]**2
    / t = p + cy
    / r[alen + alen - 2] += lo32(t)
    / cy = hi32(t)
    / r[alen + alen - 1] = cy
    movd    4(%ebx), %mm0       / mm0 = a[$-1]
    movd    8(%edx), %mm3       / mm3 = r[$-2]
    pmuludq %mm0, %mm0      / mm0 = p = a[$-1]**2
    paddq   %mm0, %mm2      / mm2 = t = p + cy
    paddq   %mm3, %mm2      / mm2 = r[$-2] + t
    movd    %mm2, 8(%edx)       / r[$-2] = lo32(r[$-2] + t)
    psrlq   $32, %mm2       / mm2 = cy = hi32(r[$-2] + t)
    movd    12(%edx), %mm3
    paddq   %mm3, %mm2
    movd    %mm2, 12(%edx)      / r[$-1] += cy

.L35:
    emms
    popl    %esi
    popl    %edi
    popl    %ebx

#if defined(MMX_MANAGE)
    ret
    SET_SIZE(big_sqr_vec_sse2_fc)
#else
    leave
    ret
    SET_SIZE(big_sqr_vec_sse2)
#endif


#if defined(MMX_MANAGE)
    ENTRY(big_sqr_vec_sse2)
    pushl   %ebp
    movl    %esp, %ebp
    KPREEMPT_DISABLE
    TEST_TS(%ebx)
    pushl   %ebx
    jnz .sqr_no_save
    pushl   %edi
    SAVE_MMX_0TO4(%edi)
    call    big_sqr_vec_sse2_fc
    RSTOR_MMX_0TO4(%edi)
    popl    %edi
    jmp .sqr_rtn

.sqr_no_save:
    call    big_sqr_vec_sse2_fc

.sqr_rtn:
    popl    %ebx
    movl    %ebx, %cr0
    KPREEMPT_ENABLE
    leave
    ret
    SET_SIZE(big_sqr_vec_sse2)

#endif  /* MMX_MANAGE */

/ ------------------------------------------------------------------------
/       UMUL Implementations
/ ------------------------------------------------------------------------


/ r = a * digit, r and a are vectors of length len
/ returns the carry digit
/ Does not use any MMX, SSE, or SSE2 instructions.
/ Uses x86 unsigned 32 X 32 -> 64 multiply instruction, MUL.
/ This is a fall-back implementation for x86 models that do not support
/ the PMULUDQ instruction.
/
/ uint32_t
/ big_mul_set_vec_umul(uint32_t *r, uint32_t *a, int len, uint32_t digit)
/
/ r      8(%ebp)    %edx    %edi
/ a     12(%ebp)    %ebx    %esi
/ len       16(%ebp)    %ecx
/ digit     20(%ebp)    %esi

    ENTRY(big_mul_set_vec_umul)
    pushl   %ebp
    movl    %esp, %ebp
    pushl   %esi
    pushl   %edi
    pushl   %ebx
    movl    16(%ebp), %ecx
    xorl    %ebx, %ebx  / cy = 0
    testl   %ecx, %ecx
    movl    8(%ebp), %edi
    movl    12(%ebp), %esi
    je  .L57

.L55:
    movl    (%esi), %eax    / eax = a[i]
    leal    4(%esi), %esi   / ++a
    mull    20(%ebp)    / edx:eax = a[i] * digit
    addl    %ebx, %eax
    adcl    $0, %edx    / edx:eax = a[i] * digit + cy
    movl    %eax, (%edi)    / r[i] = product[31..0]
    movl    %edx, %ebx  / cy = product[63..32]
    leal    4(%edi), %edi   / ++r
    decl    %ecx        / --len
    jnz .L55        / while (len != 0)
.L57:
    movl    %ebx, %eax
    popl    %ebx
    popl    %edi
    popl    %esi
    leave
    ret
    SET_SIZE(big_mul_set_vec_umul)


/ r = r + a * digit, r and a are vectors of length len
/ returns the carry digit
/ Does not use any MMX, SSE, or SSE2 instructions.
/ Uses x86 unsigned 32 X 32 -> 64 multiply instruction, MUL.
/ This is a fall-back implementation for x86 models that do not support
/ the PMULUDQ instruction.
/
/ uint32_t
/ big_mul_add_vec_umul(uint32_t *r, uint32_t *a, int len, uint32_t digit)
/
/ r      8(%ebp)    %edx    %edi
/ a     12(%ebp)    %ebx    %esi
/ len       16(%ebp)    %ecx
/ digit     20(%ebp)    %esi

    ENTRY(big_mul_add_vec_umul)
    pushl   %ebp
    movl    %esp, %ebp
    pushl   %esi
    pushl   %edi
    pushl   %ebx
    movl    16(%ebp), %ecx
    xorl    %ebx, %ebx  / cy = 0
    testl   %ecx, %ecx
    movl    8(%ebp), %edi
    movl    12(%ebp), %esi
    je  .L67
    .align 4
.L65:
    movl    (%esi), %eax    / eax = a[i]
    leal    4(%esi), %esi   / ++a
    mull    20(%ebp)    / edx:eax = a[i] * digit
    addl    (%edi), %eax
    adcl    $0, %edx    / edx:eax = a[i] * digit + r[i]
    addl    %ebx, %eax
    adcl    $0, %edx    / edx:eax = a[i] * digit + r[i] + cy
    movl    %eax, (%edi)    / r[i] = product[31..0]
    movl    %edx, %ebx  / cy = product[63..32]
    leal    4(%edi), %edi   / ++r
    decl    %ecx        / --len
    jnz .L65        / while (len != 0)
.L67:
    movl    %ebx, %eax
    popl    %ebx
    popl    %edi
    popl    %esi
    leave
    ret
    SET_SIZE(big_mul_add_vec_umul)

#endif  /* __lint */