From 9c5f57fc3c282322440b761056dc9db92e19725e Mon Sep 17 00:00:00 2001 From: s-kybound Date: Sun, 14 Apr 2024 00:23:43 +0800 Subject: [PATCH] Extend base math functions --- src/stdlib/base.ts | 2 + src/stdlib/core-math.ts | 275 +++++++++++++++++++++++++++++++++++++++- 2 files changed, 274 insertions(+), 3 deletions(-) diff --git a/src/stdlib/base.ts b/src/stdlib/base.ts index 40a76c3..23f8e7c 100644 --- a/src/stdlib/base.ts +++ b/src/stdlib/base.ts @@ -408,6 +408,8 @@ export const lcm: Function = (...vals: core.SchemeInteger[]) => { return vals.reduce(atomic_lcm); }; +export const square: Function = (n: core.SchemeNumber) => $42$(n, n); + // pair operations export const cons: Function = core.pair; diff --git a/src/stdlib/core-math.ts b/src/stdlib/core-math.ts index b13f79c..6918d55 100644 --- a/src/stdlib/core-math.ts +++ b/src/stdlib/core-math.ts @@ -853,10 +853,10 @@ export class SchemeComplex { // current scm-slang will force any polar complex numbers to be made // inexact, hence we opt to limit the use of polar form as much as possible. class SchemePolar { - private readonly magnitude: SchemeReal; - private readonly angle: SchemeReal; + readonly magnitude: SchemeReal; + readonly angle: SchemeReal; - constructor(magnitude: SchemeReal, angle: SchemeReal) { + private constructor(magnitude: SchemeReal, angle: SchemeReal) { this.magnitude = magnitude; this.angle = angle; } @@ -1046,3 +1046,272 @@ export const LOG10E = SchemeReal.build(Math.LOG10E); export const SQRT1_2 = SchemeReal.build(Math.SQRT1_2); // other important functions +// for now, exponentials, square roots and the like will be treated as +// inexact functions, and will return inexact results. this allows us to +// leverage on the inbuilt javascript Math library. +// additional logic is required to handle complex numbers, which we can do with +// our polar form representation. + +export const expt = (n: SchemeNumber, e: SchemeNumber): SchemeNumber => { + if (!is_number(n) || !is_number(e)) { + throw new Error("expt: expected numbers"); + } + if (!is_real(n) || !is_real(e)) { + // complex number case + // we can convert both parts to polar form and use the + // polar form exponentiation formula. + + // given a * e^(bi) and c * e^(di), + // (a * e^(bi)) ^ (c * e^(di)) can be represented by + // the general formula for complex exponentiation: + // (a^c * e^(-bd)) * e^(i(bc * ln(a) + ad)) + + // convert both numbers to polar form + const nPolar = (n.promote(NumberType.COMPLEX) as SchemeComplex).toPolar(); + const ePolar = (e.promote(NumberType.COMPLEX) as SchemeComplex).toPolar(); + + const a = nPolar.magnitude.coerce(); + const b = nPolar.angle.coerce(); + const c = ePolar.magnitude.coerce(); + const d = ePolar.angle.coerce(); + + // we can construct a new polar form following the formula above + const mag = SchemeReal.build(a ** c * Math.E ** (-b * d)); + const angle = SchemeReal.build(b * c * Math.log(a) + a * d); + + return SchemePolar.build(mag, angle).toCartesian(); + } + // coerce both numbers to javascript numbers + const base = n.coerce(); + const exponent = e.coerce(); + + // there are probably cases here i am not considering yet. + // for now, we will just use the javascript Math library and hope for the best. + return SchemeReal.build(Math.pow(base, exponent)); +}; + +export const exp = (n: SchemeNumber): SchemeNumber => { + if (!is_number(n)) { + throw new Error("exp: expected number"); + } + if (!is_real(n)) { + // complex number case + throw new Error("exp: expected real number"); + } + return SchemeReal.build(Math.exp(n.coerce())); +}; + +export const log = (n: SchemeNumber, base: SchemeNumber = E): SchemeNumber => { + if (!is_number(n) || !is_number(base)) { + throw new Error("log: expected numbers"); + } + if (!is_real(n) || !is_real(base)) { + // complex number case + // we can convert both parts to polar form and use the + // polar form logarithm formula. + // where log(a * e^(bi)) = log(a) + bi + // and log(c * e^(di)) = log(c) + di + // and so result is log(a) + bi / log(c) + di + // which is just (log(a) - log(c)) + (b / d) i + + // convert both numbers to polar form + const nPolar = (n.promote(NumberType.COMPLEX) as SchemeComplex).toPolar(); + const basePolar = ( + base.promote(NumberType.COMPLEX) as SchemeComplex + ).toPolar(); + const a = nPolar.magnitude.coerce(); + const b = nPolar.angle.coerce(); + const c = basePolar.magnitude.coerce(); + const d = basePolar.angle.coerce(); + + return SchemeComplex.build( + SchemeReal.build(Math.log(a) - Math.log(c)), + SchemeReal.build(b / d), + ); + } + return SchemeReal.build(Math.log(n.coerce()) / Math.log(base.coerce())); +}; + +export const sqrt = (n: SchemeNumber): SchemeNumber => { + if (!is_number(n)) { + throw new Error("sqrt: expected number"); + } + if (!is_real(n)) { + // complex number case + const polar = (n.promote(NumberType.COMPLEX) as SchemeComplex).toPolar(); + const mag = polar.magnitude; + const angle = polar.angle; + + // the square root of a complex number is given by + // the square root of the magnitude and half the angle + const newMag = sqrt(mag) as SchemeReal; + const newAngle = SchemeReal.build(angle.coerce() / 2); + + return SchemePolar.build(newMag, newAngle).toCartesian(); + } + let value = n.coerce(); + + if (value < 0) { + return SchemeComplex.build( + SchemeReal.INEXACT_ZERO, + SchemeReal.build(Math.sqrt(-value)), + ); + } + + return SchemeReal.build(Math.sqrt(n.coerce())); +}; + +export const sin = (n: SchemeNumber): SchemeNumber => { + if (!is_number(n)) { + throw new Error("sin: expected number"); + } + if (!is_real(n)) { + // complex number case + + // we can use euler's formula to find sin(x) for a complex number x = a + bi + // e^(ix) = cos(x) + i * sin(x) + // that can be rearranged into + // sin(x) = (e^(ix) - e^(-ix)) / 2i + // and finally into + // sin(x) = (sin(a) * (e^(-b) + e^(b)) / 2) + i * (cos(a) * (e^(-b) - e^(b)) / 2) + const complex = n.promote(NumberType.COMPLEX) as SchemeComplex; + const real = complex.getReal(); + const imaginary = complex.getImaginary(); + const a = real.coerce(); + const b = imaginary.coerce(); + return SchemeComplex.build( + SchemeReal.build((Math.sin(a) * (Math.exp(-b) + Math.exp(b))) / 2), + SchemeReal.build((Math.cos(a) * (Math.exp(-b) - Math.exp(b))) / 2), + ); + } + return SchemeReal.build(Math.sin(n.coerce())); +}; + +export const cos = (n: SchemeNumber): SchemeNumber => { + if (!is_number(n)) { + throw new Error("cos: expected number"); + } + if (!is_real(n)) { + // complex number case + + // we can use euler's formula to find cos(x) for a complex number x = a + bi + // e^(ix) = cos(x) + i * sin(x) + // that can be rearranged into + // cos(x) = (e^(ix) + e^(-ix)) / 2 + // and finally into + // cos(x) = (cos(a) * (e^(-b) + e^(b)) / 2) - i * (sin(a) * (e^(-b) - e^(b)) / 2) + const complex = n.promote(NumberType.COMPLEX) as SchemeComplex; + const real = complex.getReal(); + const imaginary = complex.getImaginary(); + const a = real.coerce(); + const b = imaginary.coerce(); + return SchemeComplex.build( + SchemeReal.build((Math.cos(a) * (Math.exp(-b) + Math.exp(b))) / 2), + SchemeReal.build((-Math.sin(a) * (Math.exp(-b) - Math.exp(b))) / 2), + ); + } + return SchemeReal.build(Math.cos(n.coerce())); +}; + +export const tan = (n: SchemeNumber): SchemeNumber => { + if (!is_number(n)) { + throw new Error("tan: expected number"); + } + if (!is_real(n)) { + // complex number case + const sinValue = sin(n); + const cosValue = cos(n); + return atomic_divide(sinValue, cosValue); + } + return SchemeReal.build(Math.tan(n.coerce())); +}; + +export const asin = (n: SchemeNumber): SchemeNumber => { + if (!is_number(n)) { + throw new Error("asin: expected number"); + } + if (!is_real(n)) { + // complex number case + // asin(n) = -i * ln(i * n + sqrt(1 - n^2)) + // we already have the building blocks needed to compute this + const i = SchemeComplex.build( + SchemeInteger.EXACT_ZERO, + SchemeInteger.build(1), + ); + return atomic_multiply( + atomic_negate(i), + log( + atomic_add( + atomic_multiply(i, n), + sqrt(atomic_subtract(SchemeInteger.build(1), atomic_multiply(n, n))), + ), + ), + ); + } + return SchemeReal.build(Math.asin(n.coerce())); +}; + +export const acos = (n: SchemeNumber): SchemeNumber => { + if (!is_number(n)) { + throw new Error("acos: expected number"); + } + if (!is_real(n)) { + // complex number case + // acos(n) = -i * ln(n + sqrt(n^2 - 1)) + // again, we have the building blocks needed to compute this + const i = SchemeComplex.build( + SchemeInteger.EXACT_ZERO, + SchemeInteger.build(1), + ); + return atomic_multiply( + atomic_negate(i), + log( + atomic_add( + n, + sqrt(atomic_subtract(atomic_multiply(n, n), SchemeInteger.build(1))), + ), + ), + ); + } + return SchemeReal.build(Math.acos(n.coerce())); +}; + +export const atan = (n: SchemeNumber, m?: SchemeNumber): SchemeNumber => { + if (!is_number(n)) { + throw new Error("atan: expected number"); + } + + if (m !== undefined) { + // two argument case, we construct a complex number with n + mi + // if neither n nor m are real, it's an error + if (!is_real(n) || !is_real(m)) { + throw new Error("atan: expected real numbers"); + } + return atan( + SchemeComplex.build( + n as SchemeInteger | SchemeRational | SchemeReal, + m as SchemeInteger | SchemeRational | SchemeReal, + ), + ); + } + + if (!is_real(n)) { + // complex number case + // atan(n) = 1/2 * i * ln((1 - i * n) / (1 + i * n)) + const i = SchemeComplex.build( + SchemeInteger.EXACT_ZERO, + SchemeInteger.build(1), + ); + return atomic_multiply( + // multiply is associative so the order here doesn't matter + atomic_multiply(SchemeRational.build(1, 2), i), + log( + atomic_divide( + atomic_subtract(SchemeInteger.build(1), atomic_multiply(i, n)), + atomic_add(SchemeInteger.build(1), atomic_multiply(i, n)), + ), + ), + ); + } + return SchemeReal.build(Math.atan(n.coerce())); +};