@ Copyright (c) 2009-2013 The Linux Foundation. All rights reserved.
@
@ Redistribution and use in source and binary forms, with or without
@ modification, are permitted provided that the following conditions are met:
@     * Redistributions of source code must retain the above copyright
@       notice, this list of conditions and the following disclaimer.
@     * Redistributions in binary form must reproduce the above copyright
@       notice, this list of conditions and the following disclaimer in the
@       documentation and/or other materials provided with the distribution.
@     * Neither the name of The Linux Foundation nor the names of its contributors may
@       be used to endorse or promote products derived from this software
@       without specific prior written permission.
@
@ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
@ AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
@ IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
@ ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
@ LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
@ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
@ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
@ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
@ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
@ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
@ POSSIBILITY OF SUCH DAMAGE.


#include <private/bionic_asm.h>
#include <private/libc_events.h>

@ Values which exist the program lifetime:
#define HIGH_WORD_MASK      d31
#define EXPONENT_MASK       d30
#define int_1               d29
#define double_1            d28
@ sign and 2^int_n fixup:
#define maxrange            r12
#define expadjustment       d7
#define literals            r10
@ Values which exist within both polynomial implementations:
#define int_n               d2
#define int_n_low           s4
#define int_n_high          s5
#define double_n            d3
#define k1                  d27
#define k2                  d26
#define k3                  d25
#define k4                  d24
@ Values which cross the boundaries between polynomial implementations:
#define ss                  d16
#define ss2                 d17
#define ss4                 d18
#define Result              d0
#define Return_hw           r1
#define Return_lw           r0
#define ylg2x               d0
@ Intermediate values only needed sometimes:
@ initial (sorted in approximate order of availability for overwriting):
#define x_hw                r1
#define x_lw                r0
#define y_hw                r3
#define y_lw                r2
#define x                   d0
#define bp                  d4
#define y                   d1
@ log series:
#define u                   d19
#define v                   d20
#define lg2coeff            d21
#define bpa                 d5
#define bpb                 d3
#define lg2const            d6
#define xmantissa           r8
#define twoto1o5            r4
#define twoto3o5            r5
#define ix                  r6
#define iEXP_MASK           r7
@ exp input setup:
#define twoto1o8mask        d3
#define twoto1o4mask        d4
#define twoto1o2mask        d1
#define ylg2x_round_offset  d16
#define ylg2x_temp          d17
#define yn_temp             d18
#define yn_round_offset     d19
#define ln2                 d5
@ Careful, overwriting HIGH_WORD_MASK, reset it if you need it again ...
#define rounded_exponent    d31
@ exp series:
#define k5                  d23
#define k6                  d22
#define k7                  d21
#define k8                  d20
#define ss3                 d19
@ overwrite double_1 (we're done with it by now)
#define k0                  d28
#define twoto1o4            d6

@instructions that gas doesn't like to encode correctly:
#define vmov_f64            fconstd
#define vmov_f32            fconsts
#define vmovne_f64          fconstdne

#define KRAIT_NO_AAPCS_VFP_MODE

ENTRY(pow)
#if defined(KRAIT_NO_AAPCS_VFP_MODE)
     @ ARM ABI has inputs coming in via r registers, lets move to a d register
    vmov            x, x_lw, x_hw
#endif
    push            {r4, r5, r6, r7, r8, r9, r10, lr}

    movw            maxrange, #0x0000
    movt            maxrange, #0x4010

    @ pre-staged bp values
    vldr            bpa, .LbpA
    vldr            bpb, .LbpB
    @ load two fifths into constant term in case we need it due to offsets
    vldr            lg2const, .Ltwofifths

    @ bp is initially 1.0, may adjust later based on x value
    vmov_f64        bp,  #0x70

    @ extract the mantissa from x for scaled value comparisons
    lsl             xmantissa, x_hw, #12

    @ twoto1o5 = 2^(1/5) (input bracketing)
    movw            twoto1o5, #0x186c
    movt            twoto1o5, #0x2611
    @ twoto3o5 = 2^(3/5) (input bracketing)
    movw            twoto3o5, #0x003b
    movt            twoto3o5, #0x8406

    @ finish extracting xmantissa
    orr             xmantissa, xmantissa, x_lw, lsr #20

    @ begin preparing a mask for normalization
    vmov.i64        HIGH_WORD_MASK, #0xffffffff00000000

    @ double_1 = (double) 1.0
    vmov_f64        double_1, #0x70

#if defined(KRAIT_NO_AAPCS_VFP_MODE)
     @ move y from r registers to a d register
    vmov            y, y_lw, y_hw
#endif

    cmp             xmantissa, twoto1o5

    vshl.i64        EXPONENT_MASK, HIGH_WORD_MASK, #20
    vshr.u64        int_1, HIGH_WORD_MASK, #63

    adr             literals, .LliteralTable

    bhi             .Lxgt2to1over5
    @ zero out lg2 constant term if don't offset our input
    vsub.f64        lg2const, lg2const, lg2const
    b               .Lxle2to1over5

.Lxgt2to1over5:
    @ if normalized x > 2^(1/5), bp = 1 + (2^(2/5)-1) = 2^(2/5)
    vadd.f64        bp, bp, bpa

.Lxle2to1over5:
    @ will need ln2 for various things
    vldr            ln2, .Lln2

    cmp             xmantissa, twoto3o5
@@@@ X Value Normalization @@@@

    @ ss = abs(x) 2^(-1024)
    vbic.i64        ss, x, EXPONENT_MASK

    @ N = (floor(log2(x)) + 0x3ff) * 2^52
    vand.i64        int_n, x, EXPONENT_MASK

    bls             .Lxle2to3over5
    @ if normalized x > 2^(3/5), bp = 2^(2/5) + (2^(4/5) - 2^(2/5) = 2^(4/5)
    vadd.f64      bp, bp, bpb
    vadd.f64      lg2const, lg2const, lg2const

.Lxle2to3over5:

    cmp             x_hw, maxrange
    cmpls           y_hw, maxrange
    movt            maxrange, #0x3f00
    cmpls           maxrange, x_hw

    @ load log2 polynomial series constants
    vldm            literals!, {k4, k3, k2, k1}

    @ s = abs(x) 2^(-floor(log2(x))) (normalize abs(x) to around 1)
    vorr.i64        ss, ss, double_1

@@@@ 3/2 (Log(bp(1+s)/(1-s))) input computation (s = (x-bp)/(x+bp)) @@@@

    vsub.f64        u, ss, bp
    vadd.f64        v, ss, bp

    bhi             .LuseFullImpl

    @ s = (x-1)/(x+1)
    vdiv.f64        ss, u, v

    @ load 2/(3log2) into lg2coeff
    vldr            lg2coeff, .Ltwooverthreeln2

    @ N = floor(log2(x)) * 2^52
    vsub.i64        int_n, int_n, double_1

@@@@ 3/2 (Log(bp(1+s)/(1-s))) polynomial series @@@@

    @ ss2 = ((x-dp)/(x+dp))^2
    vmul.f64        ss2, ss, ss
    @ ylg2x = 3.0
    vmov_f64        ylg2x, #8
    vmul.f64        ss4, ss2, ss2

    @ todo: useful later for two-way clamp
    vmul.f64        lg2coeff, lg2coeff, y

    @ N = floor(log2(x))
    vshr.s64        int_n, int_n, #52

    @ k3 = ss^2 * L4 + L3
    vmla.f64        k3, ss2, k4

    @ k1 = ss^2 * L2 + L1
    vmla.f64        k1, ss2, k2

    @ scale ss by 2/(3 ln 2)
    vmul.f64        lg2coeff, ss, lg2coeff

    @ ylg2x = 3.0 + s^2
    vadd.f64        ylg2x, ylg2x, ss2

    vcvt.f64.s32    double_n, int_n_low

    @ k1 = s^4 (s^2 L4 + L3) + s^2 L2 + L1
    vmla.f64        k1, ss4, k3

    @ add in constant term
    vadd.f64        double_n, lg2const

    @ ylg2x = 3.0 + s^2 + s^4 (s^4 (s^2 L4 + L3) + s^2 L2 + L1)
    vmla.f64        ylg2x, ss4, k1

    @ ylg2x = y 2 s / (3 ln(2)) (3.0 + s^2 + s^4 (s^4(s^2 L4 + L3) + s^2 L2 + L1)
    vmul.f64        ylg2x, lg2coeff, ylg2x

@@@@ Compute input to Exp(s) (s = y(n + log2(x)) - (floor(8 yn + 1)/8 + floor(8 ylog2(x) + 1)/8) @@@@@

    @ mask to extract bit 1 (2^-2 from our fixed-point representation)
    vshl.u64        twoto1o4mask, int_1, #1

    @ double_n = y * n
    vmul.f64        double_n, double_n, y

    @ Load 2^(1/4) for later computations
    vldr            twoto1o4, .Ltwoto1o4

    @ either add or subtract one based on the sign of double_n and ylg2x
    vshr.s64        ylg2x_round_offset, ylg2x, #62
    vshr.s64        yn_round_offset, double_n, #62

    @ move unmodified y*lg2x into temp space
    vmov            ylg2x_temp, ylg2x
    @ compute floor(8 y * n + 1)/8
    @ and floor(8 y (log2(x)) + 1)/8
    vcvt.s32.f64    ylg2x, ylg2x, #3
    @ move unmodified y*n into temp space
    vmov            yn_temp, double_n
    vcvt.s32.f64    double_n, double_n, #3

    @ load exp polynomial series constants
    vldm            literals!, {k8, k7, k6, k5, k4, k3, k2, k1}

    @ mask to extract bit 2 (2^-1 from our fixed-point representation)
    vshl.u64        twoto1o2mask, int_1, #2

    @ make rounding offsets either 1 or -1 instead of 0 or -2
    vorr.u64        ylg2x_round_offset, ylg2x_round_offset, int_1
    vorr.u64        yn_round_offset, yn_round_offset, int_1

    @ round up to the nearest 1/8th
    vadd.s32        ylg2x, ylg2x, ylg2x_round_offset
    vadd.s32        double_n, double_n, yn_round_offset

    @ clear out round-up bit for y log2(x)
    vbic.s32        ylg2x, ylg2x, int_1
    @ clear out round-up bit for yn
    vbic.s32        double_n, double_n, int_1
    @ add together the (fixed precision) rounded parts
    vadd.s64        rounded_exponent, double_n, ylg2x
    @ turn int_n into a double with value 2^int_n
    vshl.i64        int_n, rounded_exponent, #49
    @ compute masks for 2^(1/4) and 2^(1/2) fixups for fractional part of fixed-precision rounded values:
    vand.u64        twoto1o4mask, twoto1o4mask, rounded_exponent
    vand.u64        twoto1o2mask, twoto1o2mask, rounded_exponent

    @ convert back into floating point, double_n now holds (double) floor(8 y * n + 1)/8
    @                                   ylg2x now holds (double) floor(8 y * log2(x) + 1)/8
    vcvt.f64.s32    ylg2x, ylg2x, #3
    vcvt.f64.s32    double_n, double_n, #3

    @ put the 2 bit (0.5) through the roof of twoto1o2mask (make it 0x0 or 0xffffffffffffffff)
    vqshl.u64        twoto1o2mask, twoto1o2mask, #62
    @ put the 1 bit (0.25) through the roof of twoto1o4mask (make it 0x0 or 0xffffffffffffffff)
    vqshl.u64        twoto1o4mask, twoto1o4mask, #63

    @ center y*log2(x) fractional part between -0.125 and 0.125 by subtracting (double) floor(8 y * log2(x) + 1)/8
    vsub.f64        ylg2x_temp, ylg2x_temp, ylg2x
    @ center y*n fractional part between -0.125 and 0.125 by subtracting (double) floor(8 y * n + 1)/8
    vsub.f64        yn_temp, yn_temp, double_n

    @ Add fractional parts of yn and y log2(x) together
    vadd.f64        ss, ylg2x_temp, yn_temp

    @ Result = 1.0 (offset for exp(s) series)
    vmov_f64        Result, #0x70

    @ multiply fractional part of y * log2(x) by ln(2)
    vmul.f64        ss, ln2, ss

@@@@ 10th order polynomial series for Exp(s) @@@@

    @ ss2 = (ss)^2
    vmul.f64        ss2, ss, ss

    @ twoto1o2mask = twoto1o2mask & twoto1o4
    vand.u64        twoto1o2mask, twoto1o2mask, twoto1o4
    @ twoto1o2mask = twoto1o2mask & twoto1o4
    vand.u64        twoto1o4mask, twoto1o4mask, twoto1o4

    @ Result = 1.0 + ss
    vadd.f64        Result, Result, ss

    @ k7 = ss k8 + k7
    vmla.f64        k7, ss, k8

    @ ss4 = (ss*ss) * (ss*ss)
    vmul.f64        ss4, ss2, ss2

    @ twoto1o2mask = twoto1o2mask | (double) 1.0 - results in either 1.0 or 2^(1/4) in twoto1o2mask
    vorr.u64        twoto1o2mask, twoto1o2mask, double_1
    @ twoto1o2mask = twoto1o4mask | (double) 1.0 - results in either 1.0 or 2^(1/4) in twoto1o4mask
    vorr.u64        twoto1o4mask, twoto1o4mask, double_1

    @ TODO: should setup sign here, expadjustment = 1.0
    vmov_f64        expadjustment, #0x70

    @ ss3 = (ss*ss) * ss
    vmul.f64        ss3, ss2, ss

    @ k0 = 1/2 (first non-unity coefficient)
    vmov_f64        k0, #0x60

    @ Mask out non-exponent bits to make sure we have just 2^int_n
    vand.i64        int_n, int_n, EXPONENT_MASK

    @ square twoto1o2mask to get 1.0 or 2^(1/2)
    vmul.f64        twoto1o2mask, twoto1o2mask, twoto1o2mask
    @ multiply twoto2o4mask into the exponent output adjustment value
    vmul.f64        expadjustment, expadjustment, twoto1o4mask

    @ k5 = ss k6 + k5
    vmla.f64        k5, ss, k6

    @ k3 = ss k4 + k3
    vmla.f64        k3, ss, k4

    @ k1 = ss k2 + k1
    vmla.f64        k1, ss, k2

    @ multiply twoto1o2mask into exponent output adjustment value
    vmul.f64        expadjustment, expadjustment, twoto1o2mask

    @ k5 = ss^2 ( ss k8 + k7 ) + ss k6 + k5
    vmla.f64        k5, ss2, k7

    @ k1 = ss^2 ( ss k4 + k3 ) + ss k2 + k1
    vmla.f64        k1, ss2, k3

    @ Result = 1.0 + ss + 1/2 ss^2
    vmla.f64      Result, ss2, k0

    @ Adjust int_n so that it's a double precision value that can be multiplied by Result
    vadd.i64        expadjustment, int_n, expadjustment

    @ k1 = ss^4 ( ss^2 ( ss k8 + k7 ) + ss k6 + k5 ) + ss^2 ( ss k4 + k3 ) + ss k2 + k1
    vmla.f64        k1, ss4, k5

    @ Result = 1.0 + ss + 1/2 ss^2 + ss^3 ( ss^4 ( ss^2 ( ss k8 + k7 ) + ss k6 + k5 ) + ss^2 ( ss k4 + k3 ) + ss k2 + k1 )
    vmla.f64        Result, ss3, k1

    @ multiply by adjustment (sign*(rounding ? sqrt(2) : 1) * 2^int_n)
    vmul.f64        Result, expadjustment, Result

.LleavePow:
#if defined(KRAIT_NO_AAPCS_VFP_MODE)
    @ return Result (FP)
    vmov            Return_lw, Return_hw, Result
#endif
.LleavePowDirect:
    @ leave directly returning whatever is in Return_lw and Return_hw
    pop             {r4, r5, r6, r7, r8, r9, r10, pc}

.LuseFullImpl:
    pop             {r4, r5, r6, r7, r8, r9, r10, lr}
    b               __full_ieee754_pow

.align 6
.LliteralTable:
@ Least-sqares tuned constants for 11th order (log2((1+s)/(1-s)):
.LL4: @ ~3/11
    .long       0x53a79915, 0x3fd1b108
.LL3: @ ~1/3
    .long       0x9ca0567a, 0x3fd554fa
.LL2: @ ~3/7
    .long       0x1408e660, 0x3fdb6db7
.LL1: @ ~3/5
    .long       0x332D4313, 0x3fe33333

@ Least-squares tuned constants for 10th order exp(s):
.LE10: @ ~1/3628800
    .long       0x25c7ba0a, 0x3e92819b
.LE9: @ ~1/362880
    .long       0x9499b49c, 0x3ec72294
.LE8: @ ~1/40320
    .long       0xabb79d95, 0x3efa019f
.LE7: @ ~1/5040
    .long       0x8723aeaa, 0x3f2a019f
.LE6: @ ~1/720
    .long       0x16c76a94, 0x3f56c16c
.LE5: @ ~1/120
    .long       0x11185da8, 0x3f811111
.LE4: @ ~1/24
    .long       0x5555551c, 0x3fa55555
.LE3: @ ~1/6
    .long       0x555554db, 0x3fc55555

.LbpA: @ (2^(2/5) - 1)
    .long       0x4ee54db1, 0x3fd472d1

.LbpB: @ (2^(4/5) - 2^(2/5))
    .long       0x1c8a36cf, 0x3fdafb62

.Ltwofifths: @
    .long       0x9999999a, 0x3fd99999

.Ltwooverthreeln2:
    .long       0xDC3A03FD, 0x3FEEC709

.Lln2: @ ln(2)
    .long       0xFEFA39EF, 0x3FE62E42

.Ltwoto1o4: @ 2^1/4
    .long       0x0a31b715, 0x3ff306fe
END(pow)