andre/cosmopolitan
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Make pow() conform to standard definition

Browse Source
This commit is contained in:
Justine Tunney
2021-03-03 09:05:21 -08:00
parent 754974faaa
commit 8af91bcbe7
15 changed files with 253 additions and 85 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
libc/tinymath/copysign.S
View File

@ -36,9 +36,9 @@ copysign:
.rodata.cst16
.Lnan: .long 0xffffffff
.long 0x7fffffff
.long 0
.long 0
.Lneg0: .long 0
.long -2147483648
.long 0
.long 0
.long 0x00000000
.long 0x00000000
.Lneg0: .long 0x00000000
.long 0x80000000
.long 0x00000000
.long 0x00000000

1
libc/tinymath/cos.S
View File

@ -17,7 +17,6 @@
│ PERFORMANCE OF THIS SOFTWARE. │
╚─────────────────────────────────────────────────────────────────────────────*/
#include "libc/macros.internal.h"
.source __FILE__
// Returns cosine of 𝑥.
//

3
libc/tinymath/exp10.S
View File

@ -23,8 +23,7 @@
// @param 𝑥 is double scalar in low half of %xmm0
// @return double scalar in low half of %xmm0
// @see pow(), exp()
exp10:
ezlea exp10l,ax
exp10: ezlea exp10l,ax
jmp _d2ld2
.endfn exp10,globl
.alias exp10,pow10

1
libc/tinymath/fabsl.S
View File

@ -17,7 +17,6 @@
│ PERFORMANCE OF THIS SOFTWARE. │
╚─────────────────────────────────────────────────────────────────────────────*/
#include "libc/macros.internal.h"
.source __FILE__
// Returns absolute value of 𝑥.
//

67
libc/tinymath/powl.c Normal file
View File

@ -0,0 +1,67 @@
/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi│
╞══════════════════════════════════════════════════════════════════════════════╡
│ Copyright 2021 Justine Alexandra Roberts Tunney │
│ │
│ 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. │
╚─────────────────────────────────────────────────────────────────────────────*/
#include "libc/math.h"
/**
* Returns 𝑥^𝑦.
*/
long double powl(long double x, long double y) {
long double t, u;
if (!isunordered(x, y)) {
if (!isinf(y)) {
if (!isinf(x)) {
if (x) {
asm("fyl2x" : "=t"(u) : "0"(fabsl(x)), "u"(y) : "st(1)");
asm("fprem" : "=t"(t) : "0"(u), "u"(1.L));
asm("f2xm1" : "=t"(t) : "0"(t));
asm("fscale" : "=t"(t) : "0"(t + 1), "u"(u));
return copysignl(t, x);
} else if (y > 0) {
return 0;
} else if (!y) {
return 1;
} else if (y == truncl(y) && ((int64_t)y & 1)) {
return copysignl(INFINITY, x);
} else {
return INFINITY;
}
} else if (signbit(x)) {
if (!y) return 1;
x = y < 0 ? 0 : INFINITY;
if (y == truncl(y) && ((int64_t)y & 1)) x = -x;
return x;
} else if (y < 0) {
return 0;
} else if (y > 0) {
return INFINITY;
} else {
return 1;
}
} else {
x = fabsl(x);
if (x < 1) return signbit(y) ? INFINITY : 0;
if (x > 1) return signbit(y) ? 0 : INFINITY;
return 1;
}
} else if (!y || x == 1) {
return 1;
} else {
return NAN;
}
}

43
libc/tinymath/powl.S → libc/tinymath/powlfinite.c
View File

@ -1,7 +1,7 @@
/*-*- mode:unix-assembly; indent-tabs-mode:t; tab-width:8; coding:utf-8 -*-│
│vi: set et ft=asm ts=8 tw=8 fenc=utf-8 :vi│
/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi│
╞══════════════════════════════════════════════════════════════════════════════╡
│ Copyright 2020 Justine Alexandra Roberts Tunney │
│ Copyright 2021 Justine Alexandra Roberts Tunney │
│ │
│ Permission to use, copy, modify, and/or distribute this software for │
│ any purpose with or without fee is hereby granted, provided that the │
@ -16,37 +16,8 @@
│ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR │
│ PERFORMANCE OF THIS SOFTWARE. │
╚─────────────────────────────────────────────────────────────────────────────*/
#include "libc/runtime/pc.internal.h"
#include "libc/macros.internal.h"
#include "libc/math.h"
// Returns 𝑥^𝑦.
//
// @param 𝑥 is an 80-bit long double passed on stack in 16-bytes
// @param 𝑦 is the power, also pushed on stack, in reverse order
// @return result of exponentiation on FPU stack in %st
// @note Sun's fdlibm needs 2kLOC to do this for RISC lool
// @define z=y*log2(fabs(x)),copysign(trunc(exp2(fmod(z,1)))*exp2(z),x)
powl: push %rbp
mov %rsp,%rbp
.profilable
fldt 32(%rbp)
fldt 16(%rbp)
fxam
fstsw
fabs
fyl2x
fld1
fld %st(1)
fprem
f2xm1
faddp
fscale
fxch
fstp %st
test $FPU_C1>>8,%ah
jz 1f
fchs
1: pop %rbp
ret
.endfn powl,globl
.alias powl,__powl_finite
long double __powl_finite(long double x, long double y) {
return powl(x, y);
}

25
libc/tinymath/truncl.S
View File

@ -17,18 +17,23 @@
│ PERFORMANCE OF THIS SOFTWARE. │
╚─────────────────────────────────────────────────────────────────────────────*/
#include "libc/macros.internal.h"
.source __FILE__
truncl: .profilable
sub $24,%rsp
fldt 32(%rsp)
fnstcw 14(%rsp)
movzwl 14(%rsp),%eax
// Rounds to integer, toward zero.
//
// @param 𝑥 is long double passed on stack
// @return long double in %st
truncl: pushq %rbp
mov %rsp,%rbp
.profilable
sub $16,%rsp
fnstcw -2(%rbp)
fldt 16(%rbp)
movzwl -2(%rbp),%eax
or $0b1100,%ah # round to zero
mov %ax,12(%rsp)
fldcw 12(%rsp)
mov %ax,-4(%rbp)
fldcw -4(%rbp)
frndint
fldcw 14(%rsp)
add $24,%rsp
fldcw -2(%rbp)
leave
ret
.endfn truncl,globl