Approaches to the implementation of generalized complex numbers in the Julia language
AApproaches to the implementation of generalized complex numbers in the Julialanguage
Migran N. Gevorkyan, * Anna V. Korolkova, † and Dmitry S. Kulyabov
1, 2, ‡ Department of Applied Probability and Informatics,Peoples’ Friendship University of Russia (RUDN University),6 Miklukho-Maklaya St, Moscow, 117198, Russian Federation Laboratory of Information TechnologiesJoint Institute for Nuclear Research6 Joliot-Curie, Dubna, Moscow region, 141980, Russian Federation
In problems of mathematical physics, to study the structures of spaces using the Cayley–Kleinmodels in theoretical calculations, the use of generalized complex numbers is required. In thecase of computational experiments, such tasks require their high-quality implementation in aprogramming language. The proposed small implementation of generalized complex numbers inmodern programming languages have several disadvantages. In this article, we propose using theJulia language as the language for implementing generalized complex numbers, not least because itsupports the multiple dispatch mechanism.The paper demonstrates the approach to the implementation of one of the types of generalizedcomplex numbers, namely dual numbers. We place particular emphasis on the description of the use ofthe multiple dispatch mechanism to implement a new numerical type. The resulting implementationof dual numbers can be considered as a prototype for a complete software module for supportinggeneralized complex numbers.
Keywords: complex numbers, parabolic complex numbers, dual numbers, multiple dispatch, Julia * [email protected] † [email protected] ‡ [email protected] a r X i v : . [ c s . M S ] J u l I. INTRODUCTION
There is an approach that allows one to generalize the complex numbers and get three different classes of generalizedcomplex numbers [1–4]. In this approach, we set square equation 𝑧 + 𝑝𝑧 + 𝑞 = 0 . with a determinant Δ = 𝑝 − 𝑝𝑞. Depending on the sign of the determinant, one gets the following results are obtained complex numbers:• Δ < —elliptic (normal complex numbers);• Δ = 0 —parabolic (dual complex numbers);• Δ > —hyperbolic (split complex numbers or double numbers).Initially, these systems of complex numbers were introduced to describe Cayley–Klein models [5, 6]. However, theyhave other uses. In particular, dual complex numbers can be used for problems of automatic differentiation [7].Existing implementations of generalized complex numbers are designed for specific applications. For example, dualnumbers are used for automatic differentiation. The goal of this work is to create a pure implementation of dualnumbers that is as close as possible to their mathematical definition. We use Julia programming language for this task. A. Article structure
The II section provides sufficient description of dual complex numbers. We use constructive definition, affecting onlyoperations with dual complex numbers. We try to avoid unnecessary mathematical abstractions. In the section III wegive a General idea of the multiple dispatch mechanism — the fundamental idea of Julia language. Also, the multipledispatch is the basis for implementation of any user-defined type in Julia. In the section IV we describe in detail ourimplementation of dual complex numbers in Julia.
B. Notations and conventions
To denote a dual complex unit, we will use the symbol ε . II. DUAL NUMBERS
The dual number 𝑧 is defined algebraically as follows: 𝑧 = 𝑎 + ε 𝑏, 𝑎, 𝑏 ∈ R , ε = 0 , ε ̸ = 0 . The value of 𝑎 will be called the real part, and 𝑏 — imaginary. There is no confusion in the terminology, because inthis article we not use the ordinary complex numbers. A. The algebraic form
Algebraic properties for two numbers 𝑧 = 𝑎 + ε 𝑏 и 𝑧 = 𝑎 + ε 𝑏 . Addition: 𝑧 + 𝑧 = ( 𝑎 + 𝑎 ) + ε ( 𝑏 + 𝑏 ) . Subtraction: 𝑧 − 𝑧 = ( 𝑎 − 𝑎 ) + ε ( 𝑏 − 𝑏 ) . Multiplication: 𝑧 · 𝑧 = ( 𝑎 + ε 𝑏 ) · ( 𝑎 + ε 𝑏 ) = 𝑎 𝑎 + ε ( 𝑎 𝑏 + 𝑏 𝑎 ) . Conjunction: ¯ 𝑧 = 𝑎 − ε 𝑏 , 𝑧 ¯ 𝑧 = ( 𝑎 + ε 𝑏 ) · ( 𝑎 − ε 𝑏 ) = 𝑎 . Absolute value: | 𝑧 | = 𝑎 . The modulus of a dual number can be negative, and it can also be calculated by thefollowing formula: | 𝑧 | = 12 ( 𝑧 + ¯ 𝑧 ) = 12 ( 𝑎 + ε 𝑏 + 𝑎 − ε 𝑏 ) = 𝑎. Devision:
The division of two dual numbers 𝑧 and 𝑧 is defined for all 𝑧 such that | 𝑧 | ̸ = 0 : 𝑧 𝑧 = 𝑎 + ε 𝑏 𝑎 + ε 𝑏 = ( 𝑎 + ε 𝑏 )( 𝑎 − ε 𝑏 )( 𝑎 + ε 𝑏 )( 𝑎 − ε 𝑏 ) = 𝑎 𝑎 + ε 𝑏 𝑎 − 𝑏 𝑎 𝑎 . From the above properties, it can be seen that in imaginary parts of the numbers don’t contribute to the real partof the result.Consider a dual number of the form 𝑐 ε , where 𝑐 ∈ R . This is a dual number with zero absolute value. Such numbershave the following property: ( 𝑐 ε ) · ( 𝑑 ε ) = 𝑐𝑑 ε = 0 , from which it follows that for any number 𝑐 ε there is a number 𝑑 ε such that the product of these numbers is . Suchnumbers are called zero divisors . B. Trigonometric form
A dual number 𝑧 such that | 𝑧 | ̸ = 0 can be written as follows: 𝑎 + ε 𝑏 = 𝑎 (︁ 𝑏𝑎 ε )︁ = 𝑟 (1 + 𝜙 ε ) , where the values 𝑟 = | 𝑧 | = 𝑎 and ε = 𝐴𝑟𝑔𝑧 = 𝑏/𝑎 are called module and argument of the dual number 𝑧 , respectively.This form of a dual number is some analog of the trigonometric form of an ordinary complex number and we willcontinue to call it trigonometric form of the dual number.For the conjugate number ¯ 𝑧 , the trigonometric form is: ¯ 𝑧 = 𝑟 (1 − 𝜙 ε ) , In this case, | ¯ 𝑧 | = | 𝑧 | is an absolute value, and 𝐴𝑟𝑔 ¯ 𝑧 = − 𝑏/𝑎 is an argument of a dual number in trigonometric form.The analogy with the trigonometric form of a complex number continues further when considering multiplicationand division. When multiplying two dual numbers 𝑧 and 𝑧 , we get 𝑧 = 𝑟 (1 + 𝜙 ε ) · 𝑟 (1 + 𝜙 ε ) = 𝑟 𝑟 (︁ 𝜙 + 𝜙 ) ε )︁ , in other words, when multiplying, the arguments are added and the modules are multiplied.In the case of division, we get 𝑧 𝑧 = 𝑟 (1 + 𝜙 ε ) 𝑟 (1 + 𝜙 ε ) = 𝑟 (1 + 𝜙 ε ) 𝑟 (1 − 𝜙 ε ) 𝑟 (1 + 𝜙 ε ) 𝑟 (1 − 𝜙 ε ) = 𝑟 (1 + ( 𝜙 − 𝜙 ) ε ) 𝑟 (1 − 𝜙 ε + 𝜙 ε ) = 𝑟 𝑟 (1 + ( 𝜙 − 𝜙 ) ε ) , in other words, when dividing, arguments are subtracted, and modules are divided.Note that in the case of division, the number 𝑧 must have a non-zero module | 𝑧 | ̸ = 0 .In trigonometric form, the operation of raising to the natural power of 𝑛 looks especially simple: 𝑧 𝑛 = (︀ 𝑟 (1 + 𝜙 ε ) )︀ 𝑛 = 𝑟 · 𝑟 · . . . · 𝑟 ⏟ ⏞ 𝑛 (︀ 𝜙 + 𝜙 + . . . + 𝜙 ) ⏟ ⏞ 𝑛 )︀ = 𝑟 𝑛 (1 + 𝑛𝜙 ε ) = 𝑎 𝑛 + ε 𝑛𝑎 𝑛 − 𝑏. For finding 𝑛 √ 𝑧 = 𝑛 √︀ 𝑟 (1 + 𝜙 ε ) , 𝑛 ∈ N , assume that 𝑛 √ 𝑧 = 𝑧 , then by definition 𝑧 𝑛 = 𝑧 , hence 𝑟 𝑛 (1 + 𝑛𝜙 ε ) = 𝑧 = 𝑟 (1 + 𝜙 ε ) , what gives 𝑟 = 𝑛 √ 𝑟, 𝜙 = 𝜙𝑛 , 𝑛 √ 𝑧 = 𝑛 √ 𝑟 (︁ 𝜙𝑛 ε )︁ . Note that for an odd 𝑛 , the root 𝑛 √ 𝑧 always exists, and for an even 𝑛 — only for the number 𝑧 with a non-negativemodule | 𝑧 | = 𝑟 ̸ = 0 . In algebraic form, the formula for the root is: 𝑛 √ 𝑎 + 𝑏 ε = 𝑛 √ 𝑎 + 𝑏𝑎 − 𝑛𝑛 𝑛 ε . C. Matrix form
All algebraic operations on dual numbers can be reduced to matrix operations if we consider: ε ↔ (︂ )︂ 𝑧 = 𝑎 + 𝑏 ε ↔ (︂ 𝑎 𝑏 𝑎 )︂ . Then, for example: 𝑧 · 𝑧 ↔ (︂ 𝑎 𝑏 𝑎 )︂ (︂ 𝑎 𝑏 𝑎 )︂ = (︂ 𝑎 𝑎 𝑎 𝑏 + 𝑎 𝑏 𝑎 𝑎 )︂ ↔ 𝑎 𝑎 + ( 𝑎 𝑏 + 𝑎 𝑏 ) ε . If the programming language supports vectorization, then in theory the matrix form of dual numbers can give you aperformance gain. However, the effect is unlikely to be significant in practice.
D. Taylor series expansion
We can use ε = ε = . . . = ε 𝑛 = 0 ∀ 𝑛 ∈ N , to get: exp( 𝑏 ε ) = ∞ ∑︁ 𝑛 =0 ( 𝑏 ε ) 𝑛 𝑛 ! = 1 + 𝑏 ε + 𝑏 ε
2! + . . . = 1 + 𝑏 ε , exp( 𝑎 + 𝑏 ε ) = 𝑒 𝑎 𝑒 𝑏 ε = 𝑒 𝑎 (1 + 𝑏 ε ) . A more general formula is derived from the Taylor series for the function 𝑓 ( 𝑧 ) at the point 𝑎 : 𝑓 ( 𝑎 + ε 𝑏 ) = ∞ ∑︁ 𝑛 =0 𝑓 ( 𝑛 ) ( 𝑎 )( 𝑎 + ε 𝑏 − 𝑎 ) 𝑛 𝑛 ! = ∞ ∑︁ 𝑛 =0 𝑓 ( 𝑛 ) ( 𝑎 ) ε 𝑛 𝑏 𝑛 𝑛 ! = 𝑓 ( 𝑎 ) + 𝑓 ′ ( 𝑎 ) 𝑏 ε + 𝑓 ′′ ( 𝑎 ) ε 𝑏
2! + . . . = 𝑓 ( 𝑎 ) + 𝑓 ′ ( 𝑎 ) 𝑏 ε . As a result, we get an extremely important equation: 𝑓 ( 𝑎 + ε 𝑏 ) = 𝑓 ( 𝑎 ) + 𝑓 ′ ( 𝑎 ) 𝑏 ε , which gives a way to calculate the values of functions from a dual number, if the value of the derivative 𝑓 ′ ( 𝑎 ) is known.On the other hand, the same formula allows one to calculate the value of the first derivative at the point 𝑎 . E. Elementary functions of dual numbers
The formula 𝑓 ( 𝑎 + ε 𝑏 ) = 𝑓 ( 𝑎 ) + 𝑓 ′ ( 𝑎 ) 𝑏 ε allows you to extend elementary functions to the set of dual numbers, sincethe right part of the formula contains only the values of the function 𝑓 from the real number 𝑎 .Here is a brief summary of basic elementary functions for illustration.Trigonometric function Inverse trigonometric functions sin( 𝑎 + ε 𝑏 ) = sin 𝑎 + 𝑏 ε cos 𝑎 arcsin( 𝑎 + ε 𝑏 ) = arcsin 𝑎 + 𝑏 ε / √ − 𝑎 cos( 𝑎 + ε 𝑏 ) = cos 𝑎 − 𝑏 ε sin 𝑎 arccos( 𝑎 + ε 𝑏 ) = arccos 𝑎 − 𝑏 ε / √ − 𝑎 tg( 𝑎 + ε 𝑏 ) = tg 𝑎 + 𝑏 ε / cos 𝑎 arctg( 𝑎 + ε 𝑏 ) = arctg 𝑎 + 𝑏 ε / (1 + 𝑎 )ctg( 𝑎 + ε 𝑏 ) = ctg 𝑎 − 𝑏 ε / sin 𝑎 arcctg( 𝑎 + ε 𝑏 ) = arctg 𝑎 − 𝑏 ε / (1 + 𝑎 ) Power functions Logarithmic functions and exponent ( 𝑎 + ε 𝑏 ) 𝑛 = 𝑎 𝑛 + 𝑛𝑎 𝑛 − 𝑏 ε exp( 𝑎 + ε 𝑏 ) = exp { 𝑎 } + 𝑏 ε exp { 𝑎 } 𝑛 √ 𝑎 + ε = 𝑛 √ 𝑎 (︁ 𝑏 ε 𝑛𝑎 )︁ log 𝑐 ( 𝑎 + ε 𝑏 ) = log 𝑐 𝑎 + 𝑏 ε /𝑎 ln 𝑎 F. Calculating the first derivative of a real value function using dual numbers
The formula 𝑓 ( 𝑎 + ε 𝑏 ) = 𝑓 ( 𝑎 ) + 𝑓 ′ ( 𝑎 ) 𝑏 ε can also be used to calculate the value of the first derivative of the realfunction 𝑓 ( 𝑥 ) at the point 𝑎 , if the value of the function 𝑓 ( 𝑎 + ε 𝑏 ) is known. From the point of view of analyticalcalculations, this method of finding the derivative does not make sense, since the value 𝑓 ( 𝑎 + ε 𝑏 ) is obtained analyticallyfrom the same formula. However, it provides a numerical method for obtaining the value of the first derivative usingnumerical calculations without any additional error. This method of finding the derivative with help of dual numbersis called automatic differentiation . It also should mensioned that automatic differentiation using dual numbers arelimited to derivatives of the first order.Let’s first look at a simple example of automatic differentiation, and then discuss the general implementationprinciple for programming languages.As an example, let’s find the derivative of the function 𝑓 ( 𝑥 ) = 𝑥 sin 𝑥 : 𝑓 ( 𝑎 + ε 𝑏 ) = ( 𝑎 + ε 𝑏 ) sin( 𝑎 + ε 𝑏 ) . Knowing the value of sin( 𝑎 + ε 𝑏 ) = sin 𝑎 + 𝑏 ε cos 𝑎 , it is easy to calculate 𝑓 ( 𝑎 + ε 𝑏 ) : 𝑓 ( 𝑎 + ε 𝑏 ) = ( 𝑎 + ε )(sin 𝑎 + 𝑏 ε cos 𝑎 ) = 𝑎 sin 𝑎 + (sin 𝑎 + 𝑎 cos 𝑎 ) 𝑏 ε . Comparing with the general formula 𝑓 ( 𝑎 + ε 𝑏 ) = 𝑓 ( 𝑎 ) + 𝑓 ′ ( 𝑎 ) 𝑏 ε , we get the value 𝑓 ′ ( 𝑎 ) = sin 𝑎 + 𝑏 ε cos 𝑎 .In general, to find 𝑓 ′ ( 𝑎 ) , it is enough to find the value of 𝑓 ′ ( 𝑎 + ε 𝑏 ) , take the imaginary part of this dual numberand divide it by 𝑏 . In addition, since the choice of the number 𝑏 is arbitrary, we can choose it equal to and get rid ofthe need to divide by 𝑏 .To apply automatic differentiation using dual numbers, we need the programming language which allows user-definedtypes, as well as overloading arithmetic operations and elementary functions for this type. This is especially easy forobject-oriented languages and languages that support function overloading or multiple dispatching. III. DYNAMICAL DISPATCH IN JULIA LANGUAGE
The language Julia [8] appeared relatively recently, but has already gained popularity as a language for scientificcomputing. We will assume that readers are already familiar with this language, and briefly focus only on the conceptof multiple dispatching [9–12], which is the basis of the language, its understanding is essential for further presentation.
Dynamic dispatch is a mechanism that allows to select the specific implementation of a polymorphic function oroperator from a set and call it in a specific case.
Multiple dispatching is based on dynamic dispatching. In this case, the choice of exact implementation of polymorphicfunction is made based on the type, number, and order of the function’s arguments. This is the runtime polymorphicdispatch. In addition to the term «multiple dispatching», the term multimethod is also used.The mechanism of multiple dispatch is similar to the mechanism of functions and operators overload, implemented,for example, in the C++ language. Function overloading, however, is performed exclusively at the compilation stage,whereas multiple dispatching must also work at the program run-time (run-time polymorphism).Julia supports overloading functions at the compilation time if all data types used in the function can be casted atthe compilation stage (so called type stable functions). The JIT compiler creates efficient implementations for eachcombination of argument types.If it is impossible to cast data types at the compilation stage (type unstable function), the dynamic dispatchingmechanism is enabled. The compiler will not be able to create a specialized version, but will create a generic versionthat runs slowly.This approach allows one to combine the speed of a compiled language with strict static typing with the flexibilityof an interpreted language with dynamic typing.
IV. DUAL NUMBERS IMPLEMENTATIONA. Existing implementations of dual numbers in Julia
The authors are familiar at least with two realizations of dual numbers in Julia language:• type
Dual in the module for automatic differentiation ForwardDiff [7];• separate module DualNumbers [13], the development of which is frozen in favor of ForwardDiff.In this section we will describe our proposed implementation of dual numbers in the Julia language. It is based on abuilt-in type of complex numbers
Complex [14], and on the DualNumbers module. The proposed implementation inthis paper serves only as a example of creating a custom data type for the Julia language.We put clarity of presentation at the first place, so many computational optimizations were deliberately omitted infavor of a larger clarity. For example, in the DualNumbers module, the dual number 𝑧 = 𝑎 + ε 𝑏 can have ordinarycomplex number components 𝑎 and 𝑏 . Also in DualNumbers for implementing elementary functions 𝑓 ( 𝑧 ) from dualnumbers 𝑧 a third-party module is used, which defines differentiation rules for functions from the real variable. Thisallows to automate the definition of the 𝑓 ( 𝑧 ) function by the formula 𝑓 ( 𝑎 + ε 𝑏 ) = 𝑓 ( 𝑎 ) + 𝑓 ′ ( 𝑎 ) 𝑏 ε , since it is notnecessary to explicitly write out derivatives of 𝑓 ′ ( 𝑎 ) .In our implementation, 𝑎 and 𝑏 are only real values and elementary functions are defined explicitly. B. Data structure
Adding any custom data type in Julia starts with defining the data structure. For dual numbers we define thestructure
Dual : struct Dual{T<:Real} <: Number"real part"x::T"imaginary part"y::Tend The structure contains only two fields: the real x and the imaginary y parts of the dual number. Both fields in thestructure must have the same parametric type T , which must be a subtype of the abstract type Real , as indicated bythe <: operator. The Dual type itself is a subtype of the abstract type
Number . In this way, the dual numbers areembedded in an existing hierarchy of types (see figure 1).
NumberDual Complex RealIntegerSignedInt8 Int16 Int32 Int64 Int128Bool UnsignedUInt8 UInt16 UInt32 UInt64 UInt128Rational AbstractFloatFloat16 Float32 Float64
Legend:
Abstract typePrimitive typeComposite types
Figure 1. Our Dual type in Julia built-in types hierarchy
Immediately after declaring a new structure, we can create variables of the Dual type using the default constructors.We must specify both fields of the structure as arguments of the Dual function: z = Dual(1, 2)
Dual{Float64}(1, 2.2)
If the argument type is the same, we don’t need to explicitly specify the T parameter. If the arguments have differenttypes, the parametric type must be specified. In our example, both arguments are casted to the Float64 type.
C. Overloading functions in the show example
To print the values of variables of the Dual type we must to overload the show function. This function is intended forformatted printing of data to standard output and to REPL. The write , print , and println functions are used tooutput a minimal text representation of data. The last two functions print data only to the standard output stream.For simple data, one can make no difference between formatted and compact views and overload only show function.Then the other functions will call show .In terms of multiple dispatch, overloading a function means adding a new method to the function. To add a methodto the show function, we define it as follows: Base.show(io::IO, z::Dual) = show(io, z.x, " + ", z.y, " ε ") The short syntax for function definition is used here. The prefix
Base. is needed, since the standard methods of the show function are located in the
Base module. Functions from this module can be called directly without a prefix, butfor overloading the
Base prefix is essential.It is worth noting that our version of the show method is simplified, since it does not take into account special cases,for example, the negative imaginary part.Next, we will add methods for many functions from
Base , so that our
Dual type becomes sufficiently functional.
D. Additional constructors
Let us create some additional constructors. To do this, we should overload the default constructor
Dual .In total we will set 5 additional constructors:
Dual(x::Real, y::Real)= Dual(promote(x, y)...)
Dual{T}(x::Real) where {T<:Real} = Dual{T}(x, 0)
Dual(x::Real) = Dual(promote(x, 0)...)
Dual{T}(z::Dual) where {T<:Real} = Dual{T}(z.x, z.y)
Dual(z::Dual) = Dual(promote(z.x, z.y)...)
The first constructor allows us not to specify the T parameter every time if the arguments have different types. The promote function converts the type of arguments passed to it to a common type and returns the result as a tuple. ThePostfix operator ... unpacks the tuple and passes its elements as arguments to the constructor function. The languagecore defines casting rules for all subtypes of the abstract type Real , so now the constructor will work correctly for anycombination of arguments, as long as the
T<:Real is true. For example, the following code will work correctly:
Dual(1//2, π ) π * ε We passed a rational number (type
Rational ) and a built-in global constant (number 𝜋 ) of the type Float64 to theconstructor. After that, the type casting rule worked and both arguments were converted to the more general type
Float64 .The second and third additional constructors allow one to omit the imaginary part if it is equal to zero:
Dual{Float32}(1) ε , 2 constructor Dual(1//2) ε , 3 constructor Constructor number 2 is a parametric function that is declared using the where construct. The T parameter is asubtype of the abstract type Real . Constructor number 3 works similarly to first constructor. The fourth and fifthconstructors allow one to pass in an argument to the constructor other dual number.For more convenience, we can also create a separate constant for the imaginary unit ε : const ε = Dual(0, 1) After overloading arithmetic operations, this constant will allow us to create new dual numbers using an expressionas close as possible to their algebraic notation: z = 1 + 2 ε E. Access to structure fields
Structures in Julia are immutable by default, that is, once we create the variable z , we can access the fields of thestructure itself using the point operator, but we cannot modify the value of these fields: @show z.x z. y z.x = 1 To create mutable data types, the mutable struct structure is provided. However, for our example, an immutablestructure is more appropriate, since the requirement for mutability imposes performance restrictions, which is undesirablefor a numeric type.Julia doesn’t have access modifiers for fields like public or private , so structure fields are always readable. However,it is considered a good programming style to encapsulate structure fields and provide an interface for accessing them.This style is used, for example, for the built-in type of complex numbers.There are two ways to encapsulate fields. The first method is to create special interface functions for accessing fields.In our case it is sufficient to define the following functions: Base.real(z::Dual) = z.xBase.imag(z::Dual) = z.yBase.reim(z::Dual) = (z.x, z.y)
The same functions are defined in the built-in Complex module. Next, we should call only these functions everywhereto get the structure fields, instead of accessing the fields by name directly. This will allow developers to refactor thestructure in the future, for example, rename or add new fields. Backward compatibility is easy to maintain by rewritingonly the interface functions. Performance will not suffer, since the JIT compiler will replace calls to these functionsdirectly with their code, since they are extremely simple (inline functions).The second approach is to use the getproperty function, which has been available since version 0.7 of the Julialanguage. Overloading this function allows us to set additional names for accessing fields in the structure. So, if wewrite the following method for example: function Base.getproperty(z::Dual, p::Symbol)if p in (:real, :a, :r) return getfield(z, :x)elseif p in (:imag, :b, :i) return getfield(z, :y)elsereturn getfield(z, p)endend then we will be able to access the real and imaginary parts of the number 𝑧 in four different ways: z.x == z.a == z.r == z.real z.y == z.b == z.i == z.imag This approach is more flexible, since when changing the structure for backward compatibility, it is enough to modifyonly one method getproperty , not all the interface functions (if the structure is complex it can by lots of suchfunctions). In addition, the programmer can use any name from the aliases list to accessing the structure fields.For the Dual type, we applied both approaches, since the presence of the functions real — equivalent to
Re( 𝑧 ) , and imag — equivalent to Im( 𝑧 ) is mathematically justified. F. Unary functions
The
Base module defines the functions one and zero , which return a multiplication identity element and an additionidentity element, respectively. In the case of dual numbers, the multiplication identity is the real unit, and the additionidentity is the real zero.For the Dual type, we define the following single-line methods:
Base.one(::Type{Dual{T}}) where T <: Real = Dual(one(T))Base.one(z::Dual) = Dual(one(z.x))Base.zero(::Type{Dual{T}}) where T <: Real = Dual(zero(T))Base.zero(z::Dual) = Dual(zero(z.x))
The first version of the one and zero functions takes the data type as an argument and returns one or zero of thistype, respectively. Since for dual numbers it is and from R , it is sufficient to return a dual number with a zeroimaginary part. The real part will be and of the parametric type T . For all standard types T<:Real the one and zero methods are defined in the
Base module, which we used.The second variant of functions takes a specific object of the
Dual type as an argument and also returns a identityusing methods from the
Base module.For complex numbers in
Base , the conjugation functions conj , the absolute value abs , and the argument arg aredefined:
Base.conj(z::Dual) = Dual(z.x, -z.y)Base.abs(z::Dual) = abs(z.x)Base.abs2(z::Dual) = z.x^2function Base.arg(z::Dual)@assert !iszero(z.x)return z.y / z.xend
For the arg method, we have provided a check for the inequality of the real part of the dual number to zero, sinceotherwise we will get a division by zero. Using the standard function iszero allows us not to worry about accountingfor errors in the representation of real numbers using floating-point numbers. The @assert macro throws an exception
AssertionError if the actual part is equal to zero.The inv function defines the inverse number for this number z : function Base.inv(z::Dual)@assert !iszero(z.x)return Dual(1/z.x, -z.y/z.x^2)end There should also be an exception for the null real part.
G. Comparison function
The
Base module also defines a number of functions that return the truth, if a particular condition is true. Listsome of these functons are:• isreal is real number;• isinteger is integer;• isfinite the number is finite;• isnan the argument has the type
NaN ;• isinf the argument has the type Inf ;• iszero the number is an addition identity (zero);• isone the number is a multiplication identity (one).Methods for these functions for the Dual one-to-one case repeats methods for the built-in
Complex type, so here wedo not provide their source code.In addition to unary functions for dual numbers, it makes sense to implement methods for the comparison operator == . Operators in Julia are no different from functions, and adding methods for them is exactly the same. The onlydifference is that we need to add the character : after Base. , and since the == operator consists of two characters, wemust frame it with parentheses: Base.:(==)(z::Dual, u::Dual) = (real(z) == real(u)) && (imag(z) == imag(u))
It makes sense to compare dual numbers with real numbers as well, if the dual number has a zero imaginary part.In order for the operator to commute for arguments of different types, we must define two methods:
Base.:(==)(z::Dual, x::Real) = isreal(z) && real(z) == xBase.:(==)(x::Real, z::Dual) = isreal(z) && real(z) == x H. Arithmetic operations
We also need to overload arithmetic operators such as the unary operators + and - and the binary operators + , - , * and / . The implementation of the + , - , and * operators is trivial: Base.:+(z::Dual) = Dual(+z.x, +z.y)Base.:-(z::Dual) = Dual(-z.x, -z.y)Base.:+(z::Dual, u::Dual) = Dual(z.x + u.x, z.y + u.y)Base.:-(z::Dual, u::Dual) = Dual(z.x - u.x, z.y - u.y)Base.:*(z::Dual, u::Dual) = Dual(z.x * u.x, z.x * u.y + z.y * u.x)
In the case of the / operator, we should take into account the equality of the divisor to zero: function Base.:/(z::Dual, u::Dual)@assert !iszero(u.x)return Dual(z.x/u.x, (z.y*u.x-u.y*z.x)/u.x^2)end For binary operators, we do not need to implement separate cases of arguments of different types, because in thiscase the promotion mechanism to a common type will call promote function automatically.
I. Types promotion
To make type promotion mechanism to work, we must define the type promotion rules. Without these rules, forexample, the next operation fails with an error: >>Dual(1, 2) + 2ERROR: promotion of types Dual{Int64} and Int64 failed to change any arguments
This error occurs because there is no rule that can be used to determine the common type for numbers of the type
Dual{Int64} and the type
Int64 .To define such a rule, we have to add a method for the promote_rule function from the
Base module. Any numberof the
Real type can participate in the arithmetic operation with a dual number, since this number can be interpretedas a dual with a zero imaginary part:
Base.promote_rule(::Type{Dual{T}}, ::Type{S}) where {T<:Real, S<:Real} = Dual{promote_type(T, S)}
We should also provide a rule that will work for operators with two dual numbers with different parametric types T and S . In this case, we need to find a common type for the T and S types: Base.promote_rule(::Type{Dual{T}}, ::Type{Dual{S}}) where {T<:Real, S<:Real} = Dual{promote_type(T, S)}
J. Raising to a rational degree
To raise a number to a power, the ^ operator is used, which should also be overloaded for integer and rationalpowers. In the case of an integer degree, we use the formula ( 𝑎 + ε 𝑏 ) 𝑛 = 𝑎 𝑛 + 𝑛𝑏𝑎 𝑛 − ε . There are several special cases to consider:• if 𝑛 = 0 , then 𝑧 𝑛 = 1 ;• if 𝑛 = 1 , then 𝑧 𝑛 = 𝑧 ;• if 𝑛 > , the formula works for any 𝑎 and 𝑏 ;• if 𝑛 < , the condition 𝑎 ̸ = 0 must be met.1If we consider all these cases, we get the following implementation: function Base.:^(z::Dual, n::Integer)::Dualx, y = reim(z)if n == 0return one(z)elseif n == 1return zelseif n > 1return Dual(x^n, n*y*x^(n-1))else if iszero(x)throw(DomainError(:x, "negative exponentiation is only defined for real(z)!= 0"))elsereturn Dual(x^n, n*y*x^(n-1))endendend If the real part of the dual number is zero ( iszero(x) ) and 𝑛 < , the function throws an exception DomainError (out of the range of acceptable values).For a rational degree, a more complex formula is used: ( 𝑎 + ε 𝑏 ) 𝑛𝑚 = 𝑎 𝑛𝑚 (︂ ε 𝑛𝑚 𝑏𝑎 )︂ = 𝑎 𝑛𝑚 + ε 𝑛𝑚 𝑏𝑎 𝑛𝑚 − , for which the cases of even and odd 𝑚 should be provided. For an odd 𝑚 , the range of acceptable values includes 𝑎 < , and for an even 𝑚 , you should limit yourself to 𝑎 > : function Base.:^(z::Dual, q::Rational)::Dualx, y = reim(z)n, d = numerator(q), denominator(q)if n == 0return one(z)elseif isodd(d)return Dual(x^q, q*y*x^(q-1))elseif x < 0throw(DomainError(:x, "even radical for dual number z is only defined for real(z)>=0"))elsereturn Dual(x^q, q*y*x^(q-1))endendend It makes sense to separately overload the functions for square and cubic roots sqrt and cbrt : function Base.sqrt(z::Dual)::Dualx, y = reim(z)if x < 0throw(DomainError(:x, "sqrt for dual number z is only defined for real(z)>=0"))elsereturn Dual(√x, y/(2*√x))endendBase.cbrt(z::Dual) = Dual(cbrt(z.x), z.y*cbrt(z.x)/(3z.x)) K. Elementary functions
Elementary functions are calculated using the formula 𝑓 ( 𝑎 + ε 𝑏 ) = 𝑓 ( 𝑎 ) + 𝑓 ′ ( 𝑎 ) 𝑏 ε . Note that many of them are notdefined for the case of the zero real part of the number 𝑧 = 𝑎 + ε 𝑏 : function Base.exp(z::Dual)e = exp(z.x)return Dual(e, e * z.y)endBase.log(z::Dual) = Dual(log(z.x), z.y/z.x)Base.log(b, z::Dual) = Dual(log(b, z.x), (z.y/z.x) * log(z.x)) Base.sin(z::Dual) = Dual(sin(z.x), +z.y*cos(z.x))Base.cos(z::Dual) = Dual(cos(z.x), -z.y*sin(z.x))Base.tan(z::Dual) = Dual(tan(z.x), z.y/cos(z.x)^2)Base.cot(z::Dual) = Dual(cot(z.x), -z.y/sin(z.x)^2)Base.asin(z::Dual) = Dual(asin(z.x), z.y/sqrt(1-z.x^2))Base.acos(z::Dual) = Dual(acos(z.x), -z.y/sqrt(1-z.x^2))Base.atan(z::Dual) = Dual(atan(z.x), z.y/(1+z.x^2))Base.acot(z::Dual) = Dual(acot(z.x), -z.y/(1+z.x^2))
Base.sinh(z::Dual) = Dual(sinh(z.x), z.y*cosh(z.x))Base.cosh(z::Dual) = Dual(cosh(z.x), z.y*sinh(z.x))Base.tanh(z::Dual) = Dual(tanh(z.x), z.y/cosh(z.x)^2)Base.coth(z::Dual) = Dual(coth(z.x), z.y/sinh(z.x)^2)
V. RESULTS
The paper describes the preliminary implementation of dual complex numbers and basic operations on them in theJulia language.
VI. DISCUSSION
The implementation of dual complex numbers on Julia is completely based on the mechanism of multiple dispatching.Thus, we not only implemented a certain set of operations on dual numbers, but also demonstrated the power of thismechanism.It should also be noted that in contrast to the implementation of the
Dual type in the automatic differentiationpackage
ForwardDiff [7], our proposed implementation is cleaner. For example, in the above package, an ordinarycomplex number can be used as a coefficient before a dual complex unit. It is clear that this is due to the specifics ofusing dual numbers in this package for automatic differentiation. But this type of number is more likely to belong toquaternions [15, 16], rather than the proper complex numbers.
VII. CONCLUSION
In this paper, a prototype of the implementation of generalized complex numbers in the Julia language wasmade, namely, the implementation of dual complex numbers. Having a multiple dispatching mechanism in the Julia3language makes it very easy to implement new numeric types within the existing programming language infrastructure.We propose to further extend this prototype to more general implementation of generalized complex numbers andgeneralized quaternions.
ACKNOWLEDGMENTS
The publication has been prepared with the support of the Russian Foundation for Basic Research (RFBR) accordingto the research project No 19-01-00645. [1] I. M. Yaglom, B. A. Rozenfel’d, E. U. Yasinskaya, Projective Metrics, Russian Mathematical Surveys 19 (1964) 49–107.doi:10.1070/RM1964v019n05ABEH001159.[2] I. M. Yaglom, Complex Numbers in Geometry, Academic Press, 1968.[3] B. A. Rozenfel’d, I. M. Yaglom, On the geometries of the simplest algebras, Mat. Sbornik N. S. 28(70) (1951) 205–216.[4] D. S. Kulyabov, A. V. Korolkova, M. N. Gevorkyan, Hyperbolic numbers as Einstein numbers, Journal of Physics:Conference Series 1557 (2020) 012027.1–5. doi:10.1088/1742-6596/1557/1/012027.[5] A. Cayley, IV. A sixth memoir upon quantics, Philosophical Transactions of the Royal Society of London 149 (1859) 61–90.doi:10.1098/rstl.1859.0004.[6] F. Klein, Ueber die sogenannte Nicht-Euklidische Geometrie, in: Gau ß und die Anf¨ange der nicht-euklidischen Geometrie,volume 4 of Teubner-Archiv zur Mathematik , Springer-Verlag Wien, Wien, 1985, pp. 224–238. doi:10.1007/978-3-7091-9511-6_5.[7] Forward Mode Automatic Differentiation for Julia, 2020. URL: https://github.com/JuliaDiff/DualNumbers.jl .[8] The Julia Language, 2020. URL: https://julialang.org/ .[9] J. Bezanson, J. Chen, B. Chung, S. Karpinski, V. B. Shah, J. Vitek, L. Zoubritzky, Julia: dynamism and performancereconciled by design, Proceedings of the ACM on Programming Languages 2 (2018) 1–23. doi:10.1145/3276490.[10] F. Zappa Nardelli, J. Belyakova, A. Pelenitsyn, B. Chung, J. Bezanson, J. Vitek, Julia subtyping: a rational reconstruction,Proceedings of the ACM on Programming Languages 2 (2018) 1–27. doi:10.1145/3276483.[11] J. Bezanson, A. Edelman, S. Karpinski, V. B. Shah, Julia: A fresh approach to numerical computing, SIAM Review 59(2017) 65–98. doi:10.1137/141000671. arXiv:1411.1607 .[12] J. Bezanson, S. Karpinski, V. B. Shah, A. Edelman, Julia: A Fast Dynamic Language for Technical Computing (2012)1–27. arXiv:1209.5145 .[13] Julia package for representing dual numbers and for performing dual algebra, 2020. URL: https://github.com/JuliaDiff/ForwardDiff.jl .[14] Official Julia language GitHub repository, 2020. URL: https://github.com/JuliaLang/julia .[15] W. R. Hamilton, Elements of Quaternions, Cambridge University Press, Cambridge, 1866. doi:10.1017/CBO9780511707162.[16] A. P. Yefremov, Quaternions and Biquaternions: Algebra, Geometry and Physical Theories, 2005. arXiv:0501055 . одходы к реализации обобщённых комплексных чисел на языке Julia М. Н. Геворкян, * А. В. Королькова, † и Д. С. Кулябов
1, 2, ‡ Кафедра прикладной информатики и теории вероятностей,Российский университет дружбы народов,117198, Москва, ул. Миклухо-Маклая, д. 6 Лаборатория информационных технологий,Объединённый институт ядерных исследований,ул. Жолио-Кюри 6, Дубна, Московская область, Россия, 141980
В задачах математической физики для исследования структур пространств с применениеммоделей Кэли–Клейна в теоретических расчётах требуется использование обобщённых ком-плексных чисел. В случае проведения вычислительных экспериментов для такого рода задачнеобходима их качественная реализация в языке программирования. Предлагаемые малочис-ленные реализации обобщённых комплексных чисел в современных языках программированияимеют ряд недостатков. В данной статье мы предлагаем в качестве языка реализации обобщён-ных комплексных чисел использовать язык Julia, не в последнюю очередь из-за поддержки иммеханизма множественной диспетчеризации.В работе демонстрируется подход к реализации одного из типов обобщённых комплексныхчисел, а именно дуальных чисел. Особый упор мы делаем на описание использования меха-низма множественной диспетчеризации для реализации нового числового типа. Полученнуюреализацию дуальных чисел можно рассматривать как прототип для полного программногомодуля поддержки обобщённых комплексных чисел.
Ключевыеслова: комплексные числа, параболические комплексные числа, дуальные числа, множествен-ная диспетчеризация, Julia * [email protected] † [email protected] ‡ [email protected] a r X i v : . [ c s . M S ] J u l I. ВВЕДЕНИЕ
Существует подход, который позволяет обобщить обычные комплексные числа и получить три различныхкласса обобщённых комплексных чисел [1–4]. В этом подходе для определения трёх видов комплексных чиселзадаётся квадратное уравнение 𝑧 + 𝑝𝑧 + 𝑞 = 0 . с детерминантом Δ = 𝑝 − 𝑝𝑞. В зависимости от знака детерминанта получаются следующие комплексные числа:• Δ < — эллиптические (обычные комплексные числа);• Δ = 0 — параболические (дуальные комплексные числа);• Δ > — гиперболические (двойные комплексные (расщеплённые) числа).Первоначально эти системы комплексных чисел вводились для описания моделей Кэли–Клейна [5, 6]. Однакоони имеют и другие применения. В частности, дуальные комплексные числа могут использоваться в задачахавтоматического дифференцирования [7].Интерес представляет чистая реализация обобщённых комплексных чисел, не привязанная к частным при-менениям. Существующие реализации обобщённых комплексных чисел в рамках пакетов языка Julia — неудовлетворительны. В данной статье мы предлагаем свой подход по реализации дуальных комплексных чисел вJulia. A. Структура статьи
В разделе II приводится достаточно подробное описание дуальных комплексных чисел. Описание конструктив-ное, затрагивающее только операции с дуальными комплексными числами. Мы постарались избегать излишнихматематических абстракций. В разделе III даётся общее понятие о механизме множественной диспетчеризации вJulia. На использовании этого механизма и базируется предлагаемая нами реализация дуальных комплексныхчисел. В разделе IV собственно и описывается наша реализация дуальных комплексных чисел на Julia.
B. Обозначения и соглашения
Для обозначения дуальной комплексной единицы мы будем использовать символ ε . II. ДУАЛЬНЫЕ ЧИСЛА
Дуальное число 𝑧 алгебраически определяется следующим образом: 𝑧 = 𝑎 + ε 𝑏, 𝑎, 𝑏 ∈ R , ε = 0 , ε ̸ = 0 . Величину 𝑎 будем называть действительной частью, а 𝑏 — мнимой. Так как в рамках данной статьи не будетвстречаться обычных комплексных чисел, то путаницы в терминологии не возникнет. A. Алгебраическая форма
Алгебраические свойства для двух чисел 𝑧 = 𝑎 + ε 𝑏 и 𝑧 = 𝑎 + ε 𝑏 . Сложение: 𝑧 + 𝑧 = ( 𝑎 + 𝑎 ) + ε ( 𝑏 + 𝑏 ) . Вычитание: 𝑧 − 𝑧 = ( 𝑎 − 𝑎 ) + ε ( 𝑏 − 𝑏 ) . Умножение: 𝑧 · 𝑧 = ( 𝑎 + ε 𝑏 ) · ( 𝑎 + ε 𝑏 ) = 𝑎 𝑎 + ε ( 𝑎 𝑏 + 𝑏 𝑎 ) . Сопряжённое число: ¯ 𝑧 = 𝑎 − ε 𝑏 , 𝑧 ¯ 𝑧 = ( 𝑎 + ε 𝑏 ) · ( 𝑎 − ε 𝑏 ) = 𝑎 . Абсолютное значение числа: | 𝑧 | = 𝑎 .Модуль дуального числа может быть отрицательным, кроме того его можно вычислить по следующейформуле: | 𝑧 | = 12 ( 𝑧 + ¯ 𝑧 ) = 12 ( 𝑎 + ε 𝑏 + 𝑎 − ε 𝑏 ) = 𝑎. Деление:
Деление двух дуальных чисел 𝑧 и 𝑧 определено для всех 𝑧 таких, что | 𝑧 | ̸ = 0 : 𝑧 𝑧 = 𝑎 + ε 𝑏 𝑎 + ε 𝑏 = ( 𝑎 + ε 𝑏 )( 𝑎 − ε 𝑏 )( 𝑎 + ε 𝑏 )( 𝑎 − ε 𝑏 ) = 𝑎 𝑎 + ε 𝑏 𝑎 − 𝑏 𝑎 𝑎 . Из вышеперечисленных свойств видно, что ни в одном случае мнимые части чисел не вносят вклад вдействительную часть результата.Рассмотрим дуальное число вида 𝑐 ε , где 𝑐 ∈ R . Это дуальное число с нулевым абсолютным значением. Такиечисла обладают следующим свойством: ( 𝑐 ε ) · ( 𝑑 ε ) = 𝑐𝑑 ε = 0 , из которого следует, что для любого числа 𝑐 ε существует число 𝑑 ε , такое что произведение этих чисел равняется . Такие числа называются делителями нуля . B. Тригонометрическая форма
Дуальное число 𝑧 , такое что | 𝑧 | ̸ = 0 , можно записать в следующем виде: 𝑎 + ε 𝑏 = 𝑎 (︁ 𝑏𝑎 ε )︁ = 𝑟 (1 + 𝜙 ε ) , где величины 𝑟 = | 𝑧 | = 𝑎 и ε = 𝐴𝑟𝑔𝑧 = 𝑏/𝑎 называются модулем и аргументом дуального числа 𝑧 соответственно.Такая форма записи дуального числа является некоторым аналогом тригонометрической формы обычногокомплексного числа и далее будем её называть тригонометрической формой дуального числа.Для сопряжённого числа ¯ 𝑧 тригонометрическая форма имеет вид: ¯ 𝑧 = 𝑟 (1 − 𝜙 ε ) , При этом | ¯ 𝑧 | = | 𝑧 | — модуль, а 𝐴𝑟𝑔 ¯ 𝑧 = − 𝑏/𝑎 — аргумент дуального комплексного числа в тригонометрическойформе.Аналогия с тригонометрической формой комплексного числа продолжается далее при рассмотрении умноженияи деления. При умножении двух дуальных чисел 𝑧 и 𝑧 получим 𝑧 = 𝑟 (1 + 𝜙 ε ) · 𝑟 (1 + 𝜙 ε ) = 𝑟 𝑟 (︁ 𝜙 + 𝜙 ) ε )︁ , то есть при умножении аргументы складываются, а модули умножаются.В случае же деления, получим 𝑧 𝑧 = 𝑟 (1 + 𝜙 ε ) 𝑟 (1 + 𝜙 ε ) = 𝑟 (1 + 𝜙 ε ) 𝑟 (1 − 𝜙 ε ) 𝑟 (1 + 𝜙 ε ) 𝑟 (1 − 𝜙 ε ) = 𝑟 (1 + ( 𝜙 − 𝜙 ) ε ) 𝑟 (1 − 𝜙 ε + 𝜙 ε ) = 𝑟 𝑟 (1 + ( 𝜙 − 𝜙 ) ε ) , то есть при делении аргументы вычитаются, а модули делятся.Обратите внимание, что в случае деления число 𝑧 должно иметь ненулевой модуль | 𝑧 | ̸ = 0 .В тригонометрической форме особенно просто выглядит операция возведения в натуральную степень 𝑛 : 𝑧 𝑛 = (︀ 𝑟 (1 + 𝜙 ε ) )︀ 𝑛 = 𝑟 · 𝑟 · . . . · 𝑟 ⏟ ⏞ 𝑛 (︀ 𝜙 + 𝜙 + . . . + 𝜙 ) ⏟ ⏞ 𝑛 )︀ = 𝑟 𝑛 (1 + 𝑛𝜙 ε ) = 𝑎 𝑛 + ε 𝑛𝑎 𝑛 − 𝑏. Для нахождения 𝑛 √ 𝑧 = 𝑛 √︀ 𝑟 (1 + 𝜙 ε ) , 𝑛 ∈ N , предположим, что 𝑛 √ 𝑧 = 𝑧 , тогда по определению 𝑧 𝑛 = 𝑧 ,следовательно 𝑟 𝑛 (1 + 𝑛𝜙 ε ) = 𝑧 = 𝑟 (1 + 𝜙 ε ) , что даёт 𝑟 = 𝑛 √ 𝑟, 𝜙 = 𝜙𝑛 , 𝑛 √ 𝑧 = 𝑛 √ 𝑟 (︁ 𝜙𝑛 ε )︁ . Следует заметить, что при нечётном 𝑛 корень 𝑛 √ 𝑧 существует всегда, а при чётном 𝑛 — только для числа 𝑧 снеотрицательным модулем | 𝑧 | = 𝑟 ̸ = 0 . В алгебраической форме формула для извлечения корня имеет вид: 𝑛 √ 𝑎 + 𝑏 ε = 𝑛 √ 𝑎 + 𝑏𝑎 − 𝑛𝑛 𝑛 ε . C. Матричная форма
Все алгебраические операции над дуальными числами можно свести к матричным операциям, если положить ε ↔ (︂ )︂ 𝑧 = 𝑎 + 𝑏 ε ↔ (︂ 𝑎 𝑏 𝑎 )︂ . Тогда, например: 𝑧 · 𝑧 ↔ (︂ 𝑎 𝑏 𝑎 )︂ (︂ 𝑎 𝑏 𝑎 )︂ = (︂ 𝑎 𝑎 𝑎 𝑏 + 𝑎 𝑏 𝑎 𝑎 )︂ ↔ 𝑎 𝑎 + ( 𝑎 𝑏 + 𝑎 𝑏 ) ε . Теоретически с вычислительной точки зрения сводить действия над дуальными числами к матричнымоперациям может быть оправданно в случае, если язык программирования поддерживает векторизацию операций.Однако вряд ли на практике эффект будет значителен.
D. Разложение в ряд Тейлора
Используя свойство ε = ε = . . . = ε 𝑛 = 0 ∀ 𝑛 ∈ N , получим: exp( 𝑏 ε ) = ∞ ∑︁ 𝑛 =0 ( 𝑏 ε ) 𝑛 𝑛 ! = 1 + 𝑏 ε + 𝑏 ε
2! + . . . = 1 + 𝑏 ε , exp( 𝑎 + 𝑏 ε ) = 𝑒 𝑎 𝑒 𝑏 ε = 𝑒 𝑎 (1 + 𝑏 ε ) . Более общая формула выводится из ряда Тейлора для функции 𝑓 ( 𝑧 ) в точке 𝑎 : 𝑓 ( 𝑎 + ε 𝑏 ) = ∞ ∑︁ 𝑛 =0 𝑓 ( 𝑛 ) ( 𝑎 )( 𝑎 + ε 𝑏 − 𝑎 ) 𝑛 𝑛 ! = ∞ ∑︁ 𝑛 =0 𝑓 ( 𝑛 ) ( 𝑎 ) ε 𝑛 𝑏 𝑛 𝑛 ! = 𝑓 ( 𝑎 ) + 𝑓 ′ ( 𝑎 ) 𝑏 ε + 𝑓 ′′ ( 𝑎 ) ε 𝑏
2! + . . . = 𝑓 ( 𝑎 ) + 𝑓 ′ ( 𝑎 ) 𝑏 ε . В результате получаем крайне важную формулу: 𝑓 ( 𝑎 + ε 𝑏 ) = 𝑓 ( 𝑎 ) + 𝑓 ′ ( 𝑎 ) 𝑏 ε , которая даёт способ вычисления значений функций от дуального числа, если известно значение производной 𝑓 ′ ( 𝑎 ) . С другой стороны, эта же формула позволяет вычислить значение производной в точке 𝑎 . E. Элементарные функции от дуальных чисел
Формула 𝑓 ( 𝑎 + ε 𝑏 ) = 𝑓 ( 𝑎 ) + 𝑓 ′ ( 𝑎 ) 𝑏 ε позволяет распространить элементарные функции на множество дуальныхчисел, так как в правой части формулы находятся лишь значения функции 𝑓 от действительного числа 𝑎 .Приведём для иллюстрации небольшую сводку основных элементарных функций.Trigonometric function Inverse trigonometric functions sin( 𝑎 + ε 𝑏 ) = sin 𝑎 + 𝑏 ε cos 𝑎 arcsin( 𝑎 + ε 𝑏 ) = arcsin 𝑎 + 𝑏 ε / √ − 𝑎 cos( 𝑎 + ε 𝑏 ) = cos 𝑎 − 𝑏 ε sin 𝑎 arccos( 𝑎 + ε 𝑏 ) = arccos 𝑎 − 𝑏 ε / √ − 𝑎 tg( 𝑎 + ε 𝑏 ) = tg 𝑎 + 𝑏 ε / cos 𝑎 arctg( 𝑎 + ε 𝑏 ) = arctg 𝑎 + 𝑏 ε / (1 + 𝑎 )ctg( 𝑎 + ε 𝑏 ) = ctg 𝑎 − 𝑏 ε / sin 𝑎 arcctg( 𝑎 + ε 𝑏 ) = arctg 𝑎 − 𝑏 ε / (1 + 𝑎 ) Power functions Logarithmic functions and exponent ( 𝑎 + ε 𝑏 ) 𝑛 = 𝑎 𝑛 + 𝑛𝑎 𝑛 − 𝑏 ε exp( 𝑎 + ε 𝑏 ) = exp { 𝑎 } + 𝑏 ε exp { 𝑎 } 𝑛 √ 𝑎 + ε = 𝑛 √ 𝑎 (︁ 𝑏 ε 𝑛𝑎 )︁ log 𝑐 ( 𝑎 + ε 𝑏 ) = log 𝑐 𝑎 + 𝑏 ε /𝑎 ln 𝑎 F. Вычисление первой производной от вещественной функции с помощью дуальных чисел
Формулу 𝑓 ( 𝑎 + ε 𝑏 ) = 𝑓 ( 𝑎 ) + 𝑓 ′ ( 𝑎 ) 𝑏 ε можно применять и для вычисления значения первой производнойвещественной функции 𝑓 ( 𝑥 ) в точке 𝑎 , если известно значении функции 𝑓 ( 𝑎 + ε 𝑏 ) . С точки зрения аналитическихвычислений такой способ нахождения производной не имеет смысла, так как сама величина 𝑓 ( 𝑎 + ε 𝑏 ) получаетсяаналитически из известной формулы для первой производной 𝑓 ′ ( 𝑥 ) . Однако она даёт численный метод дляполучения значения первой производной с помощью компьютерных вычислений без какой-либо дополнительнойпогрешности. Такой способ нахождения производной получил название автоматического дифференцирования сиспользованием дуальных чисел.Рассмотрим вначале простой пример автоматического дифференцирования, а далее обсудим общий принципреализации для языков программирования.В качестве примера найдём производную от функции 𝑓 ( 𝑥 ) = 𝑥 sin 𝑥 : 𝑓 ( 𝑎 + ε 𝑏 ) = ( 𝑎 + ε 𝑏 ) sin( 𝑎 + ε 𝑏 ) . Зная значение sin( 𝑎 + ε 𝑏 ) = sin 𝑎 + 𝑏 ε cos 𝑎 , легко вычислить 𝑓 ( 𝑎 + ε 𝑏 ) : 𝑓 ( 𝑎 + ε 𝑏 ) = ( 𝑎 + ε )(sin 𝑎 + 𝑏 ε cos 𝑎 ) = 𝑎 sin 𝑎 + (sin 𝑎 + 𝑎 cos 𝑎 ) 𝑏 ε . Сравнивая с общей формулой 𝑓 ( 𝑎 + ε 𝑏 ) = 𝑓 ( 𝑎 ) + 𝑓 ′ ( 𝑎 ) 𝑏 ε , получаем значение 𝑓 ′ ( 𝑎 ) = sin 𝑎 + 𝑏 ε cos 𝑎 .В общем случае для нахождения 𝑓 ′ ( 𝑎 ) достаточно найти значение 𝑓 ′ ( 𝑎 + ε 𝑏 ) , взять мнимую часть этогодуального числа и поделить на 𝑏 . Кроме того, так как выбор числа 𝑏 произволен, то можно выбрать его равным и избавиться от необходимости деления на 𝑏 .Для применения автоматического дифференцирования с помощью дуальных чисел в языке программированиядолжен быть определён тип дуальных чисел, а также арифметические операции и элементарные функциидля данного типа. Особенно просто это осуществляется для объектно-ориентированных языков и языков,поддерживающих перегрузку функций или множественную диспетчеризацию. Также следует иметь в виду, чтоавтоматическое дифференцирование с помощью дуальных чисел ограничено производными первого порядка. III. ДИНАМИЧЕСКАЯ ДИСПЕТЧЕРИЗАЦИЯ В ЯЗЫКЕ JULIA
Язык Julia [8] появился сравнительно недавно, однако успел приобрести популярность как язык для научныхвычислений. Будем предполагать, что читатели уже знакомы с этим языком, и кратко остановимся лишь наконцепции множественной диспетчеризации [9–12], которая лежит в основе языка, её понимание существеннодля дальнейшего изложения.
Динамическая диспетчеризация (dynamic dispatch) — это механизм, который позволяет выбрать конкретнуюреализацию полиморфной функции или оператора из множества и вызвать в конкретном случае.
Множественная диспетчеризация основывается на динамической диспетчеризации. В этом случае выборреализации полиморфной функции делается исходя из типа, количества и порядка следования аргументовфункции. Таким образом реализуется полиморфизм времени выполнения (runtime polymorphic dispatch). Крометермина «множественная диспетчеризация» также употребляется термин мультиметод .Механизм множественной диспетчеризации похож на механизм перегрузки функций и операторов, реали-зованный, например, в языке C++. Перегрузка функций, однако, осуществляется исключительно на стадиикомпиляции, тогда как множественная диспетчеризация должна работать также и на стадии выполненияпрограммы (полиморфизм времени выполнения).Julia поддерживает перегрузку функций на стадии компиляции в случае, если все типы данных, используемые вфункции, выводимы на стадии компиляции (type stable functions). При этом JIT-компилятор создаёт эффективныереализации для каждой комбинации типов аргументов.Если же типы данных невыводимы на стадии компиляции (type unstable), то включается механизм динамиче-ской диспетчеризации. Компилятор не сможет создать специализированную версию, а создаст универсальную,которая работает медленно.Такой подход позволяет совмещать быстродействие компилируемого языка со строгой статической типизациейс гибкостью интерпретируемого языка с динамической типизацией.
IV. РЕАЛИЗАЦИЯ ДУАЛЬНЫХ ЧИСЕЛ НА ЯЗЫКЕ JULIAA. Существующие реализации дуальных чисел на Julia
Авторам известны по меньшей мере две реализации дуальных чисел на языке Julia:• тип
Dual в модуле для автоматического дифференцирования ForwardDiff [7];• отдельный модуль DualNumbers [13], развитие которого заморожено в пользу ForwardDiff.В данном разделе мы опишем предлагаемую нами реализацию дуальных чисел на языке Julia. Она основанакак на встроенном типе обычных комплексных чисел
Complex [14], так и на модуле DualNumbers. Предлагаемаяв данной работе реализация служит исключительно в качестве учебного примера создания пользовательскоготипа данных для языка Julia.Во главу угла мы поставили ясность изложения, поэтому многие вычислительные оптимизации были умыш-ленно опущены в пользу большей ясности. Так, например, в модуле DualNumbers дуальное число 𝑧 = 𝑎 + ε 𝑏 может иметь комплексные компоненты 𝑎 и 𝑏 . Также в DualNumbers для реализации элементарных функций отдуальных чисел 𝑓 ( 𝑧 ) используется сторонний модуль, где определены символьные правила дифференцированиядля функций от действительного переменного. Это позволяет автоматизировать определение функции 𝑓 ( 𝑧 ) поформуле 𝑓 ( 𝑎 + ε 𝑏 ) = 𝑓 ( 𝑎 ) + 𝑓 ′ ( 𝑎 ) 𝑏 ε , так как не нужно явно выписывать производные 𝑓 ′ ( 𝑎 ) .В нашей реализации 𝑎 и 𝑏 — исключительно вещественные, а элементарные функции определены явно. B. Объявление структуры данных
Добавление любого пользовательского типа данных в Julia начинается с определения структуры данных. Длядуальных чисел определим структуру
Dual : struct Dual{T<:Real} <: Number"real part"x::T"imaginary part"y::Tend Структура содержит всего два поля: действительную x и мнимую y части дуального числа. Оба поля структурыдолжны иметь одинаковый параметрический тип T , который должен быть подтипом абстрактного типа Real ,о чём говорит оператор <: . Сам же тип Dual является подтипом абстрактного типа
Number . Таким образомдуальные числа встраиваются в уже существующую иерархию типов (см. рис. 1).Сразу после объявления новой структуры мы можем создать переменные типа Dual, используя конструкторыпо умолчанию. Мы обязательно должны задать оба поля структуры в виде аргументов функции Dual: z = Dual(1, 2)
Dual{Float64}(1, 2.2)
NumberDual Complex RealIntegerSignedInt8 Int16 Int32 Int64 Int128Bool UnsignedUInt8 UInt16 UInt32 UInt64 UInt128Rational AbstractFloatFloat16 Float32 Float64
Легенда:
Абстрактный типПримитивный типСтруктура
Рис. 1. Our Dual type in Julia built-in types hierarchy
Если тип аргументов одинаков, то явное указание параметра T не требуется. Если же аргументы имеютразный тип, то указание параметрического типа обязательно. В нашем примере оба аргумента приводятся ктипу Float64 . C. Перегрузка функций на примере show
Для возможности распечатывать значения переменных типа Dual, следует перегрузить функцию show . Этафункция предназначена для форматированной распечатки данных в стандартный вывод и в REPL. Функции write , print и println используются для вывода минимального текстового представления данных. Из нихпоследние две распечатывают данные только в стандартный поток вывода. Для простых данных можно неделать разницы между форматированным и компактным представлением и перегрузить только show . Тогдаостальные функции будут использовать show в своей работе.В терминах множественной диспетчеризации перегрузка функции означает добавление нового метода кфункции. При добавлении метода для функции show следует определить её следующим образом: Base.show(io::IO, z::Dual) = show(io, z.x, " + ", z.y, " ε ") Здесь использовался сокращённый синтаксис задания функций. Префикс
Base. нужен, так как стандартныеметоды функции show находятся в модуле
Base . Функции из этого модуля доступны для вызова непосредственнобез префикса, но при необходимости перегрузки применение
Base. обязательно.Стоит отметить, что наш вариант метода show упрощён, так как не учитывает частные случаи, например,отрицательную мнимую часть.Далее мы добавим методы для многих функций из
Base , для того чтобы наш тип
Dual стал в достаточноймере функциональным.
D. Дополнительные конструкторы
Создадим дополнительные конструкторы. Для этого следует перегрузить конструктор по умолчанию
Dual .Всего мы зададим 5 дополнительных конструкторов:
Dual(x::Real, y::Real)= Dual(promote(x, y)...)
Dual{T}(x::Real) where {T<:Real} = Dual{T}(x, 0)
Dual(x::Real) = Dual(promote(x, 0)...)
Dual{T}(z::Dual) where {T<:Real} = Dual{T}(z.x, z.y)
Dual(z::Dual) = Dual(promote(z.x, z.y)...)
Первый конструктор позволит нам не указывать всякий раз параметр T , если аргументы имеют разный тип.Функция promote осуществляет приведение типов переданных ей аргументов к общему типу и возвращаетрезультат в виде кортежа. Постфиксный оператор ... распаковывает картеж и передаёт его элементы в видеаргументов в функцию-конструктор. В ядре языка определены правила приведения для всех подтипов абстракт-ного типа Real , поэтому теперь конструктор будет корректно работать для любой комбинации аргументов,главное чтобы выполнялось правило
T<:Real . Например, следующий код сработает корректно:
Dual(1//2, π ) π * ε Мы передали в конструктор рациональное число (тип
Rational ) и встроенную глобальную константу (число 𝜋 ) типа Float64 . После чего сработало правило приведения типов и оба аргумента были приведены к болееобщему типу
Float64 .Второй и третий дополнительные конструкторы позволят не указывать мнимую часть, в случае если онаравна нулю:
Dual{Float32}(1) ε , 2 конструктор Dual(1//2) ε , 3 конструктор Конструктор №2 является параметрической функцией, которая объявляется с использованием конструкции where . Параметр T является подтипом абстрактного типа Real . Конструктор №3 работает аналогично конструк-тору №1. Четвёртый и пятый конструкторы позволяют передавать в качестве аргумента конструктора другоедуальное число.Для большего удобства также можно создать отдельную константу для мнимой единицы ε : const ε = Dual(0, 1) После перегрузки арифметических операций эта константа позволит создавать новые дуальные числа спомощью выражения, максимально близкого к алгебраической их записи: z = 1 + 2 ε E. Доступ к полям структуры
Структуры в Julia по умолчанию неизменяемы (immutable), то есть, один раз создав переменную z , мы можемобращаться к полям самой структуры через оператор точка, однако не можем модифицировать значение этихполей: @show z.x z.y z.x = 1 Для создания изменяемых (mutable) типов данных предусмотрена структура mutable struct , однако длянашего примера более подходящим будет именно неизменяемая структура, так как требование изменяемостинакладывает ограничения в производительности, что для числового типа нежелательно.В Julia нет модификаторов доступа к полям типа public или private , поэтому поля структуры всегда доступ-ны. Однако хорошим стилем программирования считается инкапсуляция полей структуры и предоставлениеинтерфейса для доступа к ним. В таком стиле выполнен, например, встроенный тип комплексных чисел.Существуют два способа инкапсуляции полей. Первый способ заключается в создании специальных интер-фейсных функций для доступа к полям. В нашем случае достаточно определить следующие функции:
Base.real(z::Dual) = z.xBase.imag(z::Dual) = z.yBase.reim(z::Dual) = (z.x, z.y)
Такие же функции определены во встроенном модуле Complex. Далее везде для получения полей струк-туры следует вызывать только эти функции вместо обращения к полям по именам напрямую. Это позволитразработчикам в дальнейшем производить рефакторинг структуры, например, переименовать или добавлятьновые поля. Обратную совместимость при этом легко сохранить, переписав лишь интерфейсные функции.Производительность при этом страдать не будет, так как оптимизирующий JIT-компилятор заменит вызовыэтих функций непосредственно на их код, так как они крайне просты (inline functions).Второй подход заключается в использовании функции getproperty , которая появилась начиная с версии0.7 языка Julia. Перегрузка этой функции позволяет задавать дополнительные имена для обращения к полямструктуры. Так, если для нашего примера написать следующий метод: function Base.getproperty(z::Dual, p::Symbol)if p in (:real, :a, :r) return getfield(z, :x)elseif p in (:imag, :b, :i) return getfield(z, :y) elsereturn getfield(z, p)endend то мы получим возможность обращаться к действительной и мнимой частям числа 𝑧 четырьмя разнымиспособами: z.x == z.a == z.r == z.real z.y == z.b == z.i == z.imag Такой подход более гибок, так как при изменении структуры для обратной совместимости достаточно будетмодифицировать лишь один метод getproperty , а не все интерфейсные функции, которых может быть много,если структура сложная. Кроме того, пользователю будут доступны все варианты обращения к полям структуры.Для типа Dual мы применили оба подхода, так как математически обосновано наличие функций real —эквивалента
Re( 𝑧 ) , и imag — эквивалента Im( 𝑧 ) . F. Унарные функции
В модуле
Base определены функции one и zero , которые возвращают нейтральный элемент по умножениюи по сложению соответственно. В случае дуальных чисел нейтралом по умножению является вещественнаяединица, а нейтралом по сложению — вещественный ноль.Для типа Dual определим следующие однострочные методы: Base.one(::Type{Dual{T}}) where T <: Real = Dual(one(T))Base.one(z::Dual) = Dual(one(z.x))Base.zero(::Type{Dual{T}}) where T <: Real = Dual(zero(T))Base.zero(z::Dual) = Dual(zero(z.x))
Первый вариант функций one и zero принимает в качестве аргумента тип данных и возвращает единицу илиноль этого типа соответственно. Так как для дуальных чисел это и из R , то достаточно вернуть дуальное числос нулевой мнимой частью. Вещественная же часть будет и параметрического типа T . Для всех стандартныхтипов T <: Real методы one и zero определены в модуле Base , чем мы и воспользовались.Второй вариант функций в качестве аргумента принимает уже конкретный объект типа
Dual и такжевозвращает нейтрал с помощью методов из модуля
Base .Для комплексных чисел в
Base определены функции сопряжения conj , абсолютного значения abs и аргумента arg : Base.conj(z::Dual) = Dual(z.x, -z.y)Base.abs(z::Dual) = abs(z.x)Base.abs2(z::Dual) = z.x^2function Base.arg(z::Dual)@assert !iszero(z.x)return z.y / z.xend
Для метода arg мы предусмотрели проверку на неравенство нулю действительной части дуального числа, таккак иначе мы получим деление на ноль. Использование стандартной функции iszero позволяет не беспокоитьсяоб учёте погрешностей представления вещественных чисел с помощью чисел с плавающей запятой. Макрос @assert выбросит исключение
AssertionError в случае равенства нулю действительной части.Функция inv определяет обратное число для данного числа z : function Base.inv(z::Dual)@assert !iszero(z.x)return Dual(1/z.x, -z.y/z.x^2)end Здесь также следует предусмотреть исключение для нулевой действительной части.0
G. Функции сравнения
В модуле
Base также определён целый ряд функций, возвращающих истину, в случае выполнение того илииного условия. Перечислим некоторые из этих функций:• isreal — число действительное;• isinteger — число целое;• isfinite — число конечно;• isnan — аргумент имеет тип
NaN ;• isinf — аргумент имеет тип Inf ;• iszero — число является нейтралом по сложению (ноль);• isone — число является нейтралом по умножению (единица).Методы для этих функций на случай типа Dual один в один повторяет методы для встроенного типа
Complex ,поэтому здесь мы не приводим их исходный код.Кроме унарных функций для дуальных чисел имеет смысл реализовать методы для оператора сравнения == . Операторы в Julia ничем не отличаются от функций и добавление методов для них проходит точно также.Единственное отличие заключается в необходимости добавления символа : после Base. , а так как оператор == состоит из двух символов, то необходимо обрамить его круглыми скобками: Base.:(==)(z::Dual, u::Dual) = (real(z) == real(u)) && (imag(z) == imag(u))
Имеет смысл сравнивать дуальные числа также с действительными, в случае если у дуального числа естьнулевая мнимая часть. Для того чтобы оператор коммутировал для аргументов разного типа, необходимоопределить два метода:
Base.:(==)(z::Dual, x::Real) = isreal(z) && real(z) == xBase.:(==)(x::Real, z::Dual) = isreal(z) && real(z) == x
H. Арифметические операции
Также необходимо перегрузить арифметические операторы, такие как унарные операторы + и - и бинарные + , - , * и / . Реализация операторов + , - и * тривиальна: Base.:+(z::Dual) = Dual(+z.x, +z.y)Base.:-(z::Dual) = Dual(-z.x, -z.y)Base.:+(z::Dual, u::Dual) = Dual(z.x + u.x, z.y + u.y)Base.:-(z::Dual, u::Dual) = Dual(z.x - u.x, z.y - u.y)Base.:*(z::Dual, u::Dual) = Dual(z.x * u.x, z.x * u.y + z.y * u.x)
В случае же оператора / следует учесть равенство нулю делителя: function Base.:/(z::Dual, u::Dual)@assert !iszero(u.x)return Dual(z.x/u.x, (z.y*u.x-u.y*z.x)/u.x^2)end Для бинарных операторов не надо отдельно реализовывать случаи аргументов разного типа, потому что вэтом случае должен срабатывать механизм приведения к общему типу с помощью функции promote .1 I. Приведение типов
Для работы механизма приведения типов следует определить правила приведения типов . Без этих правил,например, следующая операция будет оканчиваться ошибкой: >>Dual(1, 2) + 2ERROR: promotion of types Dual{Int64} and Int64 failed to change any arguments
Эта ошибка возникает из-за того, что нет правила, по которому можно определить общий тип для чисел типа
Dual{Int64} и типа
Int64 .Для определения такого правила нужно добавить метод для функции promote_rule из модуля
Base . Варифметической операции с дуальным числом может участвовать любое число, относящееся к типу
Real , таккак это число можно интерпретировать как дуальное с нулевой мнимой частью:
Base.promote_rule(::Type{Dual{T}}, ::Type{S}) where {T<:Real, S<:Real} = Dual{promote_type(T, S)}
Также следует предусмотреть правило, которое будет работать в случае действий с двумя дуальными числамис разными параметрическими типами T и S . В этом случае нужно найти общий тип для типов T и S : Base.promote_rule(::Type{Dual{T}}, ::Type{Dual{S}}) where {T<:Real, S<:Real} = Dual{promote_type(T, S)}
J. Возведение в рациональную степень
Для возведения числа в степень служит оператор ^ , который также следует перегрузить для целой ирациональной степени. В случае целой степени мы используем формулу ( 𝑎 + ε 𝑏 ) 𝑛 = 𝑎 𝑛 + 𝑛𝑏𝑎 𝑛 − ε . Следует предусмотреть несколько особых случаев:• если 𝑛 = 0 , то 𝑧 𝑛 = 1 ;• если 𝑛 = 1 , то 𝑧 𝑛 = 𝑧 ;• если 𝑛 > , то формула работает для любых 𝑎 и 𝑏 ;• если 𝑛 < , то обязательно выполнение условия 𝑎 ̸ = 0 .Если учесть все эти случаи, то получим следующую реализацию: function Base.:^(z::Dual, n::Integer)::Dualx, y = reim(z)if n == 0return one(z)elseif n == 1return zelseif n > 1return Dual(x^n, n*y*x^(n-1))else if iszero(x)throw(DomainError(:x, "negative exponentiation is only defined for real(z)!= 0"))elsereturn Dual(x^n, n*y*x^(n-1))endendend Если действительная часть дуального числа равна нулю ( iszero(x) ) и при этом 𝑛 < , то функция выбрасываетисключение DomainError (выход за область допустимых значений).2Для рациональной степени используется более сложная формула: ( 𝑎 + ε 𝑏 ) 𝑛𝑚 = 𝑎 𝑛𝑚 (︂ ε 𝑛𝑚 𝑏𝑎 )︂ = 𝑎 𝑛𝑚 + ε 𝑛𝑚 𝑏𝑎 𝑛𝑚 − , для которой следует предусмотреть случаи чётного и нечётного 𝑚 . При нечётном 𝑚 в область допустимыхзначений входят 𝑎 < , а в случае чётного 𝑚 следует ограничиться только 𝑎 > : function Base.:^(z::Dual, q::Rational)::Dualx, y = reim(z)n, d = numerator(q), denominator(q)if n == 0return one(z)elseif isodd(d)return Dual(x^q, q*y*x^(q-1))elseif x < 0throw(DomainError(:x, "even radical for dual number z is only defined for real(z)>=0"))elsereturn Dual(x^q, q*y*x^(q-1))endendend Имеет смысл отдельно перегрузить функции для квадратного и кубических корней sqrt и cbrt : function Base.sqrt(z::Dual)::Dualx, y = reim(z)if x < 0throw(DomainError(:x, "sqrt for dual number z is only defined for real(z)>=0"))elsereturn Dual(√x, y/(2*√x))endendBase.cbrt(z::Dual) = Dual(cbrt(z.x), z.y*cbrt(z.x)/(3z.x)) K. Элементарные функции
Элементарные функции вычисляются по формуле 𝑓 ( 𝑎 + ε 𝑏 ) = 𝑓 ( 𝑎 ) + 𝑓 ′ ( 𝑎 ) 𝑏 ε . Следует учитывать, что многиеиз них не определены для случая нулевой действительной части числа 𝑧 = 𝑎 + ε 𝑏 : function Base.exp(z::Dual)e = exp(z.x)return Dual(e, e * z.y)endBase.log(z::Dual) = Dual(log(z.x), z.y/z.x)Base.log(b, z::Dual) = Dual(log(b, z.x), (z.y/z.x) * log(z.x)) Base.sin(z::Dual) = Dual(sin(z.x), +z.y*cos(z.x))Base.cos(z::Dual) = Dual(cos(z.x), -z.y*sin(z.x))Base.tan(z::Dual) = Dual(tan(z.x), z.y/cos(z.x)^2)Base.cot(z::Dual) = Dual(cot(z.x), -z.y/sin(z.x)^2)Base.asin(z::Dual) = Dual(asin(z.x), z.y/sqrt(1-z.x^2)) Base.acos(z::Dual) = Dual(acos(z.x), -z.y/sqrt(1-z.x^2))Base.atan(z::Dual) = Dual(atan(z.x), z.y/(1+z.x^2))Base.acot(z::Dual) = Dual(acot(z.x), -z.y/(1+z.x^2))
Base.sinh(z::Dual) = Dual(sinh(z.x), z.y*cosh(z.x))Base.cosh(z::Dual) = Dual(cosh(z.x), z.y*sinh(z.x))Base.tanh(z::Dual) = Dual(tanh(z.x), z.y/cosh(z.x)^2)Base.coth(z::Dual) = Dual(coth(z.x), z.y/sinh(z.x)^2)
V. РЕЗУЛЬТАТЫ
В работе описана предварительная реализация дуальных комплексных чисел и основных операций над нимина языке Julia.
VI. ОБСУЖДЕНИЕ
Реализация дуальных комплексных чисел на Julia полностью опирается на механизм множественной дис-петчеризации. Таким образом мы не только реализовали некоторый набор операций над дуальными числами,но и продемонстрировали мощность этого механизма. Также следует отметить, что в отличие от реализациитипа
Dual в пакете автоматического дифференцирования
ForwardDiff [7], предложенная нами реализация болеечистая. Например, в вышеупомянутом пакете в качестве коэффициента перед дуальной комплексной единицейможет выступать обычное комплексное число. Понятно, что это вызвано спецификой применения дуальныхчисел в этом пакете для автоматического дифференцирования. Но данного типа числа скорее относятся ккватернионам [15, 16], нежели к собственно комплексным числам.
VII. ЗАКЛЮЧЕНИЕ
В работе был сделан прототип реализации обобщённых комплексных чисел на языке Julia, а именно реализациядуальных комплексных чисел. Наличие механизма множественной диспетчеризации в языке Julia очень упрощаетреализацию новых числовых типов в рамках существующей инфраструктуры языка программирования. Авторыпредполагают в дальнейшем расширить данный прототип до полной реализации обобщённых комплексныхчисел и обобщённых кватернионов.
БЛАГОДАРНОСТИ
Публикация подготовлена при финансовой поддержке РФФИ в рамках научного проекта № 19-01-00645. [1] Яглом И. М., Розенфельд Б. А., Ясинская Е. У. Проективные метрики // Успехи математических наук. — 1964. —Т. 19, № 5 (119). — С. 51–113.[2] Яглом И. М. Комплексные числа и их применение в геометрии // Математика, ее преподавание, приложения иистория. — 1961. — Т. 6 из Математическое просвещение, сер. 2. — С. 61–106.[3] Розенфельд Б. А., Яглом И. М. О геометриях простейших алгебр // Математический сборник. — 1951. — Т. 28(70),№ 1. — С. 205–216.[4] Kulyabov D. S., Korolkova A. V., Gevorkyan M. N. Hyperbolic numbers as Einstein numbers // Journal of Physics:Conference Series. — 2020. — 5. — Vol. 1557. — P. 012027.1–5.[5] Cayley A. IV. A sixth memoir upon quantics // Philosophical Transactions of the Royal Society of London. — 1859. — 1. —Vol. 149. — P. 61–90. [6] Klein F. Ueber die sogenannte Nicht-Euklidische Geometrie // Gauß und die Anfänge der nicht-euklidischen Geometrie. —Wien : Springer-Verlag Wien, 1985. — Bd. 4 von Teubner-Archiv zur Mathematik. — S. 224–238.[7] Forward Mode Automatic Differentiation for Julia. — 2020. — Access mode: https://github.com/JuliaDiff/DualNumbers.jl .[8] The Julia Language. — 2020. — Access mode: https://julialang.org/ .[9] Bezanson J., Chen J., Chung B., Karpinski S., Shah V. B., Vitek J., Zoubritzky L. Julia: dynamism and performancereconciled by design // Proceedings of the ACM on Programming Languages. — 2018. — 10. — Vol. 2, no. OOPSLA. —P. 1–23.[10] Zappa Nardelli F., Belyakova J., Pelenitsyn A., Chung B., Bezanson J., Vitek J. Julia subtyping: a rational reconstruction //Proceedings of the ACM on Programming Languages. — 2018. — 10. — Vol. 2, no. OOPSLA. — P. 1–27.[11] Bezanson J., Edelman A., Karpinski S., Shah V. B. Julia: A fresh approach to numerical computing // SIAM Review. —2017. — 1. — Vol. 59, no. 1. — P. 65–98. — arXiv : 1411.1607.[12] Bezanson J., Karpinski S., Shah V. B., Edelman A. Julia: A Fast Dynamic Language for Technical Computing. — 2012. —P. 1–27. — arXiv : 1209.5145.[13] Julia package for representing dual numbers and for performing dual algebra. — 2020. — Access mode: https://github.com/JuliaDiff/ForwardDiff.jl .[14] Official Julia language GitHub repository. — 2020. — Access mode: https://github.com/JuliaLang/juliahttps://github.com/JuliaLang/julia