Prototype of generic complex numbers

This commit is contained in:
Iaroslav Postovalov 2021-03-16 23:17:54 +07:00
parent c377d6cda8
commit 0e89e783fd
35 changed files with 917 additions and 788 deletions

View File

@ -20,9 +20,9 @@ import java.math.BigInteger
internal class BigIntBenchmark { internal class BigIntBenchmark {
val kmNumber = BigIntField.number(Int.MAX_VALUE) val kmNumber = BigIntField.number(Int.MAX_VALUE)
val jvmNumber = JBigIntegerField.number(Int.MAX_VALUE) val jvmNumber = JBigIntegerRing.number(Int.MAX_VALUE)
val largeKmNumber = BigIntField { number(11).pow(100_000U) } val largeKmNumber = BigIntField { number(11).pow(100_000U) }
val largeJvmNumber: BigInteger = JBigIntegerField { number(11).pow(100_000) } val largeJvmNumber: BigInteger = JBigIntegerRing { number(11).pow(100_000) }
val bigExponent = 50_000 val bigExponent = 50_000
@Benchmark @Benchmark
@ -31,7 +31,7 @@ internal class BigIntBenchmark {
} }
@Benchmark @Benchmark
fun jvmAdd(blackhole: Blackhole) = JBigIntegerField { fun jvmAdd(blackhole: Blackhole) = JBigIntegerRing {
blackhole.consume(jvmNumber + jvmNumber + jvmNumber) blackhole.consume(jvmNumber + jvmNumber + jvmNumber)
} }
@ -41,7 +41,7 @@ internal class BigIntBenchmark {
} }
@Benchmark @Benchmark
fun jvmAddLarge(blackhole: Blackhole) = JBigIntegerField { fun jvmAddLarge(blackhole: Blackhole) = JBigIntegerRing {
blackhole.consume(largeJvmNumber + largeJvmNumber + largeJvmNumber) blackhole.consume(largeJvmNumber + largeJvmNumber + largeJvmNumber)
} }
@ -56,12 +56,12 @@ internal class BigIntBenchmark {
} }
@Benchmark @Benchmark
fun jvmMultiply(blackhole: Blackhole) = JBigIntegerField { fun jvmMultiply(blackhole: Blackhole) = JBigIntegerRing {
blackhole.consume(jvmNumber * jvmNumber * jvmNumber) blackhole.consume(jvmNumber * jvmNumber * jvmNumber)
} }
@Benchmark @Benchmark
fun jvmMultiplyLarge(blackhole: Blackhole) = JBigIntegerField { fun jvmMultiplyLarge(blackhole: Blackhole) = JBigIntegerRing {
blackhole.consume(largeJvmNumber*largeJvmNumber) blackhole.consume(largeJvmNumber*largeJvmNumber)
} }
@ -71,27 +71,27 @@ internal class BigIntBenchmark {
} }
@Benchmark @Benchmark
fun jvmPower(blackhole: Blackhole) = JBigIntegerField { fun jvmPower(blackhole: Blackhole) = JBigIntegerRing {
blackhole.consume(jvmNumber.pow(bigExponent)) blackhole.consume(jvmNumber.pow(bigExponent))
} }
@Benchmark @Benchmark
fun kmParsing16(blackhole: Blackhole) = JBigIntegerField { fun kmParsing16(blackhole: Blackhole) = JBigIntegerRing {
blackhole.consume("0x7f57ed8b89c29a3b9a85c7a5b84ca3929c7b7488593".parseBigInteger()) blackhole.consume("0x7f57ed8b89c29a3b9a85c7a5b84ca3929c7b7488593".parseBigInteger())
} }
@Benchmark @Benchmark
fun kmParsing10(blackhole: Blackhole) = JBigIntegerField { fun kmParsing10(blackhole: Blackhole) = JBigIntegerRing {
blackhole.consume("236656783929183747565738292847574838922010".parseBigInteger()) blackhole.consume("236656783929183747565738292847574838922010".parseBigInteger())
} }
@Benchmark @Benchmark
fun jvmParsing10(blackhole: Blackhole) = JBigIntegerField { fun jvmParsing10(blackhole: Blackhole) = JBigIntegerRing {
blackhole.consume("236656783929183747565738292847574838922010".toBigInteger(10)) blackhole.consume("236656783929183747565738292847574838922010".toBigInteger(10))
} }
@Benchmark @Benchmark
fun jvmParsing16(blackhole: Blackhole) = JBigIntegerField { fun jvmParsing16(blackhole: Blackhole) = JBigIntegerRing {
blackhole.consume("7f57ed8b89c29a3b9a85c7a5b84ca3929c7b7488593".toBigInteger(16)) blackhole.consume("7f57ed8b89c29a3b9a85c7a5b84ca3929c7b7488593".toBigInteger(16))
} }
} }

View File

@ -17,7 +17,7 @@ import space.kscience.kmath.structures.MutableBuffer
internal class BufferBenchmark { internal class BufferBenchmark {
@Benchmark @Benchmark
fun genericDoubleBufferReadWrite() { fun genericDoubleBufferReadWrite() {
val buffer = DoubleBuffer(size) { it.toDouble() } val buffer = DoubleBuffer(size, Int::toDouble)
(0 until size).forEach { (0 until size).forEach {
buffer[it] buffer[it]
@ -26,7 +26,8 @@ internal class BufferBenchmark {
@Benchmark @Benchmark
fun complexBufferReadWrite() { fun complexBufferReadWrite() {
val buffer = MutableBuffer.complex(size / 2) { Complex(it.toDouble(), -it.toDouble()) } val buffer =
MutableBuffer.complex(MutableBuffer.Companion::double, size / 2) { Complex(it.toDouble(), -it.toDouble()) }
(0 until size / 2).forEach { (0 until size / 2).forEach {
buffer[it] buffer[it]

View File

@ -47,7 +47,7 @@ internal class NDFieldBenchmark {
private const val dim = 1000 private const val dim = 1000
private const val n = 100 private const val n = 100
private val autoField = AlgebraND.auto(DoubleField, dim, dim) private val autoField = AlgebraND.auto(DoubleField, dim, dim)
private val specializedField = AlgebraND.real(dim, dim) private val specializedField = AlgebraND.double(dim, dim)
private val genericField = AlgebraND.field(DoubleField, Buffer.Companion::boxing, dim, dim) private val genericField = AlgebraND.field(DoubleField, Buffer.Companion::boxing, dim, dim)
} }
} }

View File

@ -13,7 +13,7 @@ import org.jetbrains.bio.viktor.F64Array
import space.kscience.kmath.nd.AlgebraND import space.kscience.kmath.nd.AlgebraND
import space.kscience.kmath.nd.StructureND import space.kscience.kmath.nd.StructureND
import space.kscience.kmath.nd.auto import space.kscience.kmath.nd.auto
import space.kscience.kmath.nd.real import space.kscience.kmath.nd.double
import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.operations.DoubleField
import space.kscience.kmath.viktor.ViktorNDField import space.kscience.kmath.viktor.ViktorNDField
@ -60,7 +60,7 @@ internal class ViktorBenchmark {
// automatically build context most suited for given type. // automatically build context most suited for given type.
private val autoField = AlgebraND.auto(DoubleField, dim, dim) private val autoField = AlgebraND.auto(DoubleField, dim, dim)
private val realField = AlgebraND.real(dim, dim) private val realField = AlgebraND.double(dim, dim)
private val viktorField = ViktorNDField(dim, dim) private val viktorField = ViktorNDField(dim, dim)
} }
} }

View File

@ -12,7 +12,7 @@ import kotlinx.benchmark.State
import org.jetbrains.bio.viktor.F64Array import org.jetbrains.bio.viktor.F64Array
import space.kscience.kmath.nd.AlgebraND import space.kscience.kmath.nd.AlgebraND
import space.kscience.kmath.nd.auto import space.kscience.kmath.nd.auto
import space.kscience.kmath.nd.real import space.kscience.kmath.nd.double
import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.operations.DoubleField
import space.kscience.kmath.viktor.ViktorFieldND import space.kscience.kmath.viktor.ViktorFieldND
@ -52,7 +52,7 @@ internal class ViktorLogBenchmark {
// automatically build context most suited for given type. // automatically build context most suited for given type.
private val autoField = AlgebraND.auto(DoubleField, dim, dim) private val autoField = AlgebraND.auto(DoubleField, dim, dim)
private val realNdField = AlgebraND.real(dim, dim) private val realNdField = AlgebraND.double(dim, dim)
private val viktorField = ViktorFieldND(intArrayOf(dim, dim)) private val viktorField = ViktorFieldND(intArrayOf(dim, dim))
} }
} }

View File

@ -8,16 +8,17 @@ package space.kscience.kmath.operations
import space.kscience.kmath.complex.Complex import space.kscience.kmath.complex.Complex
import space.kscience.kmath.complex.complex import space.kscience.kmath.complex.complex
import space.kscience.kmath.nd.AlgebraND import space.kscience.kmath.nd.AlgebraND
import space.kscience.kmath.nd.double
fun main() { fun main() {
// 2d element // 2d element
val element = AlgebraND.complex(2, 2).produce { (i, j) -> val element = AlgebraND.double(2, 2).complex().produce { (i, j) ->
Complex(i.toDouble() - j.toDouble(), i.toDouble() + j.toDouble()) Complex(i.toDouble() - j.toDouble(), i.toDouble() + j.toDouble())
} }
println(element) println(element)
// 1d element operation // 1d element operation
val result = with(AlgebraND.complex(8)) { val result = with(AlgebraND.double(8).complex()) {
val a = produce { (it) -> i * it - it.toDouble() } val a = produce { (it) -> i * it - it.toDouble() }
val b = 3 val b = 3
val c = Complex(1.0, 1.0) val c = Complex(1.0, 1.0)

View File

@ -12,19 +12,19 @@ import space.kscience.kmath.linear.transpose
import space.kscience.kmath.nd.AlgebraND import space.kscience.kmath.nd.AlgebraND
import space.kscience.kmath.nd.StructureND import space.kscience.kmath.nd.StructureND
import space.kscience.kmath.nd.as2D import space.kscience.kmath.nd.as2D
import space.kscience.kmath.nd.real import space.kscience.kmath.nd.double
import space.kscience.kmath.operations.invoke import space.kscience.kmath.operations.*
import kotlin.system.measureTimeMillis import kotlin.system.measureTimeMillis
fun main() { fun main() {
val dim = 1000 val dim = 1000
val n = 1000 val n = 1000
val realField = AlgebraND.real(dim, dim) val doubleField = AlgebraND.double(dim, dim)
val complexField: ComplexFieldND = AlgebraND.complex(dim, dim) val complexField = doubleField.complex()
val realTime = measureTimeMillis { val realTime = measureTimeMillis {
realField { doubleField {
var res: StructureND<Double> = one var res: StructureND<Double> = one
repeat(n) { repeat(n) {
res += 1.0 res += 1.0
@ -36,9 +36,9 @@ fun main() {
val complexTime = measureTimeMillis { val complexTime = measureTimeMillis {
complexField { complexField {
var res: StructureND<Complex> = one var res: StructureND<Complex<Double>> = one
repeat(n) { repeat(n) {
res += 1.0 res += Complex(1.0, 0.0)
} }
} }
} }
@ -48,18 +48,16 @@ fun main() {
fun complexExample() { fun complexExample() {
//Create a context for 2-d structure with complex values //Create a context for 2-d structure with complex values
ComplexField { AlgebraND.double(4, 8).complex().run {
nd(4, 8) {
//a constant real-valued structure //a constant real-valued structure
val x = one * 2.5 val x = one * 2.5
operator fun Number.plus(other: Complex) = Complex(this.toDouble() + other.re, other.im) operator fun Number.plus(other: Complex<Double>) = Complex(toDouble() + other.re, other.im)
//a structure generator specific to this context //a structure generator specific to this context
val matrix = produce { (k, l) -> k + l * i } val matrix = produce { (k, l) -> k + l * i }
//Perform sum //Perform sum
val sum = matrix + x + 1.0 val sum = matrix + x + Complex(1.0,0.0)
//Represent the sum as 2d-structure and transpose //Represent the sum as 2d-structure and transpose
sum.as2D().transpose() sum.as2D().transpose()
} }
}
} }

View File

@ -32,8 +32,8 @@ fun main() {
// automatically build context most suited for given type. // automatically build context most suited for given type.
val autoField = AlgebraND.auto(DoubleField, dim, dim) val autoField = AlgebraND.auto(DoubleField, dim, dim)
// specialized nd-field for Double. It works as generic Double field as well. // specialized nd-field for Double. It works as generic Double field as well
val realField = AlgebraND.real(dim, dim) val realField = AlgebraND.double(dim, dim)
//A generic boxing field. It should be used for objects, not primitives. //A generic boxing field. It should be used for objects, not primitives.
val boxingField = AlgebraND.field(DoubleField, Buffer.Companion::boxing, dim, dim) val boxingField = AlgebraND.field(DoubleField, Buffer.Companion::boxing, dim, dim)
// Nd4j specialized field. // Nd4j specialized field.

View File

@ -6,6 +6,7 @@
package space.kscience.kmath.ast package space.kscience.kmath.ast
import space.kscience.kmath.complex.Complex import space.kscience.kmath.complex.Complex
import space.kscience.kmath.complex.ComplexDoubleField
import space.kscience.kmath.complex.ComplexField import space.kscience.kmath.complex.ComplexField
import space.kscience.kmath.expressions.evaluate import space.kscience.kmath.expressions.evaluate
import space.kscience.kmath.operations.Algebra import space.kscience.kmath.operations.Algebra
@ -17,15 +18,15 @@ internal class TestParser {
@Test @Test
fun evaluateParsedMst() { fun evaluateParsedMst() {
val mst = "2+2*(2+2)".parseMath() val mst = "2+2*(2+2)".parseMath()
val res = ComplexField.evaluate(mst) val res = ComplexDoubleField.evaluate(mst)
assertEquals(Complex(10.0, 0.0), res) assertEquals(Complex(10.0, 0.0), res)
} }
@Test @Test
fun evaluateMstSymbol() { fun evaluateMstSymbol() {
val mst = "i".parseMath() val mst = "i".parseMath()
val res = ComplexField.evaluate(mst) val res = ComplexDoubleField.evaluate(mst)
assertEquals(ComplexField.i, res) assertEquals(ComplexDoubleField.i, res)
} }

View File

@ -15,12 +15,11 @@ import space.kscience.kmath.streaming.spread
import space.kscience.kmath.structures.* import space.kscience.kmath.structures.*
/** /**
* Streaming and buffer transformations * Streaming and buffer transformations
*/ */
public object Transformations { public object Transformations {
private fun Buffer<Complex>.toArray(): Array<org.apache.commons.math3.complex.Complex> = private fun Buffer<Complex<Double>>.toArray(): Array<org.apache.commons.math3.complex.Complex> =
Array(size) { org.apache.commons.math3.complex.Complex(get(it).re, get(it).im) } Array(size) { org.apache.commons.math3.complex.Complex(get(it).re, get(it).im) }
private fun Buffer<Double>.asArray() = if (this is DoubleBuffer) { private fun Buffer<Double>.asArray() = if (this is DoubleBuffer) {
@ -40,14 +39,14 @@ public object Transformations {
public fun fourier( public fun fourier(
normalization: DftNormalization = DftNormalization.STANDARD, normalization: DftNormalization = DftNormalization.STANDARD,
direction: TransformType = TransformType.FORWARD, direction: TransformType = TransformType.FORWARD,
): SuspendBufferTransform<Complex, Complex> = { ): SuspendBufferTransform<Complex<Double>, Complex<Double>> = {
FastFourierTransformer(normalization).transform(it.toArray(), direction).asBuffer() FastFourierTransformer(normalization).transform(it.toArray(), direction).asBuffer()
} }
public fun realFourier( public fun realFourier(
normalization: DftNormalization = DftNormalization.STANDARD, normalization: DftNormalization = DftNormalization.STANDARD,
direction: TransformType = TransformType.FORWARD, direction: TransformType = TransformType.FORWARD,
): SuspendBufferTransform<Double, Complex> = { ): SuspendBufferTransform<Double, Complex<Double>> = {
FastFourierTransformer(normalization).transform(it.asArray(), direction).asBuffer() FastFourierTransformer(normalization).transform(it.asArray(), direction).asBuffer()
} }
@ -76,10 +75,10 @@ public object Transformations {
* Process given [Flow] with commons-math fft transformation * Process given [Flow] with commons-math fft transformation
*/ */
@FlowPreview @FlowPreview
public fun Flow<Buffer<Complex>>.FFT( public fun Flow<Buffer<Complex<Double>>>.FFT(
normalization: DftNormalization = DftNormalization.STANDARD, normalization: DftNormalization = DftNormalization.STANDARD,
direction: TransformType = TransformType.FORWARD, direction: TransformType = TransformType.FORWARD,
): Flow<Buffer<Complex>> { ): Flow<Buffer<Complex<Double>>> {
val transform = Transformations.fourier(normalization, direction) val transform = Transformations.fourier(normalization, direction)
return map { transform(it) } return map { transform(it) }
} }
@ -89,7 +88,7 @@ public fun Flow<Buffer<Complex>>.FFT(
public fun Flow<Buffer<Double>>.FFT( public fun Flow<Buffer<Double>>.FFT(
normalization: DftNormalization = DftNormalization.STANDARD, normalization: DftNormalization = DftNormalization.STANDARD,
direction: TransformType = TransformType.FORWARD, direction: TransformType = TransformType.FORWARD,
): Flow<Buffer<Complex>> { ): Flow<Buffer<Complex<Double>>> {
val transform = Transformations.realFourier(normalization, direction) val transform = Transformations.realFourier(normalization, direction)
return map(transform) return map(transform)
} }
@ -103,10 +102,10 @@ public fun Flow<Double>.FFT(
bufferSize: Int = Int.MAX_VALUE, bufferSize: Int = Int.MAX_VALUE,
normalization: DftNormalization = DftNormalization.STANDARD, normalization: DftNormalization = DftNormalization.STANDARD,
direction: TransformType = TransformType.FORWARD, direction: TransformType = TransformType.FORWARD,
): Flow<Complex> = chunked(bufferSize).FFT(normalization, direction).spread() ): Flow<Complex<Double>> = chunked(bufferSize).FFT(normalization, direction).spread()
/** /**
* Map a complex flow into real flow by taking real part of each number * Map a complex flow into real flow by taking real part of each number
*/ */
@FlowPreview @FlowPreview
public fun Flow<Complex>.real(): Flow<Double> = map { it.re } public fun Flow<Complex<Double>>.real(): Flow<Double> = map { it.re }

View File

@ -5,223 +5,248 @@
package space.kscience.kmath.complex package space.kscience.kmath.complex
import space.kscience.kmath.memory.MemoryReader import space.kscience.kmath.operations.*
import space.kscience.kmath.memory.MemorySpec import kotlin.js.JsName
import space.kscience.kmath.memory.MemoryWriter import kotlin.jvm.JvmName
import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.operations.ExtendedField
import space.kscience.kmath.operations.Norm
import space.kscience.kmath.operations.NumbersAddOperations
import space.kscience.kmath.operations.ScaleOperations
import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.MemoryBuffer
import space.kscience.kmath.structures.MutableBuffer
import space.kscience.kmath.structures.MutableMemoryBuffer
import kotlin.math.*
/** /**
* This complex's conjugate. * Represents generic complex value consisting of real and imaginary part.
*
* @param T the type of components.
* @property re The real component.
* @property im The imaginary component.
*/ */
public val Complex.conjugate: Complex public data class Complex<out T : Any>(public val re: T, public val im: T) {
get() = Complex(re, -im) /**
* Converts this complex number to string formatted like `[re] + i * [im]`.
*/
override fun toString(): String = "$re + i * $im"
}
/** /**
* This complex's reciprocal. * The algebra of [Complex].
*
* @param T the type of components.
* @property algebra the algebra over [T].
*/ */
public val Complex.reciprocal: Complex public open class ComplexAlgebra<T : Any>(public open val algebra: NumericAlgebra<T>) : NumericAlgebra<Complex<T>> {
get() {
val scale = re * re + im * im
return Complex(re / scale, -im / scale)
}
/**
* Absolute value of complex number.
*/
public val Complex.r: Double
get() = sqrt(re * re + im * im)
/**
* An angle between vector represented by complex number and X axis.
*/
public val Complex.theta: Double
get() = atan(im / re)
private val PI_DIV_2 = Complex(PI / 2, 0)
/**
* A field of [Complex].
*/
@OptIn(UnstableKMathAPI::class)
public object ComplexField : ExtendedField<Complex>, Norm<Complex, Complex>, NumbersAddOperations<Complex>,
ScaleOperations<Complex> {
override val zero: Complex = 0.0.toComplex()
override val one: Complex = 1.0.toComplex()
/** /**
* The imaginary unit. * The imaginary unit.
*/ */
public val i: Complex by lazy { Complex(0.0, 1.0) } public open val i: Complex<T> by lazy {
algebra { Complex(number(0), number(1)) }
override fun Complex.unaryMinus(): Complex = Complex(-re, -im)
override fun number(value: Number): Complex = Complex(value.toDouble(), 0.0)
override fun scale(a: Complex, value: Double): Complex = Complex(a.re * value, a.im * value)
override fun add(a: Complex, b: Complex): Complex = Complex(a.re + b.re, a.im + b.im)
// override fun multiply(a: Complex, k: Number): Complex = Complex(a.re * k.toDouble(), a.im * k.toDouble())
override fun multiply(a: Complex, b: Complex): Complex =
Complex(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re)
override fun divide(a: Complex, b: Complex): Complex = when {
abs(b.im) < abs(b.re) -> {
val wr = b.im / b.re
val wd = b.re + wr * b.im
if (wd.isNaN() || wd == 0.0)
throw ArithmeticException("Division by zero or infinity")
else
Complex((a.re + a.im * wr) / wd, (a.im - a.re * wr) / wd)
} }
b.im == 0.0 -> throw ArithmeticException("Division by zero") override fun number(value: Number): Complex<T> =
algebra { Complex(algebra.number(value), algebra.number(0)) }
else -> { override fun bindSymbol(value: String): Complex<T> = if (value == "i") i else super.bindSymbol(value)
val wr = b.re / b.im }
val wd = b.im + wr * b.re
if (wd.isNaN() || wd == 0.0) /**
throw ArithmeticException("Division by zero or infinity") * The group of [Complex].
else *
Complex((a.re * wr + a.im) / wd, (a.im * wr - a.re) / wd) * @param T the type of components.
} */
public open class ComplexGroup<T : Any, out A>(override val algebra: A) : ComplexAlgebra<T>(algebra),
Group<Complex<T>> where A : NumericAlgebra<T>, A : Group<T> {
override val zero: Complex<T> by lazy {
algebra { Complex(zero, zero) }
} }
override operator fun Complex.div(k: Number): Complex = Complex(re / k.toDouble(), im / k.toDouble()) /**
* This complex's conjugate.
*/
public val Complex<T>.conjugate: Complex<T>
get() = Complex(re, algebra { -im })
override fun sin(arg: Complex): Complex = i * (exp(-i * arg) - exp(i * arg)) / 2.0 override fun add(a: Complex<T>, b: Complex<T>): Complex<T> = algebra { Complex(a.re + b.re, a.im + b.im) }
override fun cos(arg: Complex): Complex = (exp(-i * arg) + exp(i * arg)) / 2.0
override fun tan(arg: Complex): Complex { override fun Complex<T>.unaryMinus(): Complex<T> = algebra { Complex(-re, -im) }
@JsName("unaryMinus_T")
public operator fun T.unaryMinus(): Complex<T> = algebra { Complex(-this@unaryMinus, zero) }
@JsName("unaryPlus_T")
public operator fun T.unaryPlus(): Complex<T> = algebra { Complex(this@unaryPlus, zero) }
public operator fun T.plus(b: Complex<T>): Complex<T> = add(+this, b)
public operator fun Complex<T>.plus(b: T): Complex<T> = add(this, +b)
public operator fun T.minus(b: Complex<T>): Complex<T> = add(+this, -b)
public operator fun Complex<T>.minus(b: T): Complex<T> = add(this, -b)
}
/**
* The ring of [Complex].
*
* @param T the type of components.
*/
public open class ComplexRing<T : Any, out A>(override val algebra: A) : ComplexGroup<T, A>(algebra),
Ring<Complex<T>> where A : NumericAlgebra<T>, A : Ring<T> {
override val one: Complex<T> by lazy {
algebra { Complex(one, zero) }
}
override val i: Complex<T> by lazy {
algebra { Complex(zero, one) }
}
override fun multiply(a: Complex<T>, b: Complex<T>): Complex<T> =
algebra { Complex(a.re * b.re - a.im * b.im, a.im * b.re + a.re * b.im) }
public operator fun T.times(b: Complex<T>): Complex<T> = multiply(+this, b)
public operator fun Complex<T>.times(b: T): Complex<T> = multiply(this, +b)
}
/**
* [ComplexRing] instance for [ByteRing].
*/
public val ComplexByteRing: ComplexRing<Byte, ByteRing> = ComplexRing(ByteRing)
/**
* [ComplexRing] instance for [ShortRing].
*/
public val ComplexShortRing: ComplexRing<Short, ShortRing> = ComplexRing(ShortRing)
/**
* [ComplexRing] instance for [IntRing].
*/
public val ComplexIntRing: ComplexRing<Int, IntRing> = ComplexRing(IntRing)
/**
* [ComplexRing] instance for [LongRing].
*/
public val ComplexLongRing: ComplexRing<Long, LongRing> = ComplexRing(LongRing)
/**
* The field of [Complex].
*/
public open class ComplexField<T : Any, out A>(override val algebra: A) : ComplexRing<T, A>(algebra),
Field<Complex<T>> where A : Field<T> {
/**
* This complex's reciprocal.
*/
public val Complex<T>.reciprocal: Complex<T>
get() = algebra {
val scale = re * re + im * im
Complex(re / scale, -im / scale)
}
override fun divide(a: Complex<T>, b: Complex<T>): Complex<T> = a * b.reciprocal
override fun number(value: Number): Complex<T> = super<ComplexRing>.number(value)
override fun scale(a: Complex<T>, value: Double): Complex<T> =
algebra { Complex(a.re * value, a.im * value) }
override operator fun Complex<T>.div(k: Number): Complex<T> =
algebra { Complex(re / k.toDouble(), im / k.toDouble()) }
public operator fun T.div(b: Complex<T>): Complex<T> = divide(+this, b)
public operator fun Complex<T>.div(b: T): Complex<T> = divide(this, +b)
@JsName("scale_T")
public fun scale(a: T, value: Double): Complex<T> = scale(+a, value)
}
/**
* [ComplexRing] instance for [BigIntField].
*/
public val ComplexBigIntField: ComplexField<BigInt, BigIntField> = ComplexField(BigIntField)
/**
* The extended field of Complex.
*/
public open class ComplexExtendedField<T : Any, out A>(override val algebra: A) : ComplexField<T, A>(algebra),
ExtendedField<Complex<T>>, Norm<Complex<T>, T> where A : ExtendedField<T> {
private val two by lazy { one + one }
/**
* The *r* component of the polar form of this number.
*/
public val Complex<T>.r: T
get() = norm(this)
/**
* The *&theta;* component of the polar form of this number.
*/
public val Complex<T>.theta: T
get() = algebra { atan(im / re) }
override fun bindSymbol(value: String): Complex<T> =
if (value == "i") i else super<ExtendedField>.bindSymbol(value)
override fun sin(arg: Complex<T>): Complex<T> = i * (exp(-i * arg) - exp(i * arg)) / two
override fun cos(arg: Complex<T>): Complex<T> = (exp(-i * arg) + exp(i * arg)) / two
override fun tan(arg: Complex<T>): Complex<T> {
val e1 = exp(-i * arg) val e1 = exp(-i * arg)
val e2 = exp(i * arg) val e2 = exp(i * arg)
return i * (e1 - e2) / (e1 + e2) return i * (e1 - e2) / (e1 + e2)
} }
override fun asin(arg: Complex): Complex = -i * ln(sqrt(1 - (arg * arg)) + i * arg) override fun asin(arg: Complex<T>): Complex<T> = -i * ln(sqrt(one - (arg * arg)) + i * arg)
override fun acos(arg: Complex): Complex = PI_DIV_2 + i * ln(sqrt(1 - (arg * arg)) + i * arg) override fun acos(arg: Complex<T>): Complex<T> =
(pi / two) + i * ln(sqrt(one - (arg * arg)) + i * arg)
override fun atan(arg: Complex): Complex { override fun atan(arg: Complex<T>): Complex<T> = algebra {
val iArg = i * arg val iArg = i * arg
return i * (ln(1 - iArg) - ln(1 + iArg)) / 2 return i * (ln(this@ComplexExtendedField.one - iArg) - ln(this@ComplexExtendedField.one + iArg)) / 2
} }
override fun power(arg: Complex, pow: Number): Complex = if (arg.im == 0.0) override fun power(arg: Complex<T>, pow: Number): Complex<T> = algebra {
arg.re.pow(pow.toDouble()).toComplex() if (arg.im == 0.0)
Complex(arg.re.pow(pow.toDouble()), algebra.zero)
else else
exp(pow * ln(arg)) exp(pow * ln(arg))
}
override fun exp(arg: Complex): Complex = exp(arg.re) * (cos(arg.im) + i * sin(arg.im)) override fun exp(arg: Complex<T>): Complex<T> =
Complex(algebra.exp(arg.re), algebra.zero) * Complex(algebra.cos(arg.im), algebra.sin(arg.im))
override fun ln(arg: Complex): Complex = ln(arg.r) + i * atan2(arg.im, arg.re) override fun ln(arg: Complex<T>): Complex<T> = algebra { Complex(ln(norm(arg)), atan(arg.im / arg.re)) }
override fun norm(arg: Complex<T>): T = algebra { sqrt(arg.re * arg.re + arg.im * arg.im) }
/** @JsName("norm_T_3")
* Adds complex number to real one. @JvmName("norm\$T")
* public fun norm(arg: T): T = algebra { sqrt(arg * arg) }
* @receiver the augend.
* @param c the addend.
* @return the sum.
*/
public operator fun Double.plus(c: Complex): Complex = add(this.toComplex(), c)
/** @JsName("sin_T")
* Subtracts complex number from real one. public fun sin(arg: T): Complex<T> = sin(+arg)
*
* @receiver the minuend.
* @param c the subtrahend.
* @return the difference.
*/
public operator fun Double.minus(c: Complex): Complex = add(this.toComplex(), -c)
/** @JsName("cos_T")
* Adds real number to complex one. public fun cos(arg: T): Complex<T> = cos(+arg)
*
* @receiver the augend.
* @param d the addend.
* @return the sum.
*/
public operator fun Complex.plus(d: Double): Complex = d + this
/** @JsName("tan_T")
* Subtracts real number from complex one. public fun tan(arg: T): Complex<T> = tan(+arg)
*
* @receiver the minuend.
* @param d the subtrahend.
* @return the difference.
*/
public operator fun Complex.minus(d: Double): Complex = add(this, -d.toComplex())
/** @JsName("asin_T")
* Multiplies real number by complex one. public fun asin(arg: T): Complex<T> = asin(+arg)
*
* @receiver the multiplier.
* @param c the multiplicand.
* @receiver the product.
*/
public operator fun Double.times(c: Complex): Complex = Complex(c.re * this, c.im * this)
override fun norm(arg: Complex): Complex = sqrt(arg.conjugate * arg) @JsName("acos_T")
public fun acos(arg: T): Complex<T> = acos(+arg)
override fun bindSymbolOrNull(value: String): Complex? = if (value == "i") i else null @JsName("atan_T")
public fun atan(arg: T): Complex<T> = atan(+arg)
@JsName("power_T")
public fun power(arg: T, pow: Number): Complex<T> = power(+arg, pow)
@JsName("exp_T")
public fun exp(arg: T): Complex<T> = exp(+arg)
@JsName("ln_T")
public fun ln(arg: T): Complex<T> = ln(+arg)
} }
/** /**
* Represents `double`-based complex number. * [ComplexRing] instance for [DoubleField].
*
* @property re The real part.
* @property im The imaginary part.
*/ */
@OptIn(UnstableKMathAPI::class) public val ComplexDoubleField: ComplexExtendedField<Double, DoubleField> = ComplexExtendedField(DoubleField)
public data class Complex(val re: Double, val im: Double) {
public constructor(re: Number, im: Number) : this(re.toDouble(), im.toDouble())
public constructor(re: Number) : this(re.toDouble(), 0.0)
override fun toString(): String = "($re + i * $im)"
public companion object : MemorySpec<Complex> {
override val objectSize: Int
get() = 16
override fun MemoryReader.read(offset: Int): Complex =
Complex(readDouble(offset), readDouble(offset + 8))
override fun MemoryWriter.write(offset: Int, value: Complex) {
writeDouble(offset, value.re)
writeDouble(offset + 8, value.im)
}
}
}
/** /**
* Creates a complex number with real part equal to this real. * [ComplexRing] instance for [FloatField].
*
* @receiver the real part.
* @return the new complex number.
*/ */
public fun Number.toComplex(): Complex = Complex(this) public val ComplexFloatField: ComplexExtendedField<Double, DoubleField> = ComplexExtendedField(DoubleField)
/**
* Creates a new buffer of complex numbers with the specified [size], where each element is calculated by calling the
* specified [init] function.
*/
public inline fun Buffer.Companion.complex(size: Int, init: (Int) -> Complex): Buffer<Complex> =
MemoryBuffer.create(Complex, size, init)
/**
* Creates a new buffer of complex numbers with the specified [size], where each element is calculated by calling the
* specified [init] function.
*/
public inline fun MutableBuffer.Companion.complex(size: Int, init: (Int) -> Complex): MutableBuffer<Complex> =
MutableMemoryBuffer.create(Complex, size, init)

View File

@ -1,124 +0,0 @@
/*
* Copyright 2018-2021 KMath contributors.
* Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
*/
package space.kscience.kmath.complex
import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.nd.AlgebraND
import space.kscience.kmath.nd.BufferND
import space.kscience.kmath.nd.BufferedFieldND
import space.kscience.kmath.nd.StructureND
import space.kscience.kmath.operations.ExtendedField
import space.kscience.kmath.operations.NumbersAddOperations
import space.kscience.kmath.structures.Buffer
import kotlin.contracts.InvocationKind
import kotlin.contracts.contract
/**
* An optimized nd-field for complex numbers
*/
@OptIn(UnstableKMathAPI::class)
public class ComplexFieldND(
shape: IntArray,
) : BufferedFieldND<Complex, ComplexField>(shape, ComplexField, Buffer.Companion::complex),
NumbersAddOperations<StructureND<Complex>>,
ExtendedField<StructureND<Complex>> {
override val zero: BufferND<Complex> by lazy { produce { zero } }
override val one: BufferND<Complex> by lazy { produce { one } }
override fun number(value: Number): BufferND<Complex> {
val d = value.toComplex() // minimize conversions
return produce { d }
}
//
// @Suppress("OVERRIDE_BY_INLINE")
// override inline fun map(
// arg: AbstractNDBuffer<Double>,
// transform: DoubleField.(Double) -> Double,
// ): RealNDElement {
// check(arg)
// val array = RealBuffer(arg.strides.linearSize) { offset -> DoubleField.transform(arg.buffer[offset]) }
// return BufferedNDFieldElement(this, array)
// }
//
// @Suppress("OVERRIDE_BY_INLINE")
// override inline fun produce(initializer: DoubleField.(IntArray) -> Double): RealNDElement {
// val array = RealBuffer(strides.linearSize) { offset -> elementContext.initializer(strides.index(offset)) }
// return BufferedNDFieldElement(this, array)
// }
//
// @Suppress("OVERRIDE_BY_INLINE")
// override inline fun mapIndexed(
// arg: AbstractNDBuffer<Double>,
// transform: DoubleField.(index: IntArray, Double) -> Double,
// ): RealNDElement {
// check(arg)
// return BufferedNDFieldElement(
// this,
// RealBuffer(arg.strides.linearSize) { offset ->
// elementContext.transform(
// arg.strides.index(offset),
// arg.buffer[offset]
// )
// })
// }
//
// @Suppress("OVERRIDE_BY_INLINE")
// override inline fun combine(
// a: AbstractNDBuffer<Double>,
// b: AbstractNDBuffer<Double>,
// transform: DoubleField.(Double, Double) -> Double,
// ): RealNDElement {
// check(a, b)
// val buffer = RealBuffer(strides.linearSize) { offset ->
// elementContext.transform(a.buffer[offset], b.buffer[offset])
// }
// return BufferedNDFieldElement(this, buffer)
// }
override fun power(arg: StructureND<Complex>, pow: Number): BufferND<Complex> = arg.map { power(it, pow) }
override fun exp(arg: StructureND<Complex>): BufferND<Complex> = arg.map { exp(it) }
override fun ln(arg: StructureND<Complex>): BufferND<Complex> = arg.map { ln(it) }
override fun sin(arg: StructureND<Complex>): BufferND<Complex> = arg.map { sin(it) }
override fun cos(arg: StructureND<Complex>): BufferND<Complex> = arg.map { cos(it) }
override fun tan(arg: StructureND<Complex>): BufferND<Complex> = arg.map { tan(it) }
override fun asin(arg: StructureND<Complex>): BufferND<Complex> = arg.map { asin(it) }
override fun acos(arg: StructureND<Complex>): BufferND<Complex> = arg.map { acos(it) }
override fun atan(arg: StructureND<Complex>): BufferND<Complex> = arg.map { atan(it) }
override fun sinh(arg: StructureND<Complex>): BufferND<Complex> = arg.map { sinh(it) }
override fun cosh(arg: StructureND<Complex>): BufferND<Complex> = arg.map { cosh(it) }
override fun tanh(arg: StructureND<Complex>): BufferND<Complex> = arg.map { tanh(it) }
override fun asinh(arg: StructureND<Complex>): BufferND<Complex> = arg.map { asinh(it) }
override fun acosh(arg: StructureND<Complex>): BufferND<Complex> = arg.map { acosh(it) }
override fun atanh(arg: StructureND<Complex>): BufferND<Complex> = arg.map { atanh(it) }
}
/**
* Fast element production using function inlining
*/
public inline fun BufferedFieldND<Complex, ComplexField>.produceInline(initializer: ComplexField.(Int) -> Complex): BufferND<Complex> {
contract { callsInPlace(initializer, InvocationKind.EXACTLY_ONCE) }
val buffer = Buffer.complex(strides.linearSize) { offset -> ComplexField.initializer(offset) }
return BufferND(strides, buffer)
}
public fun AlgebraND.Companion.complex(vararg shape: Int): ComplexFieldND = ComplexFieldND(shape)
/**
* Produce a context for n-dimensional operations inside this real field
*/
public inline fun <R> ComplexField.nd(vararg shape: Int, action: ComplexFieldND.() -> R): R {
contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
return ComplexFieldND(shape).action()
}

View File

@ -0,0 +1,273 @@
/*
* Copyright 2018-2021 KMath contributors.
* Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
*/
package space.kscience.kmath.complex
import space.kscience.kmath.memory.MemoryReader
import space.kscience.kmath.memory.MemorySpec
import space.kscience.kmath.memory.MemoryWriter
import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.operations.*
import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.MemoryBuffer
import space.kscience.kmath.structures.MutableBuffer
import space.kscience.kmath.structures.MutableMemoryBuffer
import kotlin.math.*
/**
* This quaternion's conjugate.
*/
@UnstableKMathAPI
public val DoubleQuaternion.conjugate: DoubleQuaternion
get() = DoubleQuaternionField { z - x * i - y * j - z * k }
/**
* This quaternion's reciprocal.
*/
@UnstableKMathAPI
public val DoubleQuaternion.reciprocal: DoubleQuaternion
get() = DoubleQuaternionField {
val n = norm(this@reciprocal)
return conjugate / (n * n)
}
/**
* Absolute value of the quaternion.
*/
@UnstableKMathAPI
public val DoubleQuaternion.r: Double
get() = sqrt(w * w + x * x + y * y + z * z)
/**
* A field of [DoubleQuaternion].
*/
@UnstableKMathAPI
public object DoubleQuaternionField : Field<DoubleQuaternion>, Norm<DoubleQuaternion, DoubleQuaternion>,
PowerOperations<DoubleQuaternion>,
ExponentialOperations<DoubleQuaternion>, NumbersAddOperations<DoubleQuaternion>, ScaleOperations<DoubleQuaternion> {
override val zero: DoubleQuaternion = DoubleQuaternion(0)
override val one: DoubleQuaternion = DoubleQuaternion(1)
/**
* The `i` quaternion unit.
*/
public val i: DoubleQuaternion = DoubleQuaternion(0, 1)
/**
* The `j` quaternion unit.
*/
public val j: DoubleQuaternion = DoubleQuaternion(0, 0, 1)
/**
* The `k` quaternion unit.
*/
public val k: DoubleQuaternion = DoubleQuaternion(0, 0, 0, 1)
override fun add(a: DoubleQuaternion, b: DoubleQuaternion): DoubleQuaternion =
DoubleQuaternion(a.w + b.w, a.x + b.x, a.y + b.y, a.z + b.z)
override fun scale(a: DoubleQuaternion, value: Double): DoubleQuaternion =
DoubleQuaternion(a.w * value, a.x * value, a.y * value, a.z * value)
override fun multiply(a: DoubleQuaternion, b: DoubleQuaternion): DoubleQuaternion = DoubleQuaternion(
a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z,
a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y,
a.w * b.y - a.x * b.z + a.y * b.w + a.z * b.x,
a.w * b.z + a.x * b.y - a.y * b.x + a.z * b.w,
)
override fun divide(a: DoubleQuaternion, b: DoubleQuaternion): DoubleQuaternion {
val s = b.w * b.w + b.x * b.x + b.y * b.y + b.z * b.z
return DoubleQuaternion(
(b.w * a.w + b.x * a.x + b.y * a.y + b.z * a.z) / s,
(b.w * a.x - b.x * a.w - b.y * a.z + b.z * a.y) / s,
(b.w * a.y + b.x * a.z - b.y * a.w - b.z * a.x) / s,
(b.w * a.z - b.x * a.y + b.y * a.x - b.z * a.w) / s,
)
}
override fun power(arg: DoubleQuaternion, pow: Number): DoubleQuaternion {
if (pow is Int) return pwr(arg, pow)
if (floor(pow.toDouble()) == pow.toDouble()) return pwr(arg, pow.toInt())
return exp(pow * ln(arg))
}
private fun pwr(x: DoubleQuaternion, a: Int): DoubleQuaternion = when {
a < 0 -> -(pwr(x, -a))
a == 0 -> one
a == 1 -> x
a == 2 -> pwr2(x)
a == 3 -> pwr3(x)
a == 4 -> pwr4(x)
else -> {
val x4 = pwr4(x)
var y = x4
repeat((1 until a / 4).count()) { y *= x4 }
if (a % 4 == 3) y *= pwr3(x)
if (a % 4 == 2) y *= pwr2(x)
if (a % 4 == 1) y *= x
y
}
}
private fun pwr2(x: DoubleQuaternion): DoubleQuaternion {
val aa = 2 * x.w
return DoubleQuaternion(x.w * x.w - (x.x * x.x + x.y * x.y + x.z * x.z), aa * x.x, aa * x.y, aa * x.z)
}
private fun pwr3(x: DoubleQuaternion): DoubleQuaternion {
val a2 = x.w * x.w
val n1 = x.x * x.x + x.y * x.y + x.z * x.z
val n2 = 3.0 * a2 - n1
return DoubleQuaternion(x.w * (a2 - 3 * n1), x.x * n2, x.y * n2, x.z * n2)
}
private fun pwr4(x: DoubleQuaternion): DoubleQuaternion {
val a2 = x.w * x.w
val n1 = x.x * x.x + x.y * x.y + x.z * x.z
val n2 = 4 * x.w * (a2 - n1)
return DoubleQuaternion(a2 * a2 - 6 * a2 * n1 + n1 * n1, x.x * n2, x.y * n2, x.z * n2)
}
override fun exp(arg: DoubleQuaternion): DoubleQuaternion {
val un = arg.x * arg.x + arg.y * arg.y + arg.z * arg.z
if (un == 0.0) return DoubleQuaternion(exp(arg.w))
val n1 = sqrt(un)
val ea = exp(arg.w)
val n2 = ea * sin(n1) / n1
return DoubleQuaternion(ea * cos(n1), n2 * arg.x, n2 * arg.y, n2 * arg.z)
}
override fun ln(arg: DoubleQuaternion): DoubleQuaternion {
val nu2 = arg.x * arg.x + arg.y * arg.y + arg.z * arg.z
if (nu2 == 0.0)
return if (arg.w > 0)
DoubleQuaternion(ln(arg.w), 0, 0, 0)
else {
val l = ComplexDoubleField { ln(arg.w) }
DoubleQuaternion(l.re, l.im, 0, 0)
}
val a = arg.w
check(nu2 > 0)
val n = sqrt(a * a + nu2)
val th = acos(a / n) / sqrt(nu2)
return DoubleQuaternion(ln(n), th * arg.x, th * arg.y, th * arg.z)
}
override operator fun Number.plus(b: DoubleQuaternion): DoubleQuaternion =
DoubleQuaternion(toDouble() + b.w, b.x, b.y, b.z)
override operator fun Number.minus(b: DoubleQuaternion): DoubleQuaternion =
DoubleQuaternion(toDouble() - b.w, -b.x, -b.y, -b.z)
override operator fun DoubleQuaternion.plus(b: Number): DoubleQuaternion =
DoubleQuaternion(w + b.toDouble(), x, y, z)
override operator fun DoubleQuaternion.minus(b: Number): DoubleQuaternion =
DoubleQuaternion(w - b.toDouble(), x, y, z)
override operator fun Number.times(b: DoubleQuaternion): DoubleQuaternion =
DoubleQuaternion(toDouble() * b.w, toDouble() * b.x, toDouble() * b.y, toDouble() * b.z)
override fun DoubleQuaternion.unaryMinus(): DoubleQuaternion = DoubleQuaternion(-w, -x, -y, -z)
override fun norm(arg: DoubleQuaternion): DoubleQuaternion = sqrt(arg.conjugate * arg)
override fun bindSymbolOrNull(value: String): DoubleQuaternion? = when (value) {
"i" -> i
"j" -> j
"k" -> k
else -> null
}
override fun number(value: Number): DoubleQuaternion = DoubleQuaternion(value)
override fun sinh(arg: DoubleQuaternion): DoubleQuaternion = (exp(arg) - exp(-arg)) / 2.0
override fun cosh(arg: DoubleQuaternion): DoubleQuaternion = (exp(arg) + exp(-arg)) / 2.0
override fun tanh(arg: DoubleQuaternion): DoubleQuaternion = (exp(arg) - exp(-arg)) / (exp(-arg) + exp(arg))
override fun asinh(arg: DoubleQuaternion): DoubleQuaternion = ln(sqrt(arg * arg + one) + arg)
override fun acosh(arg: DoubleQuaternion): DoubleQuaternion = ln(arg + sqrt((arg - one) * (arg + one)))
override fun atanh(arg: DoubleQuaternion): DoubleQuaternion = (ln(arg + one) - ln(one - arg)) / 2.0
}
/**
* Represents `double`-based quaternion.
*
* @property w The first component.
* @property x The second component.
* @property y The third component.
* @property z The fourth component.
*/
@UnstableKMathAPI
public data class DoubleQuaternion(
val w: Double, val x: Double, val y: Double, val z: Double,
) {
public constructor(w: Number, x: Number, y: Number, z: Number) : this(
w.toDouble(),
x.toDouble(),
y.toDouble(),
z.toDouble(),
)
public constructor(w: Number, x: Number, y: Number) : this(w.toDouble(), x.toDouble(), y.toDouble(), 0.0)
public constructor(w: Number, x: Number) : this(w.toDouble(), x.toDouble(), 0.0, 0.0)
public constructor(w: Number) : this(w.toDouble(), 0.0, 0.0, 0.0)
public constructor(wx: Complex<Number>, yz: Complex<Number>) : this(wx.re, wx.im, yz.re, yz.im)
public constructor(wx: Complex<Number>) : this(wx.re, wx.im, 0, 0)
init {
require(!w.isNaN()) { "w-component of quaternion is not-a-number" }
require(!x.isNaN()) { "x-component of quaternion is not-a-number" }
require(!y.isNaN()) { "x-component of quaternion is not-a-number" }
require(!z.isNaN()) { "x-component of quaternion is not-a-number" }
}
/**
* Returns a string representation of this quaternion.
*/
override fun toString(): String = "($w + $x * i + $y * j + $z * k)"
public companion object : MemorySpec<DoubleQuaternion> {
override val objectSize: Int
get() = 32
override fun MemoryReader.read(offset: Int): DoubleQuaternion =
DoubleQuaternion(readDouble(offset),
readDouble(offset + 8),
readDouble(offset + 16),
readDouble(offset + 24))
override fun MemoryWriter.write(offset: Int, value: DoubleQuaternion) {
writeDouble(offset, value.w)
writeDouble(offset + 8, value.x)
writeDouble(offset + 16, value.y)
writeDouble(offset + 24, value.z)
}
}
}
/**
* Creates a new buffer of quaternions with the specified [size], where each element is calculated by calling the
* specified [init] function.
*/
@UnstableKMathAPI
public inline fun Buffer.Companion.quaternion(size: Int, init: (Int) -> DoubleQuaternion): Buffer<DoubleQuaternion> =
MemoryBuffer.create(DoubleQuaternion, size, init)
/**
* Creates a new buffer of quaternions with the specified [size], where each element is calculated by calling the
* specified [init] function.
*/
@UnstableKMathAPI
public inline fun MutableBuffer.Companion.quaternion(
size: Int,
init: (Int) -> DoubleQuaternion,
): MutableBuffer<DoubleQuaternion> =
MutableMemoryBuffer.create(DoubleQuaternion, size, init)

View File

@ -1,274 +0,0 @@
/*
* Copyright 2018-2021 KMath contributors.
* Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
*/
package space.kscience.kmath.complex
import space.kscience.kmath.memory.MemoryReader
import space.kscience.kmath.memory.MemorySpec
import space.kscience.kmath.memory.MemoryWriter
import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.operations.*
import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.MemoryBuffer
import space.kscience.kmath.structures.MutableBuffer
import space.kscience.kmath.structures.MutableMemoryBuffer
import kotlin.math.*
/**
* This quaternion's conjugate.
*/
public val Quaternion.conjugate: Quaternion
get() = QuaternionField { z - x * i - y * j - z * k }
/**
* This quaternion's reciprocal.
*/
public val Quaternion.reciprocal: Quaternion
get() {
QuaternionField {
val n = norm(this@reciprocal)
return conjugate / (n * n)
}
}
/**
* Absolute value of the quaternion.
*/
public val Quaternion.r: Double
get() = sqrt(w * w + x * x + y * y + z * z)
/**
* A field of [Quaternion].
*/
@OptIn(UnstableKMathAPI::class)
public object QuaternionField : Field<Quaternion>, Norm<Quaternion, Quaternion>, PowerOperations<Quaternion>,
ExponentialOperations<Quaternion>, NumbersAddOperations<Quaternion>, ScaleOperations<Quaternion> {
override val zero: Quaternion = 0.toQuaternion()
override val one: Quaternion = 1.toQuaternion()
/**
* The `i` quaternion unit.
*/
public val i: Quaternion = Quaternion(0, 1)
/**
* The `j` quaternion unit.
*/
public val j: Quaternion = Quaternion(0, 0, 1)
/**
* The `k` quaternion unit.
*/
public val k: Quaternion = Quaternion(0, 0, 0, 1)
override fun add(a: Quaternion, b: Quaternion): Quaternion =
Quaternion(a.w + b.w, a.x + b.x, a.y + b.y, a.z + b.z)
override fun scale(a: Quaternion, value: Double): Quaternion =
Quaternion(a.w * value, a.x * value, a.y * value, a.z * value)
override fun multiply(a: Quaternion, b: Quaternion): Quaternion = Quaternion(
a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z,
a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y,
a.w * b.y - a.x * b.z + a.y * b.w + a.z * b.x,
a.w * b.z + a.x * b.y - a.y * b.x + a.z * b.w,
)
override fun divide(a: Quaternion, b: Quaternion): Quaternion {
val s = b.w * b.w + b.x * b.x + b.y * b.y + b.z * b.z
return Quaternion(
(b.w * a.w + b.x * a.x + b.y * a.y + b.z * a.z) / s,
(b.w * a.x - b.x * a.w - b.y * a.z + b.z * a.y) / s,
(b.w * a.y + b.x * a.z - b.y * a.w - b.z * a.x) / s,
(b.w * a.z - b.x * a.y + b.y * a.x - b.z * a.w) / s,
)
}
override fun power(arg: Quaternion, pow: Number): Quaternion {
if (pow is Int) return pwr(arg, pow)
if (floor(pow.toDouble()) == pow.toDouble()) return pwr(arg, pow.toInt())
return exp(pow * ln(arg))
}
private fun pwr(x: Quaternion, a: Int): Quaternion = when {
a < 0 -> -(pwr(x, -a))
a == 0 -> one
a == 1 -> x
a == 2 -> pwr2(x)
a == 3 -> pwr3(x)
a == 4 -> pwr4(x)
else -> {
val x4 = pwr4(x)
var y = x4
repeat((1 until a / 4).count()) { y *= x4 }
if (a % 4 == 3) y *= pwr3(x)
if (a % 4 == 2) y *= pwr2(x)
if (a % 4 == 1) y *= x
y
}
}
private fun pwr2(x: Quaternion): Quaternion {
val aa = 2 * x.w
return Quaternion(x.w * x.w - (x.x * x.x + x.y * x.y + x.z * x.z), aa * x.x, aa * x.y, aa * x.z)
}
private fun pwr3(x: Quaternion): Quaternion {
val a2 = x.w * x.w
val n1 = x.x * x.x + x.y * x.y + x.z * x.z
val n2 = 3.0 * a2 - n1
return Quaternion(x.w * (a2 - 3 * n1), x.x * n2, x.y * n2, x.z * n2)
}
private fun pwr4(x: Quaternion): Quaternion {
val a2 = x.w * x.w
val n1 = x.x * x.x + x.y * x.y + x.z * x.z
val n2 = 4 * x.w * (a2 - n1)
return Quaternion(a2 * a2 - 6 * a2 * n1 + n1 * n1, x.x * n2, x.y * n2, x.z * n2)
}
override fun exp(arg: Quaternion): Quaternion {
val un = arg.x * arg.x + arg.y * arg.y + arg.z * arg.z
if (un == 0.0) return exp(arg.w).toQuaternion()
val n1 = sqrt(un)
val ea = exp(arg.w)
val n2 = ea * sin(n1) / n1
return Quaternion(ea * cos(n1), n2 * arg.x, n2 * arg.y, n2 * arg.z)
}
override fun ln(arg: Quaternion): Quaternion {
val nu2 = arg.x * arg.x + arg.y * arg.y + arg.z * arg.z
if (nu2 == 0.0)
return if (arg.w > 0)
Quaternion(ln(arg.w), 0, 0, 0)
else {
val l = ComplexField { ln(arg.w.toComplex()) }
Quaternion(l.re, l.im, 0, 0)
}
val a = arg.w
check(nu2 > 0)
val n = sqrt(a * a + nu2)
val th = acos(a / n) / sqrt(nu2)
return Quaternion(ln(n), th * arg.x, th * arg.y, th * arg.z)
}
override operator fun Number.plus(b: Quaternion): Quaternion = Quaternion(toDouble() + b.w, b.x, b.y, b.z)
override operator fun Number.minus(b: Quaternion): Quaternion =
Quaternion(toDouble() - b.w, -b.x, -b.y, -b.z)
override operator fun Quaternion.plus(b: Number): Quaternion = Quaternion(w + b.toDouble(), x, y, z)
override operator fun Quaternion.minus(b: Number): Quaternion = Quaternion(w - b.toDouble(), x, y, z)
override operator fun Number.times(b: Quaternion): Quaternion =
Quaternion(toDouble() * b.w, toDouble() * b.x, toDouble() * b.y, toDouble() * b.z)
override fun Quaternion.unaryMinus(): Quaternion = Quaternion(-w, -x, -y, -z)
override fun norm(arg: Quaternion): Quaternion = sqrt(arg.conjugate * arg)
override fun bindSymbolOrNull(value: String): Quaternion? = when (value) {
"i" -> i
"j" -> j
"k" -> k
else -> null
}
override fun number(value: Number): Quaternion = value.toQuaternion()
override fun sinh(arg: Quaternion): Quaternion = (exp(arg) - exp(-arg)) / 2.0
override fun cosh(arg: Quaternion): Quaternion = (exp(arg) + exp(-arg)) / 2.0
override fun tanh(arg: Quaternion): Quaternion = (exp(arg) - exp(-arg)) / (exp(-arg) + exp(arg))
override fun asinh(arg: Quaternion): Quaternion = ln(sqrt(arg * arg + one) + arg)
override fun acosh(arg: Quaternion): Quaternion = ln(arg + sqrt((arg - one) * (arg + one)))
override fun atanh(arg: Quaternion): Quaternion = (ln(arg + one) - ln(one - arg)) / 2.0
}
/**
* Represents `double`-based quaternion.
*
* @property w The first component.
* @property x The second component.
* @property y The third component.
* @property z The fourth component.
*/
@OptIn(UnstableKMathAPI::class)
public data class Quaternion(
val w: Double, val x: Double, val y: Double, val z: Double,
) {
public constructor(w: Number, x: Number, y: Number, z: Number) : this(
w.toDouble(),
x.toDouble(),
y.toDouble(),
z.toDouble(),
)
public constructor(w: Number, x: Number, y: Number) : this(w.toDouble(), x.toDouble(), y.toDouble(), 0.0)
public constructor(w: Number, x: Number) : this(w.toDouble(), x.toDouble(), 0.0, 0.0)
public constructor(w: Number) : this(w.toDouble(), 0.0, 0.0, 0.0)
public constructor(wx: Complex, yz: Complex) : this(wx.re, wx.im, yz.re, yz.im)
public constructor(wx: Complex) : this(wx.re, wx.im, 0, 0)
init {
require(!w.isNaN()) { "w-component of quaternion is not-a-number" }
require(!x.isNaN()) { "x-component of quaternion is not-a-number" }
require(!y.isNaN()) { "x-component of quaternion is not-a-number" }
require(!z.isNaN()) { "x-component of quaternion is not-a-number" }
}
/**
* Returns a string representation of this quaternion.
*/
override fun toString(): String = "($w + $x * i + $y * j + $z * k)"
public companion object : MemorySpec<Quaternion> {
override val objectSize: Int
get() = 32
override fun MemoryReader.read(offset: Int): Quaternion =
Quaternion(readDouble(offset), readDouble(offset + 8), readDouble(offset + 16), readDouble(offset + 24))
override fun MemoryWriter.write(offset: Int, value: Quaternion) {
writeDouble(offset, value.w)
writeDouble(offset + 8, value.x)
writeDouble(offset + 16, value.y)
writeDouble(offset + 24, value.z)
}
}
}
/**
* Creates a quaternion with real part equal to this real.
*
* @receiver the real part.
* @return a new quaternion.
*/
public fun Number.toQuaternion(): Quaternion = Quaternion(this)
/**
* Creates a quaternion with `w`-component equal to `re`-component of given complex and `x`-component equal to
* `im`-component of given complex.
*
* @receiver the complex number.
* @return a new quaternion.
*/
public fun Complex.toQuaternion(): Quaternion = Quaternion(this)
/**
* Creates a new buffer of quaternions with the specified [size], where each element is calculated by calling the
* specified [init] function.
*/
public inline fun Buffer.Companion.quaternion(size: Int, init: (Int) -> Quaternion): Buffer<Quaternion> =
MemoryBuffer.create(Quaternion, size, init)
/**
* Creates a new buffer of quaternions with the specified [size], where each element is calculated by calling the
* specified [init] function.
*/
public inline fun MutableBuffer.Companion.quaternion(size: Int, init: (Int) -> Quaternion): MutableBuffer<Quaternion> =
MutableMemoryBuffer.create(Quaternion, size, init)

View File

@ -0,0 +1,94 @@
/*
* Copyright 2018-2021 KMath contributors.
* Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
*/
package space.kscience.kmath.complex
import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.BufferFactory
import space.kscience.kmath.structures.MutableBuffer
import space.kscience.kmath.structures.MutableBufferFactory
private class ComplexBuffer<out T : Any>(factory: BufferFactory<T>, override val size: Int, init: (Int) -> Complex<T>) :
Buffer<Complex<T>> {
private val re: Buffer<T>
private val im: Buffer<T>
init {
val tmp = Array(size, init)
re = factory(size) { tmp[it].re }
im = factory(size) { tmp[it].im }
}
override fun get(index: Int): Complex<T> = Complex(re[index], im[index])
override fun iterator(): Iterator<Complex<T>> = object : AbstractIterator<Complex<T>>() {
private val a = re.iterator()
private val b = im.iterator()
override fun computeNext() = if (a.hasNext() && b.hasNext())
setNext(Complex(a.next(), b.next()))
else
done()
}
}
/**
* Creates a new buffer of complex elements with the specified [size], where each element is calculated by calling the
* specified [init] function.
*/
public fun <T : Any> Buffer.Companion.complex(
factory: BufferFactory<T>,
size: Int,
init: (Int) -> Complex<T>,
): Buffer<Complex<T>> = ComplexBuffer(factory, size, init)
private class MutableComplexBuffer<T : Any> private constructor(
override val size: Int,
private val re: MutableBuffer<T>,
private val im: MutableBuffer<T>,
) : MutableBuffer<Complex<T>> {
private constructor(
factory: MutableBufferFactory<T>,
size: Int,
tmp: Array<Complex<T>>,
) : this(size, factory(size) { tmp[it].re }, factory(size) { tmp[it].im })
constructor(
factory: MutableBufferFactory<T>,
size: Int,
init: (Int) -> Complex<T>,
) : this(factory, size, Array(size, init))
override fun get(index: Int): Complex<T> = Complex(re[index], im[index])
override fun set(index: Int, value: Complex<T>) {
re[index] = value.re
im[index] = value.im
}
override fun iterator(): Iterator<Complex<T>> = object : AbstractIterator<Complex<T>>() {
private val a = re.iterator()
private val b = im.iterator()
override fun computeNext() = if (a.hasNext() && b.hasNext())
setNext(Complex(a.next(), b.next()))
else
done()
}
override fun copy(): MutableBuffer<Complex<T>> = MutableComplexBuffer(size, re.copy(), im.copy())
}
/**
* Creates a new buffer of complex elements with the specified [size], where each element is calculated by calling the
* specified [init] function.
*/
public fun <T : Any> MutableBuffer.Companion.complex(
factory: MutableBufferFactory<T>,
size: Int,
init: (Int) -> Complex<T>,
): MutableBuffer<Complex<T>> = MutableComplexBuffer(factory, size, init)

View File

@ -0,0 +1,24 @@
/*
* Copyright 2018-2021 KMath contributors.
* Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
*/
package space.kscience.kmath.complex
import space.kscience.kmath.linear.BufferedLinearSpace
import space.kscience.kmath.operations.ExtendedField
import space.kscience.kmath.operations.NumericAlgebra
import space.kscience.kmath.operations.Ring
import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.BufferFactory
public class ComplexLinearSpace<T : Any, out A>(
elementContext: A,
bufferFactory: BufferFactory<T>,
) : BufferedLinearSpace<Complex<T>, ComplexRing<T, A>>(
ComplexRing(elementContext),
{ size, init -> Buffer.complex(bufferFactory, size, init) },
) where A : Ring<T>, A : NumericAlgebra<T>
public fun <T : Any, A> BufferedLinearSpace<T, A>.complex(): ComplexLinearSpace<T, A> where A : ExtendedField<T>, A : NumericAlgebra<T> =
ComplexLinearSpace(elementAlgebra, bufferFactory)

View File

@ -0,0 +1,68 @@
/*
* Copyright 2018-2021 KMath contributors.
* Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
*/
package space.kscience.kmath.complex
import space.kscience.kmath.nd.BufferedExtendedFieldND
import space.kscience.kmath.nd.BufferedFieldND
import space.kscience.kmath.nd.BufferedGroupND
import space.kscience.kmath.nd.BufferedRingND
import space.kscience.kmath.operations.*
import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.BufferFactory
public open class ComplexGroupND<T : Any, out A>(
shape: IntArray,
elementContext: A,
bufferFactory: BufferFactory<T>,
) : BufferedGroupND<Complex<T>, ComplexGroup<T, A>>(
shape,
ComplexGroup(elementContext),
{ size, init -> Buffer.complex(bufferFactory, size, init) },
) where A : Group<T>, A : NumericAlgebra<T>
public fun <T : Any, A> BufferedGroupND<T, A>.complex(): ComplexGroupND<T, A> where A : Group<T>, A : NumericAlgebra<T> =
ComplexGroupND(shape, elementContext, bufferFactory)
public open class ComplexRingND<T : Any, out A>(
shape: IntArray,
elementContext: A,
bufferFactory: BufferFactory<T>,
) : BufferedRingND<Complex<T>, ComplexRing<T, A>>(
shape,
ComplexRing(elementContext),
{ size, init -> Buffer.complex(bufferFactory, size, init) },
) where A : Ring<T>, A : NumericAlgebra<T>
public fun <T : Any, A> BufferedRingND<T, A>.complex(): ComplexRingND<T, A> where A : Ring<T>, A : NumericAlgebra<T> =
ComplexRingND(shape, elementContext, bufferFactory)
public open class ComplexFieldND<T : Any, out A>(
shape: IntArray,
elementContext: A,
bufferFactory: BufferFactory<T>,
) : BufferedFieldND<Complex<T>, ComplexField<T, A>>(
shape,
ComplexField(elementContext),
{ size, init -> Buffer.complex(bufferFactory, size, init) },
) where A : Field<T>, A : NumericAlgebra<T>
public fun <T : Any, A> BufferedFieldND<T, A>.complex(): ComplexFieldND<T, A> where A : Field<T>, A : NumericAlgebra<T> =
ComplexFieldND(shape, elementContext, bufferFactory)
public open class ComplexExtendedFieldND<T : Any, out A>(
shape: IntArray,
elementContext: A,
bufferFactory: BufferFactory<T>,
) : BufferedExtendedFieldND<Complex<T>, ComplexExtendedField<T, A>>(
shape,
ComplexExtendedField(elementContext),
{ size, init -> Buffer.complex(bufferFactory, size, init) },
) where A : ExtendedField<T>, A : NumericAlgebra<T>
public fun <T : Any, A> BufferedExtendedFieldND<T, A>.complex(): ComplexExtendedFieldND<T, A> where A : ExtendedField<T>, A : NumericAlgebra<T> =
ComplexExtendedFieldND(shape, elementContext, bufferFactory)

View File

@ -6,13 +6,14 @@
package space.kscience.kmath.complex package space.kscience.kmath.complex
import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.MutableBuffer
import kotlin.test.Test import kotlin.test.Test
import kotlin.test.assertEquals import kotlin.test.assertEquals
class ComplexBufferSpecTest { class ComplexBufferSpecTest {
@Test @Test
fun testComplexBuffer() { fun testComplexBuffer() {
val buffer = Buffer.complex(20) { Complex(it.toDouble(), -it.toDouble()) } val buffer = Buffer.complex(MutableBuffer.Companion::double, 20) { Complex(it.toDouble(), -it.toDouble()) }
assertEquals(Complex(5.0, -5.0), buffer[5]) assertEquals(Complex(5.0, -5.0), buffer[5])
} }
} }

View File

@ -5,12 +5,11 @@
package space.kscience.kmath.complex package space.kscience.kmath.complex
import space.kscience.kmath.operations.DoubleField
import space.kscience.kmath.operations.invoke import space.kscience.kmath.operations.invoke
import kotlin.math.PI import space.kscience.kmath.operations.pi
import kotlin.math.abs
import kotlin.test.Test import kotlin.test.Test
import kotlin.test.assertEquals import kotlin.test.assertEquals
import kotlin.test.assertTrue
internal class ComplexFieldTest { internal class ComplexFieldTest {
// TODO make verifier classes available in this source set // TODO make verifier classes available in this source set
@ -19,63 +18,84 @@ internal class ComplexFieldTest {
@Test @Test
fun testAddition() { fun testAddition() {
assertEquals(Complex(42, 42), ComplexField { Complex(16, 16) + Complex(26, 26) }) assertEquals(Complex(42, 42), ComplexIntRing { Complex(16, 16) + Complex(26, 26) })
assertEquals(Complex(42, 16), ComplexField { Complex(16, 16) + 26 }) assertEquals(Complex(42, 16), ComplexIntRing { Complex(16, 16) + 26 })
assertEquals(Complex(42, 16), ComplexField { 26 + Complex(16, 16) }) assertEquals(Complex(42, 16), ComplexIntRing { 26 + Complex(16, 16) })
} }
@Test @Test
fun testSubtraction() { fun testSubtraction() {
assertEquals(Complex(42, 42), ComplexField { Complex(86, 55) - Complex(44, 13) }) assertEquals(Complex(42, 42), ComplexIntRing { Complex(86, 55) - Complex(44, 13) })
assertEquals(Complex(42, 56), ComplexField { Complex(86, 56) - 44 }) assertEquals(Complex(42, 56), ComplexIntRing { Complex(86, 56) - 44 })
assertEquals(Complex(42, 56), ComplexField { 86 - Complex(44, -56) }) assertEquals(Complex(42, 56), ComplexIntRing { 86 - Complex(44, -56) })
} }
@Test @Test
fun testMultiplication() { fun testMultiplication() {
assertEquals(Complex(42, 42), ComplexField { Complex(4.2, 0) * Complex(10, 10) }) assertEquals(Complex(42.0, 42.0), ComplexDoubleField { Complex(4.2, 0.0) * Complex(10.0, 10.0) })
assertEquals(Complex(42, 21), ComplexField { Complex(4.2, 2.1) * 10 }) assertEquals(Complex(42.0, 21.0), ComplexDoubleField { Complex(4.2, 2.1) * 10 })
assertEquals(Complex(42, 21), ComplexField { 10 * Complex(4.2, 2.1) }) assertEquals(Complex(42.0, 21.0), ComplexDoubleField { 10 * Complex(4.2, 2.1) })
} }
@Test @Test
fun testDivision() { fun testDivision() {
assertEquals(Complex(42, 42), ComplexField { Complex(0, 168) / Complex(2, 2) }) assertEquals(Complex(42.0, 42.0), ComplexDoubleField { Complex(0.0, 168.0) / Complex(2.0, 2.0) })
assertEquals(Complex(42, 56), ComplexField { Complex(86, 56) - 44 }) assertEquals(Complex(42.0, 56.0), ComplexDoubleField { Complex(86.0, 56.0) - 44.0 })
assertEquals(Complex(42, 56), ComplexField { 86 - Complex(44, -56) }) assertEquals(Complex(42.0, 56.0), ComplexDoubleField { 86.0 - Complex(44.0, -56.0) })
} }
@Test @Test
fun testSine() { fun testSine() {
assertEquals(ComplexField { i * sinh(one) }, ComplexField { sin(i) }) assertEquals(ComplexDoubleField { i * sinh(one) }, ComplexDoubleField { sin(i) })
assertEquals(ComplexField { i * sinh(PI.toComplex()) }, ComplexField { sin(i * PI.toComplex()) }) assertEquals(ComplexDoubleField { i * sinh(pi) }, ComplexDoubleField { sin(i * pi) })
} }
@Test @Test
fun testInverseSine() { fun testInverseSine() = ComplexDoubleField {
assertEquals(Complex(0, -0.0), ComplexField { asin(zero) }) assertEquals(norm(zero), norm(asin(zero)), 1e-10)
assertTrue(abs(ComplexField { i * asinh(one) }.r - ComplexField { asin(i) }.r) < 0.000000000000001) assertEquals(norm(i * asinh(one)), norm(i * asinh(one)), 1e-10)
} }
@Test @Test
fun testInverseHyperbolicSine() { fun testInverseHyperbolicSine() {
assertEquals( assertEquals(ComplexDoubleField { i * pi / 2 }, ComplexDoubleField { asinh(i) })
ComplexField { i * PI.toComplex() / 2 },
ComplexField { asinh(i) })
} }
@Test @Test
fun testPower() { fun testPower() {
assertEquals(ComplexField.zero, ComplexField { zero pow 2 }) assertEquals(ComplexDoubleField.zero, ComplexDoubleField { zero pow 2 })
assertEquals(ComplexField.zero, ComplexField { zero pow 2 }) assertEquals(ComplexDoubleField.zero, ComplexDoubleField { zero pow 2 })
assertEquals( assertEquals(
ComplexField { i * 8 }.let { it.im.toInt() to it.re.toInt() }, ComplexDoubleField { i * 8 }.let { it.im.toInt() to it.re.toInt() },
ComplexField { Complex(2, 2) pow 2 }.let { it.im.toInt() to it.re.toInt() }) ComplexDoubleField { Complex(2.0, 2.0) pow 2 }.let { it.im.toInt() to it.re.toInt() },
)
} }
@Test @Test
fun testNorm() { fun testNorm() {
assertEquals(2.toComplex(), ComplexField { norm(2 * i) }) assertEquals(2.0, ComplexDoubleField { norm(number(2.0) * i) })
}
@Test
fun conjugate() = ComplexDoubleField {
assertEquals(Complex(0.0, 42.0), Complex(0.0, -42.0).conjugate)
}
@Test
fun reciprocal() = ComplexDoubleField {
assertEquals(norm(Complex(0.5, 0.0)), norm(Complex(2.0, 0.0).reciprocal), 1e-10)
}
@Test
fun polar() {
val num = Complex(0.5, 2.5)
val theta = ComplexDoubleField { num.theta }
val r = ComplexDoubleField { num.r }
DoubleField {
assertEquals(cos(theta), num.re / r, 1e-10)
assertEquals(sin(theta), num.im / r, 1e-10)
}
} }
} }

View File

@ -1,36 +0,0 @@
/*
* Copyright 2018-2021 KMath contributors.
* Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
*/
package space.kscience.kmath.complex
import space.kscience.kmath.operations.invoke
import kotlin.math.sqrt
import kotlin.test.Test
import kotlin.test.assertEquals
import kotlin.test.assertTrue
internal class ComplexTest {
@Test
fun conjugate() = ComplexField { assertEquals(Complex(0, 42), Complex(0, -42).conjugate) }
@Test
fun reciprocal() = ComplexField { assertTrue((Complex(0.5, -0.0) - 2.toComplex().reciprocal).r < 1e-10) }
@Test
fun r() = ComplexField { assertEquals(sqrt(2.0), (i + 1.0.toComplex()).r) }
@Test
fun theta() = assertEquals(0.0, 1.toComplex().theta)
@Test
fun toComplex() {
assertEquals(Complex(42), 42.toComplex())
assertEquals(Complex(42.0), 42.0.toComplex())
assertEquals(Complex(42f), 42f.toComplex())
assertEquals(Complex(42.0), 42.0.toComplex())
assertEquals(Complex(42.toByte()), 42.toByte().toComplex())
assertEquals(Complex(42.toShort()), 42.toShort().toComplex())
}
}

View File

@ -17,11 +17,11 @@ internal class ExpressionFieldForComplexTest {
@Test @Test
fun testComplex() { fun testComplex() {
val expression = FunctionalExpressionField(ComplexField).run { val expression = FunctionalExpressionField(ComplexDoubleField).run {
val x = bindSymbol(x) val x = bindSymbol(x)
x * x + 2 * x + one x * x + 2 * x + one
} }
assertEquals(expression(x to Complex(1.0, 0.0)), Complex(4.0, 0.0)) assertEquals(Complex(4.0, 0.0), expression(x to Complex(1.0, 0.0)))
} }
} }

View File

@ -4,47 +4,47 @@
*/ */
package space.kscience.kmath.complex package space.kscience.kmath.complex
//
import space.kscience.kmath.operations.invoke //import space.kscience.kmath.operations.invoke
import kotlin.test.Test //import kotlin.test.Test
import kotlin.test.assertEquals //import kotlin.test.assertEquals
//
internal class QuaternionFieldTest { //internal class QuaternionFieldTest {
@Test
fun testAddition() {
assertEquals(Quaternion(42, 42), QuaternionField { Quaternion(16, 16) + Quaternion(26, 26) })
assertEquals(Quaternion(42, 16), QuaternionField { Quaternion(16, 16) + 26 })
assertEquals(Quaternion(42, 16), QuaternionField { 26 + Quaternion(16, 16) })
}
// @Test // @Test
// fun testSubtraction() { // fun testAddition() {
// assertEquals(Quaternion(42, 42), QuaternionField { Quaternion(86, 55) - Quaternion(44, 13) }) // assertEquals(Quaternion(42, 42), QuaternionField { Quaternion(16, 16) + Quaternion(26, 26) })
// assertEquals(Quaternion(42, 56), QuaternionField { Quaternion(86, 56) - 44 }) // assertEquals(Quaternion(42, 16), QuaternionField { Quaternion(16, 16) + 26 })
// assertEquals(Quaternion(42, 56), QuaternionField { 86 - Quaternion(44, -56) }) // assertEquals(Quaternion(42, 16), QuaternionField { 26 + Quaternion(16, 16) })
// } // }
//
@Test //// @Test
fun testMultiplication() { //// fun testSubtraction() {
assertEquals(Quaternion(42, 42), QuaternionField { Quaternion(4.2, 0) * Quaternion(10, 10) }) //// assertEquals(Quaternion(42, 42), QuaternionField { Quaternion(86, 55) - Quaternion(44, 13) })
assertEquals(Quaternion(42, 21), QuaternionField { Quaternion(4.2, 2.1) * 10 }) //// assertEquals(Quaternion(42, 56), QuaternionField { Quaternion(86, 56) - 44 })
assertEquals(Quaternion(42, 21), QuaternionField { 10 * Quaternion(4.2, 2.1) }) //// assertEquals(Quaternion(42, 56), QuaternionField { 86 - Quaternion(44, -56) })
} //// }
//
// @Test // @Test
// fun testDivision() { // fun testMultiplication() {
// assertEquals(Quaternion(42, 42), QuaternionField { Quaternion(0, 168) / Quaternion(2, 2) }) // assertEquals(Quaternion(42, 42), QuaternionField { Quaternion(4.2, 0) * Quaternion(10, 10) })
// assertEquals(Quaternion(42, 56), QuaternionField { Quaternion(86, 56) - 44 }) // assertEquals(Quaternion(42, 21), QuaternionField { Quaternion(4.2, 2.1) * 10 })
// assertEquals(Quaternion(42, 56) , QuaternionField { 86 - Quaternion(44, -56) }) // assertEquals(Quaternion(42, 21), QuaternionField { 10 * Quaternion(4.2, 2.1) })
// } // }
//
@Test //// @Test
fun testPower() { //// fun testDivision() {
assertEquals(QuaternionField.zero, QuaternionField { zero pow 2 }) //// assertEquals(Quaternion(42, 42), QuaternionField { Quaternion(0, 168) / Quaternion(2, 2) })
assertEquals(QuaternionField.zero, QuaternionField { zero pow 2 }) //// assertEquals(Quaternion(42, 56), QuaternionField { Quaternion(86, 56) - 44 })
//// assertEquals(Quaternion(42, 56) , QuaternionField { 86 - Quaternion(44, -56) })
assertEquals( //// }
QuaternionField { i * 8 }.let { it.x.toInt() to it.w.toInt() }, //
QuaternionField { Quaternion(2, 2) pow 2 }.let { it.x.toInt() to it.w.toInt() }) // @Test
} // fun testPower() {
} // assertEquals(QuaternionField.zero, QuaternionField { zero pow 2 })
// assertEquals(QuaternionField.zero, QuaternionField { zero pow 2 })
//
// assertEquals(
// QuaternionField { i * 8 }.let { it.x.toInt() to it.w.toInt() },
// QuaternionField { Quaternion(2, 2) pow 2 }.let { it.x.toInt() to it.w.toInt() })
// }
//}

View File

@ -0,0 +1,22 @@
/*
* Copyright 2018-2021 KMath contributors.
* Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
*/
package space.kscience.kmath.complex
import space.kscience.kmath.operations.JBigDecimalField
import space.kscience.kmath.operations.JBigIntegerRing
import java.math.BigDecimal
import java.math.BigInteger
/**
* [ComplexRing] instance for [JBigIntegerRing].
*/
public val ComplexJBigIntegerRing: ComplexRing<BigInteger, JBigIntegerRing> = ComplexRing(JBigIntegerRing)
/**
* [ComplexRing] instance for [JBigDecimalField].
*/
public val ComplexJBigDecimalField: ComplexField<BigDecimal, JBigDecimalField.Companion> =
ComplexField(JBigDecimalField)

View File

@ -15,9 +15,9 @@ import space.kscience.kmath.structures.VirtualBuffer
import space.kscience.kmath.structures.indices import space.kscience.kmath.structures.indices
public class BufferedLinearSpace<T : Any, out A : Ring<T>>( public open class BufferedLinearSpace<T : Any, out A : Ring<T>>(
override val elementAlgebra: A, override val elementAlgebra: A,
private val bufferFactory: BufferFactory<T>, public val bufferFactory: BufferFactory<T>,
) : LinearSpace<T, A> { ) : LinearSpace<T, A> {
private fun ndRing( private fun ndRing(

View File

@ -286,3 +286,11 @@ public interface FieldND<T, out F : Field<T>> : Field<StructureND<T>>, RingND<T,
// } // }
// } // }
} }
/**
* Extended field of [StructureND].
*
* @param T the type of the element contained in ND structure.
* @param F the type of extended field of structure elements.
*/
public interface ExtendedFieldND<T, out F : ExtendedField<T>> : ExtendedField<StructureND<T>>, FieldND<T, F>

View File

@ -57,11 +57,11 @@ public interface BufferAlgebraND<T, out A : Algebra<T>> : AlgebraND<T, A> {
} }
} }
public open class BufferedGroupND<T, out A : Group<T>>( public open class BufferedGroupND<T, out G : Group<T>>(
final override val shape: IntArray, final override val shape: IntArray,
final override val elementContext: A, final override val elementContext: G,
final override val bufferFactory: BufferFactory<T>, final override val bufferFactory: BufferFactory<T>,
) : GroupND<T, A>, BufferAlgebraND<T, A> { ) : GroupND<T, G>, BufferAlgebraND<T, G> {
override val strides: Strides = DefaultStrides(shape) override val strides: Strides = DefaultStrides(shape)
override val zero: BufferND<T> by lazy { produce { zero } } override val zero: BufferND<T> by lazy { produce { zero } }
override fun StructureND<T>.unaryMinus(): StructureND<T> = produce { -get(it) } override fun StructureND<T>.unaryMinus(): StructureND<T> = produce { -get(it) }
@ -75,15 +75,32 @@ public open class BufferedRingND<T, out R : Ring<T>>(
override val one: BufferND<T> by lazy { produce { one } } override val one: BufferND<T> by lazy { produce { one } }
} }
public open class BufferedFieldND<T, out R : Field<T>>( public open class BufferedFieldND<T, out F : Field<T>>(
shape: IntArray, shape: IntArray,
elementContext: R, elementContext: F,
bufferFactory: BufferFactory<T>, bufferFactory: BufferFactory<T>,
) : BufferedRingND<T, R>(shape, elementContext, bufferFactory), FieldND<T, R> { ) : BufferedRingND<T, F>(shape, elementContext, bufferFactory), FieldND<T, F> {
override fun scale(a: StructureND<T>, value: Double): StructureND<T> = a.map { it * value } override fun scale(a: StructureND<T>, value: Double): StructureND<T> = a.map { it * value }
} }
public open class BufferedExtendedFieldND<T, out F : ExtendedField<T>>(
shape: IntArray,
elementContext: F,
bufferFactory: BufferFactory<T>,
) : BufferedFieldND<T, F>(shape, elementContext, bufferFactory), ExtendedFieldND<T, F> {
public override fun sin(arg: StructureND<T>): StructureND<T> = arg.map { elementContext.sin(it) }
public override fun cos(arg: StructureND<T>): StructureND<T> = arg.map { elementContext.cos(it) }
public override fun asin(arg: StructureND<T>): StructureND<T> = arg.map { elementContext.asin(it) }
public override fun acos(arg: StructureND<T>): StructureND<T> = arg.map { elementContext.acos(it) }
public override fun atan(arg: StructureND<T>): StructureND<T> = arg.map { elementContext.atan(it) }
public override fun power(arg: StructureND<T>, pow: Number): StructureND<T> =
arg.map { elementContext.power(it, pow) }
public override fun exp(arg: StructureND<T>): StructureND<T> = arg.map { elementContext.exp(it) }
public override fun ln(arg: StructureND<T>): StructureND<T> = arg.map { elementContext.ln(it) }
}
// group factories // group factories
public fun <T, A : Ring<T>> AlgebraND.Companion.group( public fun <T, A : Ring<T>> AlgebraND.Companion.group(
space: A, space: A,
@ -140,3 +157,27 @@ public inline fun <T, A : Field<T>, R> A.ndField(
contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) } contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
return AlgebraND.field(this, bufferFactory, *shape).run(action) return AlgebraND.field(this, bufferFactory, *shape).run(action)
} }
public fun <T, A : ExtendedField<T>> AlgebraND.Companion.extendedField(
field: A,
bufferFactory: BufferFactory<T>,
vararg shape: Int,
): BufferedExtendedFieldND<T, A> = BufferedExtendedFieldND(shape, field, bufferFactory)
@Suppress("UNCHECKED_CAST")
public inline fun <reified T : Any, A : ExtendedField<T>> AlgebraND.Companion.auto(
field: A,
vararg shape: Int,
): FieldND<T, A> = when (field) {
DoubleField -> DoubleFieldND(shape) as FieldND<T, A>
else -> BufferedFieldND(shape, field, Buffer.Companion::auto)
}
public inline fun <T, A : ExtendedField<T>, R> A.ndExtendedField(
noinline bufferFactory: BufferFactory<T>,
vararg shape: Int,
action: BufferedExtendedFieldND<T, A>.() -> R,
): R {
contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
return AlgebraND.extendedField(this, bufferFactory, *shape).run(action)
}

View File

@ -7,9 +7,7 @@ package space.kscience.kmath.nd
import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.operations.DoubleField
import space.kscience.kmath.operations.ExtendedField
import space.kscience.kmath.operations.NumbersAddOperations import space.kscience.kmath.operations.NumbersAddOperations
import space.kscience.kmath.operations.ScaleOperations
import space.kscience.kmath.structures.DoubleBuffer import space.kscience.kmath.structures.DoubleBuffer
import kotlin.contracts.InvocationKind import kotlin.contracts.InvocationKind
import kotlin.contracts.contract import kotlin.contracts.contract
@ -17,10 +15,8 @@ import kotlin.contracts.contract
@OptIn(UnstableKMathAPI::class) @OptIn(UnstableKMathAPI::class)
public class DoubleFieldND( public class DoubleFieldND(
shape: IntArray, shape: IntArray,
) : BufferedFieldND<Double, DoubleField>(shape, DoubleField, ::DoubleBuffer), ) : BufferedExtendedFieldND<Double, DoubleField>(shape, DoubleField, ::DoubleBuffer),
NumbersAddOperations<StructureND<Double>>, NumbersAddOperations<StructureND<Double>> {
ScaleOperations<StructureND<Double>>,
ExtendedField<StructureND<Double>> {
override val zero: BufferND<Double> by lazy { produce { zero } } override val zero: BufferND<Double> by lazy { produce { zero } }
override val one: BufferND<Double> by lazy { produce { one } } override val one: BufferND<Double> by lazy { produce { one } }
@ -103,7 +99,7 @@ public class DoubleFieldND(
override fun atanh(arg: StructureND<Double>): BufferND<Double> = arg.map { atanh(it) } override fun atanh(arg: StructureND<Double>): BufferND<Double> = arg.map { atanh(it) }
} }
public fun AlgebraND.Companion.real(vararg shape: Int): DoubleFieldND = DoubleFieldND(shape) public fun AlgebraND.Companion.double(vararg shape: Int): DoubleFieldND = DoubleFieldND(shape)
/** /**
* Produce a context for n-dimensional operations inside this real field * Produce a context for n-dimensional operations inside this real field

View File

@ -7,7 +7,7 @@ package space.kscience.kmath.structures
import space.kscience.kmath.nd.AlgebraND import space.kscience.kmath.nd.AlgebraND
import space.kscience.kmath.nd.get import space.kscience.kmath.nd.get
import space.kscience.kmath.nd.real import space.kscience.kmath.nd.double
import space.kscience.kmath.operations.invoke import space.kscience.kmath.operations.invoke
import space.kscience.kmath.testutils.FieldVerifier import space.kscience.kmath.testutils.FieldVerifier
import kotlin.test.Test import kotlin.test.Test
@ -16,12 +16,12 @@ import kotlin.test.assertEquals
internal class NDFieldTest { internal class NDFieldTest {
@Test @Test
fun verify() { fun verify() {
(AlgebraND.real(12, 32)) { FieldVerifier(this, one + 3, one - 23, one * 12, 6.66) } (AlgebraND.double(12, 32)) { FieldVerifier(this, one + 3, one - 23, one * 12, 6.66) }
} }
@Test @Test
fun testStrides() { fun testStrides() {
val ndArray = AlgebraND.real(10, 10).produce { (it[0] + it[1]).toDouble() } val ndArray = AlgebraND.double(10, 10).produce { (it[0] + it[1]).toDouble() }
assertEquals(ndArray[5, 5], 10.0) assertEquals(ndArray[5, 5], 10.0)
} }
} }

View File

@ -17,7 +17,7 @@ import kotlin.test.assertEquals
@Suppress("UNUSED_VARIABLE") @Suppress("UNUSED_VARIABLE")
class NumberNDFieldTest { class NumberNDFieldTest {
val algebra = AlgebraND.real(3, 3) val algebra = AlgebraND.double(3, 3)
val array1 = algebra.produce { (i, j) -> (i + j).toDouble() } val array1 = algebra.produce { (i, j) -> (i + j).toDouble() }
val array2 = algebra.produce { (i, j) -> (i - j).toDouble() } val array2 = algebra.produce { (i, j) -> (i - j).toDouble() }
@ -83,7 +83,7 @@ class NumberNDFieldTest {
@Test @Test
fun testInternalContext() { fun testInternalContext() {
algebra { algebra {
(AlgebraND.real(*array1.shape)) { with(L2Norm) { 1 + norm(array1) + exp(array2) } } (AlgebraND.double(*array1.shape)) { with(L2Norm) { 1 + norm(array1) + exp(array2) } }
} }
} }
} }

View File

@ -10,9 +10,9 @@ import java.math.BigInteger
import java.math.MathContext import java.math.MathContext
/** /**
* A field over [BigInteger]. * A ring over [BigInteger].
*/ */
public object JBigIntegerField : Ring<BigInteger>, NumericAlgebra<BigInteger> { public object JBigIntegerRing : Ring<BigInteger>, NumericAlgebra<BigInteger> {
override val zero: BigInteger get() = BigInteger.ZERO override val zero: BigInteger get() = BigInteger.ZERO
override val one: BigInteger get() = BigInteger.ONE override val one: BigInteger get() = BigInteger.ONE

View File

@ -28,7 +28,7 @@ public class DoubleHistogramSpace(
public val dimension: Int get() = lower.size public val dimension: Int get() = lower.size
private val shape = IntArray(binNums.size) { binNums[it] + 2 } private val shape = IntArray(binNums.size) { binNums[it] + 2 }
override val histogramValueSpace: DoubleFieldND = AlgebraND.real(*shape) override val histogramValueSpace: DoubleFieldND = AlgebraND.double(*shape)
override val strides: Strides get() = histogramValueSpace.strides override val strides: Strides get() = histogramValueSpace.strides
private val binSize = DoubleBuffer(dimension) { (upper[it] - lower[it]) / binNums[it] } private val binSize = DoubleBuffer(dimension) { (upper[it] - lower[it]) / binNums[it] }

View File

@ -1,5 +1,6 @@
plugins { plugins {
id("ru.mipt.npm.gradle.jvm") kotlin("jvm")
id("ru.mipt.npm.gradle.common")
} }
description = "Jafama integration module" description = "Jafama integration module"

View File

@ -1,6 +1,7 @@
plugins { plugins {
id("ru.mipt.npm.gradle.jvm") kotlin("jvm")
kotlin("jupyter.api") kotlin("jupyter.api")
id("ru.mipt.npm.gradle.common")
} }
dependencies { dependencies {
@ -9,13 +10,8 @@ dependencies {
api(project(":kmath-for-real")) api(project(":kmath-for-real"))
} }
kscience{ kscience.useHtml()
useHtml() readme.maturity = ru.mipt.npm.gradle.Maturity.PROTOTYPE
}
readme {
maturity = ru.mipt.npm.gradle.Maturity.PROTOTYPE
}
kotlin.sourceSets.all { kotlin.sourceSets.all {
languageSettings.useExperimentalAnnotation("space.kscience.kmath.misc.UnstableKMathAPI") languageSettings.useExperimentalAnnotation("space.kscience.kmath.misc.UnstableKMathAPI")

View File

@ -12,7 +12,7 @@ import space.kscience.kmath.ast.rendering.FeaturedMathRendererWithPostProcess
import space.kscience.kmath.ast.rendering.MathMLSyntaxRenderer import space.kscience.kmath.ast.rendering.MathMLSyntaxRenderer
import space.kscience.kmath.ast.rendering.renderWithStringBuilder import space.kscience.kmath.ast.rendering.renderWithStringBuilder
import space.kscience.kmath.complex.Complex import space.kscience.kmath.complex.Complex
import space.kscience.kmath.complex.Quaternion import space.kscience.kmath.complex.DoubleQuaternion
import space.kscience.kmath.expressions.MST import space.kscience.kmath.expressions.MST
import space.kscience.kmath.expressions.MstRing import space.kscience.kmath.expressions.MstRing
import space.kscience.kmath.misc.PerformancePitfall import space.kscience.kmath.misc.PerformancePitfall
@ -121,18 +121,16 @@ internal class KMathJupyter : JupyterIntegration() {
}) })
} }
render<Complex> { render<Complex<Number>> {
MstRing { MstRing { number(it.re) + number(it.im) * bindSymbol("i") }.toDisplayResult()
number(it.re) + number(it.im) * bindSymbol("i")
}.toDisplayResult()
} }
render<Quaternion> { render<DoubleQuaternion> {
MstRing { MstRing {
number(it.w) + number(it.w) +
number(it.x) * bindSymbol("i") + number(it.x) * bindSymbol("i") +
number(it.x) * bindSymbol("j") + number(it.y) * bindSymbol("j") +
number(it.x) * bindSymbol("k") number(it.z) * bindSymbol("k")
}.toDisplayResult() }.toDisplayResult()
} }
} }

View File

@ -5,18 +5,14 @@ pluginManagement {
gradlePluginPortal() gradlePluginPortal()
} }
val toolsVersion = "0.10.2"
val kotlinVersion = "1.5.21" val kotlinVersion = "1.5.21"
plugins { plugins {
id("ru.mipt.npm.gradle.project") version toolsVersion id("ru.mipt.npm.gradle.project") version "0.10.2"
id("ru.mipt.npm.gradle.mpp") version toolsVersion
id("ru.mipt.npm.gradle.jvm") version toolsVersion
kotlin("multiplatform") version kotlinVersion kotlin("multiplatform") version kotlinVersion
kotlin("jvm") version kotlinVersion
kotlin("plugin.allopen") version kotlinVersion kotlin("plugin.allopen") version kotlinVersion
id("org.jetbrains.kotlinx.benchmark") version "0.3.1" id("org.jetbrains.kotlinx.benchmark") version "0.3.1"
kotlin("jupyter.api") version "0.10.0-131-1" kotlin("jupyter.api") version "0.10.0-174"
} }
} }