diff --git a/packages/compiler/src/backend/link-plan-core.test.ts b/packages/compiler/src/backend/link-plan-core.test.ts new file mode 100644 index 000000000..9c33e6fd0 --- /dev/null +++ b/packages/compiler/src/backend/link-plan-core.test.ts @@ -0,0 +1,20 @@ +import { expect, test } from "vitest"; +import { systemLibraryLinkerArgs } from "./link-plan-core.js"; +import { MACOS_ARM64_TARGET, LINUX_ARM64_GNU_TARGET, LINUX_X64_GNU_TARGET } from "./targets.js"; + +test("Linux links libm as needed and other system libraries as given", () => { + for (const target of [LINUX_X64_GNU_TARGET, LINUX_ARM64_GNU_TARGET]) { + expect(systemLibraryLinkerArgs(target, ["m"])).toEqual([ + "-Wl,--as-needed", + "-lm", + "-Wl,--no-as-needed", + ]); + expect(systemLibraryLinkerArgs(target, ["z", "m"])).toEqual([ + "-lz", + "-Wl,--as-needed", + "-lm", + "-Wl,--no-as-needed", + ]); + } + expect(systemLibraryLinkerArgs(MACOS_ARM64_TARGET, ["System"])).toEqual(["-lSystem"]); +}); diff --git a/packages/compiler/src/backend/link-plan-core.ts b/packages/compiler/src/backend/link-plan-core.ts index de5da1657..fc05d2961 100644 --- a/packages/compiler/src/backend/link-plan-core.ts +++ b/packages/compiler/src/backend/link-plan-core.ts @@ -70,3 +70,19 @@ export function executableLinkInputs(options: { ], }; } + +/** The linker arguments for an executable's system libraries. On Linux, libm + * links as needed: the runtime defines the libm functions most programs + * reach (scr_glibc_math.c) and Math's transcendentals (scr_ieee754.c), so a + * program that calls no other libm function, such as pow, starts without + * loading libm.so. */ +export function systemLibraryLinkerArgs( + target: Pick, + names: readonly string[], +): string[] { + return names.flatMap((name) => + name === "m" && target.platform === "linux" + ? ["-Wl,--as-needed", "-lm", "-Wl,--no-as-needed"] + : [`-l${name}`], + ); +} diff --git a/packages/compiler/src/backend/linker.ts b/packages/compiler/src/backend/linker.ts index a54956ed9..c4f04133e 100644 --- a/packages/compiler/src/backend/linker.ts +++ b/packages/compiler/src/backend/linker.ts @@ -22,6 +22,7 @@ import { type NativeArtifactDependency } from "./native/contracts.js"; import { RuntimePackError, stageRuntimePackArtifacts } from "./runtime-pack.js"; import type { NativeLinkPlan } from "./link-plan.js"; import { lldLinkerArgs } from "./lld.js"; +import { systemLibraryLinkerArgs } from "./link-plan-core.js"; import { createDarwinDebugSymbols, installDarwinDebugSymbols, @@ -248,7 +249,7 @@ export async function linkNativeExecutable( ...lld, ...plan.driverFlags, ...plan.inputs.map((input) => staged.replacements.get(input) ?? input), - ...plan.systemLibraries.map((name) => `-l${name}`), + ...systemLibraryLinkerArgs(plan.target, plan.systemLibraries), "-o", privateOut, ]; diff --git a/packages/compiler/src/backend/llvm/expr-primitives.ts b/packages/compiler/src/backend/llvm/expr-primitives.ts index 810f23a3a..3714d8d5e 100644 --- a/packages/compiler/src/backend/llvm/expr-primitives.ts +++ b/packages/compiler/src/backend/llvm/expr-primitives.ts @@ -47,7 +47,7 @@ import { type ArrayAccess, } from "./dense-array-access.js"; import { emitWrappingSum, integerOrNaNRange, markWrappingOperands } from "./wrapping-sums.js"; -import { emitPowerOfTwoQuotient, isPowerOfTwoCall } from "./lib-filesystem.js"; +import { emitPowF64, emitPowerOfTwoQuotient, isPowerOfTwoCall } from "./lib-filesystem.js"; /** The signed integer predicates of the number compares. */ const INTEGER_COMPARES: Record = { @@ -305,7 +305,6 @@ export function emitOperatorExpr( "===": "oeq", "!==": "une", }; - const libm: Record = { "**": "pow" }; const bit: Record = { "&": "and", "|": "or", @@ -471,32 +470,8 @@ export function emitOperatorExpr( B.line(`${t} = ${e.op === ">>>" ? "uitofp" : "sitofp"} i32 ${result} to double`); return { name: t, type: e.type, uint32: result }; } else { - const fn = libm[e.op]; - if (fn === undefined) throw new LlvmUnsupportedError(`bin:${e.op}`, e.loc); - host.declare(`declare double @${fn}(double, double)`); - // libm's pow answers 1 for pow(±1, ±Infinity) and pow(1, NaN), where - // `**` answers NaN (and LLVM folds pow(1.0, x) to 1.0): an inline - // guard selects NaN for a base of ±1 with a non-finite exponent and - // keeps the direct pow call. A finite literal exponent, or a literal - // base other than ±1, cannot meet those cases and needs no guard. - if (fn === "pow" && !libmPowAgrees(e.left, e.right)) { - host.declare(`declare double @llvm.fabs.f64(double)`); - const power = B.tmp(); - B.line(`${power} = call double @pow(double ${l.name}, double ${r.name})`); - const base = B.tmp(); - const unit = B.tmp(); - const magnitude = B.tmp(); - const finite = B.tmp(); - const reject = B.tmp(); - B.line(`${base} = call double @llvm.fabs.f64(double ${l.name})`); - B.line(`${unit} = fcmp oeq double ${base}, 1.0`); - B.line(`${magnitude} = call double @llvm.fabs.f64(double ${r.name})`); - B.line(`${finite} = fcmp olt double ${magnitude}, 0x7FF0000000000000`); - B.line(`${reject} = select i1 ${finite}, i1 false, i1 ${unit}`); - B.line(`${t} = select i1 ${reject}, double 0x7FF8000000000000, double ${power}`); - } else { - B.line(`${t} = call double @${fn}(double ${l.name}, double ${r.name})`); - } + if (e.op !== "**") throw new LlvmUnsupportedError(`bin:${e.op}`, e.loc); + return { name: emitPowF64(host, l.name, r.name, e.right), type: e.type }; } return { name: t, type: e.type }; } @@ -1473,19 +1448,3 @@ export function emitRecordLitFields( host.initField(ptr, type, inline, v); } } - -/** Whether libm pow matches `**` for these operands: every case except a - * base of ±1 with a NaN or infinite exponent agrees. */ -function libmPowAgrees(base: IrExpr, exponent: IrExpr): boolean { - const power = numericLiteral(exponent); - if (power !== null && Number.isFinite(power)) return true; - const value = numericLiteral(base); - return value !== null && !Number.isNaN(value) && Math.abs(value) !== 1; -} - -function numericLiteral(expr: IrExpr): number | null { - if (expr.kind === "numLit") return expr.value; - if (expr.kind === "unary" && expr.op === "-" && expr.operand.kind === "numLit") - return -expr.operand.value; - return null; -} diff --git a/packages/compiler/src/backend/llvm/lib-dispatch.ts b/packages/compiler/src/backend/llvm/lib-dispatch.ts index 41124cb75..b1cef8eb4 100644 --- a/packages/compiler/src/backend/llvm/lib-dispatch.ts +++ b/packages/compiler/src/backend/llvm/lib-dispatch.ts @@ -17,6 +17,7 @@ import type { LlvmEmitterContext, LibCallExpr, LibCallPrefix, LlValue } from "./ import { ASSERT_VERDICT_LIB_FNS, LIB_FN_SYMS, + PURE_MATH_LIB_FNS, mockTimerSym, USES_TIMERS_LIB_FNS, } from "./lib-shared.js"; @@ -777,7 +778,8 @@ function emitSymbolLibCall( return { name: "", type: e.type }; } const t = B.tmp(); - B.line(`${t} = call ${retDecl} @${sym}(${argList})`); + const attrs = PURE_MATH_LIB_FNS.has(e.fn) ? " nounwind willreturn memory(none)" : ""; + B.line(`${t} = call ${retDecl} @${sym}(${argList})${attrs}`); // The result joins its frame BEFORE the pending check so an unwind // releases the dummy (NULL for refcounted returns) harmlessly. const out = host.own({ name: t, type: e.type }); diff --git a/packages/compiler/src/backend/llvm/lib-filesystem.ts b/packages/compiler/src/backend/llvm/lib-filesystem.ts index f0b766936..9174f4775 100644 --- a/packages/compiler/src/backend/llvm/lib-filesystem.ts +++ b/packages/compiler/src/backend/llvm/lib-filesystem.ts @@ -1027,6 +1027,54 @@ function emitPowerOfTwo(host: LlvmEmitterContext, n: string): string { return out; } +/** `base ** exponent` and `Math.pow(base, exponent)` as V8's math::pow + * answers them: a literal exponent 2 is `base * base`, a literal 0.5 is + * `sqrt(base + 0)` (+Infinity for an infinite base), and every other + * exponent calls scr_math_pow, which applies V8's other special cases + * before the platform pow. The platform pow is never called by its libm + * name, which LLVM would rewrite (pow(2, x) to exp2(x), pow(x, -1) to 1/x) + * into answers that can differ from V8's. */ +export function emitPowF64( + host: LlvmEmitterContext, + base: string, + exponent: string, + exponentExpr: IrExpr, +): string { + const B = host.B; + const literal = + exponentExpr.kind === "numLit" + ? exponentExpr.value + : exponentExpr.kind === "unary" && + exponentExpr.op === "-" && + exponentExpr.operand.kind === "numLit" + ? -exponentExpr.operand.value + : null; + const out = B.tmp(); + if (literal === 2) { + B.line(`${out} = fmul double ${base}, ${base}`); + return out; + } + if (literal === 0.5) { + host.declare(`declare double @llvm.sqrt.f64(double)`); + host.declare(`declare double @llvm.fabs.f64(double)`); + const plus = B.tmp(), + root = B.tmp(), + magnitude = B.tmp(), + infinite = B.tmp(); + B.line(`${plus} = fadd double ${base}, ${f64Lit(0)}`); + B.line(`${root} = call double @llvm.sqrt.f64(double ${plus})`); + B.line(`${magnitude} = call double @llvm.fabs.f64(double ${base})`); + B.line(`${infinite} = fcmp oeq double ${magnitude}, 0x7FF0000000000000`); + B.line(`${out} = select i1 ${infinite}, double 0x7FF0000000000000, double ${root}`); + return out; + } + host.declare(`declare double @scr_math_pow(double, double)`); + B.line( + `${out} = call double @scr_math_pow(double ${base}, double ${exponent}) nounwind willreturn`, + ); + return out; +} + /** `Math.pow(2, n)`, the divisor emitPowerOfTwoQuotient takes. */ export function isPowerOfTwoCall(e: IrExpr): boolean { return ( @@ -1099,6 +1147,11 @@ export function emitPrimitiveLibCall(host: LlvmEmitterContext, e: LibCallExpr): const v = host.emitExpr(e.args[1]!); return { name: emitPowerOfTwo(host, v.name), type: e.type }; } + if (e.fn === "math.pow") { + const base = host.emitExpr(e.args[0]!); + const exponent = host.emitExpr(e.args[1]!); + return { name: emitPowF64(host, base.name, exponent.name, e.args[1]!), type: e.type }; + } if (e.fn === "math.floor" || e.fn === "math.ceil") { const intr = e.fn === "math.floor" ? "floor" : "ceil"; const v = host.emitExpr(e.args[0]!); diff --git a/packages/compiler/src/backend/llvm/lib-shared.ts b/packages/compiler/src/backend/llvm/lib-shared.ts index 4b3f34b80..f3167788f 100644 --- a/packages/compiler/src/backend/llvm/lib-shared.ts +++ b/packages/compiler/src/backend/llvm/lib-shared.ts @@ -35,6 +35,31 @@ export const ASSERT_VERDICT_LIB_FNS: ReadonlyMap = new Ma ["assert.partialDeepResult", null], ]); +/** The Math functions that call scr_ieee754.c. They read and write no + * memory, so calls with the same argument fold and loops may hoist them. */ +export const PURE_MATH_LIB_FNS: ReadonlySet = new Set([ + "math.sin", + "math.sinh", + "math.cos", + "math.cosh", + "math.tan", + "math.tanh", + "math.asin", + "math.asinh", + "math.acos", + "math.acosh", + "math.atan", + "math.atanh", + "math.cbrt", + "math.exp", + "math.expm1", + "math.log", + "math.log1p", + "math.log2", + "math.log10", + "math.atan2", +]); + export const LIB_FN_SYMS: Record = { "assert.expectsErrDyn": "scr_assert_expects_err_dyn", "assert.throwsValidate": "scr_assert_throws_validate", @@ -118,27 +143,29 @@ export const LIB_FN_SYMS: Record = { "str.codePoints": "scr_str_code_points", "num.fromString": "scr_string_to_number", "math.round": "scr_math_round", - "math.sin": "sin", - "math.sinh": "sinh", - "math.cos": "cos", - "math.cosh": "cosh", - "math.tan": "tan", - "math.tanh": "tanh", - "math.asin": "asin", - "math.asinh": "asinh", - "math.acos": "acos", - "math.acosh": "acosh", - "math.atan": "atan", - "math.atanh": "atanh", - "math.cbrt": "cbrt", - "math.exp": "exp", - "math.expm1": "expm1", + // Math's transcendental functions are V8's (scr_ieee754.c), so their + // results are the ones Node prints. + "math.sin": "scr_ieee754_sin", + "math.sinh": "scr_ieee754_sinh", + "math.cos": "scr_ieee754_cos", + "math.cosh": "scr_ieee754_cosh", + "math.tan": "scr_ieee754_tan", + "math.tanh": "scr_ieee754_tanh", + "math.asin": "scr_ieee754_asin", + "math.asinh": "scr_ieee754_asinh", + "math.acos": "scr_ieee754_acos", + "math.acosh": "scr_ieee754_acosh", + "math.atan": "scr_ieee754_atan", + "math.atanh": "scr_ieee754_atanh", + "math.cbrt": "scr_ieee754_cbrt", + "math.exp": "scr_ieee754_exp", + "math.expm1": "scr_ieee754_expm1", "math.sqrt": "sqrt", - "math.log": "log", - "math.log1p": "log1p", - "math.log2": "log2", - "math.log10": "log10", - "math.atan2": "atan2", + "math.log": "scr_ieee754_log", + "math.log1p": "scr_ieee754_log1p", + "math.log2": "scr_ieee754_log2", + "math.log10": "scr_ieee754_log10", + "math.atan2": "scr_ieee754_atan2", "math.pow": "scr_math_pow", // encodeURIComponent throws URIError on a lone surrogate // (MAY_THROW_LIB_FNS); decodeUriComponent is with the throwing tier below. diff --git a/packages/compiler/src/backend/native/runtime-inputs.ts b/packages/compiler/src/backend/native/runtime-inputs.ts index d32867331..891fa9312 100644 --- a/packages/compiler/src/backend/native/runtime-inputs.ts +++ b/packages/compiler/src/backend/native/runtime-inputs.ts @@ -74,6 +74,8 @@ export const EXECUTABLE_RUNTIME_SOURCES = [ "scr_slots.c", "scr_alloc.c", "scr_cipher_info.c", + "scr_ieee754.c", + "scr_glibc_math.c", ] as const; export function runtimeSrcDir(): string { diff --git a/packages/compiler/src/native/compiler.ts b/packages/compiler/src/native/compiler.ts index 0b56b0a53..f2734d2ee 100644 --- a/packages/compiler/src/native/compiler.ts +++ b/packages/compiler/src/native/compiler.ts @@ -21,7 +21,7 @@ import type { CompileRequestResult, } from "../compile-types.js"; import { LlvmUnsupportedError } from "../backend/llvm/emitter.js"; -import { executableLinkInputs } from "../backend/link-plan-core.js"; +import { executableLinkInputs, systemLibraryLinkerArgs } from "../backend/link-plan-core.js"; import { formatNativeLinkInfo } from "../backend/native-link-info-core.js"; import { emitNativeObject, @@ -527,7 +527,7 @@ export class NativeCompiler { ...lld, ...plan.driverFlags, ...plan.inputs, - ...plan.systemLibraries.map((name) => `-l${name}`), + ...systemLibraryLinkerArgs(target, plan.systemLibraries), "-o", stagedOutput, ]; diff --git a/packages/compiler/test/llvm-remainder.test.ts b/packages/compiler/test/llvm-remainder.test.ts index 9e6cb02e7..e3f585144 100644 --- a/packages/compiler/test/llvm-remainder.test.ts +++ b/packages/compiler/test/llvm-remainder.test.ts @@ -90,7 +90,7 @@ console.log(remainder(5, 2.5), signed(-7.5), maybeZero(7.5, 2.5)); expect(maybeZero).toContain("rem.float"); expect(ll).not.toContain("@llvm.fptosi.sat"); expect(ll).not.toContain("@fmod"); - expect(ll).toContain("call double @pow("); + expect(ll).toContain("call double @scr_math_pow("); // Where truncation is an integer round trip, a dividend below 2^62 takes // the i64 remainder by the power of two and adds back its fraction; a diff --git a/packages/compiler/test/ts7/baselines/order-parity/tests/corpus/7717-math-fdlibm.ts.json b/packages/compiler/test/ts7/baselines/order-parity/tests/corpus/7717-math-fdlibm.ts.json new file mode 100644 index 000000000..5c790ff61 --- /dev/null +++ b/packages/compiler/test/ts7/baselines/order-parity/tests/corpus/7717-math-fdlibm.ts.json @@ -0,0 +1,6 @@ +{ + "order": [ + "/tests/corpus/7717-math-fdlibm.ts" + ], + "diags": [] +} diff --git a/packages/compiler/test/ts7/baselines/order-parity/tests/corpus/7718-math-fdlibm.cjs.json b/packages/compiler/test/ts7/baselines/order-parity/tests/corpus/7718-math-fdlibm.cjs.json new file mode 100644 index 000000000..201049b54 --- /dev/null +++ b/packages/compiler/test/ts7/baselines/order-parity/tests/corpus/7718-math-fdlibm.cjs.json @@ -0,0 +1,6 @@ +{ + "order": [ + "/tests/corpus/7718-math-fdlibm.cjs" + ], + "diags": [] +} diff --git a/packages/runtime-darwin-arm64/scripts/build.mjs b/packages/runtime-darwin-arm64/scripts/build.mjs index 3305d69de..40f0d1c49 100644 --- a/packages/runtime-darwin-arm64/scripts/build.mjs +++ b/packages/runtime-darwin-arm64/scripts/build.mjs @@ -348,6 +348,11 @@ async function build() { "artifacts/licenses/v8-number-radix.txt", "BSD-3-Clause", ], + [ + join(vendorRoot, "v8-ieee754", "LICENSE"), + "artifacts/licenses/v8-ieee754.txt", + "BSD-3-Clause", + ], [join(vendorRoot, "unicode", "LICENSE"), "artifacts/licenses/unicode.txt", "Unicode-3.0"], [join(zlib, "LICENSE"), "artifacts/licenses/zlib.txt", "Zlib"], [join(vendorRoot, "brotli", "LICENSE"), "artifacts/licenses/brotli.txt", "MIT"], diff --git a/packages/runtime-pack-common/runtime-pack-matrix.mjs b/packages/runtime-pack-common/runtime-pack-matrix.mjs index 49e5896b6..18a2a43be 100644 --- a/packages/runtime-pack-common/runtime-pack-matrix.mjs +++ b/packages/runtime-pack-common/runtime-pack-matrix.mjs @@ -52,6 +52,10 @@ export const BASE_RUNTIME_SOURCES = [ "scr_alloc.c", // crypto.getCipherInfo's table; it links only into programs that call it. "scr_cipher_info.c", + // Math's transcendental functions as V8 computes them, and on glibc the + // exact libm functions that let a program link without libm.so. + "scr_ieee754.c", + "scr_glibc_math.c", ]; export const OPTIONAL_RUNTIME_SOURCES = [ diff --git a/packages/runtime-pack-common/scripts/build.mjs b/packages/runtime-pack-common/scripts/build.mjs index 8c0863b04..fd26a565b 100644 --- a/packages/runtime-pack-common/scripts/build.mjs +++ b/packages/runtime-pack-common/scripts/build.mjs @@ -447,6 +447,11 @@ async function build() { "artifacts/licenses/v8-number-radix.txt", "BSD-3-Clause", ], + [ + join(vendorRoot, "v8-ieee754", "LICENSE"), + "artifacts/licenses/v8-ieee754.txt", + "BSD-3-Clause", + ], [join(zlib, "LICENSE"), "artifacts/licenses/zlib.txt", "Zlib"], [join(vendorRoot, "brotli", "LICENSE"), "artifacts/licenses/brotli.txt", "MIT"], [join(mbedtls, "LICENSE"), "artifacts/licenses/mbedtls.txt", "Apache-2.0"], diff --git a/packages/runtime/src/scr_bytes.c b/packages/runtime/src/scr_bytes.c index 6399683fc..6bb56195a 100644 --- a/packages/runtime/src/scr_bytes.c +++ b/packages/runtime/src/scr_bytes.c @@ -4141,8 +4141,8 @@ static bool scr_bytes_rw_check(const ScrBytes *b, double offset, size_t width) { * exactly Node — and writes zeros downstream. */ static bool scr_bytes_check_int(const ScrBytes *b, double value, double offset, size_t width, bool sign) { - double max = sign ? exp2((double)(8 * width - 1)) - 1 : exp2((double)(8 * width)) - 1; - double min = sign ? -exp2((double)(8 * width - 1)) : 0; + double max = sign ? ldexp(1.0, (int)(8 * width - 1)) - 1 : ldexp(1.0, (int)(8 * width)) - 1; + double min = sign ? -ldexp(1.0, (int)(8 * width - 1)) : 0; if (value > max || value < min) { char recv[48]; scr_num_received(value, recv); diff --git a/packages/runtime/src/scr_bytes_dyn.c b/packages/runtime/src/scr_bytes_dyn.c index e390ddfce..5ce36bdac 100644 --- a/packages/runtime/src/scr_bytes_dyn.c +++ b/packages/runtime/src/scr_bytes_dyn.c @@ -44,8 +44,8 @@ static size_t scr_bytes_dyn_width(ScrBytesNumKind kind) { /* Node's checkInt range test and message, before the offset is looked at. */ static bool scr_bytes_dyn_value_in_range(double value, size_t width, bool sign) { - double max = sign ? exp2((double)(8 * width - 1)) - 1 : exp2((double)(8 * width)) - 1; - double min = sign ? -exp2((double)(8 * width - 1)) : 0; + double max = sign ? ldexp(1.0, (int)(8 * width - 1)) - 1 : ldexp(1.0, (int)(8 * width)) - 1; + double min = sign ? -ldexp(1.0, (int)(8 * width - 1)) : 0; if (!(value > max || value < min)) return true; char recv[48]; scr_num_received(value, recv); diff --git a/packages/runtime/src/scr_glibc_math.c b/packages/runtime/src/scr_glibc_math.c new file mode 100644 index 000000000..2ef2ed2e7 --- /dev/null +++ b/packages/runtime/src/scr_glibc_math.c @@ -0,0 +1,271 @@ +/* The libm functions that glibc executables reach without calling a Math + * transcendental: floor, ceil, trunc and round (LLVM's rounding intrinsics + * on x86-64 before SSE4.1, and the runtime's own calls), fmod (the % + * operator), sqrt (Math.sqrt's call for a negative argument), and ldexp, + * scalbn and frexp, which glibc exports from both libc and libm (the link + * would otherwise take them from libm). Each is exactly specified, so these + * definitions answer what glibc's do, the NaN payloads of NaN inputs + * included. Defining them in the executable lets a program that uses no + * other libm function link without libm.so (-lm is linked as needed), + * which saves loading it at startup. Library builds keep the embedder's + * libm. */ + +#include +#include + +#if defined(__GLIBC__) && !defined(SCR_LIB) + +/* The sanitizer runtime calls some of these functions while it initializes, + * and the dynamic loader calls the resolvers before that: none of the code + * here carries sanitizer instrumentation. no_sanitize("thread") alone keeps + * ThreadSanitizer's function entry and exit calls, which crash before its + * runtime initializes; disable_sanitizer_instrumentation removes them. */ +#pragma clang attribute push(__attribute__((no_sanitize("address", "undefined", "thread"), disable_sanitizer_instrumentation)), apply_to = function) + +static inline uint64_t sgm_bits(double d) { + uint64_t bits; + memcpy(&bits, &d, sizeof bits); + return bits; +} + +static inline double sgm_double(uint64_t bits) { + double d; + memcpy(&d, &bits, sizeof d); + return d; +} + +#define SGM_SIGN (UINT64_C(1) << 63) + +/* Hidden: the executable's own references bind to these definitions, which + * are not exported (libc.so and libm.so define the same names), so section + * GC drops the ones a program does not call. */ +#define SGM_DEFINE __attribute__((visibility("hidden"))) + +/* Values of magnitude below 2^52 round through an integer conversion: the + * truncated value, stepped by one when the rounding direction needs it, with + * the value's sign so that a zero result keeps it (ceil(-0.5) is -0). */ +static inline int sgm_small(double x) { return __builtin_fabs(x) < 0x1p52; } + +/* The unbiased exponent; 1024 for the infinities and NaN. Values with an + * exponent of 52 or more are integers already. */ +static inline int sgm_exponent(uint64_t bits) { return (int)((bits >> 52) & 0x7FF) - 1023; } + +static double sgm_floor(double x) { + if (sgm_small(x)) { + double t = (double)(int64_t)x; + return __builtin_copysign(t > x ? t - 1.0 : t, x); + } + uint64_t bits = sgm_bits(x); + int e = sgm_exponent(bits); + if (e >= 52) return e == 1024 ? x + x : x; + if (e < 0) { + if ((bits << 1) == 0) return x; + return bits & SGM_SIGN ? -1.0 : 0.0; + } + uint64_t fraction = (UINT64_C(1) << (52 - e)) - 1; + if ((bits & fraction) == 0) return x; + if (bits & SGM_SIGN) bits += fraction; + return sgm_double(bits & ~fraction); +} + +static double sgm_ceil(double x) { + if (sgm_small(x)) { + double t = (double)(int64_t)x; + return __builtin_copysign(t < x ? t + 1.0 : t, x); + } + uint64_t bits = sgm_bits(x); + int e = sgm_exponent(bits); + if (e >= 52) return e == 1024 ? x + x : x; + if (e < 0) { + if ((bits << 1) == 0) return x; + return bits & SGM_SIGN ? -0.0 : 1.0; + } + uint64_t fraction = (UINT64_C(1) << (52 - e)) - 1; + if ((bits & fraction) == 0) return x; + if (!(bits & SGM_SIGN)) bits += fraction; + return sgm_double(bits & ~fraction); +} + +static double sgm_trunc(double x) { + if (sgm_small(x)) return __builtin_copysign((double)(int64_t)x, x); + uint64_t bits = sgm_bits(x); + int e = sgm_exponent(bits); + if (e >= 52) return e == 1024 ? x + x : x; + if (e < 0) return sgm_double(bits & SGM_SIGN); + uint64_t fraction = (UINT64_C(1) << (52 - e)) - 1; + return sgm_double(bits & ~fraction); +} + +#if defined(__x86_64__) +#include +#include + +/* glibc selects SSE4.1's rounding instruction for floor, ceil and trunc when + * the CPU has it, and so do these: the same results, with no call into the + * integer path. */ +__attribute__((target("sse4.1"))) static double sgm_floor_sse41(double x) { + return _mm_cvtsd_f64(_mm_floor_sd(_mm_set_sd(x), _mm_set_sd(x))); +} +__attribute__((target("sse4.1"))) static double sgm_ceil_sse41(double x) { + return _mm_cvtsd_f64(_mm_ceil_sd(_mm_set_sd(x), _mm_set_sd(x))); +} +__attribute__((target("sse4.1"))) static double sgm_trunc_sse41(double x) { + return _mm_cvtsd_f64(_mm_round_sd(_mm_set_sd(x), _mm_set_sd(x), _MM_FROUND_TO_ZERO | _MM_FROUND_NO_EXC)); +} + +static int sgm_has_sse41(void) { + unsigned eax, ebx, ecx, edx; + return __get_cpuid(1, &eax, &ebx, &ecx, &edx) && (ecx & bit_SSE4_1) != 0; +} + +typedef double (*SgmUnary)(double); +static SgmUnary sgm_resolve_floor(void) { return sgm_has_sse41() ? sgm_floor_sse41 : sgm_floor; } +static SgmUnary sgm_resolve_ceil(void) { return sgm_has_sse41() ? sgm_ceil_sse41 : sgm_ceil; } +static SgmUnary sgm_resolve_trunc(void) { return sgm_has_sse41() ? sgm_trunc_sse41 : sgm_trunc; } +SGM_DEFINE double floor(double x) __attribute__((ifunc("sgm_resolve_floor"))); +SGM_DEFINE double ceil(double x) __attribute__((ifunc("sgm_resolve_ceil"))); +SGM_DEFINE double trunc(double x) __attribute__((ifunc("sgm_resolve_trunc"))); +#else +SGM_DEFINE double floor(double x) { return sgm_floor(x); } +SGM_DEFINE double ceil(double x) { return sgm_ceil(x); } +SGM_DEFINE double trunc(double x) { return sgm_trunc(x); } +#endif + +/* Math.sqrt is LLVM's sqrt libcall: the square root instruction inline, and + * a call to sqrt for a negative argument, where libm would set errno. This + * definition answers that call with the same instruction. */ +SGM_DEFINE double sqrt(double x) { +#if defined(__x86_64__) + return _mm_cvtsd_f64(_mm_sqrt_sd(_mm_set_sd(x), _mm_set_sd(x))); +#elif defined(__aarch64__) + double r; + __asm__("fsqrt %d0, %d1" : "=w"(r) : "w"(x)); + return r; +#else + return __builtin_sqrt(x); +#endif +} + +/* Halfway cases round away from zero. */ +SGM_DEFINE double round(double x) { + uint64_t bits = sgm_bits(x); + int e = sgm_exponent(bits); + if (e >= 52) return e == 1024 ? x + x : x; + if (e < 0) { + uint64_t sign = bits & SGM_SIGN; + return e == -1 ? sgm_double(sign | sgm_bits(1.0)) : sgm_double(sign); + } + uint64_t fraction = (UINT64_C(1) << (52 - e)) - 1; + if ((bits & fraction) == 0) return x; + bits += UINT64_C(1) << (51 - e); + return sgm_double(bits & ~fraction); +} + +/* The significand of a finite nonzero double as an integer below 2^53, and + * its exponent so that the value's magnitude is significand * 2^exponent. */ +static inline uint64_t sgm_significand(uint64_t bits, int *exponent) { + int biased = (int)((bits >> 52) & 0x7FF); + uint64_t significand = bits & ((UINT64_C(1) << 52) - 1); + if (biased == 0) { + *exponent = -1074; + return significand; + } + *exponent = biased - 1075; + return significand | (UINT64_C(1) << 52); +} + +/* x - n*y for the integer n = trunc(x/y), which is exact: the remainder + * comes from shift-and-subtract over the significands. NaN, an infinite x + * and a zero y answer (x*y)/(x*y), as glibc's fmod does. */ +SGM_DEFINE double fmod(double x, double y) { + uint64_t xb = sgm_bits(x), yb = sgm_bits(y); + uint64_t xmag = xb & ~SGM_SIGN, ymag = yb & ~SGM_SIGN; + uint64_t infinity = UINT64_C(0x7FF0000000000000); + if (ymag == 0 || xmag >= infinity || ymag > infinity) return (x * y) / (x * y); + if (xmag < ymag) return x; + if (xmag == ymag) return 0.0 * x; + int ex, ey; + uint64_t mx = sgm_significand(xb, &ex); + uint64_t my = sgm_significand(yb, &ey); + /* Normalize both to 53 significant bits, so the loop below shifts one + * exponent step at a time. */ + while (mx < (UINT64_C(1) << 52)) { + mx <<= 1; + ex--; + } + while (my < (UINT64_C(1) << 52)) { + my <<= 1; + ey--; + } + for (; ex > ey; ex--) { + if (mx >= my) mx -= my; + if (mx == 0) return 0.0 * x; + mx <<= 1; + } + if (mx >= my) mx -= my; + if (mx == 0) return 0.0 * x; + /* The remainder is mx * 2^ey, below y: rebuild it as a double. */ + while (mx < (UINT64_C(1) << 52)) { + mx <<= 1; + ey--; + } + uint64_t result; + if (ey + 1075 >= 1) { + result = ((uint64_t)(ey + 1075) << 52) | (mx & ((UINT64_C(1) << 52) - 1)); + } else { + result = mx >> (1 - (ey + 1075)); + } + return sgm_double(result | (xb & SGM_SIGN)); +} + +/* x * 2^n, rounded once. Exponents beyond one double's range scale in + * steps; the step toward the subnormals stays 53 bits above the smallest + * normal so that only the last multiplication rounds. */ +SGM_DEFINE double scalbn(double x, int n) { + if (n > 1023) { + x *= 0x1p1023; + n -= 1023; + if (n > 1023) { + x *= 0x1p1023; + n -= 1023; + if (n > 1023) n = 1023; + } + } else if (n < -1022) { + x *= 0x1p-1022 * 0x1p53; + n += 1022 - 53; + if (n < -1022) { + x *= 0x1p-1022 * 0x1p53; + n += 1022 - 53; + if (n < -1022) n = -1022; + } + } + return x * sgm_double((uint64_t)(0x3FF + n) << 52); +} + +SGM_DEFINE double ldexp(double x, int n) { return scalbn(x, n); } + +/* The significand in [0.5, 1) and the power of two: x = significand * 2^e. + * Zero, the infinities and NaN answer themselves with e = 0. */ +SGM_DEFINE double frexp(double x, int *e) { + uint64_t bits = sgm_bits(x); + int biased = (int)((bits >> 52) & 0x7FF); + if (biased == 0) { + if ((bits << 1) == 0) { + *e = 0; + return x; + } + x = frexp(x * 0x1p64, e); + *e -= 64; + return x; + } + if (biased == 0x7FF) { + *e = 0; + return x; + } + *e = biased - 1022; + return sgm_double((bits & ~(UINT64_C(0x7FF) << 52)) | (UINT64_C(0x3FE) << 52)); +} + +#pragma clang attribute pop + +#endif diff --git a/packages/runtime/src/scr_ieee754.c b/packages/runtime/src/scr_ieee754.c new file mode 100644 index 000000000..571e3f99a --- /dev/null +++ b/packages/runtime/src/scr_ieee754.c @@ -0,0 +1,2723 @@ +// The following is adapted from fdlibm (http://www.netlib.org/fdlibm). +// +// ==================================================== +// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. +// +// Developed at SunSoft, a Sun Microsystems, Inc. business. +// Permission to use, copy, modify, and distribute this +// software is freely granted, provided that this notice +// is preserved. +// ==================================================== +// +// The original source code covered by the above license above has been +// modified significantly by Google Inc. +// Copyright 2016 the V8 project authors. All rights reserved. + +/* Math's transcendental functions as V8 computes them: a C port of V8's + * src/base/ieee754.cc (Node v24.22.0's deps/v8), V8's adaptation of + * fdlibm. Node answers Math.sin, Math.exp, Math.cbrt and the other + * functions here through that file on every platform, so the compiled + * Math members call scr_ieee754_* and print the values Node prints. The + * functions are transcribed unchanged; the port renames the entry points, + * gives the kernels internal linkage, and replaces V8's helpers with C + * (bit casts, wraparound arithmetic, NaN and infinity constants). V8's + * BSD license is in ../vendor/v8-ieee754/LICENSE. + * + * Math.pow is not here: V8's math::pow calls the platform pow + * (scr_math_pow, scr_lib.c). */ + +#include +#include +#include + +#include "scr_ieee754.h" + +/* Each operation rounds on its own, as in V8's build: no fused + * multiply-add, which would change the results on targets that have one. */ +#pragma STDC FP_CONTRACT OFF + +static inline uint64_t scr_fd_bits(double d) { + uint64_t bits; + memcpy(&bits, &d, sizeof bits); + return bits; +} + +static inline double scr_fd_double(uint64_t bits) { + double d; + memcpy(&d, &bits, sizeof d); + return d; +} + +/* floor for __kernel_rem_pio2's nonnegative integer-part values, which are + * below 2^53: exact, with no libm call. */ +static inline double scr_fd_floor(double x) { + double t = (double)(int64_t)x; + return t > x ? t - 1.0 : t; +} + +/* + * The original fdlibm code used statements like: + * n0 = ((*(int*)&one)>>29)^1; * index of high word * + * ix0 = *(n0+(int*)&x); * high word of x * + * ix1 = *((1-n0)+(int*)&x); * low word of x * + * to dig two 32 bit words out of the 64 bit IEEE floating point + * value. That is non-ANSI, and, moreover, the gcc instruction + * scheduler gets it wrong. We instead use the following macros. + * Unlike the original code, we determine the endianness at compile + * time, not at run time; I don't see much benefit to selecting + * endianness at run time. + */ + +/* Get two 32 bit ints from a double. */ + +#define EXTRACT_WORDS(ix0, ix1, d) \ + do { \ + uint64_t bits = scr_fd_bits(d); \ + (ix0) = bits >> 32; \ + (ix1) = bits & 0xFFFFFFFFu; \ + } while (0) + +/* Get the more significant 32 bit int from a double. */ + +#define GET_HIGH_WORD(i, d) \ + do { \ + uint64_t bits = scr_fd_bits(d); \ + (i) = bits >> 32; \ + } while (0) + +/* Get the less significant 32 bit int from a double. */ + +#define GET_LOW_WORD(i, d) \ + do { \ + uint64_t bits = scr_fd_bits(d); \ + (i) = bits & 0xFFFFFFFFu; \ + } while (0) + +/* Set a double from two 32 bit ints. */ + +#define INSERT_WORDS(d, ix0, ix1) \ + do { \ + uint64_t bits = 0; \ + bits |= (uint64_t)(ix0) << 32; \ + bits |= (uint32_t)(ix1); \ + (d) = scr_fd_double(bits); \ + } while (0) + +/* Set the more significant 32 bits of a double from an int. */ + +#define SET_HIGH_WORD(d, v) \ + do { \ + uint64_t bits = scr_fd_bits(d); \ + bits &= 0x00000000FFFFFFFF; \ + bits |= (uint64_t)(v) << 32; \ + (d) = scr_fd_double(bits); \ + } while (0) + +/* Set the less significant 32 bits of a double from an int. */ + +#define SET_LOW_WORD(d, v) \ + do { \ + uint64_t bits = scr_fd_bits(d); \ + bits &= 0xFFFFFFFF00000000; \ + bits |= (uint32_t)(v); \ + (d) = scr_fd_double(bits); \ + } while (0) + +static int32_t __ieee754_rem_pio2(double x, double* y); +static int __kernel_rem_pio2(double* x, double* y, int e0, int nx, int prec, + const int32_t* ipio2); +static double __kernel_cos(double x, double y); +static double __kernel_sin(double x, double y, int iy); + +/* __ieee754_rem_pio2(x,y) + * + * return the remainder of x rem pi/2 in y[0]+y[1] + * use __kernel_rem_pio2() + */ +static int32_t __ieee754_rem_pio2(double x, double *y) { + /* + * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi + */ + static const int32_t two_over_pi[] = { + 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C, + 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 0x424DD2, 0xE00649, + 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 0xA73EE8, 0x8235F5, 0x2EBB44, + 0x84E99C, 0x7026B4, 0x5F7E41, 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, + 0xBDF928, 0x3B1FF8, 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, + 0x367ECF, 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, + 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 0x560330, + 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 0x91615E, 0xE61B08, + 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 0x4D7327, 0x310606, 0x1556CA, + 0x73A8C9, 0x60E27B, 0xC08C6B, + }; + + static const int32_t npio2_hw[] = { + 0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, + 0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C, + 0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A, + 0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C, + 0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB, + 0x404858EB, 0x404921FB, + }; + + /* + * invpio2: 53 bits of 2/pi + * pio2_1: first 33 bit of pi/2 + * pio2_1t: pi/2 - pio2_1 + * pio2_2: second 33 bit of pi/2 + * pio2_2t: pi/2 - (pio2_1+pio2_2) + * pio2_3: third 33 bit of pi/2 + * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) + */ + + static const double + zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ + half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ + two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ + invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ + pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ + pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ + pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */ + pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */ + pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ + pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ + + double z, w, t, r, fn; + double tx[3]; + int32_t e0, i, j, nx, n, ix, hx; + uint32_t low; + + z = 0; + GET_HIGH_WORD(hx, x); /* high word of x */ + ix = hx & 0x7FFFFFFF; + if (ix <= 0x3FE921FB) { /* |x| ~<= pi/4 , no need for reduction */ + y[0] = x; + y[1] = 0; + return 0; + } + if (ix < 0x4002D97C) { /* |x| < 3pi/4, special case with n=+-1 */ + if (hx > 0) { + z = x - pio2_1; + if (ix != 0x3FF921FB) { /* 33+53 bit pi is good enough */ + y[0] = z - pio2_1t; + y[1] = (z - y[0]) - pio2_1t; + } else { /* near pi/2, use 33+33+53 bit pi */ + z -= pio2_2; + y[0] = z - pio2_2t; + y[1] = (z - y[0]) - pio2_2t; + } + return 1; + } else { /* negative x */ + z = x + pio2_1; + if (ix != 0x3FF921FB) { /* 33+53 bit pi is good enough */ + y[0] = z + pio2_1t; + y[1] = (z - y[0]) + pio2_1t; + } else { /* near pi/2, use 33+33+53 bit pi */ + z += pio2_2; + y[0] = z + pio2_2t; + y[1] = (z - y[0]) + pio2_2t; + } + return -1; + } + } + if (ix <= 0x413921FB) { /* |x| ~<= 2^19*(pi/2), medium size */ + t = __builtin_fabs(x); + n = (int32_t)(t * invpio2 + half); + fn = (double)(n); + r = t - fn * pio2_1; + w = fn * pio2_1t; /* 1st round good to 85 bit */ + if (n < 32 && ix != npio2_hw[n - 1]) { + y[0] = r - w; /* quick check no cancellation */ + } else { + uint32_t high; + j = ix >> 20; + y[0] = r - w; + GET_HIGH_WORD(high, y[0]); + i = j - ((high >> 20) & 0x7FF); + if (i > 16) { /* 2nd iteration needed, good to 118 */ + t = r; + w = fn * pio2_2; + r = t - w; + w = fn * pio2_2t - ((t - r) - w); + y[0] = r - w; + GET_HIGH_WORD(high, y[0]); + i = j - ((high >> 20) & 0x7FF); + if (i > 49) { /* 3rd iteration need, 151 bits acc */ + t = r; /* will cover all possible cases */ + w = fn * pio2_3; + r = t - w; + w = fn * pio2_3t - ((t - r) - w); + y[0] = r - w; + } + } + } + y[1] = (r - y[0]) - w; + if (hx < 0) { + y[0] = -y[0]; + y[1] = -y[1]; + return -n; + } else { + return n; + } + } + /* + * all other (large) arguments + */ + if (ix >= 0x7FF00000) { /* x is inf or NaN */ + y[0] = y[1] = x - x; + return 0; + } + /* set z = scalbn(|x|,ilogb(x)-23) */ + GET_LOW_WORD(low, x); + SET_LOW_WORD(z, low); + e0 = (ix >> 20) - 1046; /* e0 = ilogb(z)-23; */ + SET_HIGH_WORD(z, ix - (int32_t)((uint32_t)(e0) << 20)); + for (i = 0; i < 2; i++) { + tx[i] = (double)((int32_t)(z)); + z = (z - tx[i]) * two24; + } + tx[2] = z; + nx = 3; + while (tx[nx - 1] == zero) nx--; /* skip zero term */ + n = __kernel_rem_pio2(tx, y, e0, nx, 2, two_over_pi); + if (hx < 0) { + y[0] = -y[0]; + y[1] = -y[1]; + return -n; + } + return n; +} + +/* __kernel_cos( x, y ) + * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164 + * Input x is assumed to be bounded by ~pi/4 in magnitude. + * Input y is the tail of x. + * + * Algorithm + * 1. Since cos(-x) = cos(x), we need only to consider positive x. + * 2. if x < 2^-27 (hx<0x3E400000 0), return 1 with inexact if x!=0. + * 3. cos(x) is approximated by a polynomial of degree 14 on + * [0,pi/4] + * 4 14 + * cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x + * where the remez error is + * + * | 2 4 6 8 10 12 14 | -58 + * |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2 + * | | + * + * 4 6 8 10 12 14 + * 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then + * cos(x) = 1 - x*x/2 + r + * since cos(x+y) ~ cos(x) - sin(x)*y + * ~ cos(x) - x*y, + * a correction term is necessary in cos(x) and hence + * cos(x+y) = 1 - (x*x/2 - (r - x*y)) + * For better accuracy when x > 0.3, let qx = |x|/4 with + * the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125. + * Then + * cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)). + * Note that 1-qx and (x*x/2-qx) is EXACT here, and the + * magnitude of the latter is at least a quarter of x*x/2, + * thus, reducing the rounding error in the subtraction. + */ +static inline double __kernel_cos(double x, double y) { + static const double + one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ + C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */ + C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */ + C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */ + C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */ + C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */ + C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */ + + double a, iz, z, r, qx; + int32_t ix; + GET_HIGH_WORD(ix, x); + ix &= 0x7FFFFFFF; /* ix = |x|'s high word*/ + if (ix < 0x3E400000) { /* if x < 2**27 */ + if ((int)(x) == 0) return one; /* generate inexact */ + } + z = x * x; + r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6))))); + if (ix < 0x3FD33333) { /* if |x| < 0.3 */ + return one - (0.5 * z - (z * r - x * y)); + } else { + if (ix > 0x3FE90000) { /* x > 0.78125 */ + qx = 0.28125; + } else { + INSERT_WORDS(qx, ix - 0x00200000, 0); /* x/4 */ + } + iz = 0.5 * z - qx; + a = one - qx; + return a - (iz - (z * r - x * y)); + } +} + +/* __kernel_rem_pio2(x,y,e0,nx,prec,ipio2) + * double x[],y[]; int e0,nx,prec; int ipio2[]; + * + * __kernel_rem_pio2 return the last three digits of N with + * y = x - N*pi/2 + * so that |y| < pi/2. + * + * The method is to compute the integer (mod 8) and fraction parts of + * (2/pi)*x without doing the full multiplication. In general we + * skip the part of the product that are known to be a huge integer ( + * more accurately, = 0 mod 8 ). Thus the number of operations are + * independent of the exponent of the input. + * + * (2/pi) is represented by an array of 24-bit integers in ipio2[]. + * + * Input parameters: + * x[] The input value (must be positive) is broken into nx + * pieces of 24-bit integers in double precision format. + * x[i] will be the i-th 24 bit of x. The scaled exponent + * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 + * match x's up to 24 bits. + * + * Example of breaking a double positive z into x[0]+x[1]+x[2]: + * e0 = ilogb(z)-23 + * z = scalbn(z,-e0) + * for i = 0,1,2 + * x[i] = floor(z) + * z = (z-x[i])*2**24 + * + * + * y[] output result in an array of double precision numbers. + * The dimension of y[] is: + * 24-bit precision 1 + * 53-bit precision 2 + * 64-bit precision 2 + * 113-bit precision 3 + * The actual value is the sum of them. Thus for 113-bit + * precison, one may have to do something like: + * + * long double t,w,r_head, r_tail; + * t = (long double)y[2] + (long double)y[1]; + * w = (long double)y[0]; + * r_head = t+w; + * r_tail = w - (r_head - t); + * + * e0 The exponent of x[0] + * + * nx dimension of x[] + * + * prec an integer indicating the precision: + * 0 24 bits (single) + * 1 53 bits (double) + * 2 64 bits (extended) + * 3 113 bits (quad) + * + * ipio2[] + * integer array, contains the (24*i)-th to (24*i+23)-th + * bit of 2/pi after binary point. The corresponding + * floating value is + * + * ipio2[i] * 2^(-24(i+1)). + * + * External function: + * double scalbn(), floor(); + * + * + * Here is the description of some local variables: + * + * jk jk+1 is the initial number of terms of ipio2[] needed + * in the computation. The recommended value is 2,3,4, + * 6 for single, double, extended,and quad. + * + * jz local integer variable indicating the number of + * terms of ipio2[] used. + * + * jx nx - 1 + * + * jv index for pointing to the suitable ipio2[] for the + * computation. In general, we want + * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8 + * is an integer. Thus + * e0-3-24*jv >= 0 or (e0-3)/24 >= jv + * Hence jv = max(0,(e0-3)/24). + * + * jp jp+1 is the number of terms in PIo2[] needed, jp = jk. + * + * q[] double array with integral value, representing the + * 24-bits chunk of the product of x and 2/pi. + * + * q0 the corresponding exponent of q[0]. Note that the + * exponent for q[i] would be q0-24*i. + * + * PIo2[] double precision array, obtained by cutting pi/2 + * into 24 bits chunks. + * + * f[] ipio2[] in floating point + * + * iq[] integer array by breaking up q[] in 24-bits chunk. + * + * fq[] final product of x*(2/pi) in fq[0],..,fq[jk] + * + * ih integer. If >0 it indicates q[] is >= 0.5, hence + * it also indicates the *sign* of the result. + * + */ +static int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, + const int32_t *ipio2) { + /* Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ + static const int init_jk[] = {2, 3, 4, 6}; /* initial value for jk */ + + static const double PIo2[] = { + 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */ + 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */ + 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */ + 3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */ + 1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */ + 1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */ + 2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */ + 2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */ + }; + + static const double + zero = 0.0, + one = 1.0, + two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ + twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */ + + int32_t jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih; + double z, fw, f[20], fq[20], q[20]; + + /* initialize jk*/ + jk = init_jk[prec]; + jp = jk; + + /* determine jx,jv,q0, note that 3>q0 */ + jx = nx - 1; + jv = (e0 - 3) / 24; + if (jv < 0) jv = 0; + q0 = e0 - 24 * (jv + 1); + + /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ + j = jv - jx; + m = jx + jk; + for (i = 0; i <= m; i++, j++) { + f[i] = (j < 0) ? zero : (double)(ipio2[j]); + } + + /* compute q[0],q[1],...q[jk] */ + for (i = 0; i <= jk; i++) { + for (j = 0, fw = 0.0; j <= jx; j++) fw += x[j] * f[jx + i - j]; + q[i] = fw; + } + + jz = jk; +recompute: + /* distill q[] into iq[] reversingly */ + for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) { + fw = (double)((int32_t)(twon24 * z)); + iq[i] = (int32_t)(z - two24 * fw); + z = q[j - 1] + fw; + } + + /* compute n */ + z = scalbn(z, q0); /* actual value of z */ + z -= 8.0 * scr_fd_floor(z * 0.125); /* trim off integer >= 8 */ + n = (int32_t)(z); + z -= (double)(n); + ih = 0; + if (q0 > 0) { /* need iq[jz-1] to determine n */ + i = (iq[jz - 1] >> (24 - q0)); + n += i; + iq[jz - 1] -= i << (24 - q0); + ih = iq[jz - 1] >> (23 - q0); + } else if (q0 == 0) { + ih = iq[jz - 1] >> 23; + } else if (z >= 0.5) { + ih = 2; + } + + if (ih > 0) { /* q > 0.5 */ + n += 1; + carry = 0; + for (i = 0; i < jz; i++) { /* compute 1-q */ + j = iq[i]; + if (carry == 0) { + if (j != 0) { + carry = 1; + iq[i] = 0x1000000 - j; + } + } else { + iq[i] = 0xFFFFFF - j; + } + } + if (q0 > 0) { /* rare case: chance is 1 in 12 */ + switch (q0) { + case 1: + iq[jz - 1] &= 0x7FFFFF; + break; + case 2: + iq[jz - 1] &= 0x3FFFFF; + break; + } + } + if (ih == 2) { + z = one - z; + if (carry != 0) z -= scalbn(one, q0); + } + } + + /* check if recomputation is needed */ + if (z == zero) { + j = 0; + for (i = jz - 1; i >= jk; i--) j |= iq[i]; + if (j == 0) { /* need recomputation */ + for (k = 1; jk >= k && iq[jk - k] == 0; k++) { + /* k = no. of terms needed */ + } + + for (i = jz + 1; i <= jz + k; i++) { /* add q[jz+1] to q[jz+k] */ + f[jx + i] = ipio2[jv + i]; + for (j = 0, fw = 0.0; j <= jx; j++) fw += x[j] * f[jx + i - j]; + q[i] = fw; + } + jz += k; + goto recompute; + } + } + + /* chop off zero terms */ + if (z == 0.0) { + jz -= 1; + q0 -= 24; + while (iq[jz] == 0) { + jz--; + q0 -= 24; + } + } else { /* break z into 24-bit if necessary */ + z = scalbn(z, -q0); + if (z >= two24) { + fw = (double)((int32_t)(twon24 * z)); + iq[jz] = z - two24 * fw; + jz += 1; + q0 += 24; + iq[jz] = fw; + } else { + iq[jz] = z; + } + } + + /* convert integer "bit" chunk to floating-point value */ + fw = scalbn(one, q0); + for (i = jz; i >= 0; i--) { + q[i] = fw * iq[i]; + fw *= twon24; + } + + /* compute PIo2[0,...,jp]*q[jz,...,0] */ + for (i = jz; i >= 0; i--) { + for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++) fw += PIo2[k] * q[i + k]; + fq[jz - i] = fw; + } + + /* compress fq[] into y[] */ + switch (prec) { + case 0: + fw = 0.0; + for (i = jz; i >= 0; i--) fw += fq[i]; + y[0] = (ih == 0) ? fw : -fw; + break; + case 1: + case 2: + fw = 0.0; + for (i = jz; i >= 0; i--) fw += fq[i]; + y[0] = (ih == 0) ? fw : -fw; + fw = fq[0] - fw; + for (i = 1; i <= jz; i++) fw += fq[i]; + y[1] = (ih == 0) ? fw : -fw; + break; + case 3: /* painful */ + for (i = jz; i > 0; i--) { + fw = fq[i - 1] + fq[i]; + fq[i] += fq[i - 1] - fw; + fq[i - 1] = fw; + } + for (i = jz; i > 1; i--) { + fw = fq[i - 1] + fq[i]; + fq[i] += fq[i - 1] - fw; + fq[i - 1] = fw; + } + for (fw = 0.0, i = jz; i >= 2; i--) fw += fq[i]; + if (ih == 0) { + y[0] = fq[0]; + y[1] = fq[1]; + y[2] = fw; + } else { + y[0] = -fq[0]; + y[1] = -fq[1]; + y[2] = -fw; + } + } + return n & 7; +} + +/* __kernel_sin( x, y, iy) + * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854 + * Input x is assumed to be bounded by ~pi/4 in magnitude. + * Input y is the tail of x. + * Input iy indicates whether y is 0. (if iy=0, y assume to be 0). + * + * Algorithm + * 1. Since sin(-x) = -sin(x), we need only to consider positive x. + * 2. if x < 2^-27 (hx<0x3E400000 0), return x with inexact if x!=0. + * 3. sin(x) is approximated by a polynomial of degree 13 on + * [0,pi/4] + * 3 13 + * sin(x) ~ x + S1*x + ... + S6*x + * where + * + * |sin(x) 2 4 6 8 10 12 | -58 + * |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2 + * | x | + * + * 4. sin(x+y) = sin(x) + sin'(x')*y + * ~ sin(x) + (1-x*x/2)*y + * For better accuracy, let + * 3 2 2 2 2 + * r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6)))) + * then 3 2 + * sin(x) = x + (S1*x + (x *(r-y/2)+y)) + */ +static inline double __kernel_sin(double x, double y, int iy) { + static const double + half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ + S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */ + S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */ + S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */ + S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */ + S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */ + S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */ + + double z, r, v; + int32_t ix; + GET_HIGH_WORD(ix, x); + ix &= 0x7FFFFFFF; /* high word of x */ + if (ix < 0x3E400000) { /* |x| < 2**-27 */ + if ((int)(x) == 0) return x; + } /* generate inexact */ + z = x * x; + v = z * x; + r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6))); + if (iy == 0) { + return x + v * (S1 + z * r); + } else { + return x - ((z * (half * y - v * r) - y) - v * S1); + } +} + +/* __kernel_tan( x, y, k ) + * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854 + * Input x is assumed to be bounded by ~pi/4 in magnitude. + * Input y is the tail of x. + * Input k indicates whether tan (if k=1) or + * -1/tan (if k= -1) is returned. + * + * Algorithm + * 1. Since tan(-x) = -tan(x), we need only to consider positive x. + * 2. if x < 2^-28 (hx<0x3E300000 0), return x with inexact if x!=0. + * 3. tan(x) is approximated by an odd polynomial of degree 27 on + * [0,0.67434] + * 3 27 + * tan(x) ~ x + T1*x + ... + T13*x + * where + * + * |tan(x) 2 4 26 | -59.2 + * |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2 + * | x | + * + * Note: tan(x+y) = tan(x) + tan'(x)*y + * ~ tan(x) + (1+x*x)*y + * Therefore, for better accuracy in computing tan(x+y), let + * 3 2 2 2 2 + * r = x *(T2+x *(T3+x *(...+x *(T12+x *T13)))) + * then + * 3 2 + * tan(x+y) = x + (T1*x + (x *(r+y)+y)) + * + * 4. For x in [0.67434,pi/4], let y = pi/4 - x, then + * tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y)) + * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y))) + */ +static double __kernel_tan(double x, double y, int iy) { + static const double xxx[] = { + 3.33333333333334091986e-01, /* 3FD55555, 55555563 */ + 1.33333333333201242699e-01, /* 3FC11111, 1110FE7A */ + 5.39682539762260521377e-02, /* 3FABA1BA, 1BB341FE */ + 2.18694882948595424599e-02, /* 3F9664F4, 8406D637 */ + 8.86323982359930005737e-03, /* 3F8226E3, E96E8493 */ + 3.59207910759131235356e-03, /* 3F6D6D22, C9560328 */ + 1.45620945432529025516e-03, /* 3F57DBC8, FEE08315 */ + 5.88041240820264096874e-04, /* 3F4344D8, F2F26501 */ + 2.46463134818469906812e-04, /* 3F3026F7, 1A8D1068 */ + 7.81794442939557092300e-05, /* 3F147E88, A03792A6 */ + 7.14072491382608190305e-05, /* 3F12B80F, 32F0A7E9 */ + -1.85586374855275456654e-05, /* BEF375CB, DB605373 */ + 2.59073051863633712884e-05, /* 3EFB2A70, 74BF7AD4 */ + /* one */ 1.00000000000000000000e+00, /* 3FF00000, 00000000 */ + /* pio4 */ 7.85398163397448278999e-01, /* 3FE921FB, 54442D18 */ + /* pio4lo */ 3.06161699786838301793e-17 /* 3C81A626, 33145C07 */ + }; +#define one xxx[13] +#define pio4 xxx[14] +#define pio4lo xxx[15] +#define T xxx + + double z, r, v, w, s; + int32_t ix, hx; + + GET_HIGH_WORD(hx, x); /* high word of x */ + ix = hx & 0x7FFFFFFF; /* high word of |x| */ + if (ix < 0x3E300000) { /* x < 2**-28 */ + if ((int)(x) == 0) { /* generate inexact */ + uint32_t low; + GET_LOW_WORD(low, x); + if (((ix | low) | (iy + 1)) == 0) { + return one / __builtin_fabs(x); + } else { + if (iy == 1) { + return x; + } else { /* compute -1 / (x+y) carefully */ + double a, t; + + z = w = x + y; + SET_LOW_WORD(z, 0); + v = y - (z - x); + t = a = -one / w; + SET_LOW_WORD(t, 0); + s = one + t * z; + return t + a * (s + t * v); + } + } + } + } + if (ix >= 0x3FE59428) { /* |x| >= 0.6744 */ + if (hx < 0) { + x = -x; + y = -y; + } + z = pio4 - x; + w = pio4lo - y; + x = z + w; + y = 0.0; + } + z = x * x; + w = z * z; + /* + * Break x^5*(T[1]+x^2*T[2]+...) into + * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) + + * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12])) + */ + r = T[1] + w * (T[3] + w * (T[5] + w * (T[7] + w * (T[9] + w * T[11])))); + v = z * + (T[2] + w * (T[4] + w * (T[6] + w * (T[8] + w * (T[10] + w * T[12]))))); + s = z * x; + r = y + z * (s * (r + v) + y); + r += T[0] * s; + w = x + r; + if (ix >= 0x3FE59428) { + v = iy; + return (1 - ((hx >> 30) & 2)) * (v - 2.0 * (x - (w * w / (w + v) - r))); + } + if (iy == 1) { + return w; + } else { + /* + * if allow error up to 2 ulp, simply return + * -1.0 / (x+r) here + */ + /* compute -1.0 / (x+r) accurately */ + double a, t; + z = w; + SET_LOW_WORD(z, 0); + v = r - (z - x); /* z+v = r+x */ + t = a = -1.0 / w; /* a = -1.0/w */ + SET_LOW_WORD(t, 0); + s = 1.0 + t * z; + return t + a * (s + t * v); + } + +#undef one +#undef pio4 +#undef pio4lo +#undef T +} + +/* acos(x) + * Method : + * acos(x) = pi/2 - asin(x) + * acos(-x) = pi/2 + asin(x) + * For |x|<=0.5 + * acos(x) = pi/2 - (x + x*x^2*R(x^2)) (see asin.c) + * For x>0.5 + * acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2))) + * = 2asin(sqrt((1-x)/2)) + * = 2s + 2s*z*R(z) ...z=(1-x)/2, s=sqrt(z) + * = 2f + (2c + 2s*z*R(z)) + * where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term + * for f so that f+c ~ sqrt(z). + * For x<-0.5 + * acos(x) = pi - 2asin(sqrt((1-|x|)/2)) + * = pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z) + * + * Special cases: + * if x is NaN, return x itself; + * if |x|>1, return NaN with invalid signal. + * + * Function needed: sqrt + */ +double scr_ieee754_acos(double x) { + static const double + one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ + pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ + pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ + pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */ + pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */ + pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */ + pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */ + pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */ + pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */ + pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */ + qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */ + qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */ + qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ + qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ + + double z, p, q, r, w, s, c, df; + int32_t hx, ix; + GET_HIGH_WORD(hx, x); + ix = hx & 0x7FFFFFFF; + if (ix >= 0x3FF00000) { /* |x| >= 1 */ + uint32_t lx; + GET_LOW_WORD(lx, x); + if (((ix - 0x3FF00000) | lx) == 0) { /* |x|==1 */ + if (hx > 0) + return 0.0; /* acos(1) = 0 */ + else + return pi + 2.0 * pio2_lo; /* acos(-1)= pi */ + } + return __builtin_nans(""); // acos(|x|>1) is NaN + } + if (ix < 0x3FE00000) { /* |x| < 0.5 */ + if (ix <= 0x3C600000) return pio2_hi + pio2_lo; /*if|x|<2**-57*/ + z = x * x; + p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5))))); + q = one + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4))); + r = p / q; + return pio2_hi - (x - (pio2_lo - x * r)); + } else if (hx < 0) { /* x < -0.5 */ + z = (one + x) * 0.5; + p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5))))); + q = one + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4))); + s = __builtin_sqrt(z); + r = p / q; + w = r * s - pio2_lo; + return pi - 2.0 * (s + w); + } else { /* x > 0.5 */ + z = (one - x) * 0.5; + s = __builtin_sqrt(z); + df = s; + SET_LOW_WORD(df, 0); + c = (z - df * df) / (s + df); + p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5))))); + q = one + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4))); + r = p / q; + w = r * s + c; + return 2.0 * (df + w); + } +} + +/* acosh(x) + * Method : + * Based on + * acosh(x) = log [ x + sqrt(x*x-1) ] + * we have + * acosh(x) := log(x)+ln2, if x is large; else + * acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else + * acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1. + * + * Special cases: + * acosh(x) is NaN with signal if x<1. + * acosh(NaN) is NaN without signal. + */ +double scr_ieee754_acosh(double x) { + static const double + one = 1.0, + ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */ + double t; + int32_t hx; + uint32_t lx; + EXTRACT_WORDS(hx, lx, x); + if (hx < 0x3FF00000) { /* x < 1 */ + return __builtin_nans(""); + } else if (hx >= 0x41B00000) { /* x > 2**28 */ + if (hx >= 0x7FF00000) { /* x is inf of NaN */ + return x + x; + } else { + return scr_ieee754_log(x) + ln2; /* acosh(huge)=log(2x) */ + } + } else if (((hx - 0x3FF00000) | lx) == 0) { + return 0.0; /* acosh(1) = 0 */ + } else if (hx > 0x40000000) { /* 2**28 > x > 2 */ + t = x * x; + return scr_ieee754_log(2.0 * x - one / (x + __builtin_sqrt(t - one))); + } else { /* 10.98 + * asin(x) = pi/2 - 2*(s+s*z*R(z)) + * = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo) + * For x<=0.98, let pio4_hi = pio2_hi/2, then + * f = hi part of s; + * c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z) + * and + * asin(x) = pi/2 - 2*(s+s*z*R(z)) + * = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo) + * = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c)) + * + * Special cases: + * if x is NaN, return x itself; + * if |x|>1, return NaN with invalid signal. + */ +double scr_ieee754_asin(double x) { + static const double + one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ + huge = 1.000e+300, + pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ + pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */ + pio4_hi = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */ + /* coefficient for R(x^2) */ + pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */ + pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */ + pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */ + pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */ + pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */ + pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */ + qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */ + qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */ + qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ + qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ + + double t, w, p, q, c, r, s; + int32_t hx, ix; + + t = 0; + GET_HIGH_WORD(hx, x); + ix = hx & 0x7FFFFFFF; + if (ix >= 0x3FF00000) { /* |x|>= 1 */ + uint32_t lx; + GET_LOW_WORD(lx, x); + if (((ix - 0x3FF00000) | lx) == 0) { /* asin(1)=+-pi/2 with inexact */ + return x * pio2_hi + x * pio2_lo; + } + return __builtin_nans(""); // asin(|x|>1) is NaN + } else if (ix < 0x3FE00000) { /* |x|<0.5 */ + if (ix < 0x3E400000) { /* if |x| < 2**-27 */ + if (huge + x > one) return x; /* return x with inexact if x!=0*/ + } else { + t = x * x; + } + p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5))))); + q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4))); + w = p / q; + return x + x * w; + } + /* 1> |x|>= 0.5 */ + w = one - __builtin_fabs(x); + t = w * 0.5; + p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5))))); + q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4))); + s = __builtin_sqrt(t); + if (ix >= 0x3FEF3333) { /* if |x| > 0.975 */ + w = p / q; + t = pio2_hi - (2.0 * (s + s * w) - pio2_lo); + } else { + w = s; + SET_LOW_WORD(w, 0); + c = (t - w * w) / (s + w); + r = p / q; + p = 2.0 * s * r - (pio2_lo - 2.0 * c); + q = pio4_hi - 2.0 * w; + t = pio4_hi - (p - q); + } + if (hx > 0) + return t; + else + return -t; +} +/* asinh(x) + * Method : + * Based on + * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ] + * we have + * asinh(x) := x if 1+x*x=1, + * := sign(x)*(log(x)+ln2)) for large |x|, else + * := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else + * := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2))) + */ +double scr_ieee754_asinh(double x) { + static const double + one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ + ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ + huge = 1.00000000000000000000e+300; + + double t, w; + int32_t hx, ix; + GET_HIGH_WORD(hx, x); + ix = hx & 0x7FFFFFFF; + if (ix >= 0x7FF00000) return x + x; /* x is inf or NaN */ + if (ix < 0x3E300000) { /* |x|<2**-28 */ + if (huge + x > one) return x; /* return x inexact except 0 */ + } + if (ix > 0x41B00000) { /* |x| > 2**28 */ + w = scr_ieee754_log(__builtin_fabs(x)) + ln2; + } else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */ + t = __builtin_fabs(x); + w = scr_ieee754_log(2.0 * t + one / (__builtin_sqrt(x * x + one) + t)); + } else { /* 2.0 > |x| > 2**-28 */ + t = x * x; + w = scr_ieee754_log1p(__builtin_fabs(x) + t / (one + __builtin_sqrt(one + t))); + } + if (hx > 0) { + return w; + } else { + return -w; + } +} + +/* atan(x) + * Method + * 1. Reduce x to positive by atan(x) = -atan(-x). + * 2. According to the integer k=4t+0.25 chopped, t=x, the argument + * is further reduced to one of the following intervals and the + * arctangent of t is evaluated by the corresponding formula: + * + * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) + * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) ) + * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) ) + * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) ) + * [39/16,INF] atan(x) = atan(INF) + atan( -1/t ) + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ +double scr_ieee754_atan(double x) { + static const double atanhi[] = { + 4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */ + 7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */ + 9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */ + 1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */ + }; + + static const double atanlo[] = { + 2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */ + 3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */ + 1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */ + 6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */ + }; + + static const double aT[] = { + 3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */ + -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */ + 1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */ + -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */ + 9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */ + -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */ + 6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */ + -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */ + 4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */ + -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */ + 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */ + }; + + static const double one = 1.0, huge = 1.0e300; + + double w, s1, s2, z; + int32_t ix, hx, id; + + GET_HIGH_WORD(hx, x); + ix = hx & 0x7FFFFFFF; + if (ix >= 0x44100000) { /* if |x| >= 2^66 */ + uint32_t low; + GET_LOW_WORD(low, x); + if (ix > 0x7FF00000 || (ix == 0x7FF00000 && (low != 0))) + return x + x; /* NaN */ + if (hx > 0) + return atanhi[3] + *(volatile double *)(&atanlo[3]); + else + return -atanhi[3] - *(volatile double *)(&atanlo[3]); + } + if (ix < 0x3FDC0000) { /* |x| < 0.4375 */ + if (ix < 0x3E400000) { /* |x| < 2^-27 */ + if (huge + x > one) return x; /* raise inexact */ + } + id = -1; + } else { + x = __builtin_fabs(x); + if (ix < 0x3FF30000) { /* |x| < 1.1875 */ + if (ix < 0x3FE60000) { /* 7/16 <=|x|<11/16 */ + id = 0; + x = (2.0 * x - one) / (2.0 + x); + } else { /* 11/16<=|x|< 19/16 */ + id = 1; + x = (x - one) / (x + one); + } + } else { + if (ix < 0x40038000) { /* |x| < 2.4375 */ + id = 2; + x = (x - 1.5) / (one + 1.5 * x); + } else { /* 2.4375 <= |x| < 2^66 */ + id = 3; + x = -1.0 / x; + } + } + } + /* end of argument reduction */ + z = x * x; + w = z * z; + /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */ + s1 = z * (aT[0] + + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * aT[10]))))); + s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * aT[9])))); + if (id < 0) { + return x - x * (s1 + s2); + } else { + z = atanhi[id] - ((x * (s1 + s2) - atanlo[id]) - x); + return (hx < 0) ? -z : z; + } +} + +/* atan2(y,x) + * Method : + * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). + * 2. Reduce x to positive by (if x and y are unexceptional): + * ARG (x+iy) = arctan(y/x) ... if x > 0, + * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, + * + * Special cases: + * + * ATAN2((anything), NaN ) is NaN; + * ATAN2(NAN , (anything) ) is NaN; + * ATAN2(+-0, +(anything but NaN)) is +-0 ; + * ATAN2(+-0, -(anything but NaN)) is +-pi ; + * ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2; + * ATAN2(+-(anything but INF and NaN), +INF) is +-0 ; + * ATAN2(+-(anything but INF and NaN), -INF) is +-pi; + * ATAN2(+-INF,+INF ) is +-pi/4 ; + * ATAN2(+-INF,-INF ) is +-3pi/4; + * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2; + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ +double scr_ieee754_atan2(double y, double x) { + static volatile double tiny = 1.0e-300; + static const double + zero = 0.0, + pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */ + pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */ + pi = 3.1415926535897931160E+00; /* 0x400921FB, 0x54442D18 */ + static volatile double pi_lo = + 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */ + + double z; + int32_t k, m, hx, hy, ix, iy; + uint32_t lx, ly; + + EXTRACT_WORDS(hx, lx, x); + ix = hx & 0x7FFFFFFF; + EXTRACT_WORDS(hy, ly, y); + iy = hy & 0x7FFFFFFF; + if (((ix | ((lx | ((int32_t)(0u - (uint32_t)(lx)))) >> 31)) > 0x7FF00000) || + ((iy | ((ly | ((int32_t)(0u - (uint32_t)(ly)))) >> 31)) > 0x7FF00000)) { + return x + y; /* x or y is NaN */ + } + if ((((int32_t)((uint32_t)(hx) - (uint32_t)(0x3FF00000))) | lx) == 0) { + return scr_ieee754_atan(y); /* x=1.0 */ + } + m = ((hy >> 31) & 1) | ((hx >> 30) & 2); /* 2*sign(x)+sign(y) */ + + /* when y = 0 */ + if ((iy | ly) == 0) { + switch (m) { + case 0: + case 1: + return y; /* atan(+-0,+anything)=+-0 */ + case 2: + return pi + tiny; /* atan(+0,-anything) = pi */ + case 3: + return -pi - tiny; /* atan(-0,-anything) =-pi */ + } + } + /* when x = 0 */ + if ((ix | lx) == 0) return (hy < 0) ? -pi_o_2 - tiny : pi_o_2 + tiny; + + /* when x is INF */ + if (ix == 0x7FF00000) { + if (iy == 0x7FF00000) { + switch (m) { + case 0: + return pi_o_4 + tiny; /* atan(+INF,+INF) */ + case 1: + return -pi_o_4 - tiny; /* atan(-INF,+INF) */ + case 2: + return 3.0 * pi_o_4 + tiny; /*atan(+INF,-INF)*/ + case 3: + return -3.0 * pi_o_4 - tiny; /*atan(-INF,-INF)*/ + } + } else { + switch (m) { + case 0: + return zero; /* atan(+...,+INF) */ + case 1: + return -zero; /* atan(-...,+INF) */ + case 2: + return pi + tiny; /* atan(+...,-INF) */ + case 3: + return -pi - tiny; /* atan(-...,-INF) */ + } + } + } + /* when y is INF */ + if (iy == 0x7FF00000) return (hy < 0) ? -pi_o_2 - tiny : pi_o_2 + tiny; + + /* compute y/x */ + k = (iy - ix) >> 20; + if (k > 60) { /* |y/x| > 2**60 */ + z = pi_o_2 + 0.5 * pi_lo; + m &= 1; + } else if (hx < 0 && k < -60) { + z = 0.0; /* 0 > |y|/x > -2**-60 */ + } else { + z = scr_ieee754_atan(__builtin_fabs(y / x)); /* safe to do y/x */ + } + switch (m) { + case 0: + return z; /* atan(+,+) */ + case 1: + return -z; /* atan(-,+) */ + case 2: + return pi - (z - pi_lo); /* atan(+,-) */ + default: /* case 3 */ + return (z - pi_lo) - pi; /* atan(-,-) */ + } +} + +/* cos(x) + * Return cosine function of x. + * + * kernel function: + * __kernel_sin ... sine function on [-pi/4,pi/4] + * __kernel_cos ... cosine function on [-pi/4,pi/4] + * __ieee754_rem_pio2 ... argument reduction routine + * + * Method. + * Let S,C and T denote the sin, cos and tan respectively on + * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 + * in [-pi/4 , +pi/4], and let n = k mod 4. + * We have + * + * n sin(x) cos(x) tan(x) + * ---------------------------------------------------------- + * 0 S C T + * 1 C -S -1/T + * 2 -S -C T + * 3 -C S -1/T + * ---------------------------------------------------------- + * + * Special cases: + * Let trig be any of sin, cos, or tan. + * trig(+-INF) is NaN, with signals; + * trig(NaN) is that NaN; + * + * Accuracy: + * TRIG(x) returns trig(x) nearly rounded + */ +double scr_ieee754_cos(double x) { + double y[2], z = 0.0; + int32_t n, ix; + + /* High word of x. */ + GET_HIGH_WORD(ix, x); + + /* |x| ~< pi/4 */ + ix &= 0x7FFFFFFF; + if (ix <= 0x3FE921FB) { + return __kernel_cos(x, z); + } else if (ix >= 0x7FF00000) { + /* cos(Inf or NaN) is NaN */ + return x - x; + } else { + /* argument reduction needed */ + n = __ieee754_rem_pio2(x, y); + switch (n & 3) { + case 0: + return __kernel_cos(y[0], y[1]); + case 1: + return -__kernel_sin(y[0], y[1], 1); + case 2: + return -__kernel_cos(y[0], y[1]); + default: + return __kernel_sin(y[0], y[1], 1); + } + } +} + +/* exp(x) + * Returns the exponential of x. + * + * Method + * 1. Argument reduction: + * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. + * Given x, find r and integer k such that + * + * x = k*ln2 + r, |r| <= 0.5*ln2. + * + * Here r will be represented as r = hi-lo for better + * accuracy. + * + * 2. Approximation of exp(r) by a special rational function on + * the interval [0,0.34658]: + * Write + * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... + * We use a special Remes algorithm on [0,0.34658] to generate + * a polynomial of degree 5 to approximate R. The maximum error + * of this polynomial approximation is bounded by 2**-59. In + * other words, + * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5 + * (where z=r*r, and the values of P1 to P5 are listed below) + * and + * | 5 | -59 + * | 2.0+P1*z+...+P5*z - R(z) | <= 2 + * | | + * The computation of exp(r) thus becomes + * 2*r + * exp(r) = 1 + ------- + * R - r + * r*R1(r) + * = 1 + r + ----------- (for better accuracy) + * 2 - R1(r) + * where + * 2 4 10 + * R1(r) = r - (P1*r + P2*r + ... + P5*r ). + * + * 3. Scale back to obtain exp(x): + * From step 1, we have + * exp(x) = 2^k * exp(r) + * + * Special cases: + * exp(INF) is INF, exp(NaN) is NaN; + * exp(-INF) is 0, and + * for finite argument, only exp(0)=1 is exact. + * + * Accuracy: + * according to an error analysis, the error is always less than + * 1 ulp (unit in the last place). + * + * Misc. info. + * For IEEE double + * if x > 7.09782712893383973096e+02 then exp(x) overflow + * if x < -7.45133219101941108420e+02 then exp(x) underflow + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ +double scr_ieee754_exp(double x) { + static const double + one = 1.0, + halF[2] = {0.5, -0.5}, + o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ + u_threshold = -7.45133219101941108420e+02, /* 0xC0874910, 0xD52D3051 */ + ln2HI[2] = {6.93147180369123816490e-01, /* 0x3FE62E42, 0xFEE00000 */ + -6.93147180369123816490e-01}, /* 0xBFE62E42, 0xFEE00000 */ + ln2LO[2] = {1.90821492927058770002e-10, /* 0x3DEA39EF, 0x35793C76 */ + -1.90821492927058770002e-10}, /* 0xBDEA39EF, 0x35793C76 */ + invln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE */ + P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ + P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ + P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ + P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ + P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */ + E = 2.718281828459045; /* 0x4005BF0A, 0x8B145769 */ + + static volatile double + huge = 1.0e+300, + twom1000 = 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0*/ + two1023 = 8.988465674311579539e307; /* 0x1p1023 */ + + double y, hi = 0.0, lo = 0.0, c, t, twopk; + int32_t k = 0, xsb; + uint32_t hx; + + GET_HIGH_WORD(hx, x); + xsb = (hx >> 31) & 1; /* sign bit of x */ + hx &= 0x7FFFFFFF; /* high word of |x| */ + + /* filter out non-finite argument */ + if (hx >= 0x40862E42) { /* if |x|>=709.78... */ + if (hx >= 0x7FF00000) { + uint32_t lx; + GET_LOW_WORD(lx, x); + if (((hx & 0xFFFFF) | lx) != 0) + return x + x; /* NaN */ + else + return (xsb == 0) ? x : 0.0; /* exp(+-inf)={inf,0} */ + } + if (x > o_threshold) return huge * huge; /* overflow */ + if (x < u_threshold) return twom1000 * twom1000; /* underflow */ + } + + /* argument reduction */ + if (hx > 0x3FD62E42) { /* if |x| > 0.5 ln2 */ + if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ + /* TODO(rtoy): We special case exp(1) here to return the correct + * value of E, as the computation below would get the last bit + * wrong. We should probably fix the algorithm instead. + */ + if (x == 1.0) return E; + hi = x - ln2HI[xsb]; + lo = ln2LO[xsb]; + k = 1 - xsb - xsb; + } else { + k = (int)(invln2 * x + halF[xsb]); + t = k; + hi = x - t * ln2HI[0]; /* t*ln2HI is exact here */ + lo = t * ln2LO[0]; + } + x = hi - lo; + } else if (hx < 0x3E300000) { /* when |x|<2**-28 */ + if (huge + x > one) return one + x; /* trigger inexact */ + } else { + k = 0; + } + + /* x is now in primary range */ + t = x * x; + if (k >= -1021) { + INSERT_WORDS( + twopk, + 0x3FF00000 + (int32_t)((uint32_t)(k) << 20), 0); + } else { + INSERT_WORDS(twopk, 0x3FF00000 + ((uint32_t)(k + 1000) << 20), + 0); + } + c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); + if (k == 0) { + return one - ((x * c) / (c - 2.0) - x); + } else { + y = one - ((lo - (x * c) / (2.0 - c)) - hi); + } + if (k >= -1021) { + if (k == 1024) return y * 2.0 * two1023; + return y * twopk; + } else { + return y * twopk * twom1000; + } +} + +/* + * Method : + * 1.Reduced x to positive by atanh(-x) = -atanh(x) + * 2.For x>=0.5 + * 1 2x x + * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) + * 2 1 - x 1 - x + * + * For x<0.5 + * atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) + * + * Special cases: + * atanh(x) is NaN if |x| > 1 with signal; + * atanh(NaN) is that NaN with no signal; + * atanh(+-1) is +-INF with signal. + * + */ +double scr_ieee754_atanh(double x) { + static const double one = 1.0, huge = 1e300; + static const double zero = 0.0; + + double t; + int32_t hx, ix; + uint32_t lx; + EXTRACT_WORDS(hx, lx, x); + ix = hx & 0x7FFFFFFF; + if ((ix | ((lx | ((int32_t)(0u - (uint32_t)(lx)))) >> 31)) > 0x3FF00000) { + /* |x|>1 */ + return __builtin_nans(""); + } + if (ix == 0x3FF00000) { + return x > 0 ? __builtin_inf() + : -__builtin_inf(); + } + if (ix < 0x3E300000 && (huge + x) > zero) return x; /* x<2**-28 */ + SET_HIGH_WORD(x, ix); + if (ix < 0x3FE00000) { /* x < 0.5 */ + t = x + x; + t = 0.5 * scr_ieee754_log1p(t + t * x / (one - x)); + } else { + t = 0.5 * scr_ieee754_log1p((x + x) / (one - x)); + } + if (hx >= 0) + return t; + else + return -t; +} + +/* log(x) + * Return the logrithm of x + * + * Method : + * 1. Argument Reduction: find k and f such that + * x = 2^k * (1+f), + * where sqrt(2)/2 < 1+f < sqrt(2) . + * + * 2. Approximation of log(1+f). + * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) + * = 2s + 2/3 s**3 + 2/5 s**5 + ....., + * = 2s + s*R + * We use a special Reme algorithm on [0,0.1716] to generate + * a polynomial of degree 14 to approximate R The maximum error + * of this polynomial approximation is bounded by 2**-58.45. In + * other words, + * 2 4 6 8 10 12 14 + * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s + * (the values of Lg1 to Lg7 are listed in the program) + * and + * | 2 14 | -58.45 + * | Lg1*s +...+Lg7*s - R(z) | <= 2 + * | | + * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. + * In order to guarantee error in log below 1ulp, we compute log + * by + * log(1+f) = f - s*(f - R) (if f is not too large) + * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) + * + * 3. Finally, log(x) = k*ln2 + log(1+f). + * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) + * Here ln2 is split into two floating point number: + * ln2_hi + ln2_lo, + * where n*ln2_hi is always exact for |n| < 2000. + * + * Special cases: + * log(x) is NaN with signal if x < 0 (including -INF) ; + * log(+INF) is +INF; log(0) is -INF with signal; + * log(NaN) is that NaN with no signal. + * + * Accuracy: + * according to an error analysis, the error is always less than + * 1 ulp (unit in the last place). + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ +double scr_ieee754_log(double x) { + static const double /* -- */ + ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ + ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ + two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ + Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ + Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ + Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ + Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ + Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ + Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ + Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ + + static const double zero = 0.0; + + double hfsq, f, s, z, R, w, t1, t2, dk; + int32_t k, hx, i, j; + uint32_t lx; + + EXTRACT_WORDS(hx, lx, x); + + k = 0; + if (hx < 0x00100000) { /* x < 2**-1022 */ + if (((hx & 0x7FFFFFFF) | lx) == 0) { + return -__builtin_inf(); /* log(+-0)=-inf */ + } + if (hx < 0) { + return __builtin_nans(""); /* log(-#) = NaN */ + } + k -= 54; + x *= two54; /* subnormal number, scale up x */ + GET_HIGH_WORD(hx, x); + } + if (hx >= 0x7FF00000) return x + x; + k += (hx >> 20) - 1023; + hx &= 0x000FFFFF; + i = (hx + 0x95F64) & 0x100000; + SET_HIGH_WORD(x, hx | (i ^ 0x3FF00000)); /* normalize x or x/2 */ + k += (i >> 20); + f = x - 1.0; + if ((0x000FFFFF & (2 + hx)) < 3) { /* -2**-20 <= f < 2**-20 */ + if (f == zero) { + if (k == 0) { + return zero; + } else { + dk = (double)(k); + return dk * ln2_hi + dk * ln2_lo; + } + } + R = f * f * (0.5 - 0.33333333333333333 * f); + if (k == 0) { + return f - R; + } else { + dk = (double)(k); + return dk * ln2_hi - ((R - dk * ln2_lo) - f); + } + } + s = f / (2.0 + f); + dk = (double)(k); + z = s * s; + i = hx - 0x6147A; + w = z * z; + j = 0x6B851 - hx; + t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); + t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); + i |= j; + R = t2 + t1; + if (i > 0) { + hfsq = 0.5 * f * f; + if (k == 0) + return f - (hfsq - s * (hfsq + R)); + else + return dk * ln2_hi - ((hfsq - (s * (hfsq + R) + dk * ln2_lo)) - f); + } else { + if (k == 0) + return f - s * (f - R); + else + return dk * ln2_hi - ((s * (f - R) - dk * ln2_lo) - f); + } +} + +/* double log1p(double x) + * + * Method : + * 1. Argument Reduction: find k and f such that + * 1+x = 2^k * (1+f), + * where sqrt(2)/2 < 1+f < sqrt(2) . + * + * Note. If k=0, then f=x is exact. However, if k!=0, then f + * may not be representable exactly. In that case, a correction + * term is need. Let u=1+x rounded. Let c = (1+x)-u, then + * log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u), + * and add back the correction term c/u. + * (Note: when x > 2**53, one can simply return log(x)) + * + * 2. Approximation of log1p(f). + * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) + * = 2s + 2/3 s**3 + 2/5 s**5 + ....., + * = 2s + s*R + * We use a special Reme algorithm on [0,0.1716] to generate + * a polynomial of degree 14 to approximate R The maximum error + * of this polynomial approximation is bounded by 2**-58.45. In + * other words, + * 2 4 6 8 10 12 14 + * R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s +Lp6*s +Lp7*s + * (the values of Lp1 to Lp7 are listed in the program) + * and + * | 2 14 | -58.45 + * | Lp1*s +...+Lp7*s - R(z) | <= 2 + * | | + * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. + * In order to guarantee error in log below 1ulp, we compute log + * by + * log1p(f) = f - (hfsq - s*(hfsq+R)). + * + * 3. Finally, log1p(x) = k*ln2 + log1p(f). + * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) + * Here ln2 is split into two floating point number: + * ln2_hi + ln2_lo, + * where n*ln2_hi is always exact for |n| < 2000. + * + * Special cases: + * log1p(x) is NaN with signal if x < -1 (including -INF) ; + * log1p(+INF) is +INF; log1p(-1) is -INF with signal; + * log1p(NaN) is that NaN with no signal. + * + * Accuracy: + * according to an error analysis, the error is always less than + * 1 ulp (unit in the last place). + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + * + * Note: Assuming log() return accurate answer, the following + * algorithm can be used to compute log1p(x) to within a few ULP: + * + * u = 1+x; + * if(u==1.0) return x ; else + * return log(u)*(x/(u-1.0)); + * + * See HP-15C Advanced Functions Handbook, p.193. + */ +double scr_ieee754_log1p(double x) { + static const double /* -- */ + ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ + ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ + two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ + Lp1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ + Lp2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ + Lp3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ + Lp4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ + Lp5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ + Lp6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ + Lp7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ + + static const double zero = 0.0; + + double hfsq, f, c, s, z, R, u; + int32_t k, hx, hu, ax; + + GET_HIGH_WORD(hx, x); + ax = hx & 0x7FFFFFFF; + + k = 1; + if (hx < 0x3FDA827A) { /* 1+x < sqrt(2)+ */ + if (ax >= 0x3FF00000) { /* x <= -1.0 */ + if (x == -1.0) + return -__builtin_inf(); /* log1p(-1)=+inf */ + else + return __builtin_nans(""); // log1p(x<-1)=NaN + } + if (ax < 0x3E200000) { /* |x| < 2**-29 */ + if (two54 + x > zero /* raise inexact */ + && ax < 0x3C900000) /* |x| < 2**-54 */ + return x; + else + return x - x * x * 0.5; + } + if (hx > 0 || hx <= (int32_t)(0xBFD2BEC4)) { + k = 0; + f = x; + hu = 1; + } /* sqrt(2)/2- <= 1+x < sqrt(2)+ */ + } + if (hx >= 0x7FF00000) return x + x; + if (k != 0) { + if (hx < 0x43400000) { + u = 1.0 + x; + GET_HIGH_WORD(hu, u); + k = (hu >> 20) - 1023; + c = (k > 0) ? 1.0 - (u - x) : x - (u - 1.0); /* correction term */ + c /= u; + } else { + u = x; + GET_HIGH_WORD(hu, u); + k = (hu >> 20) - 1023; + c = 0; + } + hu &= 0x000FFFFF; + /* + * The approximation to sqrt(2) used in thresholds is not + * critical. However, the ones used above must give less + * strict bounds than the one here so that the k==0 case is + * never reached from here, since here we have committed to + * using the correction term but don't use it if k==0. + */ + if (hu < 0x6A09E) { /* u ~< sqrt(2) */ + SET_HIGH_WORD(u, hu | 0x3FF00000); /* normalize u */ + } else { + k += 1; + SET_HIGH_WORD(u, hu | 0x3FE00000); /* normalize u/2 */ + hu = (0x00100000 - hu) >> 2; + } + f = u - 1.0; + } + hfsq = 0.5 * f * f; + if (hu == 0) { /* |f| < 2**-20 */ + if (f == zero) { + if (k == 0) { + return zero; + } else { + c += k * ln2_lo; + return k * ln2_hi + c; + } + } + R = hfsq * (1.0 - 0.66666666666666666 * f); + if (k == 0) + return f - R; + else + return k * ln2_hi - ((R - (k * ln2_lo + c)) - f); + } + s = f / (2.0 + f); + z = s * s; + R = z * (Lp1 + + z * (Lp2 + z * (Lp3 + z * (Lp4 + z * (Lp5 + z * (Lp6 + z * Lp7)))))); + if (k == 0) + return f - (hfsq - s * (hfsq + R)); + else + return k * ln2_hi - ((hfsq - (s * (hfsq + R) + (k * ln2_lo + c))) - f); +} + +/* + * k_log1p(f): + * Return log(1+f) - f for 1+f in ~[sqrt(2)/2, sqrt(2)]. + * + * The following describes the overall strategy for computing + * logarithms in base e. The argument reduction and adding the final + * term of the polynomial are done by the caller for increased accuracy + * when different bases are used. + * + * Method : + * 1. Argument Reduction: find k and f such that + * x = 2^k * (1+f), + * where sqrt(2)/2 < 1+f < sqrt(2) . + * + * 2. Approximation of log(1+f). + * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) + * = 2s + 2/3 s**3 + 2/5 s**5 + ....., + * = 2s + s*R + * We use a special Reme algorithm on [0,0.1716] to generate + * a polynomial of degree 14 to approximate R The maximum error + * of this polynomial approximation is bounded by 2**-58.45. In + * other words, + * 2 4 6 8 10 12 14 + * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s + * (the values of Lg1 to Lg7 are listed in the program) + * and + * | 2 14 | -58.45 + * | Lg1*s +...+Lg7*s - R(z) | <= 2 + * | | + * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. + * In order to guarantee error in log below 1ulp, we compute log + * by + * log(1+f) = f - s*(f - R) (if f is not too large) + * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) + * + * 3. Finally, log(x) = k*ln2 + log(1+f). + * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) + * Here ln2 is split into two floating point number: + * ln2_hi + ln2_lo, + * where n*ln2_hi is always exact for |n| < 2000. + * + * Special cases: + * log(x) is NaN with signal if x < 0 (including -INF) ; + * log(+INF) is +INF; log(0) is -INF with signal; + * log(NaN) is that NaN with no signal. + * + * Accuracy: + * according to an error analysis, the error is always less than + * 1 ulp (unit in the last place). + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ + +static const double Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ + Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ + Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ + Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ + Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ + Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ + Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ + +/* + * We always inline k_log1p(), since doing so produces a + * substantial performance improvement (~40% on amd64). + */ +static inline double k_log1p(double f) { + double hfsq, s, z, R, w, t1, t2; + + s = f / (2.0 + f); + z = s * s; + w = z * z; + t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); + t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); + R = t2 + t1; + hfsq = 0.5 * f * f; + return s * (hfsq + R); +} + +/* + * Return the base 2 logarithm of x. See e_log.c and k_log.h for most + * comments. + * + * This reduces x to {k, 1+f} exactly as in e_log.c, then calls the kernel, + * then does the combining and scaling steps + * log2(x) = (f - 0.5*f*f + k_log1p(f)) / ln2 + k + * in not-quite-routine extra precision. + */ +double scr_ieee754_log2(double x) { + static const double + two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ + ivln2hi = 1.44269504072144627571e+00, /* 0x3FF71547, 0x65200000 */ + ivln2lo = 1.67517131648865118353e-10; /* 0x3DE705FC, 0x2EEFA200 */ + + double f, hfsq, hi, lo, r, val_hi, val_lo, w, y; + int32_t i, k, hx; + uint32_t lx; + + EXTRACT_WORDS(hx, lx, x); + + k = 0; + if (hx < 0x00100000) { /* x < 2**-1022 */ + if (((hx & 0x7FFFFFFF) | lx) == 0) { + return -__builtin_inf(); /* log(+-0)=-inf */ + } + if (hx < 0) { + return __builtin_nans(""); /* log(-#) = NaN */ + } + k -= 54; + x *= two54; /* subnormal number, scale up x */ + GET_HIGH_WORD(hx, x); + } + if (hx >= 0x7FF00000) return x + x; + if (hx == 0x3FF00000 && lx == 0) return 0.0; /* log(1) = +0 */ + k += (hx >> 20) - 1023; + hx &= 0x000FFFFF; + i = (hx + 0x95F64) & 0x100000; + SET_HIGH_WORD(x, hx | (i ^ 0x3FF00000)); /* normalize x or x/2 */ + k += (i >> 20); + y = (double)(k); + f = x - 1.0; + hfsq = 0.5 * f * f; + r = k_log1p(f); + + /* + * f-hfsq must (for args near 1) be evaluated in extra precision + * to avoid a large cancellation when x is near sqrt(2) or 1/sqrt(2). + * This is fairly efficient since f-hfsq only depends on f, so can + * be evaluated in parallel with R. Not combining hfsq with R also + * keeps R small (though not as small as a true `lo' term would be), + * so that extra precision is not needed for terms involving R. + * + * Compiler bugs involving extra precision used to break Dekker's + * theorem for spitting f-hfsq as hi+lo, unless double_t was used + * or the multi-precision calculations were avoided when double_t + * has extra precision. These problems are now automatically + * avoided as a side effect of the optimization of combining the + * Dekker splitting step with the clear-low-bits step. + * + * y must (for args near sqrt(2) and 1/sqrt(2)) be added in extra + * precision to avoid a very large cancellation when x is very near + * these values. Unlike the above cancellations, this problem is + * specific to base 2. It is strange that adding +-1 is so much + * harder than adding +-ln2 or +-log10_2. + * + * This uses Dekker's theorem to normalize y+val_hi, so the + * compiler bugs are back in some configurations, sigh. And I + * don't want to used double_t to avoid them, since that gives a + * pessimization and the support for avoiding the pessimization + * is not yet available. + * + * The multi-precision calculations for the multiplications are + * routine. + */ + hi = f - hfsq; + SET_LOW_WORD(hi, 0); + lo = (f - hi) - hfsq + r; + val_hi = hi * ivln2hi; + val_lo = (lo + hi) * ivln2lo + lo * ivln2hi; + + /* spadd(val_hi, val_lo, y), except for not using double_t: */ + w = y + val_hi; + val_lo += (y - w) + val_hi; + val_hi = w; + + return val_lo + val_hi; +} + +/* + * Return the base 10 logarithm of x + * + * Method : + * Let log10_2hi = leading 40 bits of log10(2) and + * log10_2lo = log10(2) - log10_2hi, + * ivln10 = 1/log(10) rounded. + * Then + * n = ilogb(x), + * if(n<0) n = n+1; + * x = scalbn(x,-n); + * log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x)) + * + * Note 1: + * To guarantee log10(10**n)=n, where 10**n is normal, the rounding + * mode must set to Round-to-Nearest. + * Note 2: + * [1/log(10)] rounded to 53 bits has error .198 ulps; + * log10 is monotonic at all binary break points. + * + * Special cases: + * log10(x) is NaN if x < 0; + * log10(+INF) is +INF; log10(0) is -INF; + * log10(NaN) is that NaN; + * log10(10**N) = N for N=0,1,...,22. + */ +double scr_ieee754_log10(double x) { + static const double + two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ + ivln10 = 4.34294481903251816668e-01, + log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */ + log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */ + + double y; + int32_t i, k, hx; + uint32_t lx; + + EXTRACT_WORDS(hx, lx, x); + + k = 0; + if (hx < 0x00100000) { /* x < 2**-1022 */ + if (((hx & 0x7FFFFFFF) | lx) == 0) { + return -__builtin_inf(); /* log(+-0)=-inf */ + } + if (hx < 0) { + return __builtin_nan(""); /* log(-#) = NaN */ + } + k -= 54; + x *= two54; /* subnormal number, scale up x */ + GET_HIGH_WORD(hx, x); + GET_LOW_WORD(lx, x); + } + if (hx >= 0x7FF00000) return x + x; + if (hx == 0x3FF00000 && lx == 0) return 0.0; /* log(1) = +0 */ + k += (hx >> 20) - 1023; + + i = (k & 0x80000000) >> 31; + hx = (hx & 0x000FFFFF) | ((0x3FF - i) << 20); + y = k + i; + SET_HIGH_WORD(x, hx); + SET_LOW_WORD(x, lx); + + double z = y * log10_2lo + ivln10 * scr_ieee754_log(x); + return z + y * log10_2hi; +} + +/* expm1(x) + * Returns exp(x)-1, the exponential of x minus 1. + * + * Method + * 1. Argument reduction: + * Given x, find r and integer k such that + * + * x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658 + * + * Here a correction term c will be computed to compensate + * the error in r when rounded to a floating-point number. + * + * 2. Approximating expm1(r) by a special rational function on + * the interval [0,0.34658]: + * Since + * r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ... + * we define R1(r*r) by + * r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r) + * That is, + * R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r) + * = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r)) + * = 1 - r^2/60 + r^4/2520 - r^6/100800 + ... + * We use a special Reme algorithm on [0,0.347] to generate + * a polynomial of degree 5 in r*r to approximate R1. The + * maximum error of this polynomial approximation is bounded + * by 2**-61. In other words, + * R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5 + * where Q1 = -1.6666666666666567384E-2, + * Q2 = 3.9682539681370365873E-4, + * Q3 = -9.9206344733435987357E-6, + * Q4 = 2.5051361420808517002E-7, + * Q5 = -6.2843505682382617102E-9; + * z = r*r, + * with error bounded by + * | 5 | -61 + * | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2 + * | | + * + * expm1(r) = exp(r)-1 is then computed by the following + * specific way which minimize the accumulation rounding error: + * 2 3 + * r r [ 3 - (R1 + R1*r/2) ] + * expm1(r) = r + --- + --- * [--------------------] + * 2 2 [ 6 - r*(3 - R1*r/2) ] + * + * To compensate the error in the argument reduction, we use + * expm1(r+c) = expm1(r) + c + expm1(r)*c + * ~ expm1(r) + c + r*c + * Thus c+r*c will be added in as the correction terms for + * expm1(r+c). Now rearrange the term to avoid optimization + * screw up: + * ( 2 2 ) + * ({ ( r [ R1 - (3 - R1*r/2) ] ) } r ) + * expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- ) + * ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 ) + * ( ) + * + * = r - E + * 3. Scale back to obtain expm1(x): + * From step 1, we have + * expm1(x) = either 2^k*[expm1(r)+1] - 1 + * = or 2^k*[expm1(r) + (1-2^-k)] + * 4. Implementation notes: + * (A). To save one multiplication, we scale the coefficient Qi + * to Qi*2^i, and replace z by (x^2)/2. + * (B). To achieve maximum accuracy, we compute expm1(x) by + * (i) if x < -56*ln2, return -1.0, (raise inexact if x!=inf) + * (ii) if k=0, return r-E + * (iii) if k=-1, return 0.5*(r-E)-0.5 + * (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E) + * else return 1.0+2.0*(r-E); + * (v) if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1) + * (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else + * (vii) return 2^k(1-((E+2^-k)-r)) + * + * Special cases: + * expm1(INF) is INF, expm1(NaN) is NaN; + * expm1(-INF) is -1, and + * for finite argument, only expm1(0)=0 is exact. + * + * Accuracy: + * according to an error analysis, the error is always less than + * 1 ulp (unit in the last place). + * + * Misc. info. + * For IEEE double + * if x > 7.09782712893383973096e+02 then expm1(x) overflow + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ +double scr_ieee754_expm1(double x) { + static const double + one = 1.0, + tiny = 1.0e-300, + o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ + ln2_hi = 6.93147180369123816490e-01, /* 0x3FE62E42, 0xFEE00000 */ + ln2_lo = 1.90821492927058770002e-10, /* 0x3DEA39EF, 0x35793C76 */ + invln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE */ + /* Scaled Q's: Qn_here = 2**n * Qn_above, for R(2*z) where z = hxs = + x*x/2: */ + Q1 = -3.33333333333331316428e-02, /* BFA11111 111110F4 */ + Q2 = 1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */ + Q3 = -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */ + Q4 = 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */ + Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */ + + static volatile double huge = 1.0e+300; + + double y, hi, lo, c, t, e, hxs, hfx, r1, twopk; + int32_t k, xsb; + uint32_t hx; + + GET_HIGH_WORD(hx, x); + xsb = hx & 0x80000000; /* sign bit of x */ + hx &= 0x7FFFFFFF; /* high word of |x| */ + + /* filter out huge and non-finite argument */ + if (hx >= 0x4043687A) { /* if |x|>=56*ln2 */ + if (hx >= 0x40862E42) { /* if |x|>=709.78... */ + if (hx >= 0x7FF00000) { + uint32_t low; + GET_LOW_WORD(low, x); + if (((hx & 0xFFFFF) | low) != 0) + return x + x; /* NaN */ + else + return (xsb == 0) ? x : -1.0; /* exp(+-inf)={inf,-1} */ + } + if (x > o_threshold) return huge * huge; /* overflow */ + } + if (xsb != 0) { /* x < -56*ln2, return -1.0 with inexact */ + if (x + tiny < 0.0) /* raise inexact */ + return tiny - one; /* return -1 */ + } + } + + /* argument reduction */ + if (hx > 0x3FD62E42) { /* if |x| > 0.5 ln2 */ + if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ + if (xsb == 0) { + hi = x - ln2_hi; + lo = ln2_lo; + k = 1; + } else { + hi = x + ln2_hi; + lo = -ln2_lo; + k = -1; + } + } else { + k = invln2 * x + ((xsb == 0) ? 0.5 : -0.5); + t = k; + hi = x - t * ln2_hi; /* t*ln2_hi is exact here */ + lo = t * ln2_lo; + } + x = hi - lo; + c = (hi - x) - lo; + } else if (hx < 0x3C900000) { /* when |x|<2**-54, return x */ + t = huge + x; /* return x with inexact flags when x!=0 */ + return x - (t - (huge + x)); + } else { + k = 0; + } + + /* x is now in primary range */ + hfx = 0.5 * x; + hxs = x * hfx; + r1 = one + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5)))); + t = 3.0 - r1 * hfx; + e = hxs * ((r1 - t) / (6.0 - x * t)); + if (k == 0) { + return x - (x * e - hxs); /* c is 0 */ + } else { + INSERT_WORDS( + twopk, + 0x3FF00000 + (int32_t)((uint32_t)(k) << 20), + 0); /* 2^k */ + e = (x * (e - c) - c); + e -= hxs; + if (k == -1) return 0.5 * (x - e) - 0.5; + if (k == 1) { + if (x < -0.25) + return -2.0 * (e - (x + 0.5)); + else + return one + 2.0 * (x - e); + } + if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */ + y = one - (e - x); + // TODO(mvstanton): is this replacement for the hex float + // sufficient? + // if (k == 1024) y = y*2.0*0x1p1023; + if (k == 1024) + y = y * 2.0 * 8.98846567431158e+307; + else + y = y * twopk; + return y - one; + } + t = one; + if (k < 20) { + SET_HIGH_WORD(t, 0x3FF00000 - (0x200000 >> k)); /* t=1-2^-k */ + y = t - (e - x); + y = y * twopk; + } else { + SET_HIGH_WORD(t, ((0x3FF - k) << 20)); /* 2^-k */ + y = x - (e + t); + y += one; + y = y * twopk; + } + } + return y; +} + +double scr_ieee754_cbrt(double x) { + static const uint32_t + B1 = 715094163, /* B1 = (1023-1023/3-0.03306235651)*2**20 */ + B2 = 696219795; /* B2 = (1023-1023/3-54/3-0.03306235651)*2**20 */ + + /* |1/cbrt(x) - p(x)| < 2**-23.5 (~[-7.93e-8, 7.929e-8]). */ + static const double P0 = 1.87595182427177009643, /* 0x3FFE03E6, 0x0F61E692 */ + P1 = -1.88497979543377169875, /* 0xBFFE28E0, 0x92F02420 */ + P2 = 1.621429720105354466140, /* 0x3FF9F160, 0x4A49D6C2 */ + P3 = -0.758397934778766047437, /* 0xBFE844CB, 0xBEE751D9 */ + P4 = 0.145996192886612446982; /* 0x3FC2B000, 0xD4E4EDD7 */ + + int32_t hx; + double r, s, t = 0.0, w; + uint32_t sign; + uint32_t high, low; + + EXTRACT_WORDS(hx, low, x); + sign = hx & 0x80000000; /* sign= sign(x) */ + hx ^= sign; + if (hx >= 0x7FF00000) return (x + x); /* cbrt(NaN,INF) is itself */ + + /* + * Rough cbrt to 5 bits: + * cbrt(2**e*(1+m) ~= 2**(e/3)*(1+(e%3+m)/3) + * where e is integral and >= 0, m is real and in [0, 1), and "/" and + * "%" are integer division and modulus with rounding towards minus + * infinity. The RHS is always >= the LHS and has a maximum relative + * error of about 1 in 16. Adding a bias of -0.03306235651 to the + * (e%3+m)/3 term reduces the error to about 1 in 32. With the IEEE + * floating point representation, for finite positive normal values, + * ordinary integer division of the value in bits magically gives + * almost exactly the RHS of the above provided we first subtract the + * exponent bias (1023 for doubles) and later add it back. We do the + * subtraction virtually to keep e >= 0 so that ordinary integer + * division rounds towards minus infinity; this is also efficient. + */ + if (hx < 0x00100000) { /* zero or subnormal? */ + if ((hx | low) == 0) return (x); /* cbrt(0) is itself */ + SET_HIGH_WORD(t, 0x43500000); /* set t= 2**54 */ + t *= x; + GET_HIGH_WORD(high, t); + INSERT_WORDS(t, sign | ((high & 0x7FFFFFFF) / 3 + B2), 0); + } else { + INSERT_WORDS(t, sign | (hx / 3 + B1), 0); + } + + /* + * New cbrt to 23 bits: + * cbrt(x) = t*cbrt(x/t**3) ~= t*P(t**3/x) + * where P(r) is a polynomial of degree 4 that approximates 1/cbrt(r) + * to within 2**-23.5 when |r - 1| < 1/10. The rough approximation + * has produced t such than |t/cbrt(x) - 1| ~< 1/32, and cubing this + * gives us bounds for r = t**3/x. + * + * Try to optimize for parallel evaluation as in k_tanf.c. + */ + r = (t * t) * (t / x); + t = t * ((P0 + r * (P1 + r * P2)) + ((r * r) * r) * (P3 + r * P4)); + + /* + * Round t away from zero to 23 bits (sloppily except for ensuring that + * the result is larger in magnitude than cbrt(x) but not much more than + * 2 23-bit ulps larger). With rounding towards zero, the error bound + * would be ~5/6 instead of ~4/6. With a maximum error of 2 23-bit ulps + * in the rounded t, the infinite-precision error in the Newton + * approximation barely affects third digit in the final error + * 0.667; the error in the rounded t can be up to about 3 23-bit ulps + * before the final error is larger than 0.667 ulps. + */ + uint64_t bits = scr_fd_bits(t); + bits = (bits + 0x80000000) & 0xFFFFFFFFC0000000ULL; + t = scr_fd_double(bits); + + /* one step Newton iteration to 53 bits with error < 0.667 ulps */ + s = t * t; /* t*t is exact */ + r = x / s; /* error <= 0.5 ulps; |r| < |t| */ + w = t + t; /* t+t is exact */ + r = (r - t) / (w + r); /* r-t is exact; w+r ~= 3*t */ + t = t + t * r; /* error <= 0.5 + 0.5/3 + epsilon */ + + return (t); +} + +/* sin(x) + * Return sine function of x. + * + * kernel function: + * __kernel_sin ... sine function on [-pi/4,pi/4] + * __kernel_cos ... cose function on [-pi/4,pi/4] + * __ieee754_rem_pio2 ... argument reduction routine + * + * Method. + * Let S,C and T denote the sin, cos and tan respectively on + * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 + * in [-pi/4 , +pi/4], and let n = k mod 4. + * We have + * + * n sin(x) cos(x) tan(x) + * ---------------------------------------------------------- + * 0 S C T + * 1 C -S -1/T + * 2 -S -C T + * 3 -C S -1/T + * ---------------------------------------------------------- + * + * Special cases: + * Let trig be any of sin, cos, or tan. + * trig(+-INF) is NaN, with signals; + * trig(NaN) is that NaN; + * + * Accuracy: + * TRIG(x) returns trig(x) nearly rounded + */ +double scr_ieee754_sin(double x) { + double y[2], z = 0.0; + int32_t n, ix; + + /* High word of x. */ + GET_HIGH_WORD(ix, x); + + /* |x| ~< pi/4 */ + ix &= 0x7FFFFFFF; + if (ix <= 0x3FE921FB) { + return __kernel_sin(x, z, 0); + } else if (ix >= 0x7FF00000) { + /* sin(Inf or NaN) is NaN */ + return x - x; + } else { + /* argument reduction needed */ + n = __ieee754_rem_pio2(x, y); + switch (n & 3) { + case 0: + return __kernel_sin(y[0], y[1], 1); + case 1: + return __kernel_cos(y[0], y[1]); + case 2: + return -__kernel_sin(y[0], y[1], 1); + default: + return -__kernel_cos(y[0], y[1]); + } + } +} + +/* tan(x) + * Return tangent function of x. + * + * kernel function: + * __kernel_tan ... tangent function on [-pi/4,pi/4] + * __ieee754_rem_pio2 ... argument reduction routine + * + * Method. + * Let S,C and T denote the sin, cos and tan respectively on + * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 + * in [-pi/4 , +pi/4], and let n = k mod 4. + * We have + * + * n sin(x) cos(x) tan(x) + * ---------------------------------------------------------- + * 0 S C T + * 1 C -S -1/T + * 2 -S -C T + * 3 -C S -1/T + * ---------------------------------------------------------- + * + * Special cases: + * Let trig be any of sin, cos, or tan. + * trig(+-INF) is NaN, with signals; + * trig(NaN) is that NaN; + * + * Accuracy: + * TRIG(x) returns trig(x) nearly rounded + */ +double scr_ieee754_tan(double x) { + double y[2], z = 0.0; + int32_t n, ix; + + /* High word of x. */ + GET_HIGH_WORD(ix, x); + + /* |x| ~< pi/4 */ + ix &= 0x7FFFFFFF; + if (ix <= 0x3FE921FB) { + return __kernel_tan(x, z, 1); + } else if (ix >= 0x7FF00000) { + /* tan(Inf or NaN) is NaN */ + return x - x; /* NaN */ + } else { + /* argument reduction needed */ + n = __ieee754_rem_pio2(x, y); + /* 1 -> n even, -1 -> n odd */ + return __kernel_tan(y[0], y[1], 1 - ((n & 1) << 1)); + } +} + +/* + * ES6 draft 09-27-13, section 20.2.2.12. + * Math.cosh + * Method : + * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2 + * 1. Replace x by |x| (cosh(x) = cosh(-x)). + * 2. + * [ exp(x) - 1 ]^2 + * 0 <= x <= ln2/2 : cosh(x) := 1 + ------------------- + * 2*exp(x) + * + * exp(x) + 1/exp(x) + * ln2/2 <= x <= 22 : cosh(x) := ------------------- + * 2 + * 22 <= x <= lnovft : cosh(x) := exp(x)/2 + * lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2) + * ln2ovft < x : cosh(x) := huge*huge (overflow) + * + * Special cases: + * cosh(x) is |x| if x is +INF, -INF, or NaN. + * only cosh(0)=1 is exact for finite x. + */ +double scr_ieee754_cosh(double x) { + static const double KCOSH_OVERFLOW = 710.4758600739439; + static const double one = 1.0, half = 0.5; + static volatile double huge = 1.0e+300; + + int32_t ix; + + /* High word of |x|. */ + GET_HIGH_WORD(ix, x); + ix &= 0x7FFFFFFF; + + // |x| in [0,0.5*log2], return 1+expm1(|x|)^2/(2*exp(|x|)) + if (ix < 0x3FD62E43) { + double t = scr_ieee754_expm1(__builtin_fabs(x)); + double w = one + t; + // For |x| < 2^-55, cosh(x) = 1 + if (ix < 0x3C800000) return w; + return one + (t * t) / (w + w); + } + + // |x| in [0.5*log2, 22], return (exp(|x|)+1/exp(|x|)/2 + if (ix < 0x40360000) { + double t = scr_ieee754_exp(__builtin_fabs(x)); + return half * t + half / t; + } + + // |x| in [22, log(maxdouble)], return half*exp(|x|) + if (ix < 0x40862E42) return half * scr_ieee754_exp(__builtin_fabs(x)); + + // |x| in [log(maxdouble), overflowthreshold] + if (__builtin_fabs(x) <= KCOSH_OVERFLOW) { + double w = scr_ieee754_exp(half * __builtin_fabs(x)); + double t = half * w; + return t * w; + } + + /* x is INF or NaN */ + if (ix >= 0x7FF00000) return x * x; + + // |x| > overflowthreshold. + return huge * huge; +} + +/* + * ES6 draft 09-27-13, section 20.2.2.30. + * Math.sinh + * Method : + * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2 + * 1. Replace x by |x| (sinh(-x) = -sinh(x)). + * 2. + * E + E/(E+1) + * 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x) + * 2 + * + * 22 <= x <= lnovft : sinh(x) := exp(x)/2 + * lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2) + * ln2ovft < x : sinh(x) := x*shuge (overflow) + * + * Special cases: + * sinh(x) is |x| if x is +Infinity, -Infinity, or NaN. + * only sinh(0)=0 is exact for finite x. + */ +double scr_ieee754_sinh(double x) { + static const double KSINH_OVERFLOW = 710.4758600739439, + TWO_M28 = + 3.725290298461914e-9, // 2^-28, empty lower half + LOG_MAXD = 709.7822265625; // 0x40862E42 00000000, empty lower half + static const double shuge = 1.0e307; + + double h = (x < 0) ? -0.5 : 0.5; + // |x| in [0, 22]. return sign(x)*0.5*(E+E/(E+1)) + double ax = __builtin_fabs(x); + if (ax < 22) { + // For |x| < 2^-28, sinh(x) = x + if (ax < TWO_M28) return x; + double t = scr_ieee754_expm1(ax); + if (ax < 1) { + return h * (2 * t - t * t / (t + 1)); + } + return h * (t + t / (t + 1)); + } + // |x| in [22, log(maxdouble)], return 0.5 * exp(|x|) + if (ax < LOG_MAXD) return h * scr_ieee754_exp(ax); + // |x| in [log(maxdouble), overflowthreshold] + // overflowthreshold = 710.4758600739426 + if (ax <= KSINH_OVERFLOW) { + double w = scr_ieee754_exp(0.5 * ax); + double t = h * w; + return t * w; + } + // |x| > overflowthreshold or is NaN. + // Return Infinity of the appropriate sign or NaN. + return x * shuge; +} + +/* Tanh(x) + * Return the Hyperbolic Tangent of x + * + * Method : + * x -x + * e - e + * 0. tanh(x) is defined to be ----------- + * x -x + * e + e + * 1. reduce x to non-negative by tanh(-x) = -tanh(x). + * 2. 0 <= x < 2**-28 : tanh(x) := x with inexact if x != 0 + * -t + * 2**-28 <= x < 1 : tanh(x) := -----; t = expm1(-2x) + * t + 2 + * 2 + * 1 <= x < 22 : tanh(x) := 1 - -----; t = expm1(2x) + * t + 2 + * 22 <= x <= INF : tanh(x) := 1. + * + * Special cases: + * tanh(NaN) is NaN; + * only tanh(0)=0 is exact for finite argument. + */ +double scr_ieee754_tanh(double x) { + static const volatile double tiny = 1.0e-300; + static const double one = 1.0, two = 2.0, huge = 1.0e300; + double t, z; + int32_t jx, ix; + + GET_HIGH_WORD(jx, x); + ix = jx & 0x7FFFFFFF; + + /* x is INF or NaN */ + if (ix >= 0x7FF00000) { + if (jx >= 0) + return one / x + one; /* tanh(+-inf)=+-1 */ + else + return one / x - one; /* tanh(NaN) = NaN */ + } + + /* |x| < 22 */ + if (ix < 0x40360000) { /* |x|<22 */ + if (ix < 0x3E300000) { /* |x|<2**-28 */ + if (huge + x > one) return x; /* tanh(tiny) = tiny with inexact */ + } + if (ix >= 0x3FF00000) { /* |x|>=1 */ + t = scr_ieee754_expm1(two * __builtin_fabs(x)); + z = one - two / (t + two); + } else { + t = scr_ieee754_expm1(-two * __builtin_fabs(x)); + z = -t / (t + two); + } + /* |x| >= 22, return +-1 */ + } else { + z = one - tiny; /* raise inexact flag */ + } + return (jx >= 0) ? z : -z; +} + +#undef EXTRACT_WORDS +#undef GET_HIGH_WORD +#undef GET_LOW_WORD +#undef INSERT_WORDS +#undef SET_HIGH_WORD +#undef SET_LOW_WORD + diff --git a/packages/runtime/src/scr_ieee754.h b/packages/runtime/src/scr_ieee754.h new file mode 100644 index 000000000..7dddf9889 --- /dev/null +++ b/packages/runtime/src/scr_ieee754.h @@ -0,0 +1,26 @@ +#ifndef SCR_IEEE754_H +#define SCR_IEEE754_H + +/* Math's transcendental functions as V8 computes them (scr_ieee754.c). */ +double scr_ieee754_acos(double x); +double scr_ieee754_acosh(double x); +double scr_ieee754_asin(double x); +double scr_ieee754_asinh(double x); +double scr_ieee754_atan(double x); +double scr_ieee754_atan2(double y, double x); +double scr_ieee754_atanh(double x); +double scr_ieee754_cbrt(double x); +double scr_ieee754_cos(double x); +double scr_ieee754_cosh(double x); +double scr_ieee754_exp(double x); +double scr_ieee754_expm1(double x); +double scr_ieee754_log(double x); +double scr_ieee754_log10(double x); +double scr_ieee754_log1p(double x); +double scr_ieee754_log2(double x); +double scr_ieee754_sin(double x); +double scr_ieee754_sinh(double x); +double scr_ieee754_tan(double x); +double scr_ieee754_tanh(double x); + +#endif diff --git a/packages/runtime/src/scr_js_value.c b/packages/runtime/src/scr_js_value.c index 0313fb0cd..1e2021de1 100644 --- a/packages/runtime/src/scr_js_value.c +++ b/packages/runtime/src/scr_js_value.c @@ -209,7 +209,7 @@ static double scr_js_number_arithmetic(double a, double b, const ScrStr *operati case '/': return a / b; case '%': return fmod(a, b); default: - if (operation->len == 2) return (fabs(a) == 1 && !isfinite(b)) ? NAN : pow(a, b); + if (operation->len == 2) return scr_math_pow(a, b); return a * b; } } diff --git a/packages/runtime/src/scr_json.c b/packages/runtime/src/scr_json.c index fd201cb74..d4222c5e9 100644 --- a/packages/runtime/src/scr_json.c +++ b/packages/runtime/src/scr_json.c @@ -11160,9 +11160,8 @@ ScrDyn *scr_dyn_arithmetic(const ScrDyn *left, const ScrDyn *right, const ScrStr if (result) { out = scr_dyn_new_bigint(result); scr_bigint_release(result); } } else { double a = l->v.num, b = r->v.num; - // ECMAScript differs from libm for ±1 ** ±Infinity and 1 ** NaN. double result = op == '-' ? a - b : op == '/' ? a / b : op == '%' ? fmod(a, b) - : operation->len == 2 ? ((fabs(a) == 1 && !isfinite(b)) ? NAN : pow(a, b)) : a * b; + : operation->len == 2 ? scr_math_pow(a, b) : a * b; out = scr_dyn_new_num(result); } scr_dyn_release(l); diff --git a/packages/runtime/src/scr_lib.c b/packages/runtime/src/scr_lib.c index 1a625cdf1..4cfda127b 100644 --- a/packages/runtime/src/scr_lib.c +++ b/packages/runtime/src/scr_lib.c @@ -6828,11 +6828,17 @@ double scr_math_round(double x) { return (r == 0.0 && x < 0.0) ? -0.0 : r; } -/* IEEE pow(±1, ±Infinity) and pow(1, NaN) return 1; ECMAScript requires - * NaN. The other special cases, including signed zero and odd integer - * powers, follow libm. */ +/* V8's math::pow: the platform pow, as Node's Math.pow and ** call it, + * after the cases V8 answers itself. IEEE pow(±1, ±Infinity) and + * pow(1, NaN) return 1; ECMAScript requires NaN. A NaN base becomes the + * quiet NaN, and the exponents 2 and 0.5 answer base * base and + * sqrt(base + 0) (+Infinity for an infinite base), as V8's optimizing + * compilers do; pow can differ from those in the last bit. */ double scr_math_pow(double base, double exponent) { if (isnan(exponent) || (isinf(exponent) && fabs(base) == 1.0)) return NAN; + if (isnan(base)) base = NAN; + if (exponent == 2.0) return base * base; + if (exponent == 0.5) return isinf(base) ? INFINITY : sqrt(base + 0.0); /* 2 ** n for an integer n in the double range is an exact power of two, * the value pow returns; its bits are the exponent alone. */ if (base == 2.0 && exponent >= -1074.0 && exponent <= 1023.0) { @@ -12815,8 +12821,15 @@ ScrStr *scr_num_to_string_radix_f64(double value, double requested) { bool negative = value < 0; value = fabs(value); double integer = floor(value), fraction = value - integer; - double tolerance = 0.5 * (nextafter(value, INFINITY) - value); - if (tolerance <= 0) tolerance = nextafter(0, INFINITY); + /* The next double above the positive finite value: its bits plus one + * (Infinity above the largest double). */ + uint64_t value_bits; + memcpy(&value_bits, &value, sizeof value_bits); + value_bits += 1; + double next; + memcpy(&next, &value_bits, sizeof next); + double tolerance = 0.5 * (next - value); + if (tolerance <= 0) tolerance = 0x1p-1074; if (fraction >= tolerance) { buffer[end++] = '.'; do { diff --git a/packages/runtime/src/scr_runtime.h b/packages/runtime/src/scr_runtime.h index bf01a4a9b..f19f3918b 100644 --- a/packages/runtime/src/scr_runtime.h +++ b/packages/runtime/src/scr_runtime.h @@ -31,6 +31,7 @@ #include "scr_string_methods.h" #include "scr_qs.h" #include "scr_array.h" +#include "scr_ieee754.h" #include "scr_regex.h" #include "scr_map.h" #include "scr_closure.h" diff --git a/packages/runtime/vendor/README.md b/packages/runtime/vendor/README.md index 5f75b683b..53238410e 100644 --- a/packages/runtime/vendor/README.md +++ b/packages/runtime/vendor/README.md @@ -16,6 +16,13 @@ The x86_64 and AArch64 register-save/restore and `makecontext` logic is adapted The double-to-shortest-decimal core of Ryū (d2s and its headers only — the float and fixed/exponent variants are not vendored). scr_number.c textually `#include`s d2s.c and calls its internal `d2d` digit generator; ECMA-262 digit *placement* stays in scr_number.c. The only modification to upstream files: `#include "ryu/x.h"` flattened to `#include "x.h"` in d2s.c and d2s_intrinsics.h so the directory is self-contained. To update, re-fetch the listed files at the new commit and re-apply that include flattening. +## v8-ieee754/ (derived source) + +- Project: https://github.com/v8/v8, `src/base/ieee754.cc` as Node v24.22.0 bundles it (`deps/v8`, V8 13.6.233.17) +- License: BSD-3-Clause (V8), with fdlibm's notice (see v8-ieee754/LICENSE) + +V8's adaptation of fdlibm, the code behind Node's Math.sin, Math.exp, Math.cbrt and the other transcendental functions on every platform. `src/scr_ieee754.c` is a C port of the file, so compiled Math members print the values Node prints. The function bodies are unchanged; the port renames the entry points (`scr_ieee754_*`), gives the kernels internal linkage, and replaces V8's C++ helpers (bit casts, wraparound arithmetic, NaN and infinity constants) with C. `legacy::pow` is not ported: Node's Math.pow calls the platform pow. To update, re-fetch `deps/v8/src/base/ieee754.cc` at the new Node version and repeat the same conversion. + ## quickjs-ng/ - Project: https://github.com/quickjs-ng/quickjs diff --git a/packages/runtime/vendor/v8-ieee754/LICENSE b/packages/runtime/vendor/v8-ieee754/LICENSE new file mode 100644 index 000000000..9258b08c2 --- /dev/null +++ b/packages/runtime/vendor/v8-ieee754/LICENSE @@ -0,0 +1,75 @@ +This license applies to all parts of V8 that are not externally +maintained libraries. The externally maintained libraries used by V8 +are: + + - PCRE test suite, located in + test/mjsunit/third_party/regexp-pcre/regexp-pcre.js. This is based on the + test suite from PCRE-7.3, which is copyrighted by the University + of Cambridge and Google, Inc. The copyright notice and license + are embedded in regexp-pcre.js. + + - Layout tests, located in test/mjsunit/third_party/object-keys. These are + based on layout tests from webkit.org which are copyrighted by + Apple Computer, Inc. and released under a 3-clause BSD license. + + - Strongtalk assembler, the basis of the files assembler-arm-inl.h, + assembler-arm.cc, assembler-arm.h, assembler-ia32-inl.h, + assembler-ia32.cc, assembler-ia32.h, assembler-x64-inl.h, + assembler-x64.cc, assembler-x64.h, assembler.cc and assembler.h. + This code is copyrighted by Sun Microsystems Inc. and released + under a 3-clause BSD license. + + - Valgrind client API header, located at third_party/valgrind/valgrind.h + This is released under the BSD license. + + - The Wasm C/C++ API headers, located at third_party/wasm-api/wasm.{h,hh} + This is released under the Apache license. The API's upstream prototype + implementation also formed the basis of V8's implementation in + src/wasm/c-api.cc. + +These libraries have their own licenses; we recommend you read them, +as their terms may differ from the terms below. + +Further license information can be found in LICENSE files located in +sub-directories. + +Copyright 2014, the V8 project authors. 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 Google Inc. 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 +OWNER 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. + +--- + +src/scr_ieee754.c is a port of V8's src/base/ieee754.cc, which V8 adapted +from fdlibm. fdlibm's notice: + +==================================================== +Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + +Developed at SunSoft, a Sun Microsystems, Inc. business. +Permission to use, copy, modify, and distribute this +software is freely granted, provided that this notice +is preserved. +==================================================== diff --git a/tests/corpus/1110-math-members.ts b/tests/corpus/1110-math-members.ts index 80dd78209..5f07821c5 100644 --- a/tests/corpus/1110-math-members.ts +++ b/tests/corpus/1110-math-members.ts @@ -1,20 +1,17 @@ -// Math members execute in the island engine; results exit as ordinary -// static numbers. Members with exact IEEE results print directly. The -// transcendentals round through toFixed(9): the engine uses the system -// libm and Node uses V8's fdlibm port, which can differ in the last ulp -// (tan(1), asin(0.5), acos(0.5) do on macOS). +// Math members print the values Node prints, the transcendentals +// included: they are V8's fdlibm functions. console.log(Math.floor(7.6), Math.ceil(7.2), Math.round(7.5), Math.round(-7.5), Math.trunc(-7.9)); console.log(Math.abs(-3.5), Math.sign(-12), Math.sign(0), Math.sign(4.2)); console.log(Math.sqrt(144), Math.sqrt(2), Math.pow(2, 10), Math.pow(2, 0.5)); console.log(Math.hypot(3, 4), Math.min(2, -9), Math.max(2, -9)); console.log(Math.log2(1024), Math.log10(100000), Math.exp(0), Math.log(1)); console.log(Math.PI, Math.E); -console.log(Math.cbrt(27).toFixed(9), Math.exp(1).toFixed(9), Math.log(10).toFixed(9)); -console.log(Math.sin(1).toFixed(9), Math.cos(1).toFixed(9), Math.tan(1).toFixed(9)); -console.log(Math.asin(0.5).toFixed(9), Math.acos(0.5).toFixed(9)); -console.log(Math.atan(1).toFixed(9), Math.atan2(1, 1).toFixed(9), Math.atan2(-1, -1).toFixed(9)); -// Island results chain into static arithmetic like any number. +console.log(Math.cbrt(27), Math.exp(1), Math.log(10)); +console.log(Math.sin(1), Math.cos(1), Math.tan(1)); +console.log(Math.asin(0.5), Math.acos(0.5)); +console.log(Math.atan(1), Math.atan2(1, 1), Math.atan2(-1, -1)); +// Math results chain into static arithmetic like any number. const area = Math.PI * Math.pow(2, 2); console.log(area.toFixed(6), Math.floor(area) + 1); -// Non-finite edges: strict exits accept every engine number. +// Non-finite results. console.log(Math.log(0), Math.log(-1), Math.sqrt(-1), Math.pow(0, -1)); diff --git a/tests/corpus/7717-math-fdlibm.ts b/tests/corpus/7717-math-fdlibm.ts new file mode 100644 index 000000000..91ce02adf --- /dev/null +++ b/tests/corpus/7717-math-fdlibm.ts @@ -0,0 +1,97 @@ +// Math's transcendental functions print the values Node prints: each line +// digests the result bits of one function over a few thousand inputs, so +// a difference in the last bit of any result changes the line. Sample +// values print in full, and pow's special cases follow V8's. +const bits = new Float64Array(1); +const words = new Uint32Array(bits.buffer); + +function digest(values: number[]): string { + let h1 = 0x811c9dc5; + let h2 = 0x01000193; + for (const v of values) { + // NaN prints as NaN; its bits are not part of the result. + bits[0] = Number.isNaN(v) ? 0 : v; + h1 = Math.imul(h1 ^ words[0], 0x01000193) >>> 0; + h2 = Math.imul(h2 ^ words[1], 0x01000193) >>> 0; + } + return `${h1}:${h2}`; +} + +let seed = 12345; +function next(): number { + seed = (seed * 1103515245 + 12345) % 2147483648; + return seed / 2147483648; +} + +const wide: number[] = []; +const unit: number[] = []; +const small: number[] = []; +for (let i = 0; i < 3000; i++) { + wide.push((next() - 0.5) * Math.pow(10, Math.floor(next() * 12) - 4)); + unit.push(next() * 2 - 1); + small.push((next() - 0.5) * 40); +} +const positive = wide.map((x) => Math.abs(x)); +const edges = [0, -0, 1, -1, 0.5, -0.5, 2, Math.PI, Math.E, 1e-300, 5e-324, 1e300, Infinity, -Infinity, NaN]; + +const unary: [string, (x: number) => number, number[]][] = [ + ["sin", Math.sin, wide], + ["cos", Math.cos, wide], + ["tan", Math.tan, wide], + ["asin", Math.asin, unit], + ["acos", Math.acos, unit], + ["atan", Math.atan, wide], + ["sinh", Math.sinh, small], + ["cosh", Math.cosh, small], + ["tanh", Math.tanh, small], + ["asinh", Math.asinh, wide], + ["acosh", Math.acosh, positive.map((x) => x + 1)], + ["atanh", Math.atanh, unit], + ["exp", Math.exp, small], + ["expm1", Math.expm1, small], + ["log", Math.log, positive], + ["log1p", Math.log1p, positive], + ["log2", Math.log2, positive], + ["log10", Math.log10, positive], + ["cbrt", Math.cbrt, wide], + ["sqrt", Math.sqrt, positive], +]; +for (const [name, fn, inputs] of unary) { + console.log(name, digest(inputs.map(fn)), digest(edges.map(fn)), fn(0.3), fn(0.75)); +} + +// Math.sin(i) and Math.cos(i * 0.37) for i below 4096: main's libm calls +// printed 284 of these 8192 values differently from Node. +const sines: number[] = []; +const cosines: number[] = []; +for (let i = 0; i < 4096; i++) { + sines.push(Math.sin(i)); + cosines.push(Math.cos(i * 0.37)); +} +console.log("sin(i)", digest(sines), "cos(i*0.37)", digest(cosines), Math.cos(2.59)); + +const atan2: number[] = []; +const pow: number[] = []; +for (let i = 0; i < wide.length; i++) { + atan2.push(Math.atan2(wide[i], small[i])); + pow.push(Math.pow(positive[i] % 50, small[i] / 4)); +} +console.log("atan2", digest(atan2), Math.atan2(2.5, 3), Math.atan2(-0, -1)); +console.log("pow", digest(pow), Math.pow(1.1, 0.3), 1.0000001 ** 1e7); + +// V8 answers these itself: x * x, the square root, and NaN for ±1 ** ±Infinity. +const squares = wide.map((x) => x ** 2); +const roots = positive.map((x) => Math.pow(x, 0.5)); +console.log("pow 2", digest(squares), "pow 0.5", digest(roots)); +// The same exponents when they are only known at run time. +const two = Number(process.argv[2] ?? "2"); +const half = two / 4; +console.log("pow runtime 2", digest(wide.map((x) => x ** two)), "0.5", digest(positive.map((x) => Math.pow(x, half)))); +console.log(Math.pow(-Infinity, 0.5), Math.pow(-0, 0.5), (-8) ** 0.5, 1 ** Infinity, (-1) ** -Infinity, NaN ** 0); +// Rounding is exact: Math.floor, Math.ceil and Math.trunc over the same +// inputs and over halves, signed zeros included. +const halves = wide.map((x, i) => Math.round(x * 4) / 4 - (i % 2 === 0 ? 0 : 0.5)); +console.log("floor", digest(wide.map(Math.floor)), digest(halves.map(Math.floor)), digest(edges.map(Math.floor))); +console.log("ceil", digest(wide.map(Math.ceil)), digest(halves.map(Math.ceil)), digest(edges.map(Math.ceil))); +console.log("trunc", digest(wide.map(Math.trunc)), digest(halves.map(Math.trunc)), Object.is(Math.ceil(-0.5), -0), Object.is(Math.floor(-0), -0)); +console.log(Math.hypot(3, 4), Math.hypot(1e200, 1e200), 7 % 2.5, -7.5 % 2, Math.floor(-0.5), Math.ceil(-0.5), Math.round(2.5), Math.trunc(-4.7)); diff --git a/tests/corpus/7718-math-fdlibm.cjs b/tests/corpus/7718-math-fdlibm.cjs new file mode 100644 index 000000000..35d26f7e8 --- /dev/null +++ b/tests/corpus/7718-math-fdlibm.cjs @@ -0,0 +1,12 @@ +// Math functions and ** over untyped values print the values Node prints +// (V8's fdlibm results, and V8's pow special cases). +const values = JSON.parse("[0.3, 0.75, 2.5, -7.25, 1e-8, 123456.789, 1e21]"); +const rows = []; +for (const x of values) { + rows.push([Math.sin(x), Math.cos(x), Math.tan(x), Math.exp(x / 100), Math.cbrt(x), Math.atan2(x, 3)].join(" ")); + rows.push([Math.log(Math.abs(x)), Math.log10(Math.abs(x)), Math.sinh(x / 1e5), Math.cosh(x / 1e5), Math.expm1(x / 1e6)].join(" ")); + rows.push([x ** 2, x ** 0.5, Math.abs(x) ** 0.3, x % 1.5].join(" ")); +} +console.log(rows.join("\n")); +const edges = JSON.parse('[1, -1, null]'); +console.log(edges[0] ** Infinity, edges[1] ** -Infinity, Number(edges[2]) ** 0.5, (-Infinity) ** 0.5); diff --git a/tests/harness/corpus-costs/7717-math-fdlibm.ts.json b/tests/harness/corpus-costs/7717-math-fdlibm.ts.json new file mode 100644 index 000000000..63543b972 --- /dev/null +++ b/tests/harness/corpus-costs/7717-math-fdlibm.ts.json @@ -0,0 +1 @@ +{ "plain": 3357, "san": 4900 } diff --git a/tests/harness/corpus-costs/7718-math-fdlibm.cjs.json b/tests/harness/corpus-costs/7718-math-fdlibm.cjs.json new file mode 100644 index 000000000..281ed9b93 --- /dev/null +++ b/tests/harness/corpus-costs/7718-math-fdlibm.cjs.json @@ -0,0 +1 @@ +{ "plain": 3519, "san": 4613 }