From dd418b067a127c10ca44f89f697954ab2de4b826 Mon Sep 17 00:00:00 2001 From: Chris Tate Date: Sun, 25 Oct 2026 19:48:38 -0500 Subject: [PATCH 1/2] Compute Math functions as V8 does and link libm only when a program needs it --- .../src/backend/link-plan-core.test.ts | 20 + .../compiler/src/backend/link-plan-core.ts | 16 + packages/compiler/src/backend/linker.ts | 3 +- .../src/backend/llvm/expr-primitives.ts | 47 +- .../compiler/src/backend/llvm/lib-dispatch.ts | 4 +- .../src/backend/llvm/lib-filesystem.ts | 53 + .../compiler/src/backend/llvm/lib-shared.ts | 67 +- .../src/backend/native/runtime-inputs.ts | 2 + packages/compiler/src/native/compiler.ts | 4 +- packages/compiler/test/llvm-remainder.test.ts | 2 +- .../tests/corpus/7717-math-fdlibm.ts.json | 6 + .../tests/corpus/7718-math-fdlibm.cjs.json | 6 + .../runtime-darwin-arm64/scripts/build.mjs | 5 + .../runtime-pack-matrix.mjs | 4 + .../runtime-pack-common/scripts/build.mjs | 5 + packages/runtime/src/scr_bytes.c | 4 +- packages/runtime/src/scr_bytes_dyn.c | 4 +- packages/runtime/src/scr_glibc_math.c | 271 ++ packages/runtime/src/scr_ieee754.c | 2723 +++++++++++++++++ packages/runtime/src/scr_ieee754.h | 26 + packages/runtime/src/scr_js_value.c | 2 +- packages/runtime/src/scr_json.c | 3 +- packages/runtime/src/scr_lib.c | 23 +- packages/runtime/src/scr_runtime.h | 1 + packages/runtime/vendor/README.md | 7 + packages/runtime/vendor/v8-ieee754/LICENSE | 75 + tests/corpus/1110-math-members.ts | 19 +- tests/corpus/7717-math-fdlibm.ts | 97 + tests/corpus/7718-math-fdlibm.cjs | 12 + .../corpus-costs/7717-math-fdlibm.ts.json | 1 + .../corpus-costs/7718-math-fdlibm.cjs.json | 1 + 31 files changed, 3421 insertions(+), 92 deletions(-) create mode 100644 packages/compiler/src/backend/link-plan-core.test.ts create mode 100644 packages/compiler/test/ts7/baselines/order-parity/tests/corpus/7717-math-fdlibm.ts.json create mode 100644 packages/compiler/test/ts7/baselines/order-parity/tests/corpus/7718-math-fdlibm.cjs.json create mode 100644 packages/runtime/src/scr_glibc_math.c create mode 100644 packages/runtime/src/scr_ieee754.c create mode 100644 packages/runtime/src/scr_ieee754.h create mode 100644 packages/runtime/vendor/v8-ieee754/LICENSE create mode 100644 tests/corpus/7717-math-fdlibm.ts create mode 100644 tests/corpus/7718-math-fdlibm.cjs create mode 100644 tests/harness/corpus-costs/7717-math-fdlibm.ts.json create mode 100644 tests/harness/corpus-costs/7718-math-fdlibm.cjs.json 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 } From b6a1058f24940b59c367b7618bec3405109ccb45 Mon Sep 17 00:00:00 2001 From: Chris Tate Date: Sun, 25 Oct 2026 19:55:14 -0500 Subject: [PATCH 2/2] Fuse multiply-adds in the Math port on macOS arm64, as Node's build does Node for macOS on arm64 is compiled by clang, which contracts a multiply and an add within one expression into a fused multiply-add. Its sin, cos, tan, asin, acos, atan, sinh, asinh, acosh, atanh, exp, expm1, log, log1p and log10 differ from the unfused results in the last bit for some inputs, so corpus 7717 and 7718 failed on macOS. scr_ieee754.c now sets FP_CONTRACT ON for Apple arm64 targets. Other targets keep it off and compile to the same code as before. --- packages/runtime/src/scr_ieee754.c | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/packages/runtime/src/scr_ieee754.c b/packages/runtime/src/scr_ieee754.c index 571e3f99a..798bcbcbc 100644 --- a/packages/runtime/src/scr_ieee754.c +++ b/packages/runtime/src/scr_ieee754.c @@ -32,9 +32,16 @@ #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. */ +/* Floating-point contraction follows Node's build of this file. Node for + * macOS on arm64 is compiled by clang, which fuses a multiply and an add + * within one expression into a fused multiply-add; the results differ from + * the unfused ones in the last bit for some inputs. Other builds round each + * operation on its own, as Node's x86-64 builds do. */ +#if defined(__APPLE__) && defined(__aarch64__) +#pragma STDC FP_CONTRACT ON +#else #pragma STDC FP_CONTRACT OFF +#endif static inline uint64_t scr_fd_bits(double d) { uint64_t bits;